Gradient Technique Theory: Tracing magnetic field and obtaining magnetic field strength

The Gradient Technique (GT) is a technique that employs the properties of MHD turbulence to study magnetic field from observational data. These properties include the anisotropy of velocity and magnetic field fluctuations in the MHD cascade as well as the alignment of these fluctuations with the local direction of magnetic field. The latter property is absolutely crucial for the ability of the GT to trace dynamically important magnetic field.¹¹ 1 If magnetic field is not dynamically important, it can be advected by velocity motions and passively reflect the statistics of velocity fluctuations. Thus the velocity gradients can trace the magnetic field even in this case. However, we do not consider this situation in the present paper.

The well-accepted hydrodynamic Kolmogorov turbulence [70] can be understood using a mental picture of an eddy cascade. Traditionally, MHD turbulence was considered a wave-type phenomenon, similar to acoustic turbulence [64, 73] with Alfvénic waves having isotropic distribution. This was shown to be an incorrect picture (Montgomery & Turner 104, 39). Note that further, we focus on the Alfvénic turbulence cascade only, as this cascade determines the properties of gradients we deal with. The separation of this cascade from the other two MHD motions is possible because the back-reaction of slow and fast modes of the Alfvénic cascade is marginal in the strong Alfvénic turbulence regime (39, Cho et al. 23, [122], see also appendix of [107]). Thus, the scaling properties of Alfvénic cascade stay unaltered in the compressible media in good agreement with numerical calculations [22].

The anisotropic distribution in 39 picture corresponds to Alfvénic wave vectors getting perpendicular to magnetic field. The corresponding motions are equivalent to eddies according to 87. The eddy motions were considered impossible in the presence of dynamically important magnetic field. However, the theory of turbulent reconnection (87, see Lazarian et al. 77 for a review) predicts that the magnetic field reconnection happens within one eddy turnover time. This enables unconstrained eddy motions in the direction perpendicular to the magnetic field. As a result, the Kolmogorov picture of turbulence is revived, but with the restriction that the eddy cascade is expected to involve only in the eddies with rotation axes aligned with the direction of magnetic field that surrounds the eddies.

Naturally, the magnetic field in question is those involved in the Alfvénic motion, i.e. the field that percolates the eddies. This is the local direction of magnetic field, which is an important advance in understanding Alfvénic cascade. In other words, due to turbulent reconnection, the turbulent motions perpendicular to the magnetic field are not constrained by the back-reaction of magnetic field. This presents the favorable way of turbulent cascading with most energy concentrated in eddies rotating perpendicular to the local direction of the magnetic field.²² 2 A common misconception about the MHD theory is related to the fact that the concept of the local direction of the magnetic field is not a part of the original 39 idea. This concept of local magnetic field direction was introduced and proven numerically in [24] and subsequent studies, e.g., [100, 21]. Dealing exclusively with Alfvénic turbulence simplifies our further treatment as velocities and magnetic field are related through the Alfvénic relations. Further, in particular in §9 we discuss complicated situations involving non-Alfvénic turbulence.

Accepting the natural constraint that eddies perform the motions that minimize the magnetic field bending, i.e. that due to fast reconnection they easily mix up magnetic field around them in the perpendicular direction, one can get the Kolmogorov-type condition for such perpendicular Alfvénic motions

where $v_{l}$ is the velocity of perpendicular to magnetic field eddies and the scale $l_{\bot}$ and $t_{casc,l}$ is the energy cascading time, i.e. time of the energy transfer from a scale $l_{\bot}$ scales to a smaller eddy. We use the symbol $\bot$ to explicitly reflect the fact that the cascade involves ”perpendicular” eddies. For unconstrained energy transfer from scale to scale, the cascading time is similar to the Kolmogorov turbulence

The eddies, however, can change their parallel scale $l_{\|}$. Indeed, the magnetic field mixing with the period $l_{\bot}/V_{l}$ induces a wave with a period $l_{\|}/V_{A}$, where $V_{A}$ is the Alfvén velocity. The latter condition, i.e.

was termed critical balance in the 39 theory of MHD turbulence.³³ 3 The critical balance assumes that in Eq. (3) it is $l_{\|}$ rather than the injection scale $L_{inj}$. If this is not true, the cascade is modified. The eddies cannot freely mix and the turbulence in the weak MHD turbulence regime with $l_{\|}=L_{inj}$ and perpendicular to magnetic field motions scaling $v_{l}\sim l^{1/2}$ (39, Galtier et al. 36. See also Mallet et al. 99, Meyrand et al. 103 for proofs of critical balance.). This does not change the picture for the gradient studies as the motions are still perpendicular to magnetic field. However, in the latter study the critical balance was derived for the corresponding parallel and perpendicular scales obtain in relation to the mean magnetic field rather than to the local magnetic field. The subsequent studies [24, 100, 21] testified that the critical balance and the relations between the $l_{\bot}$ and $l_{\|}$ that follow from combining Eqs. (1,2) and (3), i.e.

are valid only if $\|$ and $\bot$ scales are calculated in terms of the local direction of magnetic field. It follows naturally from our presentation of Alfvénic cascade that the scaling of eddy velocities is Kolmogorov-type (see Eq. (1), i.e.

where $b_{l}$ is the fluctuation of magnetic field.

Apart from Alfvénic motions, slow and fast modes are present in MHD turbulence. The slow modes are present even in the limit of incompressible turbulence and they follow the scaling of Alfvén modes (39,Cho & Lazarian 21, Cho & Lazarian 22, Lithwick & Goldreich 94). The fast modes are present in compressible MHD turbulence and it was demonstrated in [21, 22] that that provide an isotropic cascade similar to the acoustic turbulence. For magnetic field tracing with the GT, the properties of Alfvén and slow modes are employed. The energy of those usually dominates the energy in the fast modes (see Cho & Lazarian 21). For subsonic turbulence slow modes dominate the formation of density fluctuations and therefore $\rho_{l}\sim l_{\bot}^{1/3}$ and are elongated with the axis ratio scaling given by Eq. (4). The situation is more complicated, however, for supersonic turbulence where the shocks distort significantly the density structure within MHD turbulence (see Beresnyak et al. 8, Kowal et al. 72).

All the considerations above apply to sub-Alfvénic, i.e., the magnetic field fluctuations $\delta B$ that are less than the underlying regular magnetic field $B_{0}$. For Alfvénic motions, one can define the Mach number

where $V_{inj}$ is the injection velocity of the Alfvénic cascade. One identifies $M_{A}<1$ with sub-Alfvénic turbulence and $M_{A}=1$ with the trans-Alfvénic turbulence. If turbulence is super-Alfvénic, i.e. $M_{A}>1$, at the injection scale $L_{inj}$ the motions are hydrodynamic-like down to the scale

where the turbulence transfers to the magnetic field-dominated regime. All our arguments about trans-Alfvénic turbulence are applicable provided that use $l_{A}$ as the effective injection scale of turbulence.

For scales larger than $l_{A}$, the magnetic fields are moved by hydrodynamic flows and the gradients of velocities perpendicular to magnetic fields (see Lazarian et al. 92, Hu et al. 2024). For scales less than $l_{A}$, the velocity gradients are similar to the case of trans-Alfvénic turbulence. For this type of turbulence and $M_{A}<1$ turbulence, the gradients are also perpendicular to magnetic field [91, see]. In the latter case, this is the consequence of the Alfvénic turbulence anisotropy (39; Kandel et al. 2016, 2017a, henceforth KLP16, KLP17).

The properties of MHD turbulence depend on magnetization that is measured using the Alfvén Mach number $M_{A}$ which is the ratio of the turbulent injection velocity $V_{inj}$ to the Alfvén velocity $V_{A}$, i.e. $M_{A}=V_{inj}/V_{A}$. For sub-Alfvénic turbulence (i.e. $M_{A}<1$) at scales smaller than

the velocity scaling is given by 87:

where $l_{\bot}$is measured in the local system of reference. ⁴⁴ 4 At separations larger than $L_{inj}M_{A}^{2}$ as we mentioned earlier, the turbulence scaling follows the so-called weak cascade (See Santos-Lima et al. 2021), $v_{l}\approx V_{inj}(l_{\bot}/L_{inj})^{1/2}$. The gradient studies, however, usually employ small scales as weak turbulence often occupies only a small portion in MHD simulations. For this case the eddies get more elongated as it is obvious from the exact expression for $M_{A}<1$ turbulence scaling at small separations [87]:

For larger magnetization, $M_{A}$ gets smaller and the eddies get more and more elongated. The distribution of the direction of eddies also changes and the corresponding changes can be traced by studying the distribution of the gradient directions. For instance, the empirical study of gradients distribution in [91] revealed that the dispersion of velocity gradients measured within the sub-block employed for averaging changes as a power law of $M_{A}$. In this paper we relate this and other measures of velocity and magnetic field gradients with the analytical expectations.

Later in the paper we derive rigorously the properties of gradients according to the statistical theory of MHD turbulence. For readers’ reference, it is advantageous to discuss the gradient properties that directly follow from the MHD turbulence scaling.

For Alfvénic motions of the magnetic and velocity fluctuations are symmetric. The corresponding eddies in magnetized turbulence, as we discussed above, are elongated along the local direction of magnetic field and rotate about this direction. The fact that both velocity and magnetic fluctuations exhibit Kolmogorov-type scaling given by Eq. (5) ensures that the gradients of the amplitudes of fluctuations are largest at the smallest scales resolved. For instance, for velocities $\nabla|v_{l}|\sim v_{l}/l_{\bot}\sim l_{\bot}^{-2/3}$. The fact that both velocity and magnetic field fluctuations take place and are elongated in respect to the local direction of magnetic field means that gradients of individual eddies reflect the direction of the magnetic field at the eddy location. Taken together this means that the Alfvénic gradients are tracing magnetic fields similar to the way that aligned grains trace interstellar magnetic fields.

If one measures gradients within a single magnetized eddy, a distribution of gradient directions is obtained. The maximal amplitude of gradients is in the direction perpendicular to the local magnetic field, while the minimal value corresponding to the direction parallel to the magnetic field. The width of the gradient distribution depends on the media magnetization. Indeed, the eddy elongation increases with the media magnetization, i.e. with the decrease of the Alfvén Mach number $M_{A}$ (see Eq. 10). As a result, the distribution of the gradients measured for an eddy is getting narrower as $M_{A}$ decreases. In addition, the ratio of the maximal to minimal values of the measured gradients increases with $M_{A}$ decrease.

The properties of gradients arising from Alfvénic turbulence play the major role in the Gradient Technique (GT). Through sub-block averaging in [118] one identifies the direction of the maximal gradient amplitude, which is perpendicular to the magnetic field. Measuring the properties of the distribution of the gradients within the sub-block, e.g. the width of the distribution, provides a way to obtain $M_{A}$ [91].

As Alfvén modes impose their scaling on the slow modes, the properties of gradients that correspond to slow modes are similar to those of Alfvénic eddies. On the contrary, the fast modes have different type of scaling and the anisotropy that arises from those modes was shown in [91] to be orthogonal to that of Alfvén and slow modes. In most cases, the fast mode contribution to gradients is subdominant[89, 88] and therefore we focus in this paper on exploring of the gradients arising from Alfvén and slow modes.

The properties of density gradients vary with the sonic Mach number of the media $M_{s}$. At low $M_{s}$ the density can be passively moved by Alfvén velocity fluctuations and therefore the gradients of density will be similar to the velocity gradients. At the same time, at higher $M_{s}$ shocks are being produced. Compressions arising from shocks tend to change the directions of gradients and complicate the interpretation of the density gradients in terms of underlying magnetic field⁵⁵ 5 It was shown in [8] that even in high $M_{s}$ media the low contrast density fluctuations follow the Alfvénic 39 scaling.

The analytical study of the relation between the underlying properties of turbulence and the observables was started with the statistics of velocity available via channel maps (Lazarian & Pogosyan 80, henceforth LP00). The analytical description of the relation of $M_{A}$ and the anisotropy of the observable statistics was obtained in Lazarian & Pogosyan (2012, henceforth LP12) for the case of synchrotron emission fluctuations. Later studies in 66 and 67 provided the relation of the anisotropies in spectroscopic channel maps and velocity centroids with $M_{A}$. These theoretical studies provide the foundations for obtaining the analytical expressions relating $M_{A}$ with the statistics of the gradients of the observable fluctuations. Importantly, they relate both the isotropy degree and quadropole-to-monopole ratio to the Alfvénic Mach number based on the analytical calculations of 80, suggesting that the Alfvénic Mach number can be obtained from observations by studying the statistics of velocity observables.

Let us discuss the mathematical foundations of the gradient methods in application to study of the direction of the magnetic field. We also refer to Appendix A for details.

Observing emission from the turbulent media, one constructs the sky maps of different observables that describe the emission. First of all, this is the intensity of the emission in PPV space, where the position-position $\mathbf{X}$ POS 2D vector and the velocity $v$ coordinate determine position of individual intensities $\rho(\mathbf{X},v)$, and the related integrated quantities, such as the total intensity and velocity centroids which for optically thin lines are $I_{c}(\mathbf{X})\propto\int dvI(\mathbf{X},v)$ and $C(\mathbf{X})\propto\int dvv\rho(\mathbf{X},v)$ respectively (See §IV). The maps represent random fluctuating fields. For general exposition, let us denote this 2D signal map as $\Phi(\mathbf{X})$.

The main local statistical measure of the gradient of a 2D field $\Phi(\mathbf{X})$ is the gradient covariance tensor

where ${\bf R}$ is the separation between the point at the POS plane. Notice that $\sigma_{ij}$ is the zero separation limit of the second derivatives of the field structure function

For a statistically isotropic field, the covariance of the gradients is isotropic, $\sigma_{\nabla_{i}\nabla_{j}}=\frac{1}{2}\delta_{ij}\Delta D(R)|_{R\to 0}$. In the presence of the magnetic field, $D(\mathbf{R})$ of the signal becomes orientation dependent, depending on the angle between $\mathbf{R}$ and the projected direction of the magnetic field. This anisotropy is retained in the limit $\mathbf{R}\to 0$ and results in non-vanishing traceless part of the gradient covariance tensor

As shown in [95] and Appendix A.1, the eigen-directions of this covariance tensor coincide with the direction of local statistical anisotropy of the motions which in MHD turbulence is expected to be given by the local projected direction $\theta_{H}$ of the magnetic field.

where the magnitude of the anisotropic part $I_{2}$ in relation to the isotropic trace $I_{0}=(\nabla_{x}^{2}+\nabla_{y}^{2})D(\mathbf{R})$ is given by the ratio of the quadrupole to monopole moments of the structure function

with $D_{m}(R)=\frac{1}{2\pi}\int d\phi D(R,\phi)e^{-im\phi}$. Thus, the information about magnetic field direction can be recovered by measuring the distributions of the gradients over local patches of 2D sky and finding the eigendirections of $\sigma_{ij}$ in each of them.

The square of $I_{2}/I_{0}$ in Eq. (24) is related to the ratio of the two rotational invariants of the covariance tensor, that we will denote $J_{2}$,

The parameter $J_{2}$ determines the level of dispersion of the directions of the gradients relative to the eigen-directions of $\sigma_{ij}$. This is the quantity that carries information about the level of anisotropy of turbulence processes, and, therefore, on the parameters of the turbulence.

A similar to $J_{2}$ parameter is the ratio of the quadrupole $D_{2}$ to monopole $D_{0}$ moments of the structure function

This ratio has appeared in 84 as a convenient measure of turbulence anisotropy. It was analytically examined for Alfvén, slow, fast modes of MHD turbulence as well for some mixtures of MHD modes for both the synchrotron intensity fluctuations and the measures of fluctuations of intensities within spectroscopic data cube statistics [80, 66, 67]. These studies demonstrated that both $M_{A}$ and the composition of turbulent modes can be studied using the ratio $F$.

The difference with gradients is that the quadrupole to monopole ratio is evaluated at vanishing lag, i.e. $R\rightarrow 0$. We note that for both Alfvén and slow MHD modes $D_{2}(R)$ and $D_{0}(R)$ have the same power-law dependence on the separation $R$ (see Eq. (5) and [22]). Therefore, if these modes dominate, $J_{2}$ dependence on $R$ is very weak.⁶⁶ 6 The exact scaling of fast modes is still the issue of debates. It is suggestive that it follows the acoustic cascade $E(k)\sim k^{3/2}$ for low $M_{s}$ while the spectrum gets steeper, i.e. $\sim k^{-2}$ as $M_{s}$ increases and shocks are formed [22, 71, 96]. This makes it possible to relate $J_{2}$ to factor $F$ at finite lag

This equation relates earlier approach of 84 based on the structure function multipole ratio and the new one based on measuring the gradients.

Strictly local statistical measures, such as $\sigma_{ij}$ at a given sky position, are not accessible in observations, since one needs to average over some region of the sky to simulate the ensemble averaging. In case of a statistically homogeneous measure, the ergodicity principle shows that increasing the volume for averaging allows for ever better determination of statistical values. However, in our case the direction of the magnetic field varies on the sky, therefore statistics of the gradients is not homogeneous. In this situation, determining gradient covariance over the finite patch gives the direction of the projected magnetic field locally averaged over the patch. When choosing the patch size, one needs to balance the loss of resolution and information when increasing the patch size versus the increase of errors due to insufficient statistics if the patch is too small. The choice of the optimal patch size becomes an important parameter of the analysis.

The situation can be modelled by a patch-background split, where within a local patch of scale $L_{p}$ centered at position $\mathbf{X}_{p}$, termed sub-block in [118], the signal can be decomposed into short-wave Fourier modes with $K>K_{p}\approx 1/L_{p}$ that have approximately homogeneous but anisotropic distribution with the power spectrum $P(\mathbf{K},\hat{\Lambda}(\mathbf{X}_{p}))$, while the direction of the magnetic field $\hat{\Lambda}(\mathbf{X}_{p})$ has long-wave variations between the patches. Within a patch the formalism of the homogeneous axis-symmetric turbulence can be applied (Oughton et al. 106, 84) to the gradients, giving the local direction of the magnetic field via the eigen-direction of the gradient covariance tensor $\theta_{H}(\mathbf{X}_{p})$ and the now local parameter $\tilde{J}_{2}(\mathbf{X}_{p})$ expressed via the short-wave part of the power spectrum as

where $P_{2}(K)$ and $P_{0}(K)$ are the quadrupole and the monopole moments of the power spectrum determined within the patch at position $\mathbf{X}_{p}$, and $\hat{\Lambda}_{\mathbf{X}_{p}}$ is the unit vector in the direction of the mean magnetic field within the sub-block. The range of scales $K_{p}<K<K_{b}$ is bounded by the patch size, $K_{p}\approx 1/L_{p}$, and by the beam of the telescope or low-pass post-processing window, $W_{b}$, that smooths out fluctuations at very short modes $K>K_{b}$. Importantly, due to extra $K^{2}$ weight, the gradient covariance tensor is primarily determined by the shortest wavelengths $\sim 1/K_{b}$, which justifies the patch-background split.

One should stress that low pass smoothing of the signal by $W_{b}$ is critical for practical applications of the gradient techniques. It assures that the limit in Eq. (25) is convergent and, correspondingly, that the measurements of the gradients are stable. In addition, the low-pass filtering should be isotropic as not to modify $J_{2}$ by introducing artificial anisotropy. Observations inevitably have low-pass filtering due to finite beam which usually does not resolve turbulent motions on scales less than the turbulence dissipation scale. Numerical simulations, one the other hand, also have intrinsic low-pass filtering corresponding to the numerical dissipation scale $K_{d}$, where the spectrum $P(\mathbf{K},\hat{\Lambda}(\mathbf{X}_{p}))$ is exponentially suppressed for $K>K_{d}$. Despite that, it is advantageous in the analysis to deconvolve, as much as possible, this “instrumental” effects and smooth the data with a synthetic, preferably well behaved Gaussian low-pass filter, to maintain the maximum control over the scale and isotropy of the filtering procedure.

Due to different averaging along the line of sight in different patches, two dimensional projected spectra and, therefore $J_{2}(\mathbf{X}_{p})$ may vary significantly from patch to patch. We discuss this issue in the next section.

Projected 2D gradients are results of the line-of-sight (LOS) integration over depth $L$ of 3D signal

where LOS coordinate is chosen to be $z$. Analysing a patch we shall consider only modes with 2D wavelengths on the sky shorter than $L_{p}$. Correlation and spectral properties of the gradient, nevertheless, include scales extending past $L_{p}$ along the LOS. To study how signal accumulates when we integrate over LOS distance $L\geq L_{p}$, we start by considering a cubical volume of magnetized fluid $L_{p}\times L_{p}\times L_{p}$ which we will term c-block. Within c-block we assume that the turbulence is axis-symmetric with the mean direction of the magnetic field given by the vector ${\bf\lambda}$. Therefore, the analytical description in 84 applies to a c-block magnetic field. Over larger distances $\hat{\lambda}(z)$ may vary according to long-wave behaviour of the magnetic field. Let us see how the statistic of gradients is determined in this picture.

To use the MHD turbulence description in global system of reference (see 84), one should consider scales that are sufficiently small either in comparison with $L_{inj}$ or $l_{A}$. In fact, in most cases we want to get as high resolution as possible and deal with sub-blocks much less than any of the two scales (see Yuen & Lazarian 118). As a result, for our theoretical description, we take the patch size $L_{p}$ to be sufficiently smaller than either $L_{inj}$ for sub-Alfvénic case or $l_{A}$ for super-Alfvénic one. In a c-block of volume $L_{p}^{3}$ one can approximate 3D turbulent signal $\phi$ as anisotropic with a fixed 3D anisotropy direction $\hat{\lambda}$ and described by the power spectrum $P_{c}(k,\mathbf{\hat{k}}\cdot\hat{\lambda})$. Note that $P_{c}$ is defined locally to a particular c-block, thus label “c”.

First, we shall establish how the gradients are correlated along the line of sight. We are interested only in scaling, so we shall neglect the anisotropic nature of the spectrum, taking it in the model form

(See Lazarian & Pogosyan 82) where $m$ is the spectral index, e.g., $m=-2/3$ for Kolmogorov scaling. We compute the trace of the correlation for the gradients smoothed with a Gaussian beam at the beam scale $L_{b}$ in approximation $L_{p}\gg L_{b}$:

where the direct integration gives

This shows that the correlation decreases as $(L_{b}/z)^{-2+m}$ at LOS distances $z>L_{b}$. Thus, the gradients are added up in uncorrelated fashion during integration over the depth exceeding the smoothing scale $L_{b}$.

This is reflected in the variance of the gradients projected over c-block depth $L_{p}$

accumulating only linearly with the c-block depth $L_{p}$, since $L_{p}\gg L_{b}$. Here we point out that $L_{p}$ appears in the gradient variance only as the depth of a c-block along the LOS. From the point of view of orthogonal scales, gradients have sufficient small scale power to be determined predominantly by the beam scale $L_{b}$. Notably, $\langle\phi^{2}\rangle(L_{b}/L_{inj})^{m}$ can be viewed as a variance of differences of quantity $\phi$ at separations $L_{b}$, i.e. $\langle\Delta\phi(L_{b})^{2}\rangle$.

Let us now consider anisotropic aspects of the projected gradient variance tensor in a single c-block. This tensor structure in the Fourier space is

Line of sight integration through the c-block suppresses all the modes with $k_{z}>1/L_{p}$. Thus, as illustrated in Fig. 1, only $K>1/L_{p}$ and $k_{z}<1/L_{p}$ contribute, and one can effectively set $k_{z}$ to zero in the result

(degenerate case when this may not be accurate is when the anisotropy direction is close to the line of sight).

Subsequent c-blocks along z-axis can be considered as having different direction of the magnetic field $\hat{\lambda}(z)$, but uncorrelated gradient signal between the blocks, thus the total variance of projected gradients becoming sum of variances of individual c-blocks, $\sigma_{ij}\approx\sum_{c}\sigma_{ij}(L_{p})$

where overline designates averaging over the directions $\hat{\lambda}$ takes over the integration depth $L$,

i.e the effective 2D spectrum of projected gradients is the average of individual c-block spectra over the directions of the local to c-block magnetic field. Note that $P_{c}$ is defined only within a c-block where we can assume 84 description. The result in Eqs. (35,36) can be also described using polarization analogy. Covariance $\sigma_{ij}$ has the same rotational properties as the linear polarization tensor. Introducing the effective Stokes parameters for a c-block as $Q=\sigma_{xx}-\sigma_{yy}$ and $U=2\sigma_{xy}$, we see that these $Q$ and $U$ are simply additive as we sum c-blocks along the line of sight.

For sub-Alfvénic turbulence, the directions of magnetic field in individual c-blocks along the line of sight are correlated over the turbulence injection scale $L_{inj}$. Within this scale, both the regime of weak at $l>l_{tr}=L_{inj}M_{A}^{2}$ and strong at $l<l_{tr}$ sub-Alfvénic turbulence is present. On scales $l\gg L_{inj}$ the magnetic field fluctuations are not correlated and add in a random walk fashion. The rms angular difference of the magnetic field directions between c-blocks separated by more that $L_{inj}$ is of the order of the Alfvén Mach number $M_{A}$.

For super-Alfvénic turbulence $M_{A}>1$, the correlation length of magnetic field is determined by a smaller scale $l_{A}=LM_{A}^{-3}$. The magnetic field on scales larger than $l_{A}$ is uncorrelated and the variation of the magnetic field direction from one volume $l_{A}\times l_{A}\times l_{A}$ to another similar volume is of the order of unity.

For sub-Alfvénic turbulence in strong regime, the correlated rms deviations of c-block field directions $\lambda$ over the distance $l_{\perp}$ measured perpendicular to the mean magnetic field is characterized by the angle $\alpha$:

which can be seen from Eq. (10). This scaling proceeds up to the transitional scale between the weak and strong cascades $l_{tr}=L_{inj}M_{A}^{2}$ [75].⁷⁷ 7 At the scale $t_{tr}$ the scaling of turbulent motions changes and, as a result, the change of the angle behaves as $$\tan\alpha\approx\frac{v_{l}}{V_{A}}\approx\left(\frac{l_{\bot}}{L_{inj}}\right)^{1/2}M_{A},$$ (38) where we used the scaling of the weak turbulence $v_{l}\approx M_{A}V_{A}(l_{\bot}/L_{inj})^{1/2}$.

To estimate the LOS averaging of a c-block spectrum in Eq. (36) one need a specific model for the power distribution $P_{c}$ in a c-block. Following 87, in the strong turbulence regime at scales $kL_{inj}>M_{A}^{-2}$, the peak of power distribution lies at wave angles to the local direction of the magnetic field given by $\cos\theta_{k}=\hat{k}\cdot\hat{\lambda}\approx v_{k}/V_{A}$ where $v_{k}$ is the rms turbulent velocity at the scale $k=1/l$ (note that $v_{k}/V_{A}$ can be understood as a local Afvén Mach number). Relating $v_{k}$ to the velocity $v_{inj}$ at the injection scale $k\approx L_{inj}^{-1}$ using Eq. (9), we have $\cos\theta_{k}\sim(kL_{inj})^{-1/3}M_{A}^{4/3}$ where it is taken into account that strong turbulence regime is replaced by weak turbulence at $1<kL_{inj}<M_{A}^{-2}$. Thus, the power is strongly peaked for wavevector directions almost perpendicular to the magnetic field, $|\cos\theta_{k}|<M_{A}^{2}$, though strictly perpendicular waves are somewhat disfavoured. The shorter the wavelength, the more concentrated the power is in the perpendicular direction. At the level of individual c-block, the dispersion of directions will be dominated by the longest wavelengths present, i.e,. the ones of the scale of the c-block $k\sim L_{p}^{-1}$. Thus, similar to 84 we can expect that in the regime of strong turbulence,

provides an adequate description for power distribution between differently oriented Alfvén and slow modes. Note, however, that to use this scaling, the size of the c-block should be sufficiently small, $L_{p}<L_{inj}M_{A}^{2}$.

Adding contributions from different c-blocks along the line of sight amounts to convolution of the angular part of the $k_{z}=0$ sky projection c-block power $P_{c}(\mathbf{k})$ with distribution of the 3D direction of the magnetic-field $\widehat{\lambda}$ around the global direction of the magnetic field $\widehat{\lambda}_{0}$. We will model the distribution of $\lambda$ direction by assuming that the perturbations in the magnetic field $\delta\mathbf{b}$ are Gaussian and distributed in 3D isotropically around global $\lambda_{0}$ direction with $\langle\delta b^{2}\rangle/B_{0}^{2}=M_{A}^{2}$, which gives

where $\theta_{\widehat{\lambda}}\in(0,\pi/2)$ is the angle measuring the deviation of the headless vector $\widehat{\lambda}$ from $\widehat{\lambda}_{0}$ axis. Let us note that in the isotropic limit $M_{A}\to\infty$, $\langle\cos^{2}\theta_{\widehat{\lambda}}\rangle\to 1/3$ with respect to an arbitrary axis. In the opposite limit, $M_{A}\to 0$, $\mathcal{P}(\theta_{\widehat{\lambda}})\to\delta_{D}(\theta_{\widehat{\lambda}})$ and $\widehat{\lambda}\to\widehat{\lambda}_{0}$ with $\langle\sin^{2}(\theta_{\widehat{\lambda}})\rangle\sim M_{A}^{2}$.

Combining Eq. (39) and Eq. (40) gives the following expression for the angular part of the line-of-sight averaged power spectrum

For $M_{A}<1$ and $L_{p}/L_{inj}<M_{A}^{2}$ the first exponential term is narrower than the second and can be treated as a Dirac $\delta$-function. The resulting distribution of $\mathbf{K}$ directions is then largely determined by the local magnetic field’s direction variations $\hat{\lambda}$ itself. To approximate the integral over the directions of $\hat{\lambda}$ we notice that the maximum contribution comes from $\hat{\lambda}$’s that are both perpendicular to $\widehat{\mathbf{K}}$ and lie in the plane spanned by $\widehat{\mathbf{K}}$ and $\widehat{\mathbf{\lambda}}_{0}$ vectors. For such $\hat{\lambda}$’s we have $\left(\widehat{\mathbf{\lambda}}\cdot\widehat{\mathbf{\lambda}}_{0}\right)^{2}=1-\cos^{2}\psi_{k}\sin^{2}\gamma$ where $\psi_{k}$ is the angle between $\mathbf{K}$ and 2D projection of the mean field direction $\widehat{\mathbf{\Lambda}}_{0}$, so that $\cos\psi_{k}=\widehat{\mathbf{K}}\cdot\widehat{\mathbf{\Lambda}}_{0}$. Substituting this in Eq. (41) gives

in the approximation of small $\cos\psi_{k}$ that has two solutions $\cos^{2}\psi_{k}\approx(\psi_{k}\pm\pi/2)^{2}$. Our numerical integration suggests that Eq. (42) is an excellent approximation for $M_{A}/\sin\gamma<\pi/8$.

In Fig. 2 we plot the comparison of Eq. (42) with Eq. (41) as a function of $M_{A}/\sin\gamma$. We see that the approximation of Eq. (42) has not much visible differences from the exact form in Eq. (41). In Fig. 3 we show normalized angular power distribution obtained from the synthetic observational maps from pure Alfven modes extracted in numerical simulations “huge-0” to “huge-4” with different $M_{A}$. We can see that Fig. 2 and Fig. 3 has very similar trend as $M_{A}$ grows.

To summarize, we have demonstrated that the covariance matrix of projected gradients in a $L_{p}\times L_{p}$ sub-block on a sky is determined by the spectrum that has a concentration of power at wavemodes orthogonal to the mean projected direction $\widehat{\mathbf{\Lambda}}_{0}$ of the magnetic field within the line-of-sight cube based on the sub-block. Respectively, the eigendirection of the covariance matrix with the largest eigenvalue is orthogonal to $\widehat{\mathbf{\Lambda}}_{0}$ as well. Following this eigendirection change from sub-block to sub-block, one maps the magnetic field over wider region of the sky.

The variance of the directions of gradients at individual pixels is determined by the $J_{2}$ parameter of the covariance matrix. The $J_{2}$ parameter is determined by the properties of the turbulence and the LOS angle of the mean magnetic field. In approximation of Eq. (42) for strong turbulence it is easily calculated to be

while complete form of the functional dependence of $J_{2}$ on $M_{A}$ and $\sin\gamma$ can be determined by numerically integrating the $\sigma_{ij}$ components with the power spectrum given in Eq. (41). Fig. 4 shows our result from the numerical calculations as compared to the theoretical expectation in Eq. (43). We can see that the numerical data (from synthetics, blue points) and the theoretical points (red curve) have very similar trends.

In realistic observations, we will have noise, and the $J_{2}$ parameter is naturally smaller since noise introduces more gradient isotropy to the data set. As a result, we would like to evaluate how noise affects our measurement of $J_{2}$ in observations. Fig. 5 shows the variation of $J_{2}$ of Alfven mode velocity centroid as a function of the noise-to-signal ratio (N/S) for synthetic and realistic numerical data. To mimic observations, we apply an additional 2-pixel Gaussian smoothing kernel to our data. Again, we have to emphasize that synthetic data demonstrates the effect of beam size and dissipation scales in realistic data. We can see from Fig. 5 that the $J_{2}$ parameter indeed decreases dramatically as $N/S$ increases for both synthetic (red points in Fig. 5) and realistic MHD data (blue points in Fig. 5).

In § V we have developed a theoretical prediction for the one-point covariance tensor $\sigma_{ij}$ of the gradient field components. The anisotropy of the gradient map is determined by the $J_{2}$ parameter in an invariant coordinate way. At the same time, the two eigen-directions are one aligned and one perpendicular to the local direction of the LOS-averaged magnetic field. Here we present theoretical predictions for more practical statistics in the point of view of variance of the directions of the gradients, that follows from the results of Section V, assuming that 2D gradients can be described as Gaussian distributed, i.e. in low $M_{A}$. Notice that [95] points out that such assumption of Gaussian distribution on gradient orientation histogram is not correct in the case of low polarization percentage (i.e., $M_{A,LOS}$, the line of sight Alfvén Mach number).

Under Gaussian assumption for gradients, the theoretical expectation of the gradient orientation distribution is given by [95]:

where $\theta_{H}$ is the preferred direction of the distribution that follows the eigen-direction of $\sigma_{ij}$ with the largest eigenvalue. This preferred direction is orthogonal to the magnetic fields for the turbulence regime that we are considering in this current section. The distribution is defined and normalized over an arbitrary interval of angles with the span of $\pi$ radians (i.e., $\theta\in[-\pi/2,\pi/2)$). We see the fundamental role the parameter $J_{2}$ plays in this distribution: In maximally anisotropic case $J_{2}=1$ the angle distribution given by Equation (44) becomes the Dirac $\delta_{D}$-function around $\theta_{B}$ with zero variance of orientation. In the opposite limit where the turbulence is isotropic, $J_{2}=0$, i.e. the angle distribution is flat and the rms angle deviation $\sigma_{\theta}=\pi/(2\sqrt{3})\approx 52^{o}$.

Let us first determine the number $N$ of independent and uncorrelated samples of the gradient directions $\theta_{i}$ in a sub-block. We can estimate $N$ from the following considerations. Since gradients are determined by the short wave side of the spectrum, the correlation length $L_{corr}$ between gradients at neighboring pixels is determined by the largest experimental or synthetic beam scale of the map and the scale of dissipation that damps shortwave perturbations. We can estimate $L_{corr}\approx\mathrm{max}(W_{FWHM},L_{dis})$ where $W_{FWHM}$ is a full width at half maximum of the beam window. The resulting number of independent measurements in the patch is then