Studying the Local Magnetic Field and Anisotropy of Magnetic Turbulence by Synchrotron Polarization Derivative

Since the hydrodynamic type motions of eddies that mix magnetized fluid give rise to Alfvénic perturbations for the magnetic field, MHD turbulence was usually considered a dynamical process of the turbulent cascade dominated by the Alfvén wave interactions. The Alfvén Mach number, describing the strength of magnetic turbulence, is defined as M_A = V_L/V_A, where V_L is turbulent injection velocity at the scale L_in and is the Alfvén velocity. Correspondingly, the largest turbulent eddies would also reflect the Alfvén Mach number at the injection scale L_in, i.e., with B_L being the perturbation of the magnetic field B at the injection scale L_in. Then, the relative perturbations of magnetic fields and eddy velocities are related in a symmetrical way (Lazarian et al. 2020)L

where is the fluctuation of the magnetic field B at the scale l of the turbulent fluid and v_l is the corresponding velocity fluctuation.

The modern theory of MHD turbulence (GS95) can be understood as a collection of anisotropic eddies aligning with the direction of local magnetic fields around them, that is, the motions of the eddies with parallel scales l_∥ and perpendicular scales l_⊥ are elongated along the direction of the local magnetic field. Assuming that the motion perpendicular to the direction of the local magnetic field has a Kolmogorov scaling of , the anisotropic relation of eddy scales can be written as

This scaling was originally predicted in the case of strong trans-Alfvénic (M_A ∼ 1) incompressible turbulence, with the critical balance condition of . Later, the GS95 anisotropic scaling (see Equation (2)) was generalized to both M_A > 1 and M_A < 1 cases (LV99, see also Lazarian 2006). The former is called a super-Alfvénic turbulence regime with V_L > V_A. In the limiting case of , the turbulence presents an essentially hydrodynamic Kolmogorov behavior, i.e., v_l = V_L(l/L_in)^1/3, as the weak magnetic field with lower magnetic energy than kinetic energy has a marginal effect on turbulent dynamics. When v_l = V_A, the hydrodynamic-like behavior of the turbulence cascade would transition into a strong turbulence regime at the scale . In the range from l_A to the dissipation scale l_dis, the GS95 scaling can be maintained.

The latter, V_L < V_A, shows weak turbulence ranging from L_in to the transition scale l_B = L_inM_A², in which the magnetic field fluctuations are quasi-2D and perpendicular to the direction of magnetic fields (LV99, Galtier et al. 2000). At scales less than l_B, the relation between the eddy major axis and its minor axis is present again, i.e.,

Setting M_A = 1, we recapture the original GS95 relation (see Equation (2)).

In compressible MHD turbulence, the compressible turbulent motions are associated with slow and fast modes (GS95; Lithwick & Goldreich 2001). Numerical simulations (Cho & Lazarian 2002) confirmed that Alfvénic and slow modes have the same anisotropy (), while fast modes present an isotropic cascade ().

The interactions of relativistic electrons with turbulent magnetic fields emit synchrotron radiation that reveals the properties of magnetic field fluctuations. For simplicity, we consider a homogeneous and isotropic distribution of the electrons with a power-law form of N(γ) ∝ γ^2α − 1, where γ is the electron energy and α is a spectral index of the electrons.⁵ Thus, the synchrotron emission intensity as a function of radiative frequency ν can be written as (Ginzburg & Syrovatskii 1965)

where Γ (x) is a gamma function, B_⊥ is the magnetic field component perpendicular to the line of sight, and L is the emitting-region size along the line of sight. Other parameters (e, m_e, and c) have their usual meanings. With the spectral index of relativistic electrons, we would gain the fraction polarization degree and the linearly polarized intensity P = pI. The observable Stokes parameters Q and U are related to the polarized intensity by the polarization angle ψ₀, i.e., and . When it comes to a Faraday rotation effect, the polarization angle is expressed as with the wavelength λ. The Faraday rotation measure is given by

Here, represents a vector in the plane of the sky, is the parallel component of the magnetic field, and n_e is the thermal electron density.

With the Stokes parameters Q and U, the complex polarization vector is formulated as (LP16)

where indicates the intrinsic polarized intensity density that is treated as wavelength-independent; this simplified treatment would not change the numerical results stated below (see also Zhang et al. 2018 for confirmation).

The condition of decorrelation of the Faraday rotation measure is introduced as

where L_eff is an effective spatial depth of the Faraday rotation sampling. Therefore, the ratio of the depth sampled by the Faraday rotation to the emitting-region size can be written as

where , and σ_ϕ and are the root mean square and the mean of the Faraday rotation measure density () fluctuation, respectively. Using Equation (8), the strong and weak Faraday rotation would be characterized by L_eff/L < 1 and L_eff/L > 1, respectively. When L_eff/L < 1, Equation (6) can be split two parts, in which only z < L_eff suffers from the Faraday depolarization, while z > L_eff cannot contribute to the polarization measure (see Figure 15 in LY18 for an illustration).

As for a chosen wavelength λ, the complex polarization vector (see Equation (6)) in an effective Faraday depolarization region could be rewritten as

Considering two neighboring wavelengths λ₁ and λ₂, which respectively correspond to the close spatial positions L₁ and L₂, we can obtain the differences and :

which reflects the local magnetic turbulence within the region. From an observational point of view, providing the Stokes parameters Q and U at two neighboring frequencies, we can obtain the details of the magnetic field fluctuation originating from a locally spatial region. However, note that there is a nonlinear correspondence between L and λ, as expressed in Equation (8). In practice, we consider the statistics of the synchrotron polarization intensity derivative with respect to the squared wavelength,

to reveal the local magnetic field fluctuations.

The spatial correlation and structure functions of (any) physical variable have been traditionally employed for studying anisotropic properties. The former, i.e., the correlation function, is written as

where denotes an average over the entire volume of interest. As for the latter, the commonly used second-order structure function is given by

It is obvious that the second-order structure function is formally related to the correlation function. Both of them can be used to characterize the eddy anisotropy in MHD turbulence. For the structure function, considering two measurement directions perpendicular to each other, we can define their ratio of individual structure function as

which would reflect the eddy anisotropy of the turbulence structure.

As for the correlation function, it is easy to calculate the correlation function map with a periodic boundary condition through the fast Fourier transform. Studying the non-periodic locality of synchrotron polarization derivative map, we use the Hockney method (Hockney 1968) to solve the open-boundary convolution problem, which helps us decrease computational complexity (Yuen et al. 2018). The direction of the major axis of the correlation (structure) function contour determines the eddy major axis reflecting the magnetic field orientation.

The spectral correlation function (SCF) was initially introduced in Lazarian et al. (2002) to study the anisotropy of turbulence in channel maps of the position–position–velocity 3D space. Here, we adopt this method to explore the correlation of 3D data cubes (as a function of position–position–frequency f(x, y, ν)). The formula used is given by (e.g., Padoan et al. 2003)

where is a lag vector. The contour of would characterize the spatial scales that the spectral features begin to change. The 2D correlation can be used to produce the 1D spectrum of the correlation versus the lag length by an azimuthal average.

Following LP12, the normalized correlation function (NCF) of (any) physical variable can be written as

where is a separation vector between any two spatial points on the plane of the sky. Similarly, the normalized structure function is expressed as

Therefore, we can obtain the quadrupole moment ratio arising from the variable ,

to reveal the spatial anisotropy, where R is a radial separation and φ is the polar angle.

In the study of synchrotron gradient measurements, the Sobel operator will be used to compute an approximation of the gradient of the 2D image. This method adopts two 3x3 kernels that are convolved with the original 2D image to obtain approximations of the horizontal and vertical derivatives. And then the recipe of sub-block averaging (Yuen & Lazarian 2017) for the gradient map obtained is used to determine the gradients of the subregion on the 2D image. In each subregion of interest, we use a Gaussian fitting method for the gradient angles to get an optimal direction characterized by the peak of the fitting.

For the magnetic field data cubes obtained by numerical simulation, we can know the intrinsic direction of the underlying magnetic field. Therefore, the magnetic field direction obtained by correlation function anisotropy and gradient techniques can be used to make a comparison with the inherent magnetic field direction. In the case of correlation function anisotropy, the major axis direction of the contour determines the magnetic field direction, whereas in the case of gradient technique, the direction of rotated 90° gradients identifies the directions of the magnetic field.

We adopt a reduction factor to measure the correspondence between the measured magnetic field direction and the intrinsic magnetic field direction (see González-Casanova & Lazarian 2017)

which is called the alignment measure (AM), analogous to the Rayleigh reduction factor in dust alignment theory suggested by Greenberg (1968). The parameter θ in Equation (19) is the angle between the measured magnetic field direction and the intrinsic magnetic fields. AM = ±1 represents a perfect alignment, while AM = 0 no alignment.

The following equations are used to describe the interstellar magnetic turbulence environment where the synchrotron radiation is emitted:

Here, is the gas pressure, t is the evolution time of the turbulent fluid, is the current density, and is an external driving force. Additionally, an isothermal equation of state would close the above equations.

Numerically, the third-order-accurate hybrid, essentially non-oscillatory code is used to obtain data cubes reaching steady state. In our simulation, an external magnetic field with is set along the x-axis (horizontal) direction and the random turbulence is driven by a solenoidal driving force at the wavenumber k = 2.5. The resulting 3D data cubes with a numerical resolution of 512³, including the properties of three 3D magnetic fields, three 3D velocities, and one 3D density, are listed in Table 1, where they are depicted by the Alfvénic and sonic Mach numbers.

Table 1.  Data Cubes with a Numerical Resolution of 512³ Generated in the Simulation of Compressible MHD Turbulence

Model	MA	Ms		β
run1	0.65	0.48	0.164	3.668
run2	0.70	0.87	0.579	1.295
run3	0.55	4.46	0.467	0.030
run4	0.50	9.92	0.465	0.005

Note. δ B_rms denotes the root mean square of a random magnetic field, the regular magnetic field, and the plasma parameter.

Download table as:  ASCIITypeset image

Adopting a cylindrical coordinate system fixed on the eddies (following CV00), we obtain the local magnetic field via

where is parallel to the major-axis direction of eddies. In the local reference frame of the magnetic field is the corresponding second-order structure function as

with and .

As a case in point, with the run3 magnetic field data cubes listed in Table 1, Figure 1 shows the structures of scale-dependent anisotropic eddies. The eddy anisotropies in 3D space are plotted in the left panel, while the eddy anisotropies within the 100th 2D slice are plotted along the major and minor axes of eddies (middle panel) and along the x and y axes (right panel), respectively. As seen in Figure 1, the scale-dependent anisotropies of the eddies of MHD turbulence can indeed be revealed in the local reference frame, as predicted in GS95 and confirmed in CV00. The smaller the scale is, the more pronounced the degree of anisotropy is. This figure shows the result of the anisotropic structure of an eddy within a 2D slice, and compares the synchrotron observations presented below. However, even for a plane-like (a slice) geometry, from the synchrotron observational data we have no a priori information about the magnetic field and eddy. As a result, we can study their properties only in the artificially selected “local” coordinate system to observe the eddy structures, as done in the right panel of Figure 1.

Zoom In Zoom Out Reset image size

For the sake of simplicity and for the purposes of simulating real observations and the feasibility of various statistical testing techniques, we consider a scenario with spatially coincident polarized synchrotron emission and Faraday rotation regions, with the assumption of a 1 kpc spatial scale, thermal electron density 0.01 cm⁻³, and magnetic field strength 1 μG. Focusing on low-frequency observations from the Low Frequency Array for radio astronomy (LOFAR), we can synthesize data cubes of Stokes parameters I, Q, and U ranging from frequencies ν = 10 to 240 MHz in the bandwidth of Δν = 0.46 MHz, on the basis of the 3D simulation data listed in Table 1. This choice of bandwidth can ensure a sufficiently small spatial interval in the line-of-sight direction to achieve successful measurements of the local magnetic field.

Figure 2 depicts the maps at individual frequency points in service of obtaining a qualitative understanding of the synthesized observation data. The upper panels of Figure 2 are the images of synchrotron polarization intensity , while the middle panels correspond to images of the synchrotron polarization intensity derivative relative to the squared wavelength dP/dλ². These images are computed for sub-Alfvénic and subsonic turbulence, using the run1 data of Table 1. Their structures are extended along the (x-axis) horizontal direction because the mean magnetic field is set to this direction. We find that the percentage of strong noise-like structure prevents observations at the lowest frequency, 10 MHz. In general, these small-scale structures can be smoothed by a Gaussian filter technique.

Zoom In Zoom Out Reset image size

The lower panels of Figure 2 are the images of the Faraday rotation measure, using data cubes from run1 to run3 of Table 1. The values of Faraday rotation measure for sub-Alfvénic and subsonic turbulence (left and middle lower panels) are smaller than those of sub-Alfvénic and supersonic turbulence (lower right panel). In cases of supersonic turbulence, this is caused by the occurrence of high-density (clump) regions arising form shock-wave interactions.

In this section, we adopt the SCF to explore the correlation of multifrequency data cubes as a function of spatial lag separation. Numerical calculation is implemented by extracting the source code from a python package called TURBUSTAT (Koch et al. 2017). With the synchrotron simulation data cubes generated by the run3 listed in Table 1, the resulting image of SCF of dP/dλ² is shown in the left panel of Figure 3, in which the solid contour line indicates the 2D fitting by an elliptical power-law model, obtaining an optimal index of −0.191. Interestingly, we see the solid contour line extending along the x-axis, which would reflect the mean magnetic field direction. Furthermore, an azimuthal average of the SCF image produces a 1D spectrum of the power-law slope of −0.182 (left upper panel), which is fitted by a linear least-squares method with the weights of the inverse squared standard deviation. We find that the 2D fitting gives a slightly greater index than the 1D fitting.

Zoom In Zoom Out Reset image size

Similarly, the fitting results of other synchrotron simulation data are as follows: −0.194 (2D) and −0.181 (1D) for run1; −0.175 (2D) and −0.162 (1D) for run2; −0.179 (2D) and −0.172 (2D) for run4. All the above fittings, in the lag range from 1 to 40 pixels, demonstrate the correlation index around −0.2. On the other hand, we also study the possible correlation of SCF values of synchrotron polarization intensity P but find a weak correlation, i.e., a fitting index less than −0.07 for four sets of synchrotron simulation data above. The SCF of P reflects the information of the same projected mean magnetic field measured at different frequencies, while the SCF of dP/dλ² reveals links between different local magnetic fields. Accordingly, we point out that an SCF analysis of dP/dλ² can be used to reveal the scale-dependent anisotropies of MHD turbulence. Similar to studies of velocity channel maps (Lazarian et al. 2002), it is expected that the frequency channel maps will reveal more properties of the underlying magnetic field; this will be studied in detail in the future.

As mentioned in Section 4.2, since an a priori local reference system cannot be fixed on the eddies for observational data, we select the Cartesian coordinate system on the plane of the sky. Figure 4 plots the contour of the structure function of dP/dλ², which is calculated by run1 of Table 1 at frequencies of 10 MHz ((a) panel), 100 MHz ((b) panel), 200 MHz ((c) panel), and 235 MHz ((d) panel). According to Equation (8), those frequencies are related to the spatial positions 1.4, 134.8, 570.7, and 790.2 pc (a 1 kpc extent length adopted along the line of sight), respectively, where the anisotropic structures of dP/dλ² qualitatively reflect those of the eddies. Exceptionally, no correlation for the 10 MHz simulation can be seen due to the influence of extremely strong numerical noise.

Zoom In Zoom Out Reset image size

Quantitative analysis of the structure function of dP/dλ² is carried out in the left panel of Figure 5, according to Equation (14). As shown, the statistics at high frequencies present stronger anisotropy than those at low frequencies. The degree of anisotropy of dP/dλ² decreases with increasing spatial scale, which reveals the scale-dependent anisotropy of MHD turbulence (GS95). Moreover, the quadrupole moment ratio (refer to Equation (18)) modulus for dP/dλ², as a function of the spatial separation, is plotted in the right panel of Figure 5 at the corresponding frequencies. The quadrupole moment ratio moduli decrease with increasing spatial separation, which is generally consistent with the results from the structure function of dP/dλ² (left panel). In particular, the quadrupole ratio modulus is more sensitive in the anisotropy analysis, as shown in the 10 MHz curves.

Zoom In Zoom Out Reset image size

On the basis of the analysis of the synchrotron polarization derivative in Section 2.3, the correlation function anisotropy of dP/dλ² is able to reveal the magnetic turbulence anisotropy of the local spatial region. Measuring the major axis direction of the contour from the correlation function of dP/dλ², we could trace the local magnetic field direction. Similar to Yuen et al. (2018) and Yang et al. (2020), we rotate the contour map to determine the orientation of the major axis. Compared with the earlier numerical implementation, which was measured at a single spatial point with a rotation step size of 1 degree, this work makes two improvements, including the decrease of the step size to 0.5 degree and the average of the measurement direction of multiple spatial positions.

The alignment measure, between intrinsic directions of the local magnetic fields and directions predicted by the CFA, is shown in Figure 6 as a function of radiation frequency (upper panel) and of the spatial depth along the line of sight (lower panel). In cases of sub-Alfvénic and subsonic turbulence, i.e., M_A < 1 and M_s < 1, the CFA can trace the local magnetic field direction apart from the low-frequency part. This is because the presence of strong numerical noise affects the anisotropic distribution of the correlation function of dP/dλ², and the lack of anisotropy causes CFA not to work. When smoothing the small-scale noise structure in the low-frequency range, we find that AM can be improved to some extent, but this is not shown here. In cases of sub-Alfvénic and super-sonic turbulence, i.e., M_A < 1 and M_s > 1, the reverse occurs: numerical noise moves to a lower-frequency part. However, AM values decrease with increasing frequency. In this regard, the formation of a shock wave results in strong density fluctuations. The density-dominated fluctuations should alter the fluctuations of random magnetic fields, leading to the change of the anisotropic structure of the polarization derivative map, on which the CFA tracing is strongly dependent.

Zoom In Zoom Out Reset image size

As seen in Equation (8), the relationship between multiple frequencies can be mapped to the local spatial positions in the direction of the line of sight. The lower panel of Figure 6 shows the AM versus the corresponding spatial depth of Faraday rotation sampling. The sampling depth for sub-Alfvénic and subsonic turbulence is greater than that for super-Alfvénic and supersonic turbulence because the high-density fluctuations of the latter (super-) produce stronger Faraday rotation than that of the former (sub-). It is understandable that the strong Faraday rotation sampling depth is shallower. As a result, the CFA is effective for tracing the local magnetic field direction in cases of sub-Alfvénic and subsonic turbulence.

Based on the gradient measurement method stated in Section 3.4, we first compute the gradients of dP/dλ² by the Sobel operator and divide the whole gradient map obtained into 8 × 8 sub-blocks, each of which has a numerical resolution of 64 × 64 pixels. We then carry out Gaussian fitting for the sub-block gradient map, the peak of which is considered the optimal gradient orientation in each sub-block region. Finally, we average all AMs, and the alignment between the intrinsic magnetic field direction and the gradient one, to obtain an AM value at each frequency.

The resulting AM is plotted in Figure 7 as a function of the frequency (upper panel) and of the spatial depth of Faraday rotation sampling (lower panel), on the basis of four sets of data cubes listed in Table 1. In this figure, most simulations show that AM remains almost unchanged as frequency increases, an exception being the scenario in which the supersonic Mach (M_s ∼ 10) turbulence has an AM that decreases slightly. In particular, the presence of numerical noise approaching 10 MHz does not impede the tracing of the magnetic field direction. It turns out that the synchrotron polarization derivative gradient (SPDG) technique provides the ideal AM for tracing the local magnetic field. In comparison with the results provided in Figure 6 from the CFA measurement (Section 5.4), we find that the capability of the SPDG is greater than that of the CFA in the measurement of magnetic field directions, which is in agreement with studies of velocity gradient techniques (Yuen et al. 2018; Yang et al. 2020). On the other hand, our previous studies demonstrated that the CFA of polarization intensities has the advantage of distinguishing compressible turbulence modes (LP12; Herron et al. 2016; Lee et al. 2019; Wang et al. 2020). Therefore, we believe that the CFA and gradient analysis of the synchrotron polarization derivative are complementary for measurements of magnetic fields.

Zoom In Zoom Out Reset image size