The turbulence energy source driving varies with different astrophysical environments. Here, to simplify considerations, we divide the driving sources into two categories: external and internal driving mechanisms. The former has been considered as a common way for turbulence driving in the simulations of interstellar turbulence, due to the convenience of setting external forces. The most typical cases are supernova explosions in ISM (Spitzer Jr, 1978; Dib et al., 2006 ; Ostriker, 2009; Chamandy and Shukurov, 2020), merger events and active galactic nuclei (AGN) outflows in the intercluster medium (Chandran, 2005; Subramanian et al., 2006; Lazarian et al., 2020) and outflow in molecular clouds (Carroll et al., 2009), and young stellar objects (Federrath, 2016). In contrast, the latter is also accepted generally, such as magnetic reconnection driving/mediating turbulence in various regimes (Huang and Bhattacharjee, 2016; Cerri and Califano, 2017; Franci et al., 2017; Kowal et al., 2017; Dong et al., 2018), the magneto-rotational instability-generated turbulence (Sellwood and Balbus, 1999; Gardiner and Stone, 2005; Fromang, 2013; Kunz et al., 2016; Riols et al., 2017; Zhdankin et al., 2017). Indeed, for realistic astrophysical processes such as molecular clouds and supernovae, their drivings may be involved in internal and external mechanisms simultaneously.

Sometimes the complexity of turbulent driving is manifested by the coexistence of multiple driving mechanisms in the same environment. One of the most classic cases is the driving process in the ISM, where various drivers may occur simultaneously: expanding shells and gravitational instability at about 1,000 pc scales arising from galaxy mergers; expanding shells, MRI and cloud collisions at about 100 pc scales from supernova explosion; and stellar feedback at 10-pc to sub-pc scales by outflow. In the turbulence research community, for simplicity, one usually uses external driving with the solenoidal or compressible approach, which is more effective at limited numerical resolution. In recent years, internal driving mechanisms have gradually become the focus of research, but only at the cost of expensive numerical resources (Beresnyak, 2017; Kowal et al., 2017).

In the spirit of understanding the nature of turbulence from first principles, we here provide the visco-resistive equations that govern the evolution of magnetized plasma as follows ¹ : ∂ ρ ∂ t + ∇ ⋅ ρ v = 0 , (1) ρ ∂ ∂ t + u ⋅ ∇ u = j × B c − ∇ p g + ρ ν ∇ 2 u , (2) ∂ p g ∂ t + u ⋅ ∇ p g + Γ p g ∇ ⋅ u = 0 , (3) ∂ B ∂ t − ∇ × u × B = η ∇ 2 B , (4) ∇ ⋅ B = 0 , (5) describing the continuity, Eq 1, the momentum conservation, Eq 2, the adiabatic invariance, Eq 3, the magnetic induction, Eq 4, and the divergence of magnetic field, Eq 5, respectively. Here, p g = c s 2 ρ is the gas pressure, c s = Γ p g / ρ the speed of sound, t the evolution time of the fluid, ν the kinematic viscosity, η = c ²/(4πσ) the magnetic diffusivity, σ the conductivity, and Γ the adiabatic index. Other physical quantities have their usual meanings. The interaction of plasma fluids and magnetic field was formulized in the aforementioned equations, by which a series of physical quantities characterizing the properties of the MHD turbulence can be better understood.

Let us replace the operator ∇ by ∇/L and set u = u u ̃ = L / t u ̃ in Eq 2 and B = B B ̃ in Eq 4, where L is the characteristic scale length at the typical velocity u, as well as u ̃ and B ̃ are dimensionless. After a simple analytical derivation (including only the viscosity term in Eq (2)), we obtain ∂ u ̃ ∂ t + u ̃ ⋅ ∇ u ̃ = ν u L ∇ 2 u ̃ , (6) ∂ B ̃ ∂ t − ∇ × u ̃ × B ̃ = η L u ∇ 2 B ̃ . (7)

Through the coefficients of the equations, we define two dimensionless physical quantities R e = L u ν , R m = L u η , (8) where R _e is the kinetic Reynolds number characterizing the ratio of the viscous timescale to the advection timescale, and R _m is the magnetic Reynolds number denoting the ratio of resistive diffusion timescale to the advection timescale. Considering the temperature T ≃ 10⁴K related to η ≃ 0.42 × 1 0 7 T eV − 3 / 2 c m 2 / s , the scale length L ≃ 1 pc, and the velocity u ≃ 1 km/s, one can obtain a large R _m ≃ 10¹⁶ as an example of turbulence setting (e.g., Brandenburg and Subramanian, 2005). The relative dominance between the kinematic viscosity and magnetic diffusivity is described by the magnetic Prandtl number: P r = R m R e = ν η . (9)

Similarly, combining the adiabatic invariance, Eq 3 and momentum conservation, Eq 2, one can obtain the sonic Mach number (Biskamp, 2003 for a detailed derivation) M s = u / c s , (10) for the high-β plasma case. This parameter can be used to describe the information of compressibility of MHD turbulence. For the low-β plasma case, the momentum conservation Eq 2 is dominated by the magnetic pressure. On multiplying Eq 2 by u = u u ̃ , the quantity characterizing the magnetization strength is called the Alfvénic Mach number M A = u / V A , (11) where V A = B / 4 π ρ is the Alfvénic velocity.

Considering current density ( ∇ × B = 4 π c j ) related to the magnetic field by Ampere’s law, the Lorentz force can be written in the following form: j × B c = B ⋅ ∇ B 4 π − ∇ B 2 8 π , (12) where the first term on the right-hand side describes the magnetic tension force and the second term describes the magnetic pressure force. The plasma parameter β, denoting the ratio of gas pressure to magnetic one, is defined as β = p g B 2 / 8 π , (13) which can be rewritten as β = 2 M A 2 M s 2 , (14) with the sonic and Alfvén Mach numbers.

As for the description of the properties of helicity in MHD turbulence, one can introduce another two quantities termed as the magnetic helicity and cross helicity, which are, respectively, written as H M = ∫ V A ⋅ B d V , (15) H C = ∫ V u ⋅ B d V . (16)

Here, A is the vector potential with B = ∇ ×A. The former can measure the twist and linkage of the magnetic field lines, and the latter can reflect the overall correlation of magnetic and velocity fields (Falgarone and Passot, 2003).

Compressibility is an important addition to the MHD turbulence theory. In fact, in a real astrophysical environment such as a molecular cloud, the MHD turbulence is essentially compressible and is composed of three plasma modes, i.e., Alfvén, slow, and fast modes. The idea of the decomposition of MHD turbulence into individual modes was originally formulated in Dobrowolny et al. (1980), where they mainly dealt with small amplitude perturbations in incompressible MHD turbulence.

The actual decomposition in a numerical simulation has been implemented using Fourier transformation (CL02; CL03). Covering the gas pressure-dominated and magnetic pressure-dominated turbulence regimes, they revealed the spectrum and anisotropy in the compressible MHD turbulence. Typically, the Alfvén mode in the high- and low-β regimes follows the scaling slope and scale-dependent anisotropy of GS95. E A k ∝ k ⊥ − 5 / 3 , k ‖ ∝ k ⊥ 2 / 3 . (28)

The slow mode, as a passive mode, is dominated by the Alfvén mode, and has properties similar to those of the Alfvén mode: E s k ∝ k ⊥ − 5 / 3 , k ‖ ∝ k ⊥ 2 / 3 . (29)

The fast mode, different from the Alfvén and slow modes, possesses isotropic properties similar to acoustic turbulence, with the following scaling relations: E f k ∝ k − 3 / 2 , k ∥ ∝ k ⊥ . (30)

Later, Kowal and Lazarian (2010) advanced the studies on mode decomposition using a wavelet method to treat large amplitude perturbations. This method enables us to trace the direction of the local magnetic field and has significant advantages over the Fourier transform done in the frame of reference of the mean magnetic field.

To explain the numerical simulations in Maron and Goldreich (2001), several studies that have attracted community attention proposed a particular process termed as dynamical alignment to modify the GS95 spectrum from k ^−5/3 to k ^−3/2 (Boldyrev, 2005; Boldyrev, 2006; Mason et al., 2006). At present, there are numerical results that are compatible with a dynamic alignment model (e.g., Chandran et al., 2015). Moreover, solar-wind observations seem to find some scale-dependent alignment (e.g., Wicks et al., 2013) or the emergence of three-dimensional anisotropy (e.g., Chen et al., 2011) which would also point toward a scale-dependent alignment of some sort (see also Wang T et al., 2020 for similar findings from the end of the MHD cascade to the kinetic scales).

At scales close to the dissipation scale, another modification of the GS95 theory was proposed by Mallet and Schekochihin (2017) for their intermittency model and Mallet et al. (2017) for their reconnection-mediated turbulence model (see also Loureiro and Boldyrev, 2017), the testing/confirming of which is difficult due to current limitations in computational abilities.

However, the anisotropy predicted in Bolryrev’s work is inconsistent with some recent numerical simulations (Beresnyak, 2019; Beresnyak and Lazarian, 2019). Many studies claimed that the deviation from the spectral index -5/3 is transient, namely, localized in the vicinity of the injection scale, and does not extend to the whole inertial range (Beresnyak and Lazarian, 2010; Beresnyak, 2013; Beresnyak, 2014). To the best of my knowledge, this issue is still open. Anyway, the current attempts still do not change the paradigm of the modern MHD turbulence theory.

Synchrotron fluctuation techniques reviewed in this article are developed closely based on the modern understanding of MHD turbulence theory (e.g., GS95; LV99; CL02) without more concern for its controversial appearance. Fortunately, modifications attempted to the fundamental theory are not expected to seriously affect the measurement of the techniques. This is because the development of measurement techniques is independent of a specific scaling slope of turbulence and is constrained in the strongly turbulent range, while modifications to turbulence theory focus on the driving or dissipation scale.

In a completely homogeneous field, the synchrotron emissivity of relativistic electrons is expressed by (Ginzburg and Syrovatskii, 1965) j h ν , θ = ∫ d E N E F E , ν , θ . (33)

The integrand F , representing the power per frequency unit emitted by a single electron of given energy and pitch angle, is given by (see Rybicki and Lightman, 1979 for more details on the derivation) F E , ν , θ = 3 e 3 B sin θ m e c 2 ν ν c ∫ ν / ν c ∞ K 5 / 3 t d t , (34) where K _5/3(t) is the modified Bessel function of order 5/3, f ( ν ν c ) = ν ν c ∫ ν / ν c ∞ K 5 / 3 ( t ) d t is often referred to as a kernel function, and ν c = 3 E 2 e B ⁡ sin ⁡ θ / 4 π m e 3 c 5 is a characteristic frequency of synchrotron emission.

The power per frequency of a single electron, F , is divided into the perpendicular and parallel components with respect to the magnetic field B (Rybicki and Lightman, 1979). F ⊥ E , ν , θ = 3 e 3 B sin θ 2 m e c 2 f ν ν c + g ν ν c , (35) F ∥ E , ν , θ = 3 e 3 B sin θ 2 m e c 2 f ν ν c − g ν ν c , (36) where g ( ν / ν c ) = ν ν c K 2 3 ( ν ν c ) and K 2 3 is the modified Bessel function of order 2/3. As a result, the degree of linear polarization can be calculated by Π ν = F ⊥ − F ∥ F ⊥ + F ∥ = g ν / ν c f ν / ν c (37) for a single electron radiation. Replacing F in Eq 33 with F ∥ or F ⊥ , we can obtain the perpendicular (j _⊥) and parallel (j _∥) components of synchrotron emissivity for electron populations. Hence, we have the degree of linear polarization for the electron populations. Π ν = j ⊥ − j ∥ j ⊥ + j ∥ . (38)

After integration, the total emissivity (Eq 33) can be simplified as j h ν , θ ∝ B ⊥ p + 1 / 2 ν − p − 1 / 2 , (39) where the spectral index of photons is α = p − 1 2 and B _⊥ = B sin θ.

In a completely random field, the calculation of synchrotron emissivity will be more complex compared to a completely homogeneous field, due to the presence of a stochastic direction. The synchrotron emissivity is formulated as (Crusius and Schlickeiser, 1986) j r ν = 1 4 π ∫ 0 2 π d φ ∫ 0 π d θ ⁡ sin ⁡ θ j h ν , θ , (40) by averaging the total spontaneously emitted power for the homogeneous field j _h (ν, θ) over all possible values of the polar (θ) and azimuthal (φ) angles. This equation implies that randomness of the magnetic field is maintained ranging from the Larmor radii of the radiating electrons to the size of the emitting source.

Defining x = ν sin θ/ν _c, Eq 40 can be rewritten as j r ν ∝ B ∫ 0 ∞ d E N E R x , (41) where x = ν ( c 1 B E 2 ) , c 1 = 3 e 4 π m e 3 c 5 and R ( x ) = x 2 ∫ 0 π d θ ⁡ sin ⁡ θ ∫ x / sin ⁡ θ ∞ K 5 / 3 ( t ) d t is associated with the Bessel function of order 5/3. The analytical results demonstrated that j _r has the same power–law relation for frequency as that of j _h (Crusius and Schlickeiser, 1986). Different from the situation of homogeneous field, j _r is related to the total magnetic field strength B and not just its perpendicular component B _⊥ in j _h expression.

When involving in the synchrotron self-absorption, one can adopt Eq 34 to obtain the following the self-absorption coefficient: κ ν , θ = c 2 8 π ν 2 ∫ E 2 d d E N E E 2 F E , ν , θ d E . (42)

We can, thus, write total radiation intensity by the radiative transfer equation I ν = j κ 1 − e − τ , (43) where τ = κℓ is the absorption optical depth, and ℓ is the size of the emitting region.

For the electrons with a power-law distribution, the absorption coefficient can be reduced to κ ν , θ = p + 10 / 3 p + 2 μ ν , θ = G p e 3 2 π m e 3 e 2 π m e 3 c 5 p / 2 K B ⁡ sin ⁡ θ p + 2 / 2 ν − p + 4 / 2 , (44) where μ ( ν , θ ) = p + 2 p + 10 / 3 κ ( ν , θ ) , and the coefficient G ( p ) = 3 4 Γ ( 3 p + 2 12 ) Γ ( 3 p + 22 12 ) is only associated with the electron exponent.

To obtain polarization properties, we herein introduce the following coefficients in two mutually perpendicular directions: κ ⊥ ν , θ = κ ν , θ + μ ν , θ , (45) κ ∥ ν , θ = κ ν , θ − μ ν , θ . (46)

By adapting Eq 43, we have radiation intensity components I ⊥ ν = j ⊥ κ ⊥ 1 − e − κ ⊥ ℓ , (47) I ‖ ν = j ‖ κ ‖ 1 − e − κ ‖ ℓ . (48)

As for an optically thin medium, i.e., κ _⊥,∥ ℓ ≪ 1, the degree of linear polarization is consistent with Eq. 38. In the case of optically thick medium, i.e., κ _⊥,∥ ℓ ≫ 1, Π ν = j ⊥ / κ ⊥ − j ∥ / κ ∥ j ⊥ / κ ⊥ + j ∥ / κ ∥ (49) can be used to calculate the degree of linear polarization.

According to the commonly used method of correlation analysis, the two-point correlation of the polarized intensities can be expressed as (LP16) 〈 P X 1 , λ 1 2 P * X 2 , λ 2 2 〉 = ∫ 0 L d z 1 ∫ 0 L d z 2 〈 P i X 1 , z 1 P i * X 2 , z 2 e 2 i λ 1 2 Φ X 1 , z 1 − λ 2 2 Φ X 2 , z 2 〉 , (54) where P* denotes the conjugate of the polarization intensity, X ₁ and X ₂ represent two different spatial positions on the plane of the sky, and λ ₁ and λ ₂ correspond to two wavelengths. Eq 54 can not only characterize the correlation of the spatial separation but also the dispersion of frequencies. Based on this, LP16 proposed two techniques, namely the polarization spatial analysis (PSA) and the polarization frequency analysis (PFA).

In analogy with the velocity channel analysis (VCA) technique (Lazarian et al., 2002), the PSA correlates the polarization signal at different spatial points of the position–position frequency (PPF) space that has similarities with the position–position volume (PPV) space. At a fixed wavelength, the PSA can be expressed as ⟨ P X 1 P * X 2 ⟩ = ∫ 0 L d z 1 ∫ 0 L d z 2 e 2 i ϕ ̄ λ 2 z 1 − z 2 × ⟨ P i X 1 , z 1 P i * X 2 , z 2 ⟩ × e − 2 λ 4 ⟨ △ Φ X 1 , z 1 − △ Φ X 2 , z 2 2 ⟩ , (55) where ϕ ̄ is the mean Faraday rotation measure density and ΔΦ is the fluctuation of Faraday rotation measure. This represents the correlation of polarization intensity as a function of spatial separation R = X ₂ − X ₁ at the same wavelength.

The PFA that is analogous to the velocity coordinate spectrum (VCS) correlates the polarization information along a single line of sight at the same spatial position. At a fixed spatial position, the PFA for different wavelengths is expressed by 〈 P λ 1 2 P * λ 2 2 〉 = ∫ 0 L d z 1 ∫ 0 L d z 2 e 2 i ϕ ̄ λ 1 2 z 1 − λ 2 2 z 2 〈 P i z 1 P i * z 2 e 2 i λ 1 2 △ Φ z 1 − λ 2 2 △ Φ z 2 〉 , (56) which is complementary to the PSA, so this enables a more comprehensive acquisition of the properties of MHD turbulence.

Another important aspect of the development of the synchrotron fluctuation technique is the construction of the statistics of polarization intensity derivative with respect to the squared wavelength λ ². The correlation statistics of dP/dλ ² is formulated as 〈 d P X 1 d λ 2 d P * X 2 d λ 2 〉 = ∫ 0 L d z 1 ∫ 0 L d z 2 e 2 i ϕ ̄ λ 2 z 1 − z 2 × 〈 P i X 1 , z 1 P i * X 2 , z 2 Δ Φ X 1 , z 1 × Δ Φ X 2 , z 2 e − 2 λ 2 i Δ Φ X 1 , z 1 − Δ Φ X 2 , z 2 〉 , (57) which is more sensitive to the fluctuations of Faraday rotation than the polarization intensity correlation given in Eq 55. The introduction of this relationship allows for more information such as the spectrum of Faraday rotation to be obtained. In addition to considering the correlation function, one can also use the structure function of polarization intensity to reveal the anisotropy of MHD turbulence (see LP16 for detailed formulation).

With the correlation analysis of synchrotron fluctuations, the power spectrum of synchrotron polarization is described through the Fourier transform of the two-point correlation tensor as E P K = 1 2 π 2 ∫ d R e − i K ⋅ R ⟨ P X 1 P * X 2 ⟩ , (58) where R = X ₁ − X ₂ ∝ 1/K is the spatial lag vector on the plane of the sky. Similarly, the power spectrum of the polarization derivative with respect to λ ² can be written as E dP K = 1 2 π 2 ∫ d R e − i K ⋅ R ⟨ d P X 1 d λ 2 d P * X 2 d λ 2 ⟩ (59) for analyzing the multi-frequency observations. Once the power spectrum distribution of the observed information is obtained through the aforementioned two equations, we can obtain the distribution characteristics of turbulent energy in space, and also understand several properties related to turbulence, such as the scaling slope, energy dissipation, and injection scales.

Following LP12, we can define the normalized correction function for an arbitrary observable ϝ as ξ ϝ R = 〈 ϝ X ϝ X + R 〉 − 〈 ϝ X 〉 2 〈 ϝ X 2 〉 − 〈 ϝ X 〉 2 , (60) where X and R denote a spatial position and a separation between any two points on the plane of the sky, respectively. This is related to the normalized structure function of the observable ϝ D ̃ ϝ = 2 1 − ξ ϝ , (61) by which we have a multi-pole expansion M ̃ n R = 1 2 π ∫ 0 2 π e − i n φ D ̃ ϝ R , φ d φ , (62) where n represents the order of multi-pole and φ is the polar angle. The order of the multi-pole expansion has following meaning: n = 0 represents the monopole reflecting the isotropy of statistics, n = 1 denotes the dipole which is zero due to the rotational symmetry of the structure function, n = 2 is called the quadrupole moment revealing the anisotropy, and n = 4 is the octupole moment which has a negligible contribution for the anisotropy. Hence, we define the ratio of the quadrupole moment to the monopole one (termed quadrupole ratio, ref. LP16) R Q = M ̃ 2 R M ̃ 0 R = ∫ 0 2 π e − 2 i φ D ̃ ϝ R , φ d φ ∫ 0 2 π D ̃ ϝ R , φ d φ (63) to characterize the degree of anisotropy, and the absolute value of which is called the quadrupole ratio modulus.

Alternatively, the ratio of structure functions in the perpendicular direction to each other can also depict the level of anisotropy and is given by R S = 〈 | ϝ ∥ x + R ∥ − ϝ ∥ x | 2 〉 〈 | ϝ ⊥ y + R ⊥ − ϝ ⊥ y | 2 〉 , (64) where R = ( R ∥ 2 + R ⊥ 2 ) 1 / 2 is the spatial separation on the plane of the sky. When R S = 1 , statistics represents an isotropy while any deviation from R S = 1 implies an appearance of anisotropy.

Combining the local alignment theory of MHD turbulence (GS95; LV99) with the synchrotron polarization fluctuation predictions (LP16), the gradient techniques were proposed to measure the properties of MHD turbulence (Lazarian et al., 2017; Lazarian and Yuen, 2018a). This is analogous to velocity gradient techniques (González-Casanova and Lazarian, 2017; Yuen and Lazarian, 2017). The feasibility lies in that magnetic field and velocity fluctuations enter symmetrically, according to the theory of Alfvénic turbulence. In the framework of GS95 MHD turbulence theory, it is expected that the alignment of the magnetic field gradient direction is perpendicular to the local magnetic field direction with the eddies (for other different anisotropy models, the sensitivity to this alignment effect may significantly depend on the scales). Since the synchrotron fluctuations are elongated in the direction of the mean magnetic fields, the magnetic fields can be traced by the direction of the spatial gradient of the observed signal.

When involved in Faraday depolarization of synchrotron radiation fluctuations, synchrotron polarization gradients (SPGs) and synchrotron polarization derivative gradients (SPDGs) can reveal more detailed information on the turbulent magnetic field. In particular, the latter can sample the spatial regions close to the observer with the sampling depth controlled by the radiation wavelength (see Figure 1), which provides a new approach to map the 3D distribution of Galactic magnetic fields.

Specifically, the de-correlation of the Faraday rotation measure is introduced as (LP16; Zhang et al., 2020) λ 2 Φ = 0.81 λ 2 ∫ 0 L eff d z n e B ‖ = 2 π , (65) where L _eff is an effective spatial depth of the Faraday rotation sampling. We can thus write the ratio of the depth sampled by the Faraday rotation to the emitting region size as L eff L ≈ 2 π λ 2 L 1 ϕ , (66) with ϕ = m a x ( 2 σ ϕ , ϕ ̄ ) . Here, σ _ϕ denotes the root mean square and ϕ ̄ Faraday rotation density. According to Eq 66, the strong and weak depolarizations can be characterized by L _eff/L < 1 and L _eff/L > 1, respectively. In the case of the strong depolarization, the entire space is divided into two different parts as shown in Figure 1, from which we can know that only the part of z < L _eff experiences the depolarization while the part of z > L _eff cannot contribute to the polarization correlation. When observing that the frequency interval △ν = c/λ ²△λ is small enough, the difference of P (λ ₁) and P (λ ₂) in the range of [L ₂, L ₁] is written as △ P ≈ ∫ L 2 L 1 d z P i X , z e 2 i λ 2 Φ X , z , (67) from which the expression of the SPDGs is given by ∇ d | P | d λ 2 ∼ λ − 2 ∇ | P λ 1 − P λ 2 | . (68)

From another perspective, namely, based on the spatial gradient analysis of Stokes parameters Q and U on the complex plane, Gaensler et al. (2011) proposed the synchrotron polarization gradients to constrain the sonic Mach number in warm or ionized ISM. After a simple mathematical derivation, the spatial gradient of polarization intensity can be expressed as | ∇ P | = ∂ Q ∂ x 2 + ∂ U ∂ x 2 + ∂ Q ∂ y 2 + ∂ U ∂ y 2 . (69)

Subsequently, considering rotational and translational invariance of coordinates, various diagnostic quantities were developed by Herron et al. (2018).

LP12 predicted that the power spectrum, correlation, and structure functions of synchrotron radiation fluctuations, termed two-point statistics, can recover the scaling slope of turbulent magnetic fields. Numerical simulations show that the correlation function does not work in revealing scaling slopes, but it can be used to explore the correlation scale and the effect of the electron spectral index on turbulence (Herron et al., 2016; Lee et al., 2016; Zhang et al., 2016). The main reason that we think is that LP16 performed theoretical studies under the condition of the power-law correlation model, from which they predicted the possibility of obtaining the scaling slope for the corresponding model. The correlation function itself, however, has a divergent behavior for numerical simulation or observational data, because Stokes parameters Q and U contain the information of negative values arising from the change in the direction of the turbulent magnetic field.

Numerical results demonstrated that the scaling slope of the underlying MHD turbulence can be recovered approximately by the two-point structure function of synchrotron radiation fluctuations (Lee et al., 2016). From the perspective of numerical simulations, the direct calculation of the spatial structure function is more time-consuming, compared with the correlation function or power spectrum. Alternatively, one can first get the statistics of the correlation function, provided that it maintains convergence, and then calculate the structure function using the relation of SF(R) = 2 [CF(0) − CF(R)]. In general, the two-point structure functions are unable to correctly reproduce the scaling of a spectrum with a slope ∼ − 3 or steeper (sometimes problems may arise already at slopes slightly shallower than − 3, if the underlying spectrum is not a perfect power-law). For instance, Lee et al. (2016) measured the scaling slope of 5/3 with the change in a numerical resolution. It was pointed out that the scaling slope of the two-point structure function, calculated directly from 2D synthetic data, was still less than 5/3 expected for the Kolmogorov spectrum, even with a 8192² resolution. Thus, using a two-point structure function to extract the scaling slope of turbulence requires more precise statistics with a higher numerical resolution. Alternatively, the multi-point structure functions are recommended to achieve better convergence of the results (see Cho and Lazarian, 2009; Cho, 2019).

Obtaining the scaling slope through the power spectrum of synchrotron fluctuations should be one of the most efficient ways, as the calculation of the power spectrum invokes the fast Fourier transform in the wavenumber space to effectively accelerate the process. Adopting both synthetic and MHD turbulence simulation data, Lee et al. (2016) obtained the power spectrum and its derivative with respect to the wavelength of synchrotron polarization arising from synchrotron polarization radiation together with Faraday rotation fluctuations. This study succeeded in obtaining the scaling slope and demonstrated that simulations are in agreement with the theoretical prediction in LP16. In addition, using simulated interferometric observations on how to recover the scaling slope was also attempted. Notice that the study mentioned previously is based on the spatially coincident synchrotron emission and Faraday rotation regions. Furthermore, considering spatially separated and compounded synchrotron radiation and Faraday rotation fluctuations, Zhang et al. (2018) studied how to extract the scaling slope of synchrotron polarization intensities. In the short wavelength range, the power spectra revealed fluctuation statistics of the perpendicular component of turbulent magnetic fields, and the spectra reflected the fluctuation of the Faraday rotation in the long-wavelength range. As shown in Figure 2, the effects of telescope angular resolution and noise structure do not hinder the use of power spectrum methods to obtain scaling slopes of MHD turbulence in the large-scale inertial range.

Motivated by one-point statistics in LP16 (termed PFA), which is characterized by the variance of polarized emission as a function of the squared wavelength along a single line of sight, Zhang et al. (2016) studied how to extract the scaling slope of underlying MHD turbulence. To depict the level of turbulence, they defined a ratio η of the standard deviation of the line of sight turbulent magnetic field to the line of sight mean magnetic field. As shown in Figure 3, a large ratio (η ≫ 1) characterizing a turbulent field dominated regime reflects the polarization variance 〈 P 2 〉 ∝ λ − 2 , and a small ratio (η ≲ 0.2) characterizing a mean field-dominated regime reveals the polarization variance 〈 P 2 〉 ∝ λ − 2 − 2 m = λ − 10 / 3 with m = 2/3 for the Kolmogorov scaling. As a result, the turbulent spectral index was successfully recovered by the polarization variance in the case of the mean field-dominated regime, i.e., a small η. At the same time, it was pointed out that the change in cosmic ray electron spectral indices cannot affect the scaling slope measurement of magnetic turbulence. This provides a new method for recovering the scaling index of MHD turbulence and agrees with the theoretical prediction of LP16. Interestingly, this technique has been successfully applied to explain the depolarization in the optical wavebands from active galactic nuclei (Guo et al., 2017).

As described in Section 2, the prediction of the anisotropy (eddy’s elongation along the local magnetic field direction) in turbulent cascade is an important contribution of GS95. It is challenging to test the local anisotropy of MHD turbulence from an observational point of view because the observational signal is inevitably involved in the integration along the line of sight in the observer’s frame of reference. Note that integrated synchrotron fluctuations along the line of sight were predicted by LP12 to be anisotropic, with the eddy’s major axis aligned with the direction of the mean magnetic field instead of the local field within the MHD turbulence volume.

The spatial second-order structure function of radiation fluctuations is a very direct way to reveal the existence of anisotropy in MHD turbulence, although the anisotropy here does not directly correspond to the local anisotropy favored by the modern MHD turbulence theory. By visualizing a contour map for 2D structure function of observations, one can easily see the anisotropic distribution, from which the major axis direction of the elongated structure will qualitatively characterize the mean magnetic field orientation (see Figure 1 in Wang R. Y et al. (2020) for an example). Considering the ratio of two measurement directions perpendicular to each other from observed polarization intensities, called an anisotropic coefficient, we can quantitatively determine the degree of the anisotropy of the MHD turbulence in the global frame of the reference; refer to the work of Lee et al. (2019) and Wang R. Y et al. (2020) for more details.

A new more precise method to measure the anisotropy is termed as the quadrupole ratio described in Section 4. This method was proposed in LP12 that derives many analytical formulae including the predictions from the Alfvén, slow, and fast modes in the high or low β plasma regime. A numerical test of theoretical predictions on the quadrupole ratio was attempted by Herron et al. (2016) based on MHD turbulence simulation. We have to point out that their testing used a simplified version rather than a general formula. In particular, they adopted the compressible MHD turbulence data to study the isotropic version of formulae, which fails to provide proof of the usability of the analytical formula. Considering the angle change in between the mean magnetic field direction and the line of sight, recently tested the theoretical prediction of LP12 in detail and confirmed the correctness of the analytical formulae, providing convenience to study the anisotropy of MHD turbulence by using the analytical expressions.

How the anisotropic structure changes with observation frequency were studied in Lee et al. (2019) by the quadrupole ratio. Their studies showed that in the high-frequency range, the observed polarization exhibits the averaged structures of both foreground and background regions, while in the low-frequency range, the large-scale structures are wiped out by a strong Faraday depolarization and the small-scale anisotropy reflects the magnetic field structure from the foreground region.

By decomposing the compressible MHD turbulence into three plasma modes, i.e., Alfvén, slow, and fast modes, Wang R. Y et al. (2020) demonstrated that the quadrupole ratio statistics of the polarization intensity from the Alfvén and slow modes present a significant anisotropy, while the fast mode statistics show isotropic structures (see Figure 4). Those findings are in agreement with the earlier direct numerical simulations provided in CL02. As a result, the quadrupole ratio method can be applied to the measurement of the anisotropy of MHD turbulence in various astrophysical environments.

The traditional method of measuring the magnetic field direction is called the synchrotron polarization vector, which is based on the fact that the observed polarization vectors are perpendicular to the magnetic field directions within the emitting source. Note that the limitation of this method cannot be applied to observations with Faraday rotation.

Using synchrotron intensity gradients (SIGs) together with both MHD turbulence simulation and PLANCK archive data, Lazarian et al. (2017) demonstrated the practical applicability of the gradient technique described in Section 4, comparing the measured directions of magnetic fields arising from the gradient technique with that of the traditional polarization vector method. Furthermore, SIGs were advanced in Lazarian and Yuen (2018a) to SPGs and SPDGs. Subsequently, to explore the magnetic field measurement capabilities of synchrotron gradient techniques, SPGs and SPDGs were generalized and applied to a variety of possible astrophysical environments in a series of studies including different turbulence regimes, the gradient of various diagnostic quantities, realistic observations of Galactic diffuse media, as well as spatially separated synchrotron radiation and Faraday depolarization regions.

Several statistical methods, such as Genus, Tsallis, skewness, and kurtosis statistics, are sensitive to the magnetization level M _A of MHD turbulence (e.g., Burkhart et al., 2012; Herron et al., 2016). In addition, the correlation function (e.g., Lazarian et al., 2002; Burkhart et al., 2014; Esquivel et al., 2015) and structure function (Hu et al., 2021; Xu and Hu, 2021) analysis are also suggested as alternative techniques of measuring M _A. At present, these studies are based on the velocity information analysis from MHD turbulence simulations. Similarly, the adoption of synchrotron fluctuations should be on the agenda.

In analogy with velocity gradient techniques, Lazarian et al. (2018b) first discussed that synchrotron gradient techniques should be able to measure the magnetization level. The feasibility of the two methods was explored in Carmo et al. (2020) to obtain the level of magnetization from synchrotron polarization fluctuation, namely the top-base and the circular standard deviation. They claimed that the signal-to-noise ratio was more severe for the top-base method, but still reliable with the standard deviation method. In short, the finding on the power-law relation between M _A and synchrotron gradient statistics is important for determining M _A from the current or future available data cubes.