ANISOTROPY OF THIRD-ORDER STRUCTURE FUNCTIONS IN MHD TURBULENCE

Magnetohydrodynamic (MHD) turbulence in the presence of a mean field B₀ has a tendency to become two-dimensional (2D). This tendency was early recognized by inspection of the Fourier energy spectra in direct numerical simulations (DNSs). The energy distribution is indeed anisotropic, residing in wavevectors mostly perpendicular to the mean magnetic field (Montgomery & Turner 1981; Shebalin et al. 1983; Grappin 1986). Ideally one would like to quantify the two-dimensionalization as a scaling relation between parallel and perpendicular wavenumbers having the same energy density, ${{k}_{||}}\propto k_{\bot }^{p}$. If p = 1, the anisotropy is scale independent, and the aspect ratio ${{k}_{||}}/{{k}_{\bot }}$ of the isocontour of the Fourier spectrum does not change with scales. If instead $p\lt 1,$ the aspect ratio increases with wavenumber, that is, the spectrum becomes more and more 2D at smaller and smaller scales. In DNSs, the parallel spectral extent is generally very short, due to the limited achievable Reynolds numbers, and the parallel spectrum rarely shows a power law, rendering the distinction between scale-dependent and scale-independent anisotropy difficult. In contrast, the two-point correlation in real space expressed by II-order structure functions (SFs), S, shows in general nicer power-law scaling in the parallel direction, allowing one to quantify the scale-by-scale anisotropy. In analogy with Fourier spectra, the scaling relation involves parallel and perpendicular increments that have the same energy ${{\ell }_{||}}\propto \ell _{\bot }^{p}$, also known as eddy anisotropy.

The II-order SF anisotropy has been widely studied in DNSs of incompressible MHD turbulence. Using a local mean field to identify parallel and perpendicular increments, one finds a scale-dependent anisotropy (Cho & Vishniac 2000): the anisotropy grows at smaller and smaller scales, suggesting a complete 2D state at small enough scales. Furthermore, the anisotropy is controlled by ${{b}_{{\rm rms}}}/{{B}_{0}},$ where b_rms indicates the rms amplitude of turbulent fluctuations, the stronger B₀ the stronger the anisotropy (Müller et al. 2003). The anisotropy is also stronger in regions of stronger magnetic field (Milano et al. 2001). Similarly, the anisotropy in the Fourier spectrum increases for the stronger B₀ (Oughton et al. 1998). However, employing SFs to measure the anisotropy with respect to the global mean field returns a scale-independent anisotropy (Chen et al. 2011), implying that the two-dimensionalization does not increase at smaller scales but reaches an asymptotic value.

A similar dichotomy exists in solar wind measurements, in which one can compute the two-point correlation in time from time series of data collected in situ by spacecraft and then adopt the Taylor hypothesis to obtain spatial increments. The increments are thus taken along the radial direction, but the anisotropy with respect to the magnetic field is recovered thanks to its variable direction with respect to the radial. As in DNSs, the structure function S is found to be more energetic along perpendicular increments than along parallel increments. Again, a local mean-field analysis yields a scale-dependent anisotropy (Horbury et al. 2008; Wicks et al. 2010; Chen et al. 2012), while a global mean-field analysis indicates a scale-independent anisotropy (Tessein et al. 2009). Several authors (Cho et al. 2002; Chen et al. 2011; Beresnyak 2012) showed that for strong turbulence the scale-dependent anisotropy is smoothed out in a global mean-field analysis, even in the presence of a strong mean field. On the other hand, Matthaeus et al. (2012) noted that since the local mean field is a random variable, the local II-order structure functions involve higher-order statistics and cannot be the real-space equivalent of the power spectra. However, for small enough ${{b}_{{\rm rms}}}/{{B}_{0}}$ the global and local measures are expected to coincide.

In this work we investigate the process of two-dimensionalization of MHD turbulence, focusing on the anisotropy measured in the global frame. In this frame, one can obtain a dynamical equation (labeled Von Karman–Howart–Yaglom, Politano–Pouquet (KHYPP) equation) that relates II-order and III-order SFs (Politano & Pouquet 1998), extending to incompressible MHD the Von Karman–Howart–Monin equation for incompressible hydrodynamic turbulence. According to the KHYPP equation, for stationary and homogenous turbulence in the inertial range, the divergence of the III-order structure functions ${\boldsymbol{Y}}$ is proportional to the cascade rate of turbulence $4\epsilon =-\nabla \cdot {\boldsymbol{Y}}$. The divergence is negative, implying that the cascade is achieved by removing positive correlations, and thus increasing the amplitude of S (flattening its power-law index). Thus, by characterizing the anisotropy of ${\boldsymbol{Y}}$ one can get insight into the process of two-dimensionalization of MHD turbulence. Previous studies on the III-order SFs in DNSs of the MHD equations were limited to 2D (Politano et al. 1998; Sorriso-Valvo et al. 2002), thus leaving out the issue of anisotropy. In a recent work, Lamriben et al. (2011) reported for the first time the vector ${\boldsymbol{Y}}$ measured in an experiment of rotating hydrodynamic turbulence. As rotation was increased, the anisotropy of the II-order SFs also increased. They found that the two-dimensionalization can be associated with the tilting of the vector ${\boldsymbol{Y}}$ toward the plane orthogonal to the rotation axis and that the tilting begins at small scales and then propagates to larger and larger scales.

In the present work, we carry out a similar analysis on data from three-dimensional (3D) DNSs of incompressible MHD turbulence by computing for the first time the 3D III-order SFs, in the presence or absence of a mean field. We find that the degree of two-dimensionalization as measured by S is associated with the relative excitation of pseudo- and shear-Alfvén polarizations for stationary turbulence with mean field B₀. We also analyze the full KHYPP equation in developed decaying isotropic turbulence, offering a description of the cascade in real space based on SFs, analogous to the usual Kolmogorov cascade in Fourier space. Note that while the latter is based on the assumption of locality, the KHYPP equation is free from this assumption (requiring only homogeneity), thus representing a more general description of the cascade process in MHD turbulence.

The III-order structure functions ${\boldsymbol{Y}}$ has also been computed in solar wind data to obtain the cascade rate of solar wind turbulence (Sorriso-Valvo et al. 2007; Marino et al. 2008, 2012; Macbride et al. 2008; Smith et al. 2009; Stawarz et al. 2009, 2010). These rates are consistent with the heating rate estimated from proton temperature gradients (Vasquez et al. 2007; Cranmer et al. 2009), suggesting that turbulence may supply the heating required to sustain the non-adiabatic expansion of the solar wind. However, applying the KHYPP equation to the solar wind is a bit problematic since solar wind turbulence is neither stationary nor homogenous (e.g., Hellinger et al. 2013; Gogoberidze et al. 2013; Dong et al. 2014). To obtain the cascade rate, one should compute the divergence of ${\boldsymbol{Y}}$, which is quite difficult with single-spacecraft data, since increments are taken only along one direction at a time (see, however, Osman et al. 2011 for an integral form that exploits the four Cluster spacecraft with the minimal assumption of axisymmetry). The cascade rate can thus be retrieved only assuming a form for the unknown anisotropy of ${\boldsymbol{Y}}$. Although some theoretical predictions exist (e.g., Podesta et al. 2007; Galtier 2009), two methods are commonly employed in observations that assume respectively isotropy or an anisotropic model based on the geometrical slab-plus-2D turbulence that was introduced by Matthaeus et al. (1990) to describe the two-point correlation function of solar wind turbulence. We will exploit the data from our DNSs to test the two methods employed in solar wind data against known anisotropic III-order SFs and estimate the possible systematic errors.

The plan of the paper is as follows. In Section 2 we give a brief introduction to the KHYPP equation, while in Section 3 we describe the method employed to compute SFs of II-order and III-order. The results are presented in Section 4, where we first describe the anisotropy of II-order SFs of the simulations considered. The rest of the section is dedicated to the III-order SFs. We consider first a simulation of decaying turbulence without mean magnetic field, allowing us to test the soundness of our analysis method and to verify that the time-dependent KHYPP equation holds in developed decaying turbulence. Then, we analyze two simulations of turbulence with mean field that have a different strength of anisotropy. Finally, in Section 5, we test on runs with ${{B}_{0}}\ne 0$ the two methods employed to measure the cascade rate in the solar wind. We conclude with a discussion on the results and on the application to the solar wind turbulence.

The KHYPP equation for non-stationary, anisotropic, and incompressible MHD (Politano & Pouquet 1998; Podesta 2008; Carbone et al. 2009) can be obtained from the original MHD equations written in terms of the Elsässer variables ${{{\boldsymbol{z}} }^{\pm }}={\boldsymbol{u}} \mp {\boldsymbol{b}} /\sqrt{4\pi \rho }$, by subtracting the MHD equations evaluated at different positions ${\boldsymbol{x}}$ and ${\boldsymbol{x}} +\ell$ and by averaging in the volume. Under the assumptions of incompressible, homogeneous turbulence one obtains

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} {{\partial }_{t}}\langle |{\Delta }{{{\boldsymbol{z}} }^{\pm }}{{|}^{2}}\rangle +{{\nabla }_{\ell }}\cdot \left\langle {\Delta }{{{\boldsymbol{z}} }^{\mp }}|{\Delta }{{{\boldsymbol{z}} }^{\pm {{|}^{2}}}} \right\rangle = \\ \quad -{\Pi }-{\Lambda }+2\nu \nabla _{\ell }^{2}\langle |{\Delta }{{{\boldsymbol{z}} }^{\pm }}{{|}^{2}}\rangle -4{{\epsilon }^{\pm }}, \\ \end{array}\end{eqnarray} \tag{ 1 }$$

where we have defined the two-point correlation, ${\Delta }{{{\boldsymbol{z}} }^{\pm }}({\boldsymbol{x}} ,\ell )={{{\boldsymbol{z}} }^{\pm }}({\boldsymbol{x}} +\ell )-{{{\boldsymbol{z}} }^{\pm }}({\boldsymbol{x}} )$, and $\langle ...\rangle$ stands for the volume average. These equations describe the evolution of the II-order SF for each Elsässer variable, ${{S}^{\pm }}=\langle |{\Delta }{{{\boldsymbol{z}} }^{\pm }}{{|}^{2}}\rangle$. The divergence term on the left-hand side is the III-order SF, ${{{\boldsymbol{Y}} }^{\pm }}=\left\langle {\Delta }{{{\boldsymbol{z}} }^{\mp }}|{\Delta }{{{\boldsymbol{z}} }^{\pm }}{{|}^{2}} \right\rangle$, involving products of ${\Delta }{{{\boldsymbol{z}} }^{+}}$ and ${\Delta }{{{\boldsymbol{z}} }^{-}}$, which we name Yaglom flux in the following. On the right-hand side (rhs), Π and Λ represent pressure terms and sweeping terms (responsible for the Alfvén effect), respectively, both vanishing for globally homogeneous turbulence (Carbone et al. 2009). The remaining terms represent dissipation. The first one involves the Laplacian with respect to the increments ℓ; it vanishes for vanishing viscosity (for simplicity we assumed equal viscosity and resistivity $\nu =\eta$). The second one, ${{\epsilon }^{\pm }}=-{{\partial }_{t}}{{E}^{\pm }}=\nu \langle {{{\Sigma }}_{j}}({{\partial }_{j}}z_{i}^{\pm })({{\partial }_{j}}z_{i}^{\pm })\rangle$, is the dissipation rate of the Elsässer energies (${{E}^{\pm }}=\langle |{{{\boldsymbol{z}} }^{\pm }}{{|}^{2}}/2\rangle$). In the former, the derivatives of the primitive fields ${{{\boldsymbol{z}} }^{\pm }}$ do not commutate with the averaging operation, and the dissipation rate remains finite for vanishing viscosity.

Summing the contributions of both Elsässer fields, one finally gets an expression for the total energy and cascade:

$$\begin{eqnarray}&&{{\partial }_{t}}S+{{\nabla }_{\ell }}\cdot {\boldsymbol{Y}} =-4\epsilon +2\nu \nabla _{\ell }^{2}\;S,\end{eqnarray} \tag{ 2 }$$

where $S=1/2({{S}^{+}}+{{S}^{-}})$, ${\boldsymbol{Y}} =1/2({{{\boldsymbol{Y}} }^{+}}+{{{\boldsymbol{Y}} }^{-}})$, and $\epsilon =1/2({{\epsilon }^{+}}+{{\epsilon }^{-}})$. Note that because of homogeneity in the above Equations (1) and (2), all the variables depend only on the vector separation, $\ell$.

This equation is valid for decaying turbulence and describes the classical scenario of a turbulent flow in which the dissipation of energy is achieved through a cascade of energy toward smaller scales, where fluctuations are finally damped by viscosity. In forced turbulence one should add on the rhs the forcing terms ($\mathcal{F}$) that inject energy (usually) at large scales.

For stationary turbulence (${{\partial }_{t}}S=0$) forced at large scales ($\mathcal{F}\ne 0$ only at large scales) the injection, cascade, and dissipation all occur at the same rate. At high Reynolds number one expects their respective ranges to be well separated, in analogy to the Kolmogorov cascade in Fourier space, so one has that (i) at large scales the second and last terms in Equation (2) are negligible and $\mathcal{F}=4\epsilon$, and the forcing balances the dissipation rate; (ii) at small scales the second term is negligible, $2\nu \nabla _{\ell }^{2}S=4\epsilon$, and the damping rate is equal to the dissipation rate; and finally (iii) at intermediate scales the last term is negligible, yielding

$$\begin{eqnarray}&&{{\nabla }_{\ell }}\cdot {\boldsymbol{Y}} =-4\epsilon ,\end{eqnarray} \tag{ 3 }$$

that is, the cascade rate, which is equal to the dissipation rate, is given by the divergence of the Yaglom flux. Note that Equation (3) can be used as a definition of the inertial range as being the ensemble of scales for which the equation is approximately satisfied. The definition should hold for quasi-stationary forced turbulence and for developed decaying turbulence. As we will see, for developed decaying turbulence, the time-dependent term is non-negligible only at large scales where ${{\partial }_{t}}S=-4\epsilon$ (the dissipation is balanced by the decay of the II-order SF at large scales).

We can now give a more physical interpretation of the cascade process by rewriting Equation (2) in terms of the autocorrelation function, $C={{C}^{+}}+{{C}^{-}}$, with ${{C}^{\pm }}(\ell )=\langle {{{\boldsymbol{z}} }^{\pm }}({\boldsymbol{x}} +\ell )\cdot {{{\boldsymbol{z}} }^{\pm }}({\boldsymbol{x}} )\rangle$. The autocorrelation functions are related to SFs by

$$\begin{eqnarray}&&{{S}^{\pm }}(\ell )=2{{E}^{\pm }}-2{{C}^{\pm }}(\ell ),\end{eqnarray} \tag{ 4 }$$

and using ${{\partial }_{t}}{{E}^{\pm }}=-{{\epsilon }^{\pm }}$, the KHYPP equation becomes

$$\begin{eqnarray}&&{{\partial }_{t}}C-{{\nabla }_{\ell }}\cdot {\boldsymbol{Y}} =-2\nu \nabla _{\ell }^{2}\;C,\end{eqnarray} \tag{ 5 }$$

showing that ${\boldsymbol{Y}}$ is a flux of negative correlations. A permanent flux of negative correlations toward small scales is equivalent to constantly building new small-scale gradients. For instance, the formation of 2D quasi-perpendicular turbulence will be revealed by a Yaglom flux bringing negative correlation at small perpendicular scales; hence, the vector ${\boldsymbol{Y}}$ must be quasi-uniform, parallel to the ${{\ell }_{\bot }}$ axis, and pointing toward the parallel ${{\ell }_{||}}$ axis.

Coming back to Equation (3), for turbulence in stationary state or decaying turbulence, the cascade is a constant at all inertial-range scales ($-\nabla \cdot {\boldsymbol{Y}} ={\rm const}$). Thus, the inertial range anisotropy cannot appear as different cascade rates in the parallel and perpendicular directions. The anisotropy instead will show up in the shape of the domain of $\ell$ for which $\nabla \cdot {\boldsymbol{Y}} ={\rm const}$.

To illustrate such an anisotropy, one can assume some particular symmetry of the flux ${\boldsymbol{Y}}$ to characterize the cascade, with the additional advantage of obtaining a direct relation to the cascade rate, so as to avoid computing the divergence of ${\boldsymbol{Y}}$. The simplest assumption is that of isotropic turbulence, for which ${\boldsymbol{Y}}$ depends only on the scalar increment ℓ. Rewriting the divergence in spherical coordinates, and assuming stationary conditions and vanishing viscosity, Equation (3) in the inertial range becomes the isotropic KHYPP equation, yielding

$$\begin{eqnarray}&&{{\epsilon }^{{\rm iso}}}=-\frac{3}{4}\frac{{{Y}_{\ell }}(\ell )}{\ell },\end{eqnarray} \tag{ 6 }$$

in which the Yaglom flux ${{Y}_{\ell }}={\boldsymbol{Y}} \cdot \ell /|\ell |$ is projected along the increment. This form is often used in solar wind studies, since although one does not have access to the full divergence in in situ data, the cascade rate can be obtained directly from the projected Yaglom flux. The inertial range occupies a volume that is a sector of a sphere; it can be defined as isotropic since it has the same extent and location on parallel and perpendicular increments.

A strongly anisotropic case is that of 2D turbulence, obtained when the Yaglom flux in the inertial range depends only on the in-plane increments, ${{\ell }_{\bot }}$:

$$\begin{eqnarray}&&\nabla \cdot {\boldsymbol{Y}} ={{\nabla }_{\bot }}\cdot {{{\boldsymbol{Y}} }_{\bot }}=-4\epsilon ,\end{eqnarray} \tag{ 7 }$$

where ${{{\boldsymbol{Y}} }_{\bot }}$ are the in-plane components of ${\boldsymbol{Y}}$ and ${{\nabla }_{\bot }}$ denotes derivatives with respect to the in-plane increments. The Yaglom flux can have out-of-plane components, but the cascade rate is determined only by the in-plane components. Note that a Yaglom flux having only in-plane components in the inertial range is a sufficient condition to have a 2D cascade. Assuming isotropy of the in-plane increments, that is, a dependence only on the scalar separation ${{\ell }_{\bot }}$, one obtains again a direct relation between III-order SFs and the cascade rate:

$$\begin{eqnarray}&&{{\epsilon }^{2D}}=-\frac{{{Y}_{{{\ell }_{\bot }}}}({{\ell }_{\bot }})}{2{{\ell }_{\bot }}},\end{eqnarray} \tag{ 8 }$$

with ${{Y}_{{{\ell }_{\bot }}}}={\boldsymbol{Y}} \cdot {{\ell }_{\bot }}/|{{\ell }_{\bot }}|$ indicating the projection of ${{{\boldsymbol{Y}} }_{\bot }}$ on the radial direction in cylindrical coordinates. Note that the inertial-range domain is not confined to the 2D plane even in this anisotropic turbulence. As we will see, the anisotropy of the inertial-range domain shows up in its different extent and location along parallel and perpendicular increments.

For completeness, we consider finally the case of 1D turbulence, when the III-order SF depends only on one coordinate, say, ${{\ell }_{||}}$. One obtains the cascade rate as

$$\begin{eqnarray}&&{{\epsilon }^{1D}}=-\frac{{{Y}_{||}}({{\ell }_{||}})}{4{{\ell }_{||}}}.\end{eqnarray} \tag{ 9 }$$

Note that the geometrical model of slab turbulence does not correspond to the 1D turbulence, since in the former, fluctuations are assumed to be perpendicular but to depend only on parallel wavevectors (and hence parallel separations). The slab geometrical configuration would indeed have a vanishing divergence.

We consider three high-resolution simulations of MHD turbulence whose parameters are listed in Table 1. Run A is a developed decaying simulation of incompressible MHD turbulence without mean field, representing isotropic turbulence. Run B is a simulation of weakly compressible MHD turbulence,⁷  with mean field ${{B}_{0}}=5$ and anisotropic forcing. The forcing is applied only to components perpendicular to the mean field and to wavevectors mainly perpendicular to the mean field. Thus, we are dealing with strong anisotropic turbulence of fluctuations with shear-Alfvén polarizations, a configuration akin to reduced MHD.⁸ Finally, run C is a simulation of incompressible MHD, again with mean field ${{B}_{0}}=5$, but with a forcing that is isotropic in both components and wavevectors. In the latter the forcing is actually a freezing of the modes $1\leqslant k\leqslant 2$ that maintains an equal amount of pseudo-Alfvén and shear-Alfvén polarizations at large scales, along with equipartition between magnetic and kinetic energy and between Elsässer energies. This simulation can be classified as a case of weak anisotropic turbulence in terms of the strength parameter χ (see Table 1), although it does not have the properties of classical weak turbulence (Ng & Bhattacharjee 1997; Galtier et al. 2000; Meyrand et al. 2014). Indeed, the 3D spectrum has a relatively strong excitation in the parallel direction, resulting in a peculiar anisotropy $E(k,\theta )=A(\theta ){{k}^{p}}$, with an isotropic spectral index p = −2–3/2 in all directions (corresponding to a 1D spectrum with slope $-3/2$) and all the anisotropy appearing as a power anisotropy at large scales $A(\theta )$. We will not discuss the properties and the origin of such a spectrum that can be found in Müller & Grappin (2005), Grappin & Müller (2010), Grappin et al. (2013); what is mostly relevant for the present analysis is that run C has a different 3D anisotropy compared to run B, although in both runs energy resides mainly in perpendicular wavevectors.

Table 1.  Runs and Parameters for Simulations

Run	brms	B0	Rx	χ	${{N}_{x}}\cdot {{N}_{y}}\cdot {{N}_{z}}$	${{10}^{-4}}\nu$	${{10}^{3}}\;\operatorname{Re}$	Forcing
A	0.5	0	1	⋯	10243	0.9	2.6$\to 2.2$	Decaying
B	0.9	5	5	1.8	5123	1.5	1.1	$k_{||}^{f}=1/5$, $k_{\bot }^{f}\leqslant 2$
C	1	5	1	0.2	256·10242	1	2.8	$|{{k}^{f}}|\leqslant 2$, Frozen

Note. ${{b}_{{\rm rms}}}=\sqrt{2{{E}_{b}}}$ is the rms magnetic field fluctuation. B₀ is the mean magnetic field along the x-axis. ${{R}_{x}}={{L}_{x}}/{{L}_{y}}$ is the aspect ratio of the box of size ${{L}_{y}}={{L}_{z}}=2\pi$. The parameter $\chi =t_{A}^{f}/t_{{\rm NL}}^{f}=k_{\bot }^{f}{{b}_{{\rm rms}}}/k_{||}^{f}{{B}_{0}}$ controls the strength of turbulence at forcing scales. ${{N}_{x}},{{N}_{y}},$ and N_z are the number of grid points. ν is the viscosity coefficient (equal to the resistivity η). $\operatorname{Re}={{\left[ 2\pi /(k_{\bot }^{f}{{L}_{{\rm diss}}}) \right]}^{4/3}}$ is the effective Reynolds number, where the dissipation scale is defined as ${{L}_{{\rm diss}}}={{\left( {{\nu }^{3}}/\epsilon \right)}^{1/4}}$. For run A, Re decreases with time in the indicated interval. Forced wavenumbers are normalized by L_y, with ${{k}_{||}}={{k}_{x}},$ ${{k}_{\bot }}={{\left( k_{y}^{2}+k_{z}^{2} \right)}^{1/2}}$. In run B the sound speed is ${{C}_{S}}\approx 12$ and the conductivity coefficient is $\kappa =\nu$.

Download table as:  ASCIITypeset image

We will use three measures to characterize the simulations: II-order structure functions computed in the frame defined by the local mean field (local S), and II- and III-order SFs, respectively S and ${\boldsymbol{Y}}$, computed in an absolute frame attached to the x-axis, which is the direction of the global mean field when it is present. For all structure functions, computation is made calculating increments in real space. Local S are obtained following method I of (Cho & Vishniac 2000), i.e., the local field at scale $\ell$ is defined as ${\boldsymbol{B}} _{{\bf 0}}^{\ell }({\boldsymbol{x}} )=1/2[{\boldsymbol{b}} ({\boldsymbol{x}} +\ell )+{\boldsymbol{b}} ({\boldsymbol{x}} )]$. Note that the measure of anisotropy, i.e., the ratio $S({{\ell }_{\bot }},0)/S(0,{{\ell }_{||}})$ of local S, is not unambiguously defined. In a turbulent medium fluctuations have a wide range of excited scales and the definition of the mean field depends on both scale and position. Thus, higher-order statistics may be introduced in the local S to a different extent, depending on the averaging procedure. However, the two-point average employed in this work is a good working definition, at least in simulations of homogeneous turbulence since it was shown to yield the same results as line average or volume averages, provided that the averaging scale is smaller than the correlation length (Matthaeus et al. 2012). The computation of local S is made for the whole range of available increments $\ell$, but the average is made on a subset of grid points (typically ${{N}_{x}}\times {{N}_{y}}\times {{N}_{z}}={{32}^{3}}$), which was checked to be a sufficient statistic for the anisotropy to converge. On the other hand, the III-order structure functions ${\boldsymbol{Y}}$ are signed quantities, and their computation requires large statistics to converge. Thus, given the number of grid points N in a given direction, we compute ${\boldsymbol{Y}}$ either on a smaller range of increments ($\ell =(1...N/M)dx$ with typically M = 4; dx is the grid size) or on a coarser grid of increments ($\ell =(1...N/M)Mdx$, with typically M = 4), but still on all the grid points, so the averages are made on a sample of $\gtrsim {{10}^{7}}$ data. The dissipation rate () entering the KHYPP equation (2) is obtained directly from the 3D Fourier spectra of magnetic and velocity field ($\hat{b}$ and $\hat{u},$ respectively),

$$\begin{eqnarray}&&\epsilon =-{{\partial }_{t}}E=\nu {{{\Sigma }}_{k}}{{k}^{2}}\hat{u}_{k}^{2}+\eta {{{\Sigma }}_{k}}{{k}^{2}}\hat{b}_{k}^{2}.\end{eqnarray} \tag{ 10 }$$

The terms appearing in Equation (2) are evaluated for four consecutive snapshots separated by approximately a nonlinear time (for the time derivative ${{\partial }_{t}}S$ we use a simple first-order scheme). All quantities entering the KHYPP equation are normalized to the dissipation rate and then averaged over the four snapshots. The 3D data are finally reduced for purposes of representation and analysis by performing an isotropization (averaging over polar and azimuthal angle θ and ϕ in spherical coordinates),

$$\begin{eqnarray}&&\nabla \cdot {\boldsymbol{Y}} {{|}_{{\rm iso}}}=\frac{1}{4\pi }\int \left( \frac{\partial {{Y}_{x}}}{\partial {{\ell }_{x}}}+\frac{\partial {{Y}_{y}}}{\partial {{\ell }_{y}}}+\frac{\partial {{Y}_{z}}}{\partial {{\ell }_{z}}} \right){\rm sin} \theta d\phi d\theta ,\end{eqnarray} \tag{ 11 }$$

or an axisymmetrization (averaging over the azimuthal angle ϕ in cylindrical coordinate with axis along the ℓ_x),

$$\begin{eqnarray}&&\nabla \cdot {\boldsymbol{Y}} {{|}_{{\rm axis}}}=\frac{1}{2\pi }\int \left( \frac{\partial {{Y}_{x}}}{\partial {{\ell }_{x}}}+\frac{\partial {{Y}_{y}}}{\partial {{\ell }_{y}}}+\frac{\partial {{Y}_{z}}}{\partial {{\ell }_{z}}} \right)d\phi .\end{eqnarray} \tag{ 12 }$$

In the following we drop the subscripts iso and axis and eventually mention explicitly the average procedure used for representation.

4.1. Anisotropy of II-order Structure Functions

We measured the anisotropy of II-order SFs in two frames. In the global frame the increments ${{\ell }_{||}}$ and ${{\ell }_{\bot }}$ are taken parallel and perpendicular to a fixed direction x, which is the direction of the mean field ${{{\boldsymbol{B}} }_{{\bf 0}}}$ when it is present. In the local frame the parallel and perpendicular directions are relative to the scale-dependent mean-field direction ${\boldsymbol{B}} _{{\bf 0}}^{\ell }$ (see Section 2). The measure of the anisotropy is obtained by identifying the couples of increments (${{\ell }_{||}},\;{{\ell }_{\bot }}$) at which the parallel and the perpendicular SFs have the same value: ${{S}_{||}}\equiv S({{\ell }_{||}},\;0)={{S}_{\bot }}\equiv S(0,\;{{\ell }_{\bot }})$, i.e., the function ${{\ell }_{||}}({{\ell }_{\bot }})$ measures the aspect ratio of isolevels of the SFs at different scales, also known as eddy shape. Scale-independent anisotropy results in a linear relation ${{\ell }_{||}}\propto {{\ell }_{\bot }}$, that is, an aspect ratio ${{\ell }_{\bot }}/{{\ell }_{||}}$ that does not change with scale. Conversely, scale-dependent anisotropy results in a deviation from the linear scaling ${{\ell }_{||}}\propto \ell _{\bot }^{p}$, with $p\lt 1$: the aspect ratio of SFs, ${{\ell }_{\bot }}/{{\ell }_{||}},$ increases at smaller scales, that is, eddies become more and more elongated in the parallel direction. In Figure 1 we plot ${{\ell }_{||}}({{\ell }_{\bot }})$ for the two measures of anisotropy, local (gray squares) and global (black diamonds), for the three runs listed in Table 1. The dashed line, with slope −1, is a reference for scale-independent anisotropy. The dotted line is a reference for the scale-dependent anisotropy predicted by the critical balance relation (Goldreich & Sridhar 1995):

$$\begin{eqnarray}&&{{\ell }_{||}}\propto \ell _{\bot }^{2/3}.\end{eqnarray} \tag{ 13 }$$

The insets display the same plots compensated by the critical balance anisotropy in order to better appreciate the scaling relation ${{\ell }_{||}}({{\ell }_{\bot }})$.

Zoom In Zoom Out Reset image size

In the local frame, all runs have a scale-dependent anisotropy extending to a wide range of scales. The function ${{\ell }_{||}}({{\ell }_{\bot }})$ has a slope flatter than 1; thus, the anisotropy grows with decreasing scales and eddies are more and more elongated in the parallel direction. The scaling law actually follows the critical balance anisotropy Equation (13) in the range where ${{S}_{\bot }}\propto \ell _{\bot }^{2/3}$, extending for about a decade for runs A, B, and C in the intervals $0.008\lesssim {{\ell }_{\bot }}\lesssim 0.05$, $0.02\lesssim {{\ell }_{\bot }}\lesssim 0.2$, and $0.008\lesssim {{\ell }_{\bot }}\lesssim 0.08,$ respectively. This critical balance anisotropy is quite robust since in all runs the anisotropic relation holds rather well even outside the above-mentioned range of scales where ${{S}_{||}}$ and ${{S}_{\bot }}$ have a clear power-law scaling. Note that run C is a weak turbulence simulation, and it is not obvious that it should have a critical balance anisotropy (see Galtier et al. 2005 for an explanation based on a heuristic model of anisotropic turbulence).

Consider now the anisotropy measured in the global frame (black diamonds in Figure 1), which will be more relevant for the following analysis, since it is related to the III-order SF by the KHYPP equation (2). As expected, run A is perfectly isotropic, and the aspect ratio is unity at all scales. The anisotropy of strong turbulence, run B, becomes scale independent (it has a slope equal to 1) for scales ${{\ell }_{\bot }}\lesssim 0.08$, approaching a constant ratio ${{\ell }_{||}}/{{\ell }_{\bot }}=A\approx 10$ at small scales. For weak turbulence, run C, the anisotropy becomes scale independent only at very small scales (${{\ell }_{\bot }}\lesssim 0.01$), with an aspect ratio $A\approx 5$ that is smaller compared to strong turbulence. The vertical bars in the plot bound the inertial range as identified from Equation (3).⁹ While run B has a scale-independent anisotropy that develops in the inertial range, in run C scale-independent anisotropy is attained only at dissipative scales. Thus, at inertial-range scales the anisotropy is mostly scale dependent in our weak turbulence simulation (run C): it follows the critical balance scaling ${{\ell }_{||}}\propto \ell _{\bot }^{2/3}$ even when measured in the global frame, contrary to expectations (Chen et al. 2011).

Zoom In Zoom Out Reset image size

4.2. III-order Structure Functions and KHYPP Equation

4.2.1. Isotropic Case

We consider first the isotropic case (run A) for which S is isotropic, so we expect also to find an isotropic III-order SF. In Figure 3, panel (a), we plot the Yaglom flux ${\boldsymbol{Y}}$ (III-order SF), averaged along the polar angle of cylindrical coordinates with axis ${{\ell }_{x}}\equiv {{\ell }_{||}}$, and normalized by the scalar increment ℓ. We consider relatively small scales (the largest scale is $\ell =0.5$) to highlight inertial range features, as will be clearer below. The Yaglom flux is almost radial at large scales and becomes remarkably radial at smaller scales $\ell \lesssim 0.08$. The length of the arrows increases toward the origin, indicating that the intensity of the cascade increases when approaching the inertial range (the gray area) while keeping the same (radial) direction. Note also that the arrow length is constant on circles of given ℓ, meaning that ${\boldsymbol{Y}} \approx {\rm Const}\times \ell$ as expected for isotropic turbulence, Equation (6).

Zoom In Zoom Out Reset image size

In panel (b) we plot the isocontours of the divergence of the Yaglom flux, normalized by the dissipation rate. The isocontours are roughly isotropic at large scales and become perfectly isotropic at small scales ($\ell \leqslant 0.03$). In a small interval of scales around $\ell \approx 0.01$, the divergence $-\nabla \cdot {\boldsymbol{Y}}$ has a maximum reaching the value $\approx 0.9$. Thus, the dissipation rate is approximately equal to the cascade rate, and these scales can be identified as the inertial range of turbulence. However, the divergence is not strictly a constant, as expected for the inertial range. Note that although the latter is very short, it is uniformly distributed among scales.

Finally, in panel (c) we plot in logarithmic scales, after isotropization and normalization by $4\epsilon$, all the terms appearing in the KHYPP equation, Equation (2), namely, the divergence of the Yaglom flux $-\nabla \cdot {\boldsymbol{Y}}$ (thick solid line), the dissipative term $2\nu \partial _{\ell }^{2}S$ (red triple-dot-dashed line), and the time-dependent term ${{\partial }_{t}}S$ (dotted line). The dissipative term dominates at small scales, while the time-dependent term (decay) dominates at large scales. The cascade term, $-\nabla \cdot {\boldsymbol{Y}}$, is larger than the other terms for scales $0.003\lesssim \ell \lesssim 0.03$ (the gray area in panel (a)). These two extrema can be identified with the injection scale and the dissipative scale, respectively. It is worth noting that the dissipation scale defined in this way coincides with the estimate based on II-order SFs in Figure 1. From this logarithmic plot it is clear that the inertial range is quite small in this decaying simulation because $-\nabla \cdot {\boldsymbol{Y}}$ is much larger than other terms only in a small interval centered at $\ell \approx 0.01$, where it is roughly horizontal and equal to 0.9 (it should be equal to 1 in an ideally infinite inertial range).

In the same plot we also traced, as a blue long-dashed line, the sum of the three terms just discussed, which should amount at all scales to the dissipation rate $4\epsilon$ (thick solid horizontal line) for good energy conservation. The sum is an almost horizontal straight line, only a factor of 1.2 higher than the dissipation rate for scales $\ell \gtrsim 0.005$. This confirms that the statistics is large enough to ensure convergence and that the conservation of energy holds with sufficient accuracy except at very small scales where a small numerical dissipation probably kicks in: the injection at large scale (including the decay), the cascade in the inertial range, and the dissipation at small scales all occur at the same rate.

To summarize, the analysis of the isotropic turbulence confirms the theoretical expectations: inside the inertial range the Yaglom flux is radially directed and its magnitude scales linearly with the scalar increment ℓ. The divergence of the Yaglom flux is approximately constant in the inertial range, and it is uniformly distributed among scales (isotropic). Note, however, that the extent of the inertial range is very limited due to the relatively small Reynolds number that prevents the formation of a large range of scales where the divergence of the Yaglom flux is the dominant term. With the current resolution (1024³) it is at best half a decade, indicating that this is the minimal resolution required for studies of decaying turbulence (although hyperviscosity would probably alleviate the problem). This is relevant for solar wind studies (Dong et al. 2014), in which expansion induces an additional decay of magnetic and kinetic energy on top of the decay due to the turbulent dissipation.

4.2.2. Anisotropic Case, Strong Turbulence

Let us now turn to the anisotropic run B, which is a simulation of strong turbulence with guide field and forced at large scales on components perpendicular to B₀. In Figure 3, panel (a), we plot again the III-order SF, i.e., the Yaglom flux ${\boldsymbol{Y}}$, averaged over the polar angle in cylindrical coordinates with axis along increments parallel to the mean magnetic field. At variance with the isotropic case, the arrows are colored according to the angle ${{\theta }_{\bot }}$ formed with respect to the perpendicular direction (black stands for ${{\theta }_{\bot }}\leqslant 5{}^\circ$), and their length is normalized to the perpendicular increment ${{\ell }_{\bot }}={{\left( \ell _{y}^{2}+\ell _{z}^{2} \right)}^{1/2}}$. The Yaglom flux is remarkably vertical in the inertial range (the gray area), and it is proportional to the perpendicular increments ${\boldsymbol{Y}} \propto {{\ell }_{\bot }}$ (the arrow length is uniform in the whole inertial range, after normalization). This suggest that turbulence is undergoing a purely 2D cascade, Equation (8).

In panel (b) the normalized divergence of the Yaglom flux is constant over a large interval of parallel and perpendicular increments, with a value close to 1 (the light-green area at level 0.9). However, it is not uniformly distributed among scales; the isocontour of constant divergence extends to smaller perpendicular scales, suggesting that the cascade is not removing positive correlation from the parallel direction. This is consistent with Figure 1(b), which shows a clear anisotropy of S in favor of a two-dimensionalization in the perpendicular plane.

This can be better appreciated in panel (c), where we plot cuts along the parallel and perpendicular directions of the divergence of the Yaglom flux in logarithmic scales. There is a very nice constant divergence in the perpendicular direction, with a value close to the dissipation rate (the horizontal thick solid line at level 1), covering about one decade in the range $0.01\lesssim {{\ell }_{\bot }}\lesssim 0.1$. In the parallel direction, an approximate inertial range is also found, having a smaller extent and shifted to larger scales $0.1\lesssim {{\ell }_{||}}\lesssim 0.6$.

We recall that from $-\nabla \cdot {\boldsymbol{Y}}$ (panels (b) and (c)) one can identify the location of the inertial range and its distribution among parallel and perpendicular scales. On the other hand, the simple dependence of the Yaglom flux allows us to identify the cascade as a 2D cascade with ${\boldsymbol{Y}} =-2\epsilon \;{{\ell }_{\bot }}$ We anticipate that ${{{\boldsymbol{Y}} }^{\mp }}=\left\langle {\Delta }{{{\boldsymbol{z}} }^{\pm }}|{\Delta }{{{\boldsymbol{z}} }^{\mp }}{{|}^{2}} \right\rangle$ have only perpendicular components because there is a dominance of shear-Alfvén polarizations. Indeed, these polarizations have ${\Delta }{{{\boldsymbol{z}} }^{\pm }}$ lying in the perpendicular plane, and as a consequence, the cascade is 2D. This simple picture of 2D cascade does not hold anymore as soon as pseudo polarizations, which have out-of-plane components, are energetically important.

4.2.3. Anisotropic Case, Weak Turbulence

We finally consider the case of weak turbulence with guide field (run C), which is forced isotropically at small wavevectors $1\leqslant |{\boldsymbol{k}} |\leqslant 2$ by imposing at all times (freezing) the corresponding Fourier modes of the fluctuating fields ${\boldsymbol{B}} ,\;{\boldsymbol{u}}$ (note that their $x,y,z$ components have also equal energy). In Figure 4, panel (a), one can see immediately that the Yaglom flux is oblique, with an angle ${{\theta }_{\bot }}$ that changes with scales (we normalize arrow length by the perpendicular increment $|{\boldsymbol{Y}} |/{{\ell }_{\bot }}$ as in Figure 3). This means that the III-order SF has a parallel component and a non-negligible dependence on the parallel increments, thus contributing to the cascade rate through the divergence of the Yaglom flux. Such contribution seems to decrease at small parallel scales, where the Yaglom flux becomes vertical, hinting at a milder two-dimensionalization of this weak turbulence cascade.

Zoom In Zoom Out Reset image size

Isocontours of $-\nabla \cdot {\boldsymbol{Y}}$ are plotted in panel (b), with the usual normalization by $4\epsilon$. The inertial range can be identified with the light-orange area at level 1.1, extending to a wide range of parallel and perpendicular scales. Its distribution is non-uniform and more complex than that of strong turbulence (Figure 3(b)), reflecting the dependence of the III-order SF from ${{\ell }_{\bot }}$ and ${{\ell }_{||}}$.

Panel (c) shows cuts of the divergence of the Yaglom flux along the parallel and perpendicular directions. Although the divergence is not exactly constant, the inertial range covers more than one decade in the perpendicular direction, $0.005\lesssim {{\ell }_{\bot }}\lesssim 0.1$, yielding a cascade rate that is slightly higher than the dissipation rate $4\epsilon$. On the other hand, the cut in the parallel direction is less flat, making more questionable the identification of an inertial range in this direction. The approximate parallel inertial range is shorter and located at larger scales, $0.03\lesssim {{\ell }_{||}}\lesssim 0.3$.

4.2.4. The KHYPP Equation for Pseudo-Alfvén and Shear-Alfvén Polarizations

Decomposing fluctuations in pseudo-Alfvén and shear-Alfvén polarizations (hereafter "pseudo" and "shear", respectively) proves useful to analyze in some more detail the relation between the Yaglom flux and the cascade rate in run C, which is a simulation of incompressible MHD. The decomposition is made in Fourier space, where pseudo-Alfvén polarizations and shear-Alfvén polarizations are oriented along the unitary vectors (e.g., Maron & Goldreich 2001),

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} {{{\boldsymbol{ \xi }} }_{{\rm sh}}}=\frac{{\boldsymbol{k}} \times {{{\boldsymbol{B}} }_{{\bf 0}}}}{{{\left[ 1-{{\left( {\boldsymbol{k}} \cdot {{{\boldsymbol{B}} }_{{\bf 0}}} \right)}^{2}} \right]}^{1/2}}},{{{\boldsymbol{ \xi }} }_{{\rm ps}}}=\frac{{{{\boldsymbol{B}} }_{{\bf 0}}}-({\boldsymbol{k}} \cdot {{{\boldsymbol{B}} }_{{\bf 0}}}){\boldsymbol{k}} }{{{\left[ 1-{{\left( {\boldsymbol{k}} \cdot {{{\boldsymbol{B}} }_{{\bf 0}}} \right)}^{2}} \right]}^{1/2}}}. \\ \end{array}\end{eqnarray} \tag{ 14 }$$

This decomposition is completely equivalent to the decomposition into toroidal and poloidal components of the magnetic fluctuations. In incompressible MHD the same decomposition applies to the velocity field, which is also solenoidal. Shear-Alfvén polarizations are the proper Alfvén modes in full MHD. Their component is perpendicular to both the mean field ${{{\boldsymbol{B}} }_{{\bf 0}}}$ and the wavevector ${\boldsymbol{k}}$, being incompressible. The pseudo-Alfvén polarizations are the incompressible limit of slow modes in MHD. Their component lies in the plane identified by the mean field and the wavevector, and it is again perpendicular to the wavevector. For fluctuations with strong anisotropic spectra (${{k}_{\bot }}\gt \gt {{k}_{||}}$), shear polarizations have wavevectors and components lying in the plane perpendicular to B₀; thus, they represent the 2D modes in the slab-plus-2D decomposition introduced by Matthaeus et al. (1990). Instead, pseudo polarizations have wavevectors in the plane perpendicular to B₀ but components along B₀. This polarization, which is absent in reduced MHD, is instead present in 2.5D configurations with out-of-plane mean field and should not be confused with the slab component that has wavevectors parallel to B₀ and fluctuations perpendicular to it (and is thus included in reduced MHD).

After decomposing fluctuations in pseudo- and shear-Alfvén polarizations, we go back to the real space and compute separately all the contributions to the KHYPP equation, Equation (2). Note that for strictly parallel wavevectors the pseudo-shear decomposition degenerates, so we remove such modes (the slab component) in the following analysis to avoid arbitrary partition of energy into the shear and pseudo polarizations. Using Parseval's theorem $\langle {\Delta }{\boldsymbol{z}} _{{\rm ps}}^{\pm }\cdot {\Delta }{\boldsymbol{z}} _{{\rm sh}}^{\pm }\rangle =0$, the decomposed KHYPP equation can be written as

$$\begin{eqnarray}&&{{\partial }_{t}}({{S}_{{\rm sh}}}+{{S}_{{\rm ps}}})-2\nu \nabla _{\ell }^{2}({{S}_{{\rm sh}}}+{{S}_{{\rm ps}}})+4\epsilon =\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&-{{\nabla }_{\ell }}\cdot \left\langle {\Delta }{{{\boldsymbol{z}} }_{{\rm sh}}}|{\Delta }{{{\boldsymbol{z}} }_{{\rm sh}}}{{|}^{2}}+{\Delta }{{{\boldsymbol{z}} }_{{\rm sh}}}|{\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}}{{|}^{2}} \right\rangle \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&-{{\nabla }_{\ell }}\cdot \left\langle {\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}}|{\Delta }{{{\boldsymbol{z}} }_{{\rm sh}}}{{|}^{2}}+{\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}}|{\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}}{{|}^{2}} \right\rangle \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&-{{\nabla }_{\ell }}\cdot \left\langle 2{\Delta }{{{\boldsymbol{z}} }_{{\rm sh}}}\left( {\Delta }{{{\boldsymbol{z}} }_{{\rm sh}}}\cdot {\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}} \right)+2{\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}}\left( {\Delta }{{{\boldsymbol{z}} }_{{\rm sh}}}\cdot {\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}} \right) \right\rangle ,\end{eqnarray} \tag{ 18 }$$

where we summed the ±species and dropped ±superscripts, i.e., ${{S}_{{\rm sh}}}=1/2\langle |{\Delta }{\boldsymbol{z}} _{{\rm sh}}^{+}{{|}^{2}}+|{\Delta }{\boldsymbol{z}} _{{\rm sh}}^{-}{{|}^{2}}\rangle$, ${\Delta }{{{\boldsymbol{z}} }_{{\rm sh}}}|{\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}}{{|}^{2}}=$ $1/2[{\Delta }{\boldsymbol{z}} _{{\rm sh}}^{+}|{\Delta }{\boldsymbol{z}} _{{\rm ps}}^{-}{{|}^{2}}+$ ${\Delta }{\boldsymbol{z}} _{{\rm sh}}^{-}|{\Delta }{\boldsymbol{z}} _{{\rm ps}}^{+}{{|}^{2}}]$, etc.

The III-order SF on the rhs is split into three lines containing (1) the strain of the shear polarizations on the shear and pseudo energies, (2) the strain of the pseudo polarizations on the shear and pseudo energies, and (3) the mixed terms accounting for the exchange of energy during the cascade between the shear and pseudo polarizations. We anticipate that the mixed terms are negligible in the inertial range; thus, we split the above equation into a system of equations for S_ps and S_sh, with their own cascade rate ${{\epsilon }_{{\rm sh}}}$ and ${{\epsilon }_{{\rm ps}}}$,

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\partial }_{t}}{{S}_{{\rm sh}}}-2\nu \nabla _{\ell }^{2}{{S}_{{\rm sh}}}+4{{\epsilon }_{{\rm sh}}} & = & -{{\nabla }_{\ell }}\cdot \left\langle {\Delta }{{{\boldsymbol{z}} }_{{\rm sh}}}|{\Delta }{{{\boldsymbol{z}} }_{{\rm sh}}}{{|}^{2}} \right.+\left. {\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}}|{\Delta }{{{\boldsymbol{z}} }_{{\rm sh}}}{{|}^{2}} \right\rangle , \\ & & \\ \end{array}\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\partial }_{t}}{{S}_{{\rm ps}}}-2\nu \nabla _{\ell }^{2}{{S}_{{\rm ps}}}+4{{\epsilon }_{{\rm ps}}} & = & -{{\nabla }_{\ell }}\cdot \left\langle {\Delta }{{{\boldsymbol{z}} }_{{\rm sh}}}|{\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}}{{|}^{2}} \right.+\left. {\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}}|{\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}}{{|}^{2}} \right\rangle . \\ & & \\ \end{array}\end{eqnarray} \tag{ 20 }$$

Note that each equation contains only quadratic terms of a given Alfvén polarization (shear or pseudo) even in the III-order SF. On the rhs of Equations (19) and (20), the second terms represent an active contribution to the cascade of pseudo-Alfvén polarizations: if it is vanishing or negligible, the pseudo-Alfvén polarizations can be said to be passive.

This separation will appear justified by the analysis of each term, which is presented in Figure 5 in the same format as Figures 3 and 4, i.e., with the Yaglom flux normalized by the perpendicular increment, ${\boldsymbol{Y}} /{{\ell }_{\bot }}$, and the divergence normalized by the cascade rate, $-\nabla \cdot {\boldsymbol{Y}} /4\epsilon$. We consider separately the two contributions in the rhs terms of the equation for the pseudo energies, Equation (20). In the first column we have the III-order SF accounting for the strain of shear polarizations acting on the pseudo energy, ${{{\boldsymbol{Y}} }^{{\rm SP}}}=\left\langle {\Delta }{{{\boldsymbol{z}} }_{{\rm sh}}}|{\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}}{{|}^{2}} \right\rangle$ (the first term on the rhs). The Yaglom flux (top panel) is perpendicular, as in run B, since the vectorial part of ${{{\boldsymbol{Y}} }^{{\rm SP}}}$ is made of only shear polarizations, which have only perpendicular components. Its magnitude depends mostly on ${{\ell }_{\bot }}$ as can be guessed by the length of the arrows, with a small dependence on ${{\ell }_{||}}$ at large parallel scales. This dependence is better appreciated on the map of $-\nabla \cdot {{{\boldsymbol{Y}} }^{{\rm SP}}}$ (middle panel) that has slightly inclined isocontours instead of horizontal isocontours. The inertial range identified with the green area at level $\approx 0.4$ is non-uniformly distributed in the $({{\ell }_{||}},\;{{\ell }_{\bot }})$ plane, occupying preferentially perpendicular scales ${{\ell }_{\bot }}\lt 0.05$ independently of parallel scale. In the bottom panel one can better locate the inertial range, which is much more extended toward smaller scales in the perpendicular than in the parallel direction.

Zoom In Zoom Out Reset image size

In the second column, we analyze the III-order SF accounting for the strain of pseudo polarizations on the pseudo energy, ${{{\boldsymbol{Y}} }^{{\rm PP}}}=\left\langle {\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}}|{\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}}{{|}^{2}} \right\rangle$, the second term on the rhs of Equation (20). Note that the Yaglom flux is now horizontal (except at small parallel scales ${{\ell }_{||}}\lesssim 0.03$), since for a spectrum (or II-order SF) with energy residing mainly in perpendicular wavevectors, pseudo polarizations have a dominant parallel component. The contribution to the divergence (middle panel) is almost complementary to the previous term, being larger at large parallel and perpendicular scales while it is negligible in the inertial range. This is better seen in the bottom panel, in which it is also shown clearly that this III-order SF does not form an inertial range on its own, but it is responsible for a non-constant energy flux from large to small scales in both parallel and perpendicular directions. One would be tempted to model the two terms, ${{{\boldsymbol{Y}} }^{{\rm SP}}}$ and ${{{\boldsymbol{Y}} }^{{\rm PP}}}$, as 2D+1D components, respectively, according to the orientation of the Yaglom flux. However, such a decomposition does not work since the candidate for the 1D component (${{{\boldsymbol{Y}} }^{{\rm PP}}}$) has a strong dependence on ${{\ell }_{\bot }}$ (second column, top panel), and its divergence does not yield any identifiable inertial range (second column, middle and bottom panels).

Finally, in the third column we analyze the sum of the two terms just discussed, ${{{\boldsymbol{Y}} }^{*{\rm P}}}={{{\boldsymbol{Y}} }^{{\rm PP}}}+{{{\boldsymbol{Y}} }^{{\rm SP}}}=\left\langle ({\Delta }{{{\boldsymbol{z}} }_{{\rm sh}}}+{\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}})|{\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}}{{|}^{2}} \right\rangle$, the whole rhs of Equation (20). The Yaglom flux is now oblique, with dependence on both ${{\ell }_{\bot }}$ and ${{\ell }_{||}}$, the latter being mostly inherited from the strain of pseudo polarizations (${{{\boldsymbol{Y}} }^{{\rm PP}}}$), indicating a smaller degree of two-dimensionalization for turbulence of pseudo-Alfvén polarizations, compared to the strong turbulence of run B. In the middle panel the inertial range extends to parallel and perpendicular scales in a way somewhat similar to run B, although the Yaglom flux is quite different. Comparing the bottom panels in columns one, two, and three, one can see that the strain of both shear and pseudo polarizations contributes to the formation of the wide inertial range in the perpendicular direction $0.005\lesssim {{\ell }_{\bot }}\lesssim 0.1$, as well as to the shorter inertial range found at larger scales in the parallel direction $0.03\lesssim {{\ell }_{\bot }}\lesssim 0.2$ (their ranges coincide with those one obtained without separating shear and pseudo polarizations). More precisely, shear polarizations contribute with a constant cascade at small scales ($\ell \lesssim 0.04$), while pseudo polarizations control the injection of energy into the cascade at large scales ($\ell \approx 0.1$).

We do not repeat the analysis of the strain of shear and pseudo polarizations on shear energy, Equation (19), since it yields qualitatively similar results. The only noticeable difference is that the Yaglom flux ${{{\boldsymbol{Y}} }^{{\rm PS}}}=\left\langle {\Delta }{{{\boldsymbol{z}} }^{{\rm ps}}}|{\Delta }{{{\boldsymbol{z}} }^{{\rm sh}}}{{|}^{2}} \right\rangle$ is radially directed instead of being horizontal. We plot in the fourth column the whole contribution to the cascade of shear energy, ${{{\boldsymbol{Y}} }^{*{\rm S}}}=\left\langle ({\Delta }{{{\boldsymbol{z}} }_{{\rm sh}}}+{\Delta }{{{\boldsymbol{z}} }_{{\rm ps}}})|{\Delta }{{{\boldsymbol{z}} }_{{\rm sh}}}{{|}^{2}} \right\rangle$. Comparing the third and fourth columns, one can see that (i) shear and pseudo energies cascade qualitatively in a similar way (top panels), although the Yaglom flux for shear polarizations, ${{{\boldsymbol{Y}} }^{*{\rm S}}}$, is more perpendicular (more 2D); (ii) the cascade rate is slightly stronger for shear energy than for pseudo energy (middle panels); and (iii) a clear inertial range is found for both shear and pseudo energies (bottom panel) following the decomposition of Equations (19) and (20).

In the fourth column, bottom panel, we also plot the sum of the mixing terms, Equation (18), along the parallel and perpendicular increments (gray lines). They are negligible compared to all other terms in the inertial range, indicating that the shear and pseudo polarizations cascade without exchanging their energy, as already pointed out by Maron & Goldreich (2001) for decaying simulations. However, at large scales $\ell \approx 0.2$, they are of the same order of the term accounting for the strain exerted by pseudo polarizations (II column, bottom panel in Figure 5), suggesting an exchange of energy between the shear and pseudo polarizations.

To summarize, in run C the freezing of large scales (forcing) causes an exchange of energy between pseudo and shear polarization at large scales. The two polarizations cascade in a similar manner (same rate, same inertial range) under the action (strain) of both the shear and pseudo polarizations. The flux toward small scales that is triggered by pseudo polarizations is smaller in magnitude and not constant (it is not a proper inertial range on its own), but it affects considerably the total cascade rate at large inertial-range scales ($\ell \approx 0.1$). On the other hand, the shear polarizations control the constant energy flux at small inertial-range scales ($\ell \lesssim 0.04$). Thus, according to the separation made in Equations (19) and (20), the two cascades are not totally independent. Indeed, whatever the polarization considered (S_ps or S_sh), the cascade is triggered by pseudo polarizations at the interface between injection scales and inertial range scales, while it is sustained by shear polarizations at smaller inertial range scales. The presence and importance of the pseudo polarizations are revealed by the Yaglom flux; the more pseudo polarizations are energetic, the more ${\boldsymbol{Y}}$ becomes oblique.

In measuring the cascade rate of solar wind turbulence through III-order SFs, one needs to make assumptions on the unknown anisotropy of ${\boldsymbol{Y}}$. Despite the limitations due to the imposed symmetries, this brings the advantage of increasing the statistic, on the one hand, and of avoiding the computation of the divergence of ${\boldsymbol{Y}} ,$ on the other hand. Two common assumptions usually employed to reduce in situ data will be tested against our anisotropic DNSs (runs B and C) in order to highlight possible systematic errors on the measure of the cascade rate in the solar wind.

The simplest assumption is that of isotropic turbulence, Equation (6), whose expression is repeated here for convenience,

$$\begin{eqnarray}&&{{\epsilon }^{{\rm iso}}}=-\frac{3}{4}\frac{{{Y}_{\ell }}}{\ell },\end{eqnarray} \tag{ 21 }$$

and involves projecting ${\boldsymbol{Y}}$ along the direction of increments (${{Y}_{\ell }}={\boldsymbol{Y}} \cdot \ell /\ell$). As discussed in Section 2, the cascade rate is directly obtained from the III-order SF without any need to compute derivatives. This "isotropic" method has been mainly employed for fast polar wind (Sorriso-Valvo et al. 2007; Marino et al. 2008, 2012) on Ulysses measurements.

Admitting some form of anisotropy, the minimal assumption is that of axisymmetry. To avoid computing derivatives along the two increments, one again resorts to simplified geometrical models, such as the 2D + 1D turbulence employed by Macbride et al. (2008) and Stawarz et al. (2009, 2010) on Wind and ACE data in the ecliptic solar wind. This method assumes that the III-order SF has parallel and perpendicular components that depend only on the parallel and perpendicular increments, respectively (the cascade is independent in the two directions). The total cascade rate is obtained by combining the two independent equations for (isotropic) 2D-perpendicular and for 1D-parallel cascades (Equations (8) and (9), respectively). Such a "hybrid" cascade rate reads

$$\begin{eqnarray}&&{{\epsilon }^{{\rm hybr}}}={{\epsilon }^{1{\rm D}}}+{{\epsilon }^{2{\rm D}}}={{\epsilon }^{||}}+{{\epsilon }^{\bot }}=-\left( \frac{1}{4}\frac{{{Y}_{||}}}{{{\ell }_{||}}}+\frac{1}{2}\frac{{{Y}_{\bot }}}{{{\ell }_{\bot }}} \right).\end{eqnarray} \tag{ 22 }$$

The isotropic and hybrid cascades are obtained by first computing the 3D III-order SF, and then by applying the corresponding average (isotropy or axisymmetry). Similarly, the "true" cascade rate is obtained from Equation (12):

$$\begin{eqnarray}&&{{\epsilon }^{{\rm true}}}=-\frac{1}{4}\frac{1}{2\pi }\int \left( \frac{\partial {{Y}_{x}}}{\partial {{\ell }_{x}}}+\frac{\partial {{Y}_{y}}}{\partial {{\ell }_{y}}}+\frac{\partial {{Y}_{z}}}{\partial {{\ell }_{z}}} \right)d\phi .\end{eqnarray} \tag{ 23 }$$

We now compare in Figure 6 the true cascade rate (thick solid line) with the isotropic (thin solid line) and the hybrid (thin dashed line) cascade rate for runs B and C (top and bottom panels, respectively). All cascade rates are normalized by the dissipation rate obtained directly from spectra, Equation (10). Each panel represents a cut in the $({{\ell }_{||}},\;{{\ell }_{\bot }})$ plane taken along a fixed direction that forms an angle θ with the mean-field direction.

Zoom In Zoom Out Reset image size

The isotropic cascade rate is strongly angle dependent, independently of the run considered. It returns increasing cascade rates for increasing angles, as can be expected since both runs B and C have a dominant perpendicular cascade. It overestimates the true cascade rate at large angles ($\theta \gtrsim 70{}^\circ$), while it underestimates the true cascade rate by a factor of $\approx 2$ even at $\theta =45{}^\circ$. The hybrid method instead performs very well: it does not vary with angle θ and yields the correct cascade rate at oblique angles $\theta \gtrsim 20{}^\circ$.

In the nearly parallel direction, both the hybrid and isotropic methods strongly underestimate the true cascade rate (by a factor of $\gt 10$), or they completely fail in yielding any cascade rate. Indeed, neither of the two methods is able to account for the dependence of the parallel component of the III-order SF on the perpendicular increments, ${{{\boldsymbol{Y}} }_{||}}({{\ell }_{\bot }})$, which is instead fundamental in our strongly anisotropic runs (e.g., Figure 2(a)). Indeed, the 1D cascade entering the hybrid method is basically negligible at all angles (red dashed line in Figure 6). This implies that the cascade is with a good approximation a 2D cascade in both our anisotropic runs. On the other hand, a linear scaling for the 1D cascade is found in the solar wind (e.g., Macbride et al. 2008), suggesting that the III-order SF of solar wind turbulence is less anisotropic than that one of runs B and C.

We studied the anisotropy of MHD turbulence with or without guide field (B₀) by means of SFs of II and III order, in three DNSs of MHD turbulence at moderate and high resolution (see Table 1). We consider isotropic decaying turbulence (run A), forced strong turbulence with mean field (run B), and forced weak turbulence with mean field (run C).

We used II-order SFs (S) to characterize the anisotropy with respect to the local mean field and the global mean field. The anisotropy with respect to the local mean field is scale dependent in all runs (with or without mean field) and follows the critical balance prediction ${{\ell }_{||}}\propto \ell _{\bot }^{2/3}$. This confirms previous theoretical and numerical findings (Goldreich & Sridhar 1995; Cho & Vishniac 2000; Maron & Goldreich 2001); in the local frame the anisotropy grows at smaller and smaller scales. When instead we used as a reference the global mean field, we found isotropy for MHD turbulence without mean field (run A) and strong two-dimensionalization for turbulence with mean field, with run B being more anisotropic than run C (the small-scale aspect ratio is $\approx 10$ and $\approx 5,$ respectively). Surprisingly, for run C, we found a scale-dependent anisotropy consistent with critical balance even when S is computed in the frame attached to the mean field B₀, at variance with previous analysis in DNSs (Chen et al. 2011) and in solar wind measurements (Tessein et al. 2009), which reported scale-independent anisotropy. A possible explanation is that run C has a strong magnetic field (${{b}_{{\rm rms}}}/{{B}_{0}}=1/5$) and so the local and global frames almost coincide. In favor of this explanation, the anisotropy becomes scale independent at smaller scales. However, there are other "non-standard" aspects of run C that could be related to scale-dependent anisotropy in the global frame: the special 3D spectrum that has an isotropic spectral slope but anisotropic energy levels in parallel and perpendicular direction (Grappin & Müller 2010), the strong excitation of pseudo-Alfvén polarizations (see below), and the isotropic forcing (freezing) that maintains a magnetic excess at large scales (Grappin et al. 2013). These three aspects could be all related to one another and are under investigation.

We then analyzed the vectorial III-order SF (or Yaglom flux, ${\boldsymbol{Y}}$), which is related to S computed in the global frame by a generalization to incompressible MHD (Politano & Pouquet 1998) of the Von Karman–Howart–Yaglom relation in hydrodynamics (KHYPP equation in the following). The III-order SF, ${\boldsymbol{Y}}$, is expected to be anisotropic for MHD turbulence with mean field, but it has never been reported in DNSs or experiments (except for rotating hydrodynamic turbulence experiments; Lamriben et al. 2011). Despite that some anisotropic models exist (Podesta et al. 2007; Galtier 2009, 2012), they have never been verified in DNSs.

We tested the KHYPP equation in isotropic decaying turbulence (run A) and found that the Yaglom flux is almost purely radial, its divergence is isotropic, and the maximum of $-\nabla \cdot {\boldsymbol{Y}}$ is of the order of the dissipation rate (only a factor of 0.8 smaller), confirming theoretical expectations (Politano & Pouquet 1998). However, the condition defining the inertial range, $-\nabla \cdot {\boldsymbol{Y}} ={\rm const}$, is only fairly satisfied in our decaying turbulence that has an inertial range of only half a decade despite the high resolution (1024³). More interestingly, we showed that the familiar concept of Kolmogorov cascade in Fourier space also applies in the real space with SFs. Note that no assumption on locality of nonlinear interactions is needed to obtain the KHYPP equation. Despite this, we found that injection, cascade, and dissipation occur at the same rate and that the cascade is achieved via a fairly constant nonlinear transfer in the inertial range. In particular, by comparing the divergence of the Yaglom flux with the other terms appearing in the KHYPP equation (dissipation term and decaying term), one can determine the dissipation scale and the injection scale of turbulence. This has important applications in solar wind turbulence, since, due to expansion, it is an intrinsically decaying turbulence for which the injection scale is not easily identified (Hellinger et al. 2013; Dong et al. 2014).

For strong turbulence with mean field, we found that ${\boldsymbol{Y}}$ has only perpendicular components and depends only on perpendicular increments; this is a 2D turbulence as defined in Equation (7). The corresponding inertial range extends to small perpendicular scales, while it is shorter and located at larger scales in the parallel direction: the anisotropy of ${\boldsymbol{Y}}$ emerges as a non-uniform distribution of the scales having $-\nabla \cdot {\boldsymbol{Y}} ={\rm const}$. A similar non-uniform distribution is found in the forced, weak turbulence case (run C), with an important difference concerning the orientation of ${\boldsymbol{Y}}$. The latter is no more strictly perpendicular, but oblique. In particular, the Yaglom flux is almost radial at large parallel scales in the inertial range and becomes more and more perpendicular at small parallel scales. We interpret this behavior as a weaker two-dimensionalization of turbulence, since the oblique orientation of ${\boldsymbol{Y}}$ implies a dependence on parallel and perpendicular scales, in contrast to the purely 2D case in which ${\boldsymbol{Y}} ={{{\boldsymbol{Y}} }_{\bot }}({{\ell }_{\bot }})$.

By decomposing fluctuations into shear-Alfvén and pseudo-Alfvén polarizations in run C, we also found that energetically important pseudo polarizations are responsible for the non-perpendicular III-order SF. The pseudo-shear decomposition allows us to prove that the two polarizations cascade independently, that is, without exchanging their energy, in the inertial range. One can thus write separate KHYPP equations for pseudo and shear energies. However, strictly speaking, the two cascades are not independent, since in each of them the strain of pseudo and shear polarizations enters the expression for the Yaglom flux. In particular, the strain exerted by pseudo polarizations is associated with a non-constant cascade of energy at the large inertial-range scales; thus, pseudo polarizations control the injection of energy in run C (although becoming passive at smaller inertial-range scales; e.g., Maron & Goldreich 2001; Cho & Lazarian 2003). In Grappin et al. (2013), it is shown how the equipartition of kinetic and magnetic energy at the large frozen scales is associated with the weaker anisotropy of the 3D spectrum of run C, having the property of a 3D Iroshnikov–Kraichnan turbulence. Here we argue that kinetic/magnetic equipartition at the largest scales induces an exchange of energy between pseudo and shear polarizations, thus allowing pseudo polarizations to control the cascade at the interface between injection and inertial-range scales.

Finally, on the anisotropic runs B and C we tested two methods that are commonly applied to obtain the cascade rate in the solar wind. The methods both rely on assumptions of the anisotropy of the III-order SF, allowing one to obtain the cascade rate directly from ${\boldsymbol{Y}}$ (i.e., without computing the divergence). The first method assumes that ${\boldsymbol{Y}}$ is isotropic. The second method (hybrid method) assumes axisymmetry and models the anisotropy as a geometric superposition of a 1D (parallel) cascade and an isotropic 2D (perpendicular) cascade. We found that the isotropic method is strongly angle dependent, yielding correct cascade rates only at large angles between the direction of increments and the magnetic field, translating to angles between the magnetic field direction and the solar wind direction ${{\theta }_{{\rm BV}}}\gtrsim 70{}^\circ$ (needless to say, the isotropic method works very well for run A). For fast polar wind, emanating from stable coronal holes, turbulence is expected to have a strong anisotropy (Verdini et al. 2012; Perez & Chandran 2013), comparable to our runs B and C. We thus suggest that the angle dependence of the isotropic method is at the origin of the relatively small number of intervals (13%) in which the linear scaling ${{Y}_{\ell }}\propto \ell$ is found in polar solar wind in situ data (Marino et al. 2012). In contrast, the hybrid method performs very well on our anisotropic runs, yielding an angle-independent cascade rate for ${{\theta }_{{\rm BV}}}\gtrsim 20{}^\circ$. Both methods fails to obtain any cascade rate when increments are parallel to the mean field. This is due to the dependence ${{{\boldsymbol{Y}} }_{||}}({{\ell }_{\bot }})$ that is excluded in both assumptions but present in all our runs. This, together with the good performance of the hybrid method, indicates that the anisotropic runs B and C are fairly well described by a 2D cascade model, despite their different degree of anisotropy. Note that in the solar wind a 1D cascade, actually a linear scaling ${{Y}_{||}}\propto {{\ell }_{||}}$, is indeed measured with the hybrid method, indicating that solar wind turbulence in the ecliptic has a different and probably much weaker anisotropy compared to our DNSs.

We conclude by noting that the ratio ${{b}_{{\rm rms}}}/{{B}_{0}}\approx 1/5$ employed in our simulation is lower than the actual value found in the solar wind, which is $\approx 1/2$ at a scale of a few hours. The small value of fluctuation amplitude was chosen to highlight the effects of anisotropy and makes the comparison of our results with solar wind turbulence at large scales a bit weak, being more meaningful for scales shorter than a minute. In a future work we plan to apply a similar analysis to simulations with larger ratio ${{b}_{{\rm rms}}}/{{B}_{0}},$ also dropping the assumption of axisymmetry. In fact, Cluster observations (Narita et al. 2010) showed that the solar wind turbulence is not axis-symmetric (but see Turner et al. 2011 for an explanation in terms of an observational bias). DNSs of MHD turbulence in the framework of the Expanding Box Model (Grappin et al. 1993; Grappin & Velli 1996) also show that expansion causes non-axis-symmetric 3D spectra (Dong et al. 2014). Note that expansion is also responsible for a differential decay of fluctuations; thus, cross helicity and magnetic/kinetic imbalance of fluctuations also enter the KHYPP equation in the expanding solar wind (Gogoberidze et al. 2013; Hellinger et al. 2013). Another interesting direction is to explore the role of velocity shears on the KHYPP equation (Wan et al. 2009, 2010), which have been invoked to be the main driver for the turbulence cascade in both the ecliptic (Roberts et al. 1992) and the polar solar wind (Marino et al. 2012). Again, numerical simulations with the Expanding Box Model seem to be a very promising tool to understand the combined effect of shear and expansion (e.g., Roberts & Ghosh 1999) in shaping the cascade of solar wind turbulence.

This project has received funding from the European Union's Seventh Framework Programme for research, technological development, and demonstration under grant agreement no. 284515 (SHOCK). Website: project-shock.eu/home/. A.V. acknowledges the Interuniversity Attraction Poles Programme initiated by the Belgian Science Policy Office (IAP P7/08 CHARM). P.H. acknowledges grant P209/12/2023 of the Grant Agency of the Czech Republic. Computational resources were provided by CINECA (grant 2014 HP10CLF0ZB and HP10CNMQX2M).