2.1. K41

Let us recall with maximum brevity what this Kolmogorov state is. Assume that energy is being pumped into the system at large scales and at some fixed rate $\varepsilon$. Then, in the inertial range (i.e., at small enough scales so the system is locally homogeneous but not small enough for viscosity or any other microphysics to matter yet), this same $\varepsilon$ is the constant energy flux from scale to scale. Assuming that the cascade (i.e., the passing of energy from scale to scale) is local, the energy spectrum is, by dimensional analysis,

(2.1)\begin{equation} E(k)\sim \varepsilon^{2/3}k^{{-}5/3}, \end{equation}

the famous Kolmogorov spectrum (Kolmogorov Reference Kolmogorov1941b; henceforth K41), or, in terms of typical velocity increments between points separated by a distance $\lambda$:

(2.2)\begin{equation} \delta u_\lambda \sim (\varepsilon\lambda)^{1/3}. \end{equation}

This is all obvious because the dimensions of the quantities involved are

(2.3)\begin{equation} [\varepsilon] = \frac{U^3}{L},\quad \left[\int\mathrm{d} k E(k)\right] = [\delta u_\lambda^2] = U^2,\quad [k] = [\lambda^{{-}1}] = L^{{-}1}, \end{equation}

where $U$ is a unit of velocity and $L$ of length. As we will be dealing with an incompressible medium (which is always achievable by going to small enough scales and so to sufficiently subsonic motions), its density is an irrelevant constant.Footnote ³

2.2. IK

It was Kraichnan (Reference Kraichnan1965) who appears to have been the first to realise clearly the point made above about the irreducibility of the magnetic field. He therefore argued that, if the background uniform magnetic field $\boldsymbol {B}_0$, which in velocity units is called the Alfvén speed,

(2.4)\begin{equation} v_{\mathrm{A}} = \frac{B_0}{\sqrt{4\pi\rho_0}} \end{equation}

($\rho _0$ is the mass density of the conducting medium), was to have a persistent (at small scales) role in the energy transfer from scale to scale, then the energy spectrum in the inertial range must be, again by dimensional analysis,

(2.5)\begin{equation} E(k) \sim (\varepsilon v_{\mathrm{A}})^{1/2} k^{{-}3/2}\quad\Leftrightarrow\quad \delta u_\lambda \sim (\varepsilon v_{\mathrm{A}}\lambda)^{1/4}. \end{equation}

This is known as the Iroshnikov–Kraichnan spectrum (henceforth IK; figure 1).Footnote ⁴ The scaling exponent was fixed by the requirement, put forward with the trademark combination of deep insight and slightly murky argumentation that one often finds in Kraichnan's papers, that the Alfvén time $\tau _{\mathrm {A}}\sim 1/kv_{\mathrm {A}}$ was the typical time during which interactions would occur (before build-up of correlations was arrested by perturbations propagating away from each other), so the energy flux had to be proportional to $\tau _{\mathrm {A}}$ and, therefore, to $1/v_{\mathrm {A}}$ – thus requiring them to enter in the combination $\varepsilon v _{\mathrm {A}}$.Footnote ⁵

Figure 1. IK and GS (photo of R. S. Iroshnikov courtesy of N. Lipunova and K. Bychkov, Sternberg Astronomical Institute; photo of R. H. Kraichnan courtesy of the AIP Emilio Segrè Visual Archives).

Kraichnan's prediction was viewed as self-evidently correct for 30 years, then wrong for 10 years (§ 5), then correct again (in a different sense) for another 10 years (§ 6), then had to be revised again, at small enough scales (§ 7). His own interpretation of it (which was also Iroshnikov's, arrived at independently) was certainly wrong, as it was based on the assumption – natural for a true Kolmogorovian susceptible to the great man's universalist notion of ‘restoration of symmetries’ at small scales, but, in retrospect, illogical in the context of proclaiming the unwaning importance of $\boldsymbol {B}_0$ at those same small scales – that turbulence sufficiently deep in the inertial range would be isotropic, i.e., that there is only one $k$ to be used in the dimensional analysis. In fact, one both can and should argue that, a priori, there is a $k_{\parallel }$ and a $k_{\perp }$, which represent the variation of the turbulent fields along and across $\boldsymbol {B}_0$ and need not be the same. The presence of the dimensionless ratio $k_{\parallel }/k_{\perp }$ undermines the dimensional inevitability of (2.5) and opens up space for much theorising, inspired or otherwise.

2.3. GS95

Intuitively, in a strong magnetic field, perturbations with $k_{\parallel } \ll k_{\perp }$ should be more natural than isotropic ones, as the field is frozen into the motions but hard to bend. It turns out that MHD turbulence is indeed anisotropic in this way, at all scales, however small (see figure 2). This was realised quite early on, when the first, very tentative, experimental and numerical evidence started to be looked at (Robinson & Rusbridge Reference Robinson and Rusbridge1971; Montgomery & Turner Reference Montgomery and Turner1981; Shebalin, Matthaeus & Montgomery Reference Shebalin, Matthaeus and Montgomery1983), but, interestingly, it took more than a decade after that for the IK theory to be properly revised.

Figure 2. A visualisation of numerical RMHD turbulence, courtesy of A. Beresnyak (run R5 from Beresnyak Reference Beresnyak2012a, $1536^3$). The shades of grey represent the absolute value of $\boldsymbol {Z}_{\perp }^+ = \boldsymbol {u}_{\perp } + \boldsymbol {b}_{\perp }$ (see § 3).

Dynamically, the parallel variation (on scale $l_{\parallel }\sim k _{\parallel }^{-1}$) is associated with the propagation of Alfvén (Reference Alfvén1942) waves, the wave period (or ‘propagation time’) being

(2.6)\begin{equation} \tau_{\mathrm{A}} \sim \frac{l_{{\parallel}}}{v_{\mathrm{A}}}, \end{equation}

and the perpendicular variation (on scale $\lambda \sim k _{\perp }^{-1}$) with nonlinear interactions, whose characteristic time is naïvely equal to the ‘eddy-turnover time’:

(2.7)\begin{equation} \tau_{\mathrm{nl}} \sim \frac{\lambda}{\delta u_\lambda} \end{equation}

(we shall see in § 6 why this is naïve). Here and below, $\delta u_\lambda$ is used to represent the turbulent field on the grounds that, in Alfvénic perturbations, $\delta u_\lambda \sim \delta b_\lambda$, where $\delta b$ is the magnetic perturbation in velocity units (see § 3 for a discussion with equations). Declaring the two times comparable at all scales was an inspired conjecture by Goldreich & Sridhar (Reference Goldreich and Sridhar1995, Reference Goldreich and Sridhar1997) (henceforth GS95; figure 1),Footnote ⁶ which has come to be known as the critical balance (CB). I shall discuss the physical reasons for it properly in §§ 4 and 5, but here let me just postulate it. Then, naturally, the ‘cascade time’ (i.e., the typical time to transfer energy from one perpendicular scale $\lambda$ to the next) must be of the same order as either of the two other times:

(2.8)\begin{equation} \tau_{\mathrm{c}} \sim \tau_{\mathrm{A}}\sim \tau_{\mathrm{nl}}. \end{equation}

If (2.7) is used for $\tau _{\mathrm {nl}}$, then (2.8) obviates the magnetic field and returns us to the K41 scaling (2.2), viz.,

(2.9)\begin{equation} \frac{\delta u_\lambda^2}{\tau_{\mathrm{c}}} \sim \varepsilon,\quad \tau_{\mathrm{c}}\sim\tau_{\mathrm{nl}}\sim\frac{\lambda}{\delta u_\lambda} \quad\Rightarrow\quad \delta u_\lambda\sim(\varepsilon\lambda)^{1/3} \quad\Leftrightarrow\quad E(k_{{\perp}}) \sim \varepsilon^{2/3} k_{{\perp}}^{{-}5/3}. \end{equation}

This anisotropic version of K41 is known as the Goldreich–Sridhar (or GS95) spectrum. Simultaneously, along the field,Footnote ⁷ the velocity increments must satisfy

(2.10)\begin{equation} \frac{\delta u_{l_{{\parallel}}}^2}{\tau_{\mathrm{c}}} \sim \varepsilon,\quad \tau_{\mathrm{c}}\sim\tau_{\mathrm{A}}\sim\frac{l_{{\parallel}}}{v_{\mathrm{A}}} \quad\Rightarrow\quad \delta u_{l_{{\parallel}}}\sim\left(\frac{\varepsilon l_{{\parallel}}}{v_{\mathrm{A}}}\right)^{1/2}. \end{equation}

Thus, $\boldsymbol {B}_0$'s influence does persist, but its size enters only the parallel scaling relations, not the perpendicular ones. Formally speaking, (2.9) is just the K41 dimensional argument for the perpendicular scale $\lambda$, with the CB conjecture used to justify not including $v_{\mathrm {A}}$ and $l_{\parallel }$ amongst the local governing parameters. The assumption is that the sole role of $\boldsymbol {B}_0$ is to set the value of $l_{\parallel }$ for any given $\lambda$: comparing (2.9) and (2.10), we get

(2.11)\begin{equation} l_{{\parallel}} \sim v_{\mathrm{A}}\varepsilon^{{-}1/3}\lambda^{2/3}. \end{equation}

Physically, this $l_{\parallel }$ is the distance that an Alfvénic pulse travels along the field, at speed $v_{\mathrm {A}}$, over the time $\tau _{\mathrm {nl}}$, given by (2.7), that it takes a turbulent perturbation of size $\lambda$ to break up nonlinearly. It is natural to argue, by causality, that this is the maximum distance over which any perturbation can remain correlated (Boldyrev Reference Boldyrev2005; Nazarenko & Schekochihin Reference Nazarenko and Schekochihin2011).

This narrative arc brings us approximately to the state of affairs in the mid-1990s, although the GS95 theory did not really become mainstream until the early years of this century – and soon had to be revised. Before I move on to discussing this revision (§ 6) and the modern state of the subject, I would like to put the discussion of what happens dynamically and how CB is achieved on a slightly less hand-waving basis than I have done so far. Indeed, why critical balance? Pace the causality argument, which sets the maximum $l_{\parallel }$, why can $l_{\parallel }$ not be shorter? Is the nonlinear-time estimate (2.7), crucial for the scaling (2.9), justified? What happens dynamically?

From this point on, my exposition will be more sequential, I will avoid jumping ahead to the highlights and adopt a more systematic style, rederiving carefully some of the results reviewed in this section (an already well educated – or impatient – reader is welcome to skip or skim forward at her own pace).

4.1. WT is irrelevant

On a broad-brush qualitative level, one can deal with this possibility as follows. Assume that in the energy-injection range, represented by some perpendicular scale $L_{\perp }$ and some parallel scale $L_{\parallel }=2\pi /k_{\parallel }$, Alfvén waves are generated with amplitudes $Z^\pm$ so small that

(4.1)\begin{equation} \omega_{\boldsymbol{k}}^\pm = {\pm} k_{{\parallel}} v_{\mathrm{A}} = \frac{1}{\tau_{\mathrm{A}}} \gg \frac{1}{\tau_{\mathrm{nl}}^\pm} \sim \frac{Z^\mp}{L_{{\perp}}}. \end{equation}

If they are viewed as interacting quasiparticles (‘$+$’ can only interact with ‘$-$’, and vice versa), the momentum and energy conservation in a three-wave interaction require

(4.2)\begin{equation} \left.\begin{array}{ll@{}} & {\boldsymbol{p}} + {\boldsymbol{q}} = {\boldsymbol{k}},\\ & \omega_{\boldsymbol{p}}^\mp + \omega_{\boldsymbol{q}}^\pm = \omega_{\boldsymbol{k}}^\pm\quad\Rightarrow\quad -p_{{\parallel}} + q_{{\parallel}} = k_{{\parallel}} \end{array}\right\}\quad\Rightarrow\quad q_{{\parallel}} = k_{{\parallel}},\quad p_{{\parallel}} = 0. \end{equation}

Thus, three-wave interaction in fact involves a wave (${\boldsymbol {q}}$) scattering off a 2D perturbation ($p_{\parallel }=0$, not a wave) and becoming a wave (${\boldsymbol {k}}$) with the same frequency (because $k_{\parallel }=q_{\parallel }$) and a different perpendicular wavenumber (${\boldsymbol {k}}_{\perp } = {\boldsymbol {p}}_{\perp } + {\boldsymbol {q}}_{\perp }$). Intuitively, there will be a cascade of the waves to higher $k_{\perp }$. If the amplitude of the waves does not fall off with $k_{\perp }$ faster than $k_{\perp }^{-1}$, which is equivalent to their energy spectrum being less steep than $k_{\perp }^{-3}$, then the nonlinear-interaction time will become ever shorter with larger $k_{\perp }$, even as the waves’ $k_{\parallel }$ and, therefore, their frequency stay the same. Eventually, at some perpendicular scale, which I shall call $\lambda _{\mathrm {CB}}$, the condition $\tau _{\mathrm {nl}}\gg \tau _{\mathrm {A}}$ will be broken, so we end up with $\tau _{\mathrm {nl}} \sim \tau _{\mathrm {A}}$ and can return to considerations of the strong-turbulence regime, critical balance, etc. Numerically, this transition was first captured quite recently, by Meyrand et al. (Reference Meyrand, Galtier and Kiyani2016), whose result is shown in figure 3.

Figure 3. (Decaying) MHD simulation of transition from weak to strong turbulence by Meyrand, Galtier & Kiyani (Reference Meyrand, Galtier and Kiyani2016): the upper panel shows the magnetic spectrum vs. $k_{\parallel }$ and $k_{\perp }$ (where $k_{\parallel }$ is along the global mean field), the lower one the same integrated over $k_{\parallel }$ and normalised by $k_{\perp }^{3/2}$ (see § 6 for why $k_{\perp }^{-3/2}$ rather than $k_{\perp }^{-5/3}$). A transition manifestly occurs from a $k_{\perp }^{-2}$ to a $k_{\perp }^{-3/2}$ spectrum and, simultaneously, from a state with no $k_{\parallel }$ cascade (and a relatively narrow-band parallel spectrum) to one consistent with a CB cascade (2D spectra of CB turbulence are worked out in Appendix C). (Reprinted with permission from Meyrand et al. Reference Meyrand, Galtier and Kiyani2016, copyright (2016) by the American Physical Society.)

The transition scale $\lambda _{\mathrm {CB}}$ is easy to estimate without the need for a specific WT theory. In view of (4.2), weak interactions cannot increase the characteristic parallel scale of the perturbations, which therefore remains $L_{\parallel }$. Then $\lambda _{\mathrm {CB}}$ is the perpendicular scale corresponding to $l_{\parallel }=L_{\parallel }$ in (2.11), viz.,

(4.3)\begin{equation} \lambda_{\mathrm{CB}} \sim \varepsilon^{1/2}\left(\frac{L_{{\parallel}}}{v_{\mathrm{A}}}\right)^{3/2}. \end{equation}

In fact, one does not even need to invoke the GS95 CB curve (2.11), because (4.3) is the only dimensionally correct possibility if one asks for a scale that depends on $\varepsilon$ and $\tau _{\mathrm {A}}\sim L _{\parallel }/v_{\mathrm {A}}$ only; that $L_{\parallel }$ and $v_{\mathrm {A}}$ must enter in this combination follows from the fact that $\nabla _{\parallel }$ and $v_{\mathrm {A}}$ only enter multiplying each other in the RMHD equations (3.1).

A reader who is both convinced by this argument and regards it as grounds for dismissing the WT regime as asymptotically irrelevant, can at this point skip to § 5. The rest of this section is for those restless souls who insist on worrying about what happens in weakly forced systems at $\lambda \gg \lambda _{\mathrm {CB}}$.

4.2. A sketch of WT theory

A very simple heuristic WT calculation (Ng & Bhattacharjee Reference Ng and Bhattacharjee1997; Goldreich & Sridhar Reference Goldreich and Sridhar1997) – a useful and physically transparent shortcut, and a good starting point for discussion – goes as follows.

Imagine two counterpropagating Alfvénic structures of perpendicular size $\lambda$ and parallel coherence length $L_{\parallel }$ (which cannot change in WT, as per the argument in § 4.1) passing through each other and interacting weakly. Their transit time through each other is $\tau _{\mathrm {A}}\sim L _{\parallel }/v_{\mathrm {A}}$ and the change in their amplitudes during this time is

(4.4)\begin{equation} \varDelta(\delta Z_\lambda^\pm) \sim \delta Z_\lambda^\pm\frac{\tau_{\mathrm{A}}}{\tau_{\mathrm{nl}}^\pm} \sim \frac{\delta Z^+_\lambda\delta Z^-_\lambda}{\lambda}\,\tau_{\mathrm{A}}, \end{equation}

assuming $\tau _{\mathrm {nl}}^\pm \sim \lambda /\delta Z^\mp _\lambda$. By definition of the WT regime, $\tau _{\mathrm {nl}}^\pm \gg \tau _{\mathrm {A}}$, so the amplitude change in any one interaction is small, $\varDelta (\delta Z_\lambda ^\pm ) \ll \delta Z_\lambda ^\pm$, and many such interactions are needed in order to change the amplitude $\delta Z_\lambda ^\pm$ by an amount comparable to itself, i.e., to ‘cascade’ the energy associated with scale $\lambda$ to smaller scales. Suppose that interactions occur all the time and that the kicks (4.4) accumulate as a random walk. Then the cascade time is $\tau _{\mathrm {c}}^\pm = N\tau _{\mathrm {A}}$ if after $N$ interactions the amplitude change is of order $\delta Z_\lambda ^\pm$:

(4.5)\begin{equation} \varDelta(\delta Z_\lambda^\pm)\sqrt{N}\sim \delta Z_\lambda^\pm \quad\Rightarrow\quad \frac{\tau_{\mathrm{A}}}{\tau_{\mathrm{nl}}^\pm}\sqrt{\frac{\tau_{\mathrm{c}}^\pm}{\tau_{\mathrm{A}}}} \sim 1 \quad\Rightarrow\quad \tau_{\mathrm{c}}^\pm \sim \frac{(\tau_{\mathrm{nl}}^\pm)^2}{\tau_{\mathrm{A}}}. \end{equation}

The standard Kolmogorov constant-flux requirement gives

(4.6)\begin{equation} \varepsilon^\pm \sim \frac{(\delta Z_\lambda^\pm)^2}{\tau_{\mathrm{c}}^\pm} \sim \frac{(\delta Z_\lambda^+)^2(\delta Z_\lambda^-)^2\tau_{\mathrm{A}}}{\lambda^2}. \end{equation}

Assuming for the moment that $\varepsilon ^+\sim \varepsilon ^-$ and, therefore, $\delta Z_\lambda ^+\sim \delta Z_\lambda ^-$, gets us the classic WT scalingFootnote ¹⁰

(4.7)\begin{equation} \delta Z_\lambda \sim \left(\frac{\varepsilon}{\tau_{\mathrm{A}}}\right)^{1/4}\lambda^{1/2} \quad\Leftrightarrow\quad E(k_{{\perp}})\sim \left(\frac{\varepsilon}{\tau_{\mathrm{A}}}\right)^{1/2}k_{{\perp}}^{{-}2}. \end{equation}

This scaling is indeed what one finds numerically (see figures 3 and 4) – it was first confirmed in early, semidirect simulations by Ng & Bhattacharjee (Reference Ng and Bhattacharjee1997) and Bhattacharjee & Ng (Reference Bhattacharjee and Ng2001), and then definitively by Perez & Boldyrev (Reference Perez and Boldyrev2008) and Boldyrev & Perez (Reference Boldyrev and Perez2009), leading the community to tick off WT as done and dusted.

Figure 4. Kinetic ($E$, solid lines) and magnetic ($M$, dotted lines with crosses) energy spectra for $k_{\parallel }=0$ (red), $k_{\parallel }=2\pi /L_{\parallel }$ (blue) and $k_{\parallel }=4\pi /L_{\parallel }$ (green) from an unpublished weak RMHD turbulence simulation by Yousef & Schekochihin (Reference Yousef and Schekochihin2009). The box size was $(L_{\perp },L_{\parallel })$ in the perpendicular and parallel directions, respectively, and the forcing was narrow-band, at $k_{\parallel }=2\pi /L_{\parallel }$ and $k_{\perp } = (1,2)\times 2\pi /L_{\perp }$, deep in the WT regime ($L_{\perp }\gg \lambda _{\mathrm {CB}}$). WT spectra for the case of broad-band forcing can be found in Perez & Boldyrev (Reference Perez and Boldyrev2008) and Boldyrev & Perez (Reference Boldyrev and Perez2009) and are discussed in Appendix A.4.

As anticipated in § 4.1, with the scaling (4.7), the ratio of the time scales can only stay small above a certain finite scale:

(4.8)\begin{equation} \frac{\tau_{\mathrm{A}}}{\tau_{\mathrm{nl}}} \sim \frac{\tau_{\mathrm{A}}\delta Z_\lambda}{\lambda} \sim \frac{\tau_{\mathrm{A}}^{3/4}\varepsilon^{1/4}}{\lambda^{1/2}} \ll 1 \quad\Leftrightarrow\quad \lambda \gg \varepsilon^{1/2}\tau_{\mathrm{A}}^{3/2}\sim\lambda_{\mathrm{CB}}, \end{equation}

where $\lambda _{\mathrm {CB}}$ is transition scale anticipated in (4.3). For $\lambda \lesssim \lambda _{\mathrm {CB}}$, turbulence becomes strong and, presumably, critically balanced. Thus, the WT cascade, by transferring energy to smaller scales, where nonlinear times are shorter, saws the seeds of its own destruction.

4.3. Imbalanced WT

What if $\varepsilon ^+\neq \varepsilon ^-$, viz., say, $\varepsilon ^+\gg \varepsilon ^-$? (If $\varepsilon ^+>\varepsilon ^-$ but both are of the same order, arguably the results obtained for $\varepsilon ^+\sim \varepsilon ^-$ should still work, at least on the ‘twiddle’ level.) Alas, (4.6) is patently incapable of accommodating such a case, an embarrassment first noticed by Dobrowolny, Mangeney & Veltri (Reference Dobrowolny, Mangeney and Veltri1980), who were attempting an IK-style, isotropic ($L_{\parallel }\sim \lambda$), imbalanced theory – quite wrong, as we now know (§ 2.3), but they correctly identified the issue with the imbalanced regime. They concluded that no imbalanced stationary state was possible except a pure Elsasser state. This may be true for (certain types of) decaying turbulence (see § 12.7), but is certainly not a satisfactory conclusion for a forced case where $\varepsilon ^\pm$ are externally prescribed. Arguably, the inability to describe the imbalanced case casts a shadow of doubt also on the validity of the argument in § 4.2 for balanced WT, as introducing the imbalance just lifts a kind of ‘degeneracy’ and perhaps highlights a problem with the whole story.

A way out of this difficulty, various versions of which have been explored by Galtier et al. (Reference Galtier, Nazarenko, Newell and Pouquet2000), Lithwick & Goldreich (Reference Lithwick and Goldreich2003), and Chandran (Reference Chandran2008), is to accept (4.6) but notice that it allows the two Elsasser fields to have different scaling exponents, $\delta Z^\pm _\lambda \propto \lambda ^{\gamma ^\pm }$, as long as they satisfy $\gamma ^+ + \gamma ^- = 1$. The corresponding 2D spectra of the two fields are

(4.9)\begin{equation} E_{\mathrm{2D}}^\pm(k_{{\perp}},k_{{\parallel}}) = f^\pm(k_{{\parallel}}) k_{{\perp}}^{\mu^\pm}, \quad \mu^+{+} \mu^- = -4, \end{equation}

because $\mu ^\pm = -2\gamma ^\pm -1$ and, WT permitting no changes in $k_{\parallel }$, the scaling arguments of § 4.2 apply to each $k_{\parallel }$ individually. One may then declare that the difference between $\varepsilon ^+$ and $\varepsilon ^-$ is hidden in the prefactors $f^\pm (k_{\parallel })$, which are non-universal, inaccessible to ‘twiddle’ scaling arguments about local interactions in the WT inertial range, and have to be fixed from outside it. At the large-scale end, one has to decide whether the outer scales for the two Elsasser fields are the same or different (Chandran Reference Chandran2008) and whether it is the fluxes $\varepsilon ^\pm$ or the fields’ energies at the outer scale(s) that it makes better sense to consider prescribed. At the dissipation scale, one has the option of ‘pinning’ the spectra to the same value (an idea due to Grappin, Leorat & Pouquet Reference Grappin, Leorat and Pouquet1983 and revived by Lithwick & Goldreich Reference Lithwick and Goldreich2003), and it must also be decided whether the two fields are required to start feeling viscosity at the same scale or one can do so before the other (see discussion in Beresnyak & Lazarian Reference Beresnyak and Lazarian2008, and, for strong imbalanced turbulence, in § 9.6.4). If WT breaks down before the dissipation scale is reached, some other set of ad hoc arrangements is required (see, e.g., Chandran Reference Chandran2008). Typically, the outcome is that the stronger field has a steeper spectrum than the weaker field, but their scalings are non-universal, i.e., they depend on the particular set up of the problem, at both macro- and micro-scales.

Another possibility is that (4.6) is wrong. Let me observe that the balanced ($\varepsilon ^+\sim \varepsilon ^-$, $\delta Z^+_\lambda \sim \delta Z^-_\lambda$) version of this scaling, i.e., the statement that the flux $\varepsilon$ is proportional to the fourth power of the amplitude, is less likely to be wrong than any particular assignment of ‘$+$'s and ‘$-$’s to these amplitudes: all it says is that the flux $\varepsilon$ is what it would have been in the case of strong interactions, $\sim \delta Z_\lambda ^2/\tau _{\mathrm {nl}}$ (cf. (2.9)), times the first power of the expansion parameter $\tau _{\mathrm {A}}/\tau _{\mathrm {nl}}$, i.e., the lowest order that $\varepsilon$ can be in a perturbation expansion in that parameter. Thus, one may doubt the validity of (4.6) for the imbalanced regime without rejecting the numerically confirmed $k_{\perp }^{-2}$ scaling of the balanced spectra. For example, in the (heuristic) scheme proposed by Schekochihin, Nazarenko & Yousef (Reference Schekochihin, Nazarenko and Yousef2012),

(4.10)\begin{equation} \varepsilon^\pm \sim \frac{(\delta Z_\lambda^\pm)^3\delta Z_\lambda^\mp\tau_{\mathrm{A}}}{\lambda^2}, \end{equation}

which changes nothing for balanced WT, but leads to a very different situation in the imbalanced case than (4.6), allowing perfectly good $k_{\perp }^{-2}$ spectra for both fields.

I do not go through all this in detail because, the WT regime being largely irrelevant (§ 4.1), it would also, if it really were non-universal, not be very interesting. If it is universal and something like (4.10) holds, that is interesting, but I do not know how to make much progress beyond Schekochihin et al. (Reference Schekochihin, Nazarenko and Yousef2012), whose theory does not quite match simulations (see § 4.4). I also do not know how to construct a theory of imbalanced WT that would connect smoothly to any believable theory of strong imbalanced turbulence (e.g., one presented in § 9.6). An interested reader will find some further, equally unsatisfactory, observations in Appendix A.6.

4.4. 2D condensate

It follows from the discussion in § 4.1 that the WT approximation in its standard form cannot, in fact, work for the turbulence of Alfvén waves, at least not formally, because in every three-wave interaction, one of the three waves has $k_{\parallel }=0$, so is not a wave at all, but a zero-frequency 2D perturbation, for which the nonlinear interactions are the dominant influence. If such $k_{\parallel }=0$ perturbations are forbidden, i.e., if displacements vanish at infinity, one must consider four-wave interactions (i.e., go to next order in $\tau _{\mathrm {A}}/\tau _{\mathrm {nl}}$), which gives rise to an apparently legitimate WT state, different from (4.7) (Sridhar & Goldreich Reference Sridhar and Goldreich1994). There is no particular reason to think, however, that such a restriction on displacements is legitimate in a general physical situation (Ng & Bhattacharjee Reference Ng and Bhattacharjee1996) or, even if one starts with no energy at $k_{\parallel }=0$, that such a state can be maintained, except in a box with field lines nailed down at the boundaries – failing such restrictions, a 2D ‘condensate’ must emerge (and does, in numerical simulations: see Boldyrev & Perez Reference Boldyrev and Perez2009, Wang, Boldyrev & Perez Reference Wang, Boldyrev and Perez2011, Meyrand, Kiyani & Galtier Reference Meyrand, Kiyani and Galtier2015; Meyrand et al. Reference Meyrand, Galtier and Kiyani2016, and figure 4).

Mathematically, this becomes quite obvious if we represent the solutions to (3.1) as

(4.11)\begin{equation} \boldsymbol{Z}_{{\perp}}^\pm(t,\boldsymbol{r}) = \sum_{k_{{\parallel}}}\boldsymbol{Z}^\pm_{k_{{\parallel}}}(t,x,y) e^{ik_{{\parallel}}(z \pm v_{\mathrm{A}} t)} \end{equation}

and separate the $k_{\parallel }=0$ modes from the rest:

(4.12)\begin{align} \partial_t \boldsymbol{Z}_0^\pm + \hat{{\mathcal{P}}}\boldsymbol{Z}_0^\mp \cdot \boldsymbol{\nabla}_{{\perp}}\boldsymbol{Z}_0^\pm = & -\sum_{k_{{\parallel}}\neq 0}\hat{{\mathcal{P}}}\boldsymbol{Z}_{k_{{\parallel}}}^\mp \cdot \boldsymbol{\nabla}_{{\perp}}\boldsymbol{Z}_{{-}k_{{\parallel}}}^\pm e^{{\mp} i 2k_{\parallel} v_{\mathrm{A}} t}, \end{align}

(4.13)\begin{align} \partial_t \boldsymbol{Z}_{k_{{\parallel}}}^\pm + \hat{{\mathcal{P}}}\boldsymbol{Z}_0^\mp \cdot \boldsymbol{\nabla}_{{\perp}}\boldsymbol{Z}_{k_{{\parallel}}}^\pm = & - \sum_{p_{{\parallel}}\neq 0}\hat{{\mathcal{P}}}\boldsymbol{Z}_{p_{{\parallel}}}^\mp \cdot \boldsymbol{\nabla}_{{\perp}}\boldsymbol{Z}_{k_{{\parallel}}-p_{{\parallel}}}^\pm e^{{\mp} i 2p_{{\parallel}}v_{\mathrm{A}} t}, \end{align}

where $\hat {{\mathcal {P}}}$ is the projection operator that takes care of the pressure term (see (3.2)) and has been introduced for brevity; forcing and dissipation terms have been dropped. The first of these equations, (4.12), describes the condensate – two real fields $\boldsymbol {Z}_0^\pm (x,y)$ advecting each other in the 2D plane while subject to an oscillating ‘force’ due to the mutual coupling of the Alfvén waves $\boldsymbol {Z}^\pm _{k_{\parallel }}$. These Alfvén waves, described by (4.13), are advected by the 2D field and also by each other, but the latter interaction has an oscillating factor and vanishes in the WT approximation. Even if only the Alfvén waves are forced and the condensate is not, the condensate will nevertheless be built up.

Returning to three-wave interactions then (where one of the waves is not a wave), the traditional approach has been to ignore the inapplicability of the WT approximation to the $k_{\parallel }=0$ modes by conjecturing that the function $f^\pm (k_{\parallel })$ in (4.9) is flat around $k_{\parallel }=0$ – the hypothesis of ‘spectral continuity’. One can then press on with putting MHD through the WT analytical grinder, find an evolution equation for the spectra and show that it has steady-state, constant-flux solutions of the form (4.9). This is what was done in the now-classic paper by Galtier et al. (Reference Galtier, Nazarenko, Newell and Pouquet2000) (see Appendices A.2 and A.3). In balanced turbulence, obviously, $\mu ^+=\mu ^-=2$, and we are back to (4.7); in imbalanced turbulence, one needs further scheming – see references in § 4.3.

Nazarenko (Reference Nazarenko2007) argues that the hypothesis of spectral continuity is certainly false if the nonlinear broadening of the waves’ frequencies, of order $\tau _{\mathrm {nl}}^{-1}$, is smaller than the linear frequency associated with the spacing of the $k_{\parallel }$ ‘grid’ ($=2\pi /L_{\parallel }$, the inverse parallel ‘box’ size) – i.e., if the Alfvénic perturbations at the longest finite parallel scale in the system are already in the WT limit (4.1), $v_{\mathrm {A}}/L_{\parallel }\gg \tau _{\mathrm {nl}}^{-1}$. He is right. Figure 4 is taken from a (sadly, unpublished) numerical study of weak RMHD turbulence by Yousef & Schekochihin (Reference Yousef and Schekochihin2009), who forced Alfvén waves at $k_{\parallel } = 2\pi /L_{\parallel }$, where $L_{\parallel }$ was the box size. It shows that, while the $k_{\perp }^{-2}$ scaling of the $k_{\parallel }=2\pi /L_{\parallel }$ modes is undeniable, the spectra for the unforced modes ($k_{\parallel } = 0$ and $k_{\parallel } =$ multiples of $2\pi /L_{\parallel }$) are dramatically shallower. Similar spectra were reported by Bigot & Galtier (Reference Bigot and Galtier2011) and by Meyrand et al. (Reference Meyrand, Kiyani and Galtier2015). Qualitatively similar spectra (and a simple mechanism for how they might form) were also proposed by Schekochihin et al. (Reference Schekochihin, Nazarenko and Yousef2012) – but their theory fails quantitatively, with the spectra that it predicts for all unforced modes at least one power of $k_{\perp }$ steeper than the numerical ones (e.g., their $k_{\parallel }=0$ condensate has a $\propto k_{\perp }^{-1}$ spectrum, while simulations suggest $\propto k _{\perp }^0$).

Nazarenko (Reference Nazarenko2007) expects that the conventional WT theory should survive when $k_{\parallel} v _{\mathrm {A}}\gg \tau _{\mathrm {nl}}^{-1}\gg v _{\mathrm {A}}/L_{\parallel }$. This is a situation that should be realisable in a system that is weakly and randomly forced in a broad band of frequencies (and, therefore, parallel wavenumbers). In Appendix A.4, I discuss how, and in what sense, one might defend spectral continuity for such a system; I argue that the 2D condensate in this case is a strongly turbulent, critically balanced subsystem constantly fed by the weakly turbulent waves and developing a falsifiable set of scalings, which are, indeed, continuous with the WT scalings. While there are some indications (from the simulations by Wang et al. Reference Wang, Boldyrev and Perez2011; see Appendix A.5) that these scalings might be right, I have not seen spectral continuity corroborated numerically in a definitive fashion, as even Perez & Boldyrev (Reference Perez and Boldyrev2008) and Boldyrev & Perez (Reference Boldyrev and Perez2009), who took great care to force in a broad band of $k_{\parallel }$ to make sure the conventional WT theory did apply, saw a distinct dip in $f^\pm (k_{\parallel })$ at $k_{\parallel } = 0$, associated with an emergent condensate (which is magnetically dominated; see § 10.3 and Appendix A.5). The same was true in the decaying simulations of ?see the upper panel of figure 3) meyrand15,meyrand16, where an initial small-amplitude (and so WT-compliant) state had the choice to evolve towards a continuous parallel spectrum, but refused to do so, again developing a $k_{\parallel }=0$ condensate with dramatically distinct properties, including a high degree of intermittency and a spectrum quite similar to figure 4.

Thus, the conventional WT theory is at best incomplete and at worst wrong. It is discussed further in Appendix A, where I review the WT's derivation, speculate about the structure of the condensate, and discuss a number of other WT-related issues. Here, having flagged these issues, I want to halt this digression into matters that are, arguably, of little impact, and move on to the physics-rich core of the MHD-turbulence theory.

5.1. Critical balance

Section 4 can be viewed as one long protracted justificatory piece in favour of critical balance: even if an ensemble of high-frequency Alfvén waves is stirred up very gently ($\tau _{\mathrm {nl}}\gg \tau _{\mathrm {A}}$), it will, at small enough scales, get itself into the strong-turbulence regime ($\tau _{\mathrm {nl}}\sim \tau _{\mathrm {A}}$). The opposite limit, a 2D regime with $\tau _{\mathrm {nl}}\ll \tau _{\mathrm {A}}$, is unsustainable for the very simple reason of causality: as information in an RMHD system propagates along $\boldsymbol {B}_0$ at speed $v_{\mathrm {A}}$, no structure longer than $l_{\parallel }\sim v _{\mathrm {A}}\tau _{\mathrm {nl}}$ can be kept coherent and so will break up (see Boldyrev Reference Boldyrev2005, Nazarenko & Schekochihin Reference Nazarenko and Schekochihin2011 and figure 5).

Figure 5. Critical balance in a ($2+1$)D system supporting both nonlinearity and waves (RMHD).

It is worth mentioning in passing that the CB turns out to be a very robust feature of the turbulence in the following interesting sense. With a certain appropriate definition of $\tau _{\mathrm {nl}}$ (which will be explained in § 6.1), the ratio $\tau _{\mathrm {A}}/\tau _{\mathrm {nl}}$ has been found (numerically) by Mallet et al. (Reference Mallet, Schekochihin and Chandran2015) to have a scale-invariant distribution (figure 6), a property that they dubbed refined critical balance (RCB). It gives a quantitative meaning to the somewhat vague statement $\tau _{\mathrm {A}}/\tau _{\mathrm {nl}}\sim 1$ – and becomes important in the (as it turns out, unavoidable) discussion of intermittency of MHD turbulence (§ 6.4.2).

Figure 6. Refined critical balance: this figure, taken from Mallet et al. (Reference Mallet, Schekochihin and Chandran2015), shows the probability density function (PDF) of the ratio $\chi ^+=\tau _{\mathrm {A}}/\tau _{\mathrm {nl}}^+$ with $\tau _{\mathrm {nl}}^+$ defined by (6.4). In fact, 17 PDFs are plotted here, taken at different scales within an approximately decade-wide inertial range (this was a $1024^3$ RMHD simulation) – the corresponding lines are in colour shades from blue (smaller scales) to red (larger scales), but this is barely visible because the PDFs all collapse on top of each other. The inset shows that the self-similarity does not work if $\tau _{\mathrm {nl}}^+$ is defined without the alignment angle (see § 6). (Reprinted from Mallet et al. Reference Mallet, Schekochihin and Chandran2015 by permission of the Royal Astronomical Society.)

5.2. Parallel cascade

The most straightforward – and the least controversial – consequence of CB is the scaling of parallel increments. I have already derived this result in (2.10), but let me now restate it using Elsasser fields. If it is the case that the nonlinear-interaction time and, therefore, the cascade time for $\boldsymbol {Z}_{\perp }^\pm$ are approximately the same as their propagation time $\tau _{\mathrm {A}}\sim l _{\parallel }/v_{\mathrm {A}}$, then the parallel increments $\delta Z_{l_{\parallel }}^\pm$ satisfy

(5.1)\begin{equation} \frac{(\delta Z_{l_{{\parallel}}}^\pm)^2}{\tau_{\mathrm{A}}} \sim \varepsilon^\pm\quad\Rightarrow\quad \delta Z^\pm_{l_{{\parallel}}} \sim \left(\frac{\varepsilon^\pm l_{{\parallel}}}{v_{\mathrm{A}}}\right)^{1/2} \quad\Leftrightarrow\quad E^\pm(k_{{\parallel}}) \sim \frac{\varepsilon^\pm}{v_{\mathrm{A}}}k_{{\parallel}}^{{-}2}. \end{equation}

Beresnyak (Reference Beresnyak2012a, Reference Beresnyak2015) gives two rather elegant (and related) arguments in favour of the scaling (5.1), alongside robust numerical evidence presented in the latter paper.Footnote ¹¹ First, he argues that the scaling relation (5.1) can be obtained by dimensional analysis because the RMHD equations (3.1) stay invariant if $v_{\mathrm {A}}$ and $1/k_{\parallel }$ are scaled simultaneously (see (3.5)) and so these two quantities must always appear in the combination $k_{\parallel} v _{\mathrm {A}}$ in scaling relations for any physical quantities – in the case of (5.1), energy, or field increment. Secondly, Beresnyak (Reference Beresnyak2015) notes that following the structure of the fluctuating field (calculating its increments) along the field line (in the positive $\boldsymbol {B}_0$ direction) is the MHD equivalent of following its evolution forward (for $\boldsymbol {Z}_{\perp }^-$) or backward (for $\boldsymbol {Z}_{\perp }^+$) in time and it should, therefore, be possible to infer the parallel spectrum (5.1) from the Lagrangian frequency spectrum of the turbulence (as opposed to the Eulerian one, which is dominated by large-scale sweeping effects: see Lugones et al. Reference Lugones, Dmitruk, Mininni, Wan and Matthaeus2016, Reference Lugones, Dmitruk, Mininni, Pouquet and Matthaeus2019). Estimating the energy flux as the rate of change of energy in a fluid element in the Lagrangian frame (i.e., excluding sweeping by large eddies), one obtains (Landau & Lifshitz Reference Landau and Lifshitz1987; Corrsin Reference Corrsin1963)

(5.2)\begin{equation} \varepsilon^\pm \sim (\delta Z_\tau^\pm)^2\tau^{{-}1} \quad\Leftrightarrow\quad E^\pm(\omega) \sim \varepsilon^\pm\omega^{{-}2}, \end{equation}

where $\delta Z_\tau ^\pm$ is the Lagrangian field increment over time interval $\tau$. Then (5.1) is recovered from (5.2) by changing variables $\omega =k_{\parallel} v _{\mathrm {A}}$ and letting $E^\pm (\omega )\mathrm {d}\omega = E^\pm (k_{\parallel })\mathrm {d}k_{\parallel }$.

Thus, the parallel cascade and the associated scaling (5.1) appear to be a very simple and solid property of MHD turbulence. What happens in the perpendicular direction is a more complicated story.

5.3. Local, scale-dependent anisotropy

Using instead of the parallel increments the perpendicular ones $\delta Z^\pm _\lambda$ and substituting the nonlinear time

(5.3)\begin{equation} \tau_{\mathrm{nl}}^\pm \sim \frac{\lambda}{\delta Z^\mp_\lambda} \end{equation}

for the cascade time, we recover (2.9):Footnote ¹²

(5.4)\begin{align} \frac{(\delta Z_\lambda^\pm)^2}{\tau_{\mathrm{nl}}^\pm} \sim \varepsilon^\pm &\quad\Rightarrow\quad \frac{\delta Z_\lambda^+}{\delta Z_\lambda^-}\sim\frac{\varepsilon^+}{\varepsilon^-},\quad \delta Z^\pm_\lambda \sim (\tilde\varepsilon^\pm\lambda)^{1/3}, \quad \tilde\varepsilon^\pm\equiv \frac{(\varepsilon^\pm)^2}{\varepsilon^\mp} \nonumber\\ &\quad\Rightarrow\quad E^\pm(k_{{\perp}}) \sim (\tilde\varepsilon^\pm)^{2/3}k_{{\perp}}^{{-}5/3}. \end{align}

Treating $\delta Z^\pm _\lambda$ and $\delta Z^\pm _{l_{\parallel }}$ as increments for the same structure, but measured across and along the field, and setting them equal to each other, we find a relationship between the parallel and perpendicular scales – the scale-dependent anisotropy (2.11):

(5.5)\begin{equation} l_{{\parallel}}^\pm \sim v_{\mathrm{A}} (\tilde\varepsilon^\mp)^{{-}1/3}\lambda^{2/3}. \end{equation}

The fact of scale-dependent anisotropy of MHD turbulence (if, in retrospect, not with the same confidence as the scaling (5.5)) was confirmed numerically by Cho & Vishniac (Reference Cho and Vishniac2000) and Maron & Goldreich (Reference Maron and Goldreich2001) and, in a rare triumph of theory correctly anticipating measurement, observed in the solar wind by Horbury et al. (Reference Horbury, Forman and Oughton2008), followed by many others (e.g., Podesta Reference Podesta2009; Wicks et al. Reference Wicks, Horbury, Chen and Schekochihin2010; Luo & Wu Reference Luo and Wu2010; Chen et al. Reference Chen, Mallet, Yousef, Schekochihin and Horbury2011 – a complete list is impossible here as this has now become an industry, as successful ideas do; see Chen Reference Chen2016 for a review). Figure 7 shows some of the first of those results. An important nuance is that, in order to see scale-dependent anisotropy, one must measure the parallel correlations along the perturbed, rather than global, mean magnetic field.Footnote ¹³ The reason for this is as follows.

Figure 7. (a) Parallel ($P_\parallel$) and perpendicular ($P_\perp$) spectra (Fourier and wavelet) of the magnetic fluctuations in the solar wind, measured by the Ulysses spacecraft and computed by Wicks et al. (Reference Wicks, Horbury, Chen and Schekochihin2010), with frequencies $f$ converted to wavenumbers $k$ using the Taylor hypothesis (reprinted from Wicks et al. Reference Wicks, Horbury, Chen and Schekochihin2010 by permission of the Royal Astronomical Society). (b) An earlier (historic, the first ever) measurement by Horbury, Forman & Oughton (Reference Horbury, Forman and Oughton2008) of the spectral index of these spectra as a function of angle to the local mean field (reprinted with permission from Horbury et al. Reference Horbury, Forman and Oughton2008, copyright (2008) by the American Physical Society).

Both the causality argument (Boldyrev Reference Boldyrev2005; Nazarenko & Schekochihin Reference Nazarenko and Schekochihin2011) and the Lagrangian-frequency one (Beresnyak Reference Beresnyak2015) that I invoked in §§ 5.1 and 5.2 to justify long parallel coherence lengths of the MHD fluctuations rely on the ability of Alfvénic perturbations to propagate along the magnetic field. Physically, a small such perturbation on any given scale does not know the difference between a larger perturbation on, say, a few times its scale, and the ‘true’ mean field (whatever that is, outside the ideal world of periodic simulation boxes). Thus, it will propagate along the local field and so it is along the local field that the arguments based on this propagation will apply. What if we instead measure correlations along the global mean field or, more generally, along some coarse-grained version of the exact field? Let that coarse-grained field be the average over all perpendicular scales at and below some $L_{\perp }$ (to get the global mean field, make $L_{\perp }$ the outer scale). Define Elsasser-field increments between pairs of points separated by a vector ${\boldsymbol {l}}$,

(5.6)\begin{equation} \delta\boldsymbol{Z}_{\boldsymbol{l}}^\pm = \boldsymbol{Z}_{{\perp}}^\pm(\boldsymbol{r} + {\boldsymbol{l}}) - \boldsymbol{Z}_{{\perp}}^\pm(\boldsymbol{r}), \end{equation}

and consider ${\boldsymbol {l}}$ along the exact magnetic field vs. ${\boldsymbol {l}}$ along our coarse-grained field. The perpendicular distance by which the latter vector will veer off the field line (figure 8) will be dominated by the magnetic perturbation at the largest scale that was not included in the coarse-grained field:

(5.7)\begin{equation} \varDelta l_{{\perp}}\sim l\frac{\delta b_{L_{{\perp}}}}{v_{\mathrm{A}}}. \end{equation}

If we are trying to capture parallel correlations corresponding to perturbations with perpendicular scale $\lambda \ll L _{\perp }$, then, using CB, $l/v_{\mathrm {A}}\sim \tau _{\mathrm {nl}}$, and (5.3) with $\delta Z_\lambda ^\pm \sim \delta b_\lambda$, we conclude that

(5.8)\begin{equation} \varDelta l_{{\perp}} \sim \lambda\frac{\delta b_{L_{{\perp}}}}{\delta b_\lambda} \gg \lambda, \end{equation}

i.e., in such a measurement, the parallel correlations are swamped by perpendicular decorrelation, unless, in fact, $\lambda \sim L_{\perp }$ or larger (there is no such problem with measuring perpendicular correlations: small changes in a separation vector ${\boldsymbol {l}}$ taken perpendicular to the global vs. exact field make no difference).

Figure 8. Measuring correlations along local vs. global mean field. True parallel correlations cannot be captured by a measurement along the global field $\boldsymbol {B}_0$ if the distance $\varDelta l_\perp$ (see (5.7)) by which the point-separation vector ${\boldsymbol {l}}$ along $\boldsymbol {B}_0$ ‘slips’ off the exact field line ($\boldsymbol {B}_0+\boldsymbol {b}_{\perp }$) is greater than the perpendicular decorrelation length $\lambda$ between ‘neighbouring’ field lines.

Consequently, the easiest practical way to extract correlations along the local field from either observed or numerically simulated turbulence (Chen et al. Reference Chen, Mallet, Yousef, Schekochihin and Horbury2011) is to measure field increments (5.6) for many different separation vectors ${\boldsymbol {l}}$ and to calculate for each such increment the angle between ${\boldsymbol {l}}$ and the ‘local mean field’ $\boldsymbol {B}_{\mathrm {loc}}$ defined as the arithmetic mean of the magnetic field measured at the two points involved:

(5.9)\begin{equation} \cos\phi =\frac{{\boldsymbol{l}}\cdot\boldsymbol{B}_{\mathrm{loc}}}{|{\boldsymbol{l}}||\boldsymbol{B}_{\mathrm{loc}}|}, \quad \boldsymbol{B}_{\mathrm{loc}} = \boldsymbol{B}_0 + \frac{\boldsymbol{b}_{{\perp}}(\boldsymbol{r} + {\boldsymbol{l}}) + \boldsymbol{b}_{{\perp}}(\boldsymbol{r})}{2}. \end{equation}

This amounts to coarse-graining the field always at the right scale (just) for the correlations that are being probed. One can then measure (for example) perpendicular and parallel structure functions as conditional averages:

(5.10)\begin{align} \langle(\delta Z^\pm_\lambda)^n\rangle &= \langle|\delta\boldsymbol{Z}_{\boldsymbol{l}}^\pm|^n|\phi=90^{\rm o}\rangle, \end{align}

(5.11)\begin{align} \langle(\delta Z^\pm_{l_{{\parallel}}})^n\rangle &= \langle|\delta\boldsymbol{Z}_{\boldsymbol{l}}^\pm|^n|\phi=0\rangle, \end{align}

and similarly for intermediate values of $\phi$ (as explained above, the difference between $\boldsymbol {B}_0$ and $\boldsymbol {B}_{\mathrm {loc}}$ matters only for small $\phi$). Thus, in general, one measures

(5.12)\begin{equation} \langle|\delta\boldsymbol{Z}_{\boldsymbol{l}}^\pm|^n|\phi\rangle \propto l^{\zeta_n(\phi)}. \end{equation}

Alternatively, in simulations, one can simply follow field lines to get $\delta Z^\pm _{l_{\parallel }}$ (Cho & Vishniac Reference Cho and Vishniac2000; Maron & Goldreich Reference Maron and Goldreich2001) or, as was initially done in the solar wind, use local wavelet spectra (Horbury et al. Reference Horbury, Forman and Oughton2008; Podesta Reference Podesta2009; Wicks et al. Reference Wicks, Horbury, Chen and Schekochihin2010).

It turns out (see references cited above and innumerable others) that, quite robustly, $\zeta _2(0)=1$, consistent with (5.1), whereas $\zeta _2(90^\textrm {o})$ is typically between $2/3$ and $1/2$, i.e., between GS95 and IK, in the solar wind, and rather closer to $1/2$ in numerical simulations – although this, as I will discuss in §§ 6.3 and 6.4, has been hotly disputed by Beresnyak (Reference Beresnyak2011, Reference Beresnyak2012a, Reference Beresnyak2014b, Reference Beresnyak2019).

Thus, while little doubt remains about the reality of scale-dependent anisotropy (although not necessarily of the specific scaling (5.5)) and of the $k_{\parallel }^{-2}$ spectrum (5.1), both arising from the GS95 theory, the GS95 prediction for the perpendicular spectrum (5.4) has continued to be suspect and controversial.

6.1. Alignment and anisotropy in the perpendicular plane

Let me, however, deviate from Boldyrev's original case for alignmentFootnote ¹⁶ and follow instead Chandran et al. (Reference Chandran, Schekochihin and Mallet2015) in positing that, as the two Elsasser fields advect each other, they will shear each other into an increasingly parallel configuration. As this occurs, the gradient of the advectee (say, $\boldsymbol {Z}_{\perp }^-$) in the direction of the advector ($\boldsymbol {Z}_{\perp }^+$) will get smaller than the gradient of either of them in the direction transverse to both.

In general, therefore, we must allow the possibility of a local anisotropy in the 2D plane perpendicular to the mean magnetic field (figure 9). If we do that, our ‘twiddle’ estimate of the nonlinear term in (3.1) becomes

(6.1)\begin{equation} \boldsymbol{Z}_{{\perp}}^+ \cdot \boldsymbol{\nabla}_{{\perp}}\boldsymbol{Z}_{{\perp}}^- \sim \frac{\delta Z_\lambda^+\delta Z_\lambda^-}{\xi^-}, \end{equation}

where $\xi ^-$ is the scale of variation of $\boldsymbol {Z}_{\perp }^-$ (advectee) in the direction of $\boldsymbol {Z}_{\perp }^+$ (advector), taken at scale $\lambda$, which is the scale of $\boldsymbol {Z}_{\perp }^\pm$'s variation in the direction perpendicular both to $\boldsymbol {Z}_{\perp }^+$ and to $\boldsymbol {B}_0$ (all interactions are still assumed local in scale). But, by elementary vector calculus,

(6.2)\begin{equation} \boldsymbol{Z}_{{\perp}}^+ \cdot \boldsymbol{\nabla}_{{\perp}}\boldsymbol{Z}_{{\perp}}^- - \boldsymbol{Z}_{{\perp}}^- \cdot \boldsymbol{\nabla}_{{\perp}}\boldsymbol{Z}_{{\perp}}^+ = \boldsymbol{\nabla}_{{\perp}}\times\left(\boldsymbol{Z}_{{\perp}}^- \times \boldsymbol{Z}_{{\perp}}^+\right) \sim \frac{\delta Z_\lambda^+\delta Z_\lambda^-\sin\theta_\lambda}{\lambda}, \end{equation}

where $\theta _\lambda$ the angle between the two Elsasser fields taken at scale $\lambda$. Under the scheme whereby the two fields shear each other into alignment in an, on average, symmetric fashion, we have $\xi ^+ \sim \xi ^-$ and

(6.3)\begin{equation} \frac{\lambda}{\xi} \sim \sin\theta_\lambda, \end{equation}

i.e., $\sin \theta _\lambda$ is the aspect ratio of the field structures in the perpendicular 2D plane.Footnote ¹⁷

Figure 9. Cartoon of a GS95 eddy (a) vs. a Boldyrev (Reference Boldyrev2006) aligned eddy (b). The latter has three scales: $l_{\parallel }\gg \xi \gg \lambda$ (along $\boldsymbol {B}_0$, along $\boldsymbol {b}_{\perp }$, and transverse to both). This picture is adapted from Boldyrev (Reference Boldyrev2006) (reprinted with permission from Boldyrev Reference Boldyrev2006, copyright (2006) by the American Physical Society). In the context of my discussion, the fluctuation direction should, in fact, be thought of as along $\boldsymbol {Z}_{\perp }^+$ or $\boldsymbol {Z}_{\perp }^-$ (see figure 38).

In view of all this, we ought to replace the estimate (5.3) of the nonlinear time with

(6.4)\begin{equation} \tau_{\mathrm{nl}}^\pm \sim \frac{\xi}{\delta Z_\lambda^\mp} \sim \frac{\lambda}{\delta Z_\lambda^\mp\sin\theta_\lambda}. \end{equation}

If the last expression is correct, then the more aligned the two fields are the more the nonlinearity is reduced compared to the ‘naïve’GS95 estimate (5.3). The challenge is to work out the scale dependence of the reduction factor $\sin \theta _\lambda$, which is what Boldyrev (Reference Boldyrev2006) did.

6.2. Boldyrev's alignment conjecture

A version of Boldyrev's argumentFootnote ¹⁸ is to conjecture that the fields want to be misaligned as little as possible and that this minimal degree of misalignment is set by a kind of uncertainty principle: since the direction of the local magnetic field along which these perturbations propagate can itself only be defined within a small angle $\sim \delta b_\lambda /v_{\mathrm {A}}$, the two Elsasser fields (or the velocity and the magnetic field) cannot be aligned any more precisely than this and so

(6.5)\begin{equation} \sin\theta_\lambda \sim \theta_\lambda \sim \frac{\delta b_\lambda}{v_{\mathrm{A}}}\ll1. \end{equation}

Alignment and imbalance (local or global), and, therefore the relative magnitudes of magnetic, velocity, and Elsasser fields, can be linked in a non-trivial way, which is still a matter of some debate, but I do not wish be distracted and so will postpone the discussion of that topic to § 9. Till then, I shall assume everywhere that

(6.6)\begin{equation} \varepsilon^+ \sim \varepsilon^-\quad\Rightarrow\quad \delta Z_\lambda^+ \sim \delta Z_\lambda^- \sim \delta u_\lambda\sim\delta b_\lambda. \end{equation}

This allows (6.5) to be combined with (6.4) and yield

(6.7)\begin{equation} \tau_{\mathrm{nl}} \sim \frac{\lambda v_{\mathrm{A}}}{\delta Z_\lambda^2}. \end{equation}

The constancy of flux then implies immediatelyFootnote ¹⁹

(6.8)\begin{equation} \frac{\delta Z_\lambda^2}{\tau_{\mathrm{nl}}} \sim \varepsilon\quad\Rightarrow\quad \delta Z_\lambda \sim (\varepsilon v_{\mathrm{A}}\lambda)^{1/4} \quad\Leftrightarrow\quad E(k_{{\perp}}) \sim (\varepsilon v_{\mathrm{A}})^{1/2}k_{{\perp}}^{{-}3/2}. \end{equation}

In dimensional terms, this has brought us back to the IK spectrum (2.5), except the wavenumber is now the perpendicular wavenumber and both anisotropy and CB are retained, although the relationship between the parallel and perpendicular scales changes:

(6.9)\begin{equation} \tau_{\mathrm{nl}} \sim \frac{l_{{\parallel}}}{v_{\mathrm{A}}}\quad\Rightarrow\quad l_{{\parallel}} \sim v_{\mathrm{A}}^{3/2}\varepsilon^{{-}1/2}\lambda^{1/2}. \end{equation}

Since CB remains in force, the parallel cascade stays the same as described in § 5.2.

For imminent use in what follows, let us compute the extent of the inertial range that this aligned cascade is supposed to span. Comparing the nonlinear cascade time (6.7) with the Ohmic diffusion time (assuming, for convenience that the magnetic diffusivity $\eta$ is either the same or larger than the kinematic viscosity of our MHD fluid), we find

(6.10)\begin{equation} \tau_{\mathrm{nl}} \sim \left(\frac{v_{\mathrm{A}}\lambda}{\varepsilon}\right)^{1/2}\ll \tau_{\eta} \sim \frac{\lambda^2}{\eta} \quad\Leftrightarrow\quad\lambda \gg \eta^{2/3}\left(\frac{v_{\mathrm{A}}}{\varepsilon}\right)^{1/3}\equiv\lambda_{\eta}, \end{equation}

where $\lambda _{\eta }$ is the cutoff scale – the Kolmogorov scale for this turbulence. For comparison, note that the same calculation based on the GS95 scalings (5.3) and (5.4) gives

(6.11)\begin{equation} \tau_{\mathrm{nl}}^\mathrm{GS95} \sim \varepsilon^{{-}1/3}\lambda^{2/3} \ll \tau_{\eta} \sim \frac{\lambda^2}{\eta} \quad\Leftrightarrow\quad \lambda \gg \eta^{3/4}\varepsilon^{{-}1/4}\equiv\lambda_{\eta}^\mathrm{GS95}, \end{equation}

where $\lambda _{\eta }^\mathrm {GS95}$ is the classic Kolmogorov scale.

If one embraces (6.8), one could argue that Kraichnan's dimensional argument was actually right, but it should have been used with $k_{\perp }$, rather than with $k$, because $k_{\parallel }$ is not a ‘nonlinear’ dimension. This is the style of reasoning that Kraichnan himself might have found attractive. We are about to see, however, that the result (6.8) also runs into serious trouble and needs revision.

6.3. Plot thickens

This is a very appealing theory, whose main conclusions were rapidly confirmed by a programme of numerical simulations undertaken by Boldyrev's group – in particular, the angle between velocity and magnetic field, measured in a certain opportune way,Footnote ²⁰ was reported to scale according to $\theta _\lambda \propto \lambda ^{1/4}$, as implied by (6.8) and (6.5) (Mason, Cattaneo & Boldyrev Reference Mason, Cattaneo and Boldyrev2006, Reference Mason, Cattaneo and Boldyrev2008; Mason et al. Reference Mason, Perez, Cattaneo and Boldyrev2011, Reference Mason, Perez, Boldyrev and Cattaneo2012; Perez et al. Reference Perez, Mason, Boldyrev and Cattaneo2012, Reference Perez, Mason, Boldyrev and Cattaneo2014b). The same papers confirmed the earlier numerical claims that the spectrum of MHD turbulence indeed scaled as $k_{\perp }^{-3/2}$ (figure 10a). However, the legitimacy of this conclusion was contested by Beresnyak (Reference Beresnyak2011, Reference Beresnyak2012a, Reference Beresnyak2014b, Reference Beresnyak2019), who had, in fact, been first to focus on the alignment effect numerically (Beresnyak & Lazarian Reference Beresnyak and Lazarian2006), but disputed that the numerical spectra exhibiting it were converged and argued that systematic convergence tests in fact favoured a trend towards a $k_{\perp }^{-5/3}$ spectrum at asymptotically small scales. His point was that convergence of spectra with increasing resolution ought to be checked from the dissipative end of the inertial interval and that rescaling the spectra in his simulations to the Kolmogorov scale (6.11) gave a better data collapse than rescaling them to Boldyrev's cutoff scale (6.10) (figure 10(b); to be precise, in Beresnyak Reference Beresnyak2014b, he reports that the best numerical convergence obtained by this method is to $k_{\perp }^{-1.7}$). Despite the sound and fury of the ensuing debate about the quality of the two competing sets of numerics (Perez et al. Reference Perez, Mason, Boldyrev and Cattaneo2014a; Beresnyak Reference Beresnyak2013, Reference Beresnyak2014a), it would not necessarily be obvious to anyone who took a look at their papers that their raw numerical results themselves were in fact all that different – certainly not as different as the interpretation of these results by their authors. Without dwelling on either, however, let me focus instead on a conceptual wrinkle in Boldyrev's original argument that Beresnyak (Reference Beresnyak2011) spotted and that cannot be easily dismissed.

Figure 10. The best-resolved currently available spectra of RMHD turbulence. (a) From simulations by Perez et al. (Reference Perez, Mason, Boldyrev and Cattaneo2012) (their figure 1), with Laplacian viscosity and resolution up to $2048^2\times 512$. (b) From simulations by Beresnyak (Reference Beresnyak2014b) (his figure 1, © AAS, reproduced with permission), with Laplacian viscosity (a) and with fourth-order hyperviscosity (b); the resolution for the three spectra is $1024^3$, $2048^3$ and $4096^3$. His spectra are rescaled to Kolmogorov scale (6.11) (which he denotes $\eta$). He finds poorer convergence (see his figure 2) when he rescales to Boldyrev's scale (6.24). Perez et al. (Reference Perez, Mason, Boldyrev and Cattaneo2012) appear to get a somewhat better outcome (see their figure 8) if they plot their spectra vs. $k_{\perp }\lambda _{\eta }$ where $\lambda _{\eta }$ is given by (6.24) with $\lambda _{\mathrm {CB}}$ computed in each simulation as the normalisation constant in the scaling (6.22) of $\sin \theta _\lambda$ (in their analysis, however, this is the angle between velocity and magnetic perturbations, not the Elsasser fields).

In the RMHD limit (whose applicability to MHD turbulence at sufficiently small scales we have no reason to doubt), $\delta b_\lambda /v_{\mathrm {A}}$ is an arbitrarily small quantity, and so must then be, according to (6.5), the alignment angle $\sin \theta _\lambda$. Introducing such a large depletion of the nonlinearity into (3.1) would abolish it completely in the RMHD ordering and render the system linear. The only way to keep the nonlinearity while assuming a small angle $\theta _\lambda$ is to take the angle to be small but still ordered as unity in the RMHD ordering – in other words, it cannot scale with $\epsilon$ under the RMHD rescaling symmetry (3.5). The same rescaling symmetry implies that any physical scaling that involves $v_{\mathrm {A}}$ and $l_{\parallel }$ (and no other scales) must involve them in the combination $l_{\parallel }/v_{\mathrm {A}}$ (see § 5.2 and Beresnyak Reference Beresnyak2012a), which (6.9) manifestly does not. All this flies in the face of the fact that a substantial body of numerical evidence supporting aligned MHD turbulence was obtained by means of RMHD simulations (Mason et al. Reference Mason, Perez, Cattaneo and Boldyrev2011, Reference Mason, Perez, Boldyrev and Cattaneo2012; Perez et al. Reference Perez, Mason, Boldyrev and Cattaneo2012; Beresnyak Reference Beresnyak2012a; Mallet et al. Reference Mallet, Schekochihin and Chandran2015, Reference Mallet, Schekochihin, Chandran, Chen, Horbury, Wicks and Greenan2016) – complemented by explicit evidence that full MHD simulations produce quantitatively the same alignment – so the standard recourse to casting a cloud of suspicion on the validity of an asymptotic approximation is not available in this case.

In a further blow to the conjecture (6.5), it turns out that the alignment angle between the Elsasser fields at any given scale is anticorrelated with their amplitudes (Mallet et al. Reference Mallet, Schekochihin and Chandran2015), supporting the view that the dynamical alignment is indeed dynamical, being brought about by the mutual shearing of the Elsasser fields (Chandran et al. Reference Chandran, Schekochihin and Mallet2015), rather than by the uncertainty principle (6.5) (which would imply, presumably, a positive correlation between $\theta _\lambda$ and $\delta Z_\lambda$).

The numerical evidence in favour of alignment appears to be real and solid.Footnote ²¹ While numerical simulations at currently feasible resolutions cannot definitively verify or falsify Beresnyak's expectation that alignment is but a transient feature that disappears at asymptotically small scales, they certainly show aligned, locally 3D-anisotropic turbulence over a respectable inertial subrange at least one order of magnitude wide, and probably two. This is approaching the kind of scale separations that actually exist in nature, e.g., in the solar wind (where the evidence for a $k_{\perp }^{-3/2}$ spectrum has also been firming up: see, e.g., Chen et al. Reference Chen, Bale, Bonnell, Borovikov, Bowen, Burgess, Case, Chandran, Dudok de Wit and Goetz2020), and we cannot be casually dismissive of a physical regime, even if transient, that occupies most of the phase space that we are able to probe!

6.4. Revised model of aligned MHD turbulence

6.4.1. Dimensional and RMHD-symmetry constraints

Let me make the restrictions implied by Beresnyak's objection more explicit. Under the RMHD rescaling symmetry (3.5),

(6.12)\begin{equation} \delta Z_\lambda \to \epsilon\delta Z_\lambda,\quad \varepsilon \to \frac{\epsilon^3}{a}\varepsilon,\quad v_{\mathrm{A}} \to v_{\mathrm{A}},\quad \lambda\to a \lambda \end{equation}

(note that $\varepsilon$ is the energy flux whereas $\epsilon$ is the scaling factor). Therefore, the scaling relation (6.8) becomes $\epsilon \delta Z_\lambda \sim \epsilon ^{3/4}(\varepsilon v _{\mathrm {A}}\lambda )^{1/4}$, which is obviously a contradiction. Indeed, trialling

(6.13)\begin{equation} \delta Z_\lambda \sim \varepsilon^\mu v_{\mathrm{A}}^\nu\lambda^\gamma \end{equation}

and mandating both the symmetry (6.12) and dimensional consistency, one finds that the GS95 solution (5.4), $\nu =0$ and $\gamma =\mu =1/3$, is the only possibility, which was Beresnyak's point.

It seems obvious that the only way to rescue alignment is to allow another parameter – and the (almost) obvious choice is $L_{\parallel }$, the parallel outer scale, which transforms as $L_{\parallel }\to (a/\epsilon )L_{\parallel }$. Then

(6.14)\begin{equation} \delta Z_\lambda \sim \varepsilon^\mu v_{\mathrm{A}}^\nu\lambda^\gamma L_{{\parallel}}^\delta = \varepsilon^{(1+\delta)/3}\left(\frac{L_{{\parallel}}}{v_{\mathrm{A}}}\right)^\delta\lambda^{(1-2\delta)/3}, \end{equation}

where the second expression is the result of imposing on the first the RMHD symmetry (6.12) and dimensional correctness; $\delta =0$ returns us to GS95.Footnote ²² The same argument applied to the scaling of $l_{\parallel }$ with $\varepsilon$, $v_{\mathrm {A}}$, $\lambda$ and $L_{\parallel }$ gives

(6.15)\begin{equation} l_{{\parallel}} \sim \varepsilon^{(\sigma-1)/3}v_{\mathrm{A}}^{1-\sigma}L_{{\parallel}}^\sigma\lambda^{2(1-\sigma)/3}, \end{equation}

where $\sigma$ is a free parameter. Note that both (6.14) and (6.15) manifestly contain the parallel scales and $v_{\mathrm {A}}$ in the solely allowed combinations $l_{\parallel }/v_{\mathrm {A}}$ and $L_{\parallel }/v_{\mathrm {A}}$. A reassuring consistency check is to ask what perpendicular scale $\lambda =L_{\perp }$ corresponds to $l_{\parallel }=L_{\parallel }$: this turns out to be

(6.16)\begin{equation} L_{{\perp}} \sim \varepsilon^{1/2}\left(\frac{L_{{\parallel}}}{v_{\mathrm{A}}}\right)^{3/2} \sim \lambda_{\mathrm{CB}}, \end{equation}

the very same $\lambda _{\mathrm {CB}}$, given by (4.3), at which weak turbulence becomes strong – thus seamlessly connecting any strong-turbulence theory expressed by (6.14) and (6.15) with the WT cascade discussed in § 4.Footnote ²³ Notably, if we applied such a test to (6.9), we would find the price of consistency to be $L_{\perp }=L_{\parallel }$, which is allowed but does not have to be the case in MHD and certainly cannot be the case in RMHD.

Finally, since the parallel-cascade scaling (5.1) remains beyond reasonable doubt and, as can be readily checked, respects the rescaling symmetry (3.5) (Beresnyak Reference Beresnyak2015), combining it with (6.15) and (6.14) fixes

(6.17)\begin{equation} \sigma = 2\delta. \end{equation}

Alas, CB does not help with determining $\delta$ because, in aligned turbulence, the nonlinear time (6.4) contains the unknown scale $\xi$, or, equivalently, the alignment angle $\theta _\lambda \sim \lambda /\xi$. If we did know $\delta$, CB would let us determine this angle:

(6.18)\begin{equation} \frac{l_{{\parallel}}}{v_{\mathrm{A}}} \sim \tau_{\mathrm{nl}} \sim \frac{\lambda}{\delta Z_\lambda\sin\theta_\lambda}\quad\Rightarrow\quad \sin\theta_\lambda \sim \left(\frac{\lambda}{\lambda_{\mathrm{CB}}}\right)^{2\delta}, \end{equation}

where $\lambda _{\mathrm {CB}}$ is given by (6.16). The answer that we want to get – keeping Boldyrev's scalings of everything with $\lambda$ but not with $\varepsilon$ or $v_{\mathrm {A}}$ – requires

(6.19)\begin{equation} \delta = \frac{1}{8}. \end{equation}

Then, instead of (6.8), we end up with

(6.20)\begin{equation} \delta Z_\lambda \sim \varepsilon^{3/8}\left(\frac{L_{{\parallel}}}{v_{\mathrm{A}}}\right)^{1/8}\lambda^{1/4},\quad l_{{\parallel}} \sim \varepsilon^{{-}1/4}v_{\mathrm{A}}^{3/4}L_{{\parallel}}^{1/4}\lambda^{1/2},\quad \sin\theta_\lambda \sim \varepsilon^{{-}1/8}\left(\frac{v_{\mathrm{A}}}{L_{{\parallel}}}\right)^{3/8}\lambda^{1/4}, \end{equation}

and the dissipation cutoff scale (6.10) is corrected as follows:

(6.21)\begin{equation} \tau_{\mathrm{nl}} \sim \left(\frac{L_{{\parallel}}}{\varepsilon v_{\mathrm{A}}}\right)^{1/4}\lambda^{1/2} \ll \tau_{\eta} \sim \frac{\lambda^2}{\eta}\quad\Leftrightarrow\quad \lambda \gg \eta^{2/3}\left(\frac{L_{{\parallel}}}{\varepsilon v_{\mathrm{A}}}\right)^{1/6}\equiv\lambda_{\eta}. \end{equation}

Note that, since $\lambda _{\eta }\propto \eta ^{2/3}$ still, this does not address Beresnyak's numerical evidence on the convergence of the spectra (§ 6.3) – I shall come back to this problem in § 7.

For future convenience, let me recast the scalings (6.20) in a somewhat simpler form:

(6.22)\begin{equation} \delta Z_\lambda \sim \left(\frac{\varepsilon L_{{\parallel}}}{v_{\mathrm{A}}}\right)^{1/2}\left(\frac{\lambda}{\lambda_{\mathrm{CB}}}\right)^{1/4},\quad \frac{l_{{\parallel}}}{L_{{\parallel}}} \sim \left(\frac{\lambda}{\lambda_{\mathrm{CB}}}\right)^{1/2},\quad \sin\theta_\lambda \sim \left(\frac{\lambda}{\lambda_{\mathrm{CB}}}\right)^{1/4}. \end{equation}

Defining the magnetic Reynolds number based on the CB scale (6.16) and the fluctuation amplitude at this scale,

(6.23)\begin{equation} \mathrm{Rm} = \frac{\delta Z_{\lambda_{\mathrm{CB}}}\lambda_{\mathrm{CB}}}{\eta} \sim \frac{\varepsilon}{\eta}\left(\frac{L_{{\parallel}}}{v_{\mathrm{A}}}\right)^2, \end{equation}

allows the dissipation scale (6.21) to be recast as follows:

(6.24)\begin{equation} \frac{\lambda_{\eta}}{\lambda_{\mathrm{CB}}} \sim \left(\frac{\mathrm{Rm}}{1+\mathrm{Pm}}\right)^{{-}2/3} = \widetilde{\mathrm{Re}}^{{-}2/3},\quad \mathrm{Pm} = \frac{\nu}{\eta},\quad \widetilde{\mathrm{Re}} = \frac{\delta Z_{\lambda_{\mathrm{CB}}}\lambda_{\mathrm{CB}}}{\nu+\eta}. \end{equation}

I have restored the possibility that viscosity $\nu$ might be larger than the magnetic diffusivity $\eta$: if that is the case, one must replace the latter with the former in the calculation of the dissipative cutoff, whereas if $\mathrm {Pm}\lesssim 1$, it does not matter, hence the appearance of the magnetic Prandtl number $\mathrm {Pm}$ in the combination $(1+\mathrm {Pm})$.

Yet another way to write the first of the scaling relations (6.20) is

(6.25)\begin{equation} \delta Z_\lambda \sim \varepsilon^{1/3}\lambda_{\mathrm{CB}}^{1/12}\lambda^{1/4} \quad\Leftrightarrow\quad E(k_{{\perp}})\sim \varepsilon^{2/3}\lambda_{\mathrm{CB}}^{1/6}k_{{\perp}}^{{-}3/2}. \end{equation}

This is effectively the prediction for the spectrum that Perez et al. (Reference Perez, Mason, Boldyrev and Cattaneo2012, Reference Perez, Mason, Boldyrev and Cattaneo2014b) used in their numerical convergence studies. Thus, they and I are on the same page as to what the spectrum of aligned turbulence is expected to be, although the question remains why it should be that if Boldyrev's uncertainty principle (6.5) can no longer be used.

A set of RMHD-compatible scalings (6.20), or (6.25), is also effectively what was deduced by Chandran et al. (Reference Chandran, Schekochihin and Mallet2015) and by Mallet & Schekochihin (Reference Mallet and Schekochihin2017) from a set of plausible conjectures about the dynamics and statistics of RMHD turbulence (they did not explicitly discuss the issue of the RMHD rescaling symmetry, but used normalisations that enforced it automatically). The two papers differed in their strategy for determining the exponent $\delta$; my exposition here will be a ‘heuristic’ version of Mallet & Schekochihin (Reference Mallet and Schekochihin2017).

6.4.2. Intermittency matters!

The premise of both Chandran et al. (Reference Chandran, Schekochihin and Mallet2015) and Mallet & Schekochihin (Reference Mallet and Schekochihin2017) is that in order to determine the scalings of everything, including the energy spectrum, one must have a working model of intermittency, i.e., of the way in which fluctuation amplitudes and their scale lengths in all three directions – $\lambda$, $\xi$ and $l_{\parallel }$ – are distributed in a turbulent MHD system. It may be disturbing to the reader, or viewed by her as an unnecessary complication, that we must involve ‘rare’ events – as this is what the theory of intermittency is about – in the mundane business of the scaling of the energy spectra, which are usually viewed as made up from the more ‘typical’ fluctuations. These doubts might be alleviated by the following observation. The appearance of the outer scale $L_{\parallel }$ in (6.14) suggests that the self-similarity is broken – this is somewhat analogous to what happens in hydrodynamic turbulence, where corrections to the K41 scaling (2.2) come in as powers of $\lambda /L$ (Kolmogorov Reference Kolmogorov1962; Frisch Reference Frisch1995). We may view $\delta$ as just such a correction to the self-similar GS95 result, and alignment as the physical mechanism whereby this intermittency correction arises. The main difference with the hydrodynamic case is that $\delta$ is not all that small (MHD turbulence is ‘more intermittent’ than the hydrodynamic one), the mechanism responsible for it has important consequences (§ 7), and so we care.

I shall forgo a detailed discussion of the intermittency model that Mallet & Schekochihin (Reference Mallet and Schekochihin2017) proposed; for my purposes here, a vulgarised version of their argument will suffice. They consider the turbulent field as an ensemble of structures, or fluctuations, each of which has some amplitude and three scales: parallel $l_{\parallel }$, perpendicular $\lambda$ and fluctuation-direction $\xi$ (they call this the ‘RMHD ensemble’). They make certain conjectures about the joint probability distribution of these quantities, which then allow them to fix scalings. The most crucial (and perhaps also the most arbitrary) of these conjectures is, effectively, that for all structures, $l_{\parallel }\sim \lambda ^\alpha$ with the same exponent $\alpha$, i.e., that the quantity $l_{\parallel }/\lambda ^\alpha$ has a scale-invariant distribution (this appears to be confirmed by numerical evidence: see figure 11a). They then determine the exponent $\alpha$ by considering ‘the most intense structures’Footnote ²⁴ – because it is possible to work out what the probability of encountering them is as a function both of $\lambda$ and of $l_{\parallel }$.

Figure 11. (a) Probability distribution of $l_{\parallel }/\lambda ^{1/2}$ in a $1024^3$ RMHD simulation (the shades of colour from blue to red correspond to PDFs at increasing scales within the inertial range). This plot is taken from Mallet & Schekochihin (Reference Mallet and Schekochihin2017) (where the reader will find a discussion – somewhat inconclusive – of the slope of this PDF) and illustrates how good (or otherwise) is the assumption that $l_{\parallel }/\lambda ^\alpha$ has a scale-invariant distribution (the assumption is not as good as RCB, illustrated in figure 6 based on data from the same simulation). (b) Joint probability distribution for the length $l_{\parallel }$ and width $\xi$ (in my notation) of the most intense dissipative structures (adapted from Zhdankin, Boldyrev & Uzdensky Reference Zhdankin, Boldyrev and Uzdensky2016b). This shows that $\xi \propto l _{\parallel }$, in line with (6.29). (Reprinted from Zhdankin et al. Reference Zhdankin, Boldyrev and Uzdensky2016b with the permission of AIP Publishing.) Independent simulations by J. M. Stone (private communication, 2018) support this scaling.

They then conjecture that the most intense structures in the RMHD ensemble are sheets transverse to the local perpendicular direction. Therefore, if one looks for their probability (filling fraction) in any perpendicular plane as a function of the perpendicular scale $\lambda$, one expects it to scale as

(6.26)\begin{equation} P \propto \lambda. \end{equation}

If, on the other hand, one is interested in their filling fraction in the plane locally tangent to a flux sheet (i.e., defined by the local mean field and the direction of the fluctuation vector), it is

(6.27)\begin{equation} P \propto \xi l_{{\parallel}}. \end{equation}

The next conjecture is the ‘refined critical balance’ (RCB, already advertised in § 5.1), stating that not only is $\tau _{\mathrm {nl}}\sim \tau _{\mathrm {A}}$ in some vague ‘typical’ sense, but the quantity

(6.28)\begin{equation} \chi=\frac{\delta Zl_{{\parallel}}}{\xi v_{\mathrm{A}}} \sim \frac{\tau_{\mathrm{A}}}{\tau_{\mathrm{nl}}} \end{equation}

has a scale-invariant distribution in the RMHD ensemble – this was discovered by Mallet et al. (Reference Mallet, Schekochihin and Chandran2015) to be satisfied with truly remarkable accuracy in numerically simulated RMHD turbulence (figure 6).Footnote ²⁵ If this is true for all structures, it is true for the most intense ones – and a further assumption about those is that their amplitude $\delta Z_{\mathrm {max}}$ is not a function of scale but is simply equal to some typical outer-scale value (i.e., the most intense structures are formed by the largest perturbations collapsing, or being sheared, into sheets without breaking up into smaller perturbations; see Chandran et al. Reference Chandran, Schekochihin and Mallet2015). This, together with (6.27), implies that for those structures,

(6.29)\begin{equation} \xi \sim l_{{\parallel}} \frac{\delta Z_{\mathrm{max}}}{v_{\mathrm{A}}} \quad\Rightarrow\quad P \propto l_{{\parallel}}^2 \end{equation}

(Zhdankin et al. Reference Zhdankin, Boldyrev and Uzdensky2016b confirm that $\xi \propto l _{\parallel }$ for the most intense dissipative structures: see figure 11b). Comparing (6.29) with (6.26), we conclude that $l_{\parallel }\propto \lambda ^{1/2}$ for the most intense structures and, therefore, for everyone else – by the conjecture of scale invariance of $l_{\parallel }/\lambda ^\alpha$, where we now know that $\alpha =1/2$. Comparing this with (6.15), we see that $\alpha = 2(1-\sigma )/3$, whence

(6.30)\begin{equation} \sigma = \frac{1}{4}\quad\Rightarrow\quad \delta=\frac{1}{8}, \end{equation}

the latter by virtue of (6.17). Q.e.d.: we now have the scalings (6.20).

I do not know if the reader will find this quasi-intuitive argument more (or less) convincing than the formal-looking conjectures and corollaries in Mallet & Schekochihin (Reference Mallet and Schekochihin2017). There is no need to repeat their algebra here, but hopefully the above sheds some (flickering) light – if not, perhaps a better argument will be invented one day, but all I can recommend for now is reading their paper. Notably, in their more formal treatment, not just the energy spectrum but the two-point structure functions of all orders are predicted – and turn out to be a decent fit to numerical data as it currently stands.Footnote ²⁶ The same is true about the model proposed in the earlier paper by Chandran et al. (Reference Chandran, Schekochihin and Mallet2015). Their approach is based on a much more enthusiastic engagement with dynamics: a careful analysis of how aligned and non-aligned structures might form and interact. They get $\delta \approx 0.108$, which leads to $\delta Z_\lambda \propto \lambda ^{0.26}$ – not a great deal of difference with (6.20), considering that all of this is very far from being exact science. Their approach does have the distinction, however, of emphasising particularly strongly the dynamic nature of the dynamic alignment, which arises as Elsasser fields shear each other into sheet-like structures.

6.5. 3D anisotropy

Before moving on, I would like to re-emphasise the 3D anisotropy of the aligned MHD turbulence – and the fact that this anisotropy is local, associated at every point with the three directions that themselves depend on the fluctuating fields: parallel to the magnetic field ($l_{\parallel }$), along the vector direction of the perturbed field $\boldsymbol {Z}_{\perp }^\mp$ that advects the field $\boldsymbol {Z}_{\perp }^\pm$ whose correlations we are measuring ($\xi$), and the third direction perpendicular to the other two ($\lambda$). This local 3D anisotropy is measurableFootnote ²⁷ and has indeed been observed both in the solar wind (Chen et al. Reference Chen, Mallet, Schekochihin, Horbury, Wicks and Bale2012; Verdini et al. Reference Verdini, Grappin, Alexandrova and Lion2018, Reference Verdini, Grappin, Alexandrova, Franci, Landi, Matteini and Papini2019) and in numerical simulations (Verdini & Grappin Reference Verdini and Grappin2015; Mallet et al. Reference Mallet, Schekochihin, Chandran, Chen, Horbury, Wicks and Greenan2016) – both are illustrated by figure 12. The main point of discrepancy between the true and virtual reality is the scale dependence of the anisotropy: confirmed solidly in simulations but only very tentatively in the solar wind (see footnote 21). However, progress never stops, and one can hope for better missions (Chen et al. Reference Chen, Bale, Bonnell, Borovikov, Bowen, Burgess, Case, Chandran, Dudok de Wit and Goetz2020) and even more sophisticated analysis.

Figure 12. Locally 3D-anisotropic structures in the (a) solar wind (from Verdini et al. Reference Verdini, Grappin, Alexandrova and Lion2018 © AAS, reproduced with permission) and (b) numerical simulations (here $l_{\parallel }$ is normalised to $L_{\parallel }/2\pi$ and $\lambda$ and $\xi$ to $L_{\perp }/2\pi$, hence apparent isotropy at the outer scale). These are surfaces of constant second-order structure function of the magnetic field (a) or one of the Elsasser fields (b). The three images correspond to successively smaller fluctuations and so successively smaller scales (only the last of the three is firmly in the universal inertial range). In both cases, the emergence of statistics with $l_{\parallel }\gg \xi \gg \lambda$ is manifest. In the solar wind, the route that turbulence takes to this aligned state appears to depend quite strongly on the solar-wind expansion, which distorts magnetic-field component in the radial direction compared to the azimuthal ones (Verdini & Grappin Reference Verdini and Grappin2015; Vech & Chen Reference Vech and Chen2016). The data shown in panel (a) was carefully selected to minimise this effect; without such selection, one sees structures most strongly elongated in the $\xi$ direction at the larger scales ($\xi >l_{\parallel }>\lambda$), although they too tend to the universal aligned regime at smaller scales (see Chen et al. Reference Chen, Mallet, Schekochihin, Horbury, Wicks and Bale2012, where 3D-anisotropy plots like ones shown here first appeared).

The scaling of the energy spectrum in the parallel direction (§ 5.2) was arguably the most robust and uncontroversial of the results reviewed thus far. We then occupied ourselves with the scalings of the Elsasser-field increments and of $l_{\parallel }$ vs. the perpendicular scale $\lambda$, culminating in § 6.4 with a theory that one (hopefully) can believe in. The scalings with the third, fluctuation-direction coordinate $\xi$ are very easy to obtain because the nonlinear time of the aligned cascade (6.4) has the same dependence on $\xi$ as it did on $\lambda$ in the unaligned, GS95 theory: see (5.3). Therefore,

(6.31)\begin{equation} \delta Z_\xi \sim (\varepsilon\xi)^{1/3},\quad \xi \sim \varepsilon^{1/8}\left(\frac{L_{{\parallel}}}{v_{\mathrm{A}}}\right)^{3/8}\lambda^{3/4} \sim \lambda_{\mathrm{CB}}^{1/4}\lambda^{3/4}, \end{equation}

with the latter formula following from (6.3) and (6.20) or (6.22). Thus, Elsasser fields’ spectra have exponents $-2$ in the $l_{\parallel }$ direction, $-3/2$ in the $\lambda$ direction and $-5/3$ in the $\xi$ direction (see the $n=2$ exponents in figure 13a). Let me note in passing that the ‘Kolmogorov’ scaling (6.31) will play a key part in my discussion, in Appendix D.7, of why the Lazarian & Vishniac (Reference Lazarian and Vishniac1999) notion of ‘stochastic reconnection’does not automatically invalidate the aligned cascade and return us to GS95, as an erudite reader might have been worried it would – the idea is that the Richardson (Reference Richardson1926) (super)diffusion of Lagrangian trajectories in a turbulent field is always determined by the $\xi$-dependent scaling (6.31), regardless of the nature of the cascade in $\lambda$.

Figure 13. Scaling exponents of the structure functions in RMHD turbulence simulated by Mallet et al. (Reference Mallet, Schekochihin, Chandran, Chen, Horbury, Wicks and Greenan2016) (the plot is reproduced from Mallet & Schekochihin Reference Mallet and Schekochihin2017). (a) Structure functions of the Elsasser-field increments (5.6): by definition, $\langle |\delta \boldsymbol {Z}_{\boldsymbol {l}}^+|^n\rangle \propto l^{\zeta _n}$ and $\zeta _n^\perp$, $\zeta _n^\mathrm {fluc}$, $\zeta _n^\parallel$ are exponents for $l=\lambda$, $\xi$, $l_{\parallel }$, respectively (i.e., all structure functions are conditional on the point separation being in one of the three directions of local 3D anisotropy; see §§ 5.3 and 6.5). The solid lines are for a $1024^3$ simulation (with hyperviscosity), the dashed ones are for a $512^3$ simulation, indicating how converged, or otherwise, the exponents are, and the dotted lines, in both (a) and (b), are the theoretical model by Mallet & Schekochihin (Reference Mallet and Schekochihin2017). (b) Similarly defined structure functions of the velocity (solid lines) and magnetic-field (dashed lines) increments from the same $1024^3$ simulation. The magnetic field is ‘more intermittent’ than the Elsasser fields and the latter more so than velocity. An early (possibly first) numerical measurement of this kind, highlighting the differences between scalings of different fields and in different local directions, was done by Cho, Lazarian & Vishniac (Reference Cho, Lazarian and Vishniac2003).

A good way of thinking physically of the inevitability of 3D anisotropy is to note that, from (6.4) and CB,

(6.32)\begin{equation} \xi \sim l_{{\parallel}}\frac{\delta Z_\lambda}{v_{\mathrm{A}}} \sim l_{{\parallel}} \frac{\delta b_\lambda}{v_{\mathrm{A}}}, \end{equation}

i.e., $\xi$ is the typical displacement of a fluid element and also the typical perpendicular distance a field line wanders within a structure coherent on the parallel scale $l_{\parallel }$. Fluctuations must therefore preserve coherence in their own direction at least on the scale $\xi$. They are not constrained in this way in the third direction $\lambda$, and the fluctuation direction itself has an angular uncertainty of the order of the angle $\theta _\lambda$ between the two fields, so it makes sense that the aspect ratio of the structures in the perpendicular plane should satisfy (6.3).

The dependence of the anisotropy on the local direction of the fluctuating fields makes the connection between anisotropy, alignment and intermittency more obvious: when we follow perturbed field lines to extract parallel correlations or measure one Elsasser field's decorrelation along the direction of another Elsasser field, we are clearly not calculating second-order statistics in the strict sense – and so, in formal terms, local scale-dependent anisotropy always involves correlation functions of (all) higher orders.Footnote ²⁸ Thus, it makes a certain natural sense to speak of the alignment-induced departure of MHD-turbulence spectrum from the Kolmogorovian GS95 scaling and of the 3D anisotropy of the underlying fluctuation field as an intermittency effect, as I have done here.

6.6. Higher-order statistics

In several places (e.g., in §§ 5.3 and 6.4.2), I have brushed against the more formal task of the intermittency theory – the calculation of the scaling exponents of higher-order structure functions or, equivalently, of the probability distributions of field increments – and recoiled every time, opting for ‘twiddle’ algebra and statements about spectra. A fair amount of information on these matters is available from simulations and from the solar-wind measurements: what intermittency looks like in the former is illustrated by figure 13 (a survey of previous measurements of structure functions, both in simulations and in the solar wind, can be found in Chandran et al. Reference Chandran, Schekochihin and Mallet2015). Some of what is known is perhaps understood, but much remains a mystery: for example, we do not know why the higher-order scaling exponents are generally not the same for velocity, Elsasser and magnetic fields, with the latter ‘more intermittent’ than the former, as is evident in figure 13(b) (in § 10.4, I will moot a possible connexion between intermittency and negative ‘residual energy’ – an asymmetry between the magnetic and velocity spectra seen both in numerical simulations and in the solar wind – but I do not know how to translate this into anything resembling figure 13b).

Interesting as it is, I will now leave the problem of higher-order statistics alone. We know from the (ongoing) history of hydrodynamic-turbulence theory that once this becomes the unsolved problem that everyone is working on, the scope for abstract theorising expands to fill all available space (and time) while attention paid by the outside world diminishes.Footnote ²⁹ This said, I hasten to dispel any possible impression that I do not consider intermittency of MHD turbulence an important problem: in fact, as I have argued above, intermittency as a physical phenomenon appears to be so inextricably hard-wired into the structure of MHD turbulence that any workable theory of the latter has to be a theory of its intermittency.

Finally, let me jump ahead of myself and mention also that we know nothing at all of the intermittency in ‘tearing-dominated turbulence’, which is about to be introduced (§ 7), and very little of the intermittency in the variants of MHD turbulence surveyed in §§ 9–13. In particular, the relationship between intermittency and Elsasser imbalance, local or global, appears to me to be a promising object for theoreticians’ scrutiny (see § 9.1).

7.1. Disruption by tearing

The scale at which the aligned structures will be disrupted by tearing can be estimated very easily by comparing the nonlinear time (6.4) of the aligned cascade with the growth time of the fastest tearing mode that can be triggered in an MHD sheet of a given transverse scale $\lambda$. That this growth time is a good estimate for the time that reconnection needs to break up a sheet forming as a result of ideal-MHD dynamics is not quite as obvious as it might appear, but it is true and was carefully shown to be so by Uzdensky & Loureiro (Reference Uzdensky and Loureiro2016). The maximum tearing growth rate is

(7.1)\begin{equation} \gamma \sim \frac{v_{\mathrm{A}y}}{\lambda}\, S_\lambda^{{-}1/2} (1+\mathrm{Pm})^{{-}1/4},\quad S_\lambda = \frac{v_{\mathrm{A}y}\lambda}{\eta},\quad \mathrm{Pm} = \frac{\nu}{\eta}. \end{equation}

How to derive this is reviewed in Appendix D.1 (see (D32)). Here $v_{\mathrm {A}y}$ is the Alfvén speed associated with the perturbed magnetic field that reverses at scale $\lambda$, $S_\lambda$ is the corresponding Lundquist number and $\mathrm {Pm}$ is the magnetic Prandtl number, which only matters if the viscosity $\nu$ is larger than the magnetic diffusivity $\eta$. In application to our aligned Alfvénic structures, we should replace $v_{\mathrm {A}y}\sim \delta Z_\lambda$. Then, using the scalings (6.22) to work out $\tau _{\mathrm {nl}}$, we find that the aligned cascade is faster than tearing as long as

(7.2)\begin{equation} \gamma \tau_{\mathrm{nl}} \sim \frac{S_\lambda^{{-}1/2} (1+\mathrm{Pm})^{{-}1/4}}{\sin\theta_\lambda} \ll 1 \quad\Leftrightarrow\quad \lambda \gg \mathrm{Rm}^{{-}4/7}(1+\mathrm{Pm})^{{-}2/7}\lambda_{\mathrm{CB}} \equiv \lambda_{\mathrm{D}}, \end{equation}

where $\mathrm {Rm} \sim S_{\lambda _{\mathrm {CB}}}$, as defined in (6.23), and the scale $\lambda _{\mathrm {D}}$ will henceforth be referred to as the disruption scale. At scales $\lambda \lesssim \lambda _{\mathrm {D}}$, aligned sheet-like structures can no longer retain their integrity against the onslaught of tearing.Footnote ³²

7.1.1. Some reservations and limitations

Let me observe parenthetically that one can entertain legitimate reservations about the validity of (7.1) and other formulae based on laminar-tearing theory in a noisy, turbulent environment (see, e.g., the discussion and references in footnote 103). Here suffice it to say, that at scales where $\gamma \ll \tau _{\mathrm {nl}}^{-1}$, the laminar formulae are, of course, unjustified, and at scales where $\gamma \gg \tau _{\mathrm {nl}}^{-1}$, if such situations existed, they would be perfectly fine, modulo the issue of flows, which I will explain in a moment.Footnote ³³ The disruption scale $\lambda _{\mathrm {D}}$ is where $\gamma \sim \tau _{\mathrm {nl}}^{-1}$, so the laminar-tearing formulae are presumably only order-unity wrong, which I acknowledge by using ‘$\sim$’ instead of ‘$=$’ everywhere.

Another legitimate reservation concerns application to aligned structures formed by the MHD cascade of the tearing theory derived for a magnetostatic equilibrium. Indeed, these structures are not purely magnetic sheets, but rather Elsasser ones, i.e., they have local shear flows superimposed on them. These flows are likely somewhat sub-Alfvénic (the reasons and evidence for that are explained in § 10.4), but, technically speaking, it has not been shown that they are sub-Alfvénic enough to justify application of the no-flow tearing-mode theory or even that they cannot fatally stabilise tearing (some further discussion and references on tearing with flows can be found at the end of Appendix D.4.2). Since the shear rate in these flows is no larger than $\tau _{\mathrm {nl}}^{-1}$, I think that this difficulty is a technical one, rather than a deal breaker, i.e., that what I am doing here is again only order-unity wrong (cf. Tolman, Loureiro & Uzdensky Reference Tolman, Loureiro and Uzdensky2018), but I cannot prove it rigorously. If this admission has not discouraged the reader, let us proceed.

The disruption scale $\lambda _{\mathrm {D}}$, upon comparison with the putative resistive cutoff (6.24) of the aligned cascade turns out to supersede it provided $\mathrm {Pm}$ is not too large:

(7.3)\begin{equation} \frac{\lambda_{\mathrm{D}}}{\lambda_{\eta}} \sim \left[\frac{\mathrm{Rm}}{(1+\mathrm{Pm})^{10}}\right]^{2/21}\gg1. \end{equation}

In view of the ridiculous exponents involved, this means that in a system with even moderately large $\mathrm {Pm}$ and/or not a truly huge $\mathrm {Rm}$, the aligned cascade will happily make it to the dissipation cutoff (6.24) and no further chapters are necessary in this story.Footnote ³⁴ However, I do want to tell the story in full and so will focus on situations in which the condition (7.3) is satisfied.

7.1.2. Debris of disruption

I shall turn to the question of what happens at scales below $\lambda _{\mathrm {D}}$ in § 7.2, but to do that, it is necessary first to ask what becomes of the aligned structures that are disrupted at $\lambda _{\mathrm {D}}$.

The tearing instability that disrupts them, the so-called Coppi mode, or (the fastest-growing) resistive internal kink mode (Coppi et al. Reference Coppi, Galvão, Pellat, Rosenbluth and Rutherford1976), has the wavenumber (see (D32))

(7.4)\begin{equation} k_* \sim \frac{1}{\lambda}\,S_{\lambda}^{{-}1/4}(1+\mathrm{Pm})^{1/8} \sim\frac{1}{\lambda_{\mathrm{CB}}} \mathrm{Rm}^{{-}1/4}(1+\mathrm{Pm})^{1/8}\left(\frac{\lambda}{\lambda_{\mathrm{CB}}}\right)^{{-}21/16}, \end{equation}

where (6.22) has been used to substitute for $\delta Z_\lambda$ inside $S_\lambda$. Therefore, at the disruption scale ($\lambda =\lambda _{\mathrm {D}}$),

(7.5)\begin{equation} k_*\sim\frac{1}{\lambda_{\mathrm{CB}}} \mathrm{Rm}^{1/2}(1+\mathrm{Pm})^{1/2}. \end{equation}

If referred to the length of the sheet $\xi _{\mathrm {D}}$, which depends on $\lambda _{\mathrm {D}}$ via (6.31), this wavenumber gives us an estimate for the number of islands in the growing perturbation:

(7.6)\begin{equation} N \sim k_*\xi_{\mathrm{D}} \sim \mathrm{Rm}^{1/14}(1+\mathrm{Pm})^{2/7}. \end{equation}

As this is always large (or at least $\gtrsim 1$), the mode fits comfortably into the sheet that it is trying to disrupt.Footnote ³⁵

What happens to these islands? When the tearing mode enters the nonlinear regime, the island width is (see Appendix D.2)

(7.7)\begin{equation} w \sim k_*\lambda_{\mathrm{D}}^2, \end{equation}

which is smaller than $\lambda _{\mathrm {D}}$ and so, technically speaking, the aligned structures need not be destroyed by these islands. Uzdensky & Loureiro (Reference Uzdensky and Loureiro2016) (followed by Mallet et al. Reference Mallet, Schekochihin and Chandran2017b and by Loureiro & Boldyrev Reference Loureiro and Boldyrev2017b) argue that, after the tearing mode goes nonlinear, the islands will grow further while the $X$-points between them collapse into current sheets – all of that on the same time scale (7.1) as the mode grew. It seems intuitive that, in order to break up the aligned structure (ideal-MHD sheet) that spawned them, the islands would have to get to $w\sim \lambda _{\mathrm {D}}$ (Uzdensky & Loureiro Reference Uzdensky and Loureiro2016). If they did this in a cross-section-area-preserving way, then they would be circular at the point of disruption: from (7.7), $w k_*^{-1}\sim \lambda _{\mathrm {D}}^2$. Such a set of islands (flux ropes) would be isotropic in the perpendicular plane. Another simplistic model is to imagine their width grow due to further reconnection while their length remains the same, viz., $\sim k_*^{-1}$ (simplistic because the lengthening of the inter-island current sheets is ignored). The aspect ratio of such a flux rope once it reached the width $\lambda _{\mathrm {D}}$ would be

(7.8)\begin{equation} \lambda_{\mathrm{D}} k_* \sim N \frac{\lambda_{\mathrm{D}}}{\xi_{\mathrm{D}}} \sim N \sin\theta_{\lambda_{\mathrm{D}}} \sim \mathrm{Rm}^{{-}1/14}(1+\mathrm{Pm})^{3/14}, \end{equation}

where $N$ is given by (7.6). This is a degree of alignment preserved, but reduced by a (asymptotically large) factor of $N$. At the extreme end of the range of possibilities is the view of Boldyrev & Loureiro (Reference Boldyrev and Loureiro2017) that, in fact, islands of size (7.7) are already sufficient to break up the aligned structure, so alignment is not abruptly reduced by tearing disruption. Both this idea and some nuances about the $X$-point collapse are discussed further in § 7.4.1. The key point for us here is that at the disruption scale, the aligned structures that cascade down from the inertial range are broken up by reconnection into flux ropes that are underaligned and may even be isotropic in the perpendicular plane. They are a starting point for a new kind of cascade, which I shall now proceed to consider.

7.2. Tearing-mediated turbulence

If you accept the argument that the disruption by tearing of an aligned structure at the scale $\lambda _{\mathrm {D}}$ leads to its break-up into a number of flux ropes of width $\sim \lambda _{\mathrm {D}}$, then the natural conclusion is that $\lambda _{\mathrm {D}}$ now becomes a kind of ‘outer scale’ for a new cascade. There need not be anything particularly different about this cascade compared to the standard aligned cascade except the alignment angle may now be reset to a greater value. As the disruption-scale structures interact with each other and break up into smaller structures, the latter should develop the same tendency to align as their inertial-range predecessors did. For a while, the structures in this new cascade are safe from tearing as their aspect ratio is not large enough, but eventually (i.e., at small enough scales), they too will become sufficiently aligned to be broken up by tearing modes. This leads to another disruption, another iteration of an aligned ‘mini-cascade’, and so on. Thus, if we rename our critical-balance scale $\lambda _{\mathrm {CB}}=\lambda _0$, the disruption scale $\lambda _{\mathrm {D}}=\lambda _1$, and call the subsequent disruption scales $\lambda _n$, we can think of the MHD cascade as consisting of a sequence of aligned cascades interrupted by disruption episodes. I shall calculate $\lambda _n$ and the spectral exponents for this new, composite tearing-mediated cascade in § 7.2.1, the final dissipation cutoff in § 7.2.2, the cascade's alignment properties in § 7.2.3, and its parallel scalings, again determined by CB, in § 7.2.4.

7.2.1. Mini-cascades and the spectrum of tearing-mediated turbulence

Let us calculate the disruption scales $\lambda _n$, following Mallet et al. (Reference Mallet, Schekochihin and Chandran2017b). Since the ‘mini-cascades’ that connect them are just the same as the aligned cascade whose disruption we analysed in § 7.1, we can use (7.2) to deduce a recursion relation

(7.9)\begin{equation} \lambda_{n+1} \sim S_{\lambda_n}^{{-}4/7}(1+\mathrm{Pm})^{{-}2/7}\lambda_n \end{equation}

(remembering that $\mathrm {Rm}$ was defined as the Lundquist number at scale $\lambda _{\mathrm {CB}}=\lambda _0$). To work out the Lundquist number $S_{\lambda _n}$ at scale $\lambda _n$, notice that, if the alignment angle $\theta _{\lambda }$ just below $\lambda _n$ is increased, there must be a downward jump in the amplitude of the turbulent fluctuations at $\lambda _n$: indeed, the nonlinear time (6.4) shortens compared to what it was in the aligned cascade just above $\lambda _n$, and the cascade accelerates. I shall consider the extreme possibility that the alignment is reset to being order unity, as this will effectively bracket the range of outcomes for the scalings of this cascade. Since it still has to carry the same energy flux, we have, for amplitudes just below the disruption scale ($\lambda _n^-$),

(7.10)\begin{equation} \frac{\bigl(\delta Z_{\lambda_n^-}\bigr)^3}{\lambda_n}\sim\varepsilon \quad\Rightarrow\quad \delta Z_{\lambda_n^-} \sim (\varepsilon\lambda_n)^{1/3}, \end{equation}

just a Kolmogorov (or GS95) scaling. Therefore, the Lundquist number of the $n$-th mini-cascade is

(7.11)\begin{equation} S_{\lambda_n}\sim \frac{\delta Z_{\lambda_n^-}\lambda_n}{\eta} \sim \frac{\varepsilon^{1/3}\lambda_n^{4/3}}{\eta} \sim \mathrm{Rm}\left(\frac{\lambda_n}{\lambda_{\mathrm{CB}}}\right)^{4/3}. \end{equation}

In combination with (7.9), this gives us

(7.12)\begin{equation} \frac{\lambda_n}{\lambda_{\mathrm{CB}}} \sim \left[\mathrm{Rm}^{{-}4/7}(1+\mathrm{Pm})^{{-}2/7}\right]^{\frac{21}{16}\left[1-\left(\frac{5}{21}\right)^n\right]} \to \mathrm{Rm}^{{-}3/4}(1+\mathrm{Pm})^{{-}3/8},\quad n\to\infty. \end{equation}

I shall return to this obviously suggestive (Kolmogorov!) scaling in § 7.2.2.

In the picture that I have just painted, the cascade in the tearing-mediated range looks like a ladder (figure 14), with amplitude dropping at each successive disruption scale as structures become unaligned (or less aligned, in which case the the steps of the ladder are less tall). In between the disruption scales, there are aligned ‘mini-cascades’ of the same kind as the original one discussed in § 6.4, with $k_{\perp }^{-3/2}$ spectra. This means that the overall scaling of the turbulent fluctuation amplitudes can be constrained between their scaling just below each disruption scale ($\lambda _n^-$), given by (7.10), and just above it ($\lambda _n^+$). The latter scaling, for the amplitudes of the structures just before they get disrupted can be inferred from the fact that for these structures, the tearing growth rate (7.1) must be the same as the nonlinear interaction (cascade) rate: letting $v_{\mathrm {A}y}\sim \delta Z_{\lambda _n^+}$ in (7.1), we get

(7.13)\begin{equation} \tau_{\mathrm{nl}}^{{-}1}\sim\gamma \sim\bigl(\delta Z_{\lambda_n^+}\bigr)^{1/2}\lambda_n^{{-}3/2}\eta^{1/2}(1+\mathrm{Pm})^{{-}1/4} \end{equation}

and, therefore,

(7.14)\begin{align} \frac{\bigl(\delta Z_{\lambda_n^+}\bigr)^2}{\tau_{\mathrm{nl}}} \sim \varepsilon\quad\Rightarrow\quad \delta Z_{\lambda_n^+} &\sim \varepsilon^{2/5}\eta^{{-}1/5}(1+\mathrm{Pm})^{1/10}\lambda_n^{3/5} \nonumber\\ &\sim \left(\frac{\varepsilon L_{{\parallel}}}{v_{\mathrm{A}}}\right)^{1/2}\left(\frac{\lambda_{\mathrm{D}}}{\lambda_{\mathrm{CB}}}\right)^{1/4} \left(\frac{\lambda_n}{\lambda_{\mathrm{D}}}\right)^{3/5}. \end{align}

The last expression puts this result explicitly in contact with the inertial-range scaling (6.22). Thus, the tearing-mediated-range spectrum is (Mallet et al. Reference Mallet, Schekochihin and Chandran2017b)Footnote ³⁶

(7.15)\begin{equation} \varepsilon^{2/3}k_{{\perp}}^{{-}5/3} \lesssim E(k_{{\perp}}) \lesssim \varepsilon^{4/5}\eta^{{-}2/5}(1+\mathrm{Pm})^{1/5}k_{{\perp}}^{{-}11/5}. \end{equation}

Since the $-11/5$ upper envelope is steeper than the $-5/3$ lower one, the two converge and eventually meet at

(7.16)\begin{equation} \lambda_\infty \sim \eta^{3/4}\varepsilon^{{-}1/4}(1+\mathrm{Pm})^{{-}3/8} \sim \lambda_{\mathrm{D}}^{21/16}\lambda_{\mathrm{CB}}^{{-}5/16}, \end{equation}

which is, of course, the scale (7.12) in the limit $n\to \infty$. For $\mathrm {Pm}\gg 1$, this will, in fact, be superseded by the Kolmogorov cutoff to be derived in § 7.2.2.

Figure 14. Spectrum of MHD turbulence and transition to tearing-mediated cascade (see (7.15)) (adapted from Mallet et al. Reference Mallet, Schekochihin and Chandran2017b). The width of the tearing-mediated range is, of course, exaggerated in this cartoon. The spectral slopes of the ‘mini-cascades’ between $\lambda _n^{-1}$ and $\lambda _{n+1}^{-1}$ are all $k_{\perp }^{-3/2}$, but the overall envelope is $k_{\perp }^{-11/5}$. Note that, modulo $\mathrm {Pm}$ dependence, the disrupted aligned cascade and a putative unaligned GS95 $k_{\perp }^{-5/3}$ spectrum, starting at $\lambda _{\mathrm {CB}}$, terminate at the same, Kolmogorov, scale (7.19).

While in the above construction, the tearing-mediated-range spectrum is pictured as a succession of ‘steps’ representing the ‘mini-cascades’ that connect the successive disruption scales (figure 14), the reality will almost certainly look more like some overall power-law spectrum with a slope for which the upper $-11/5$ bound (7.15) seems to be a good estimate. Indeed, the tearing disruptions will be happening within intermittently distributed aligned structures of different amplitudes and sizes, on which the disruption scales will depend (Mallet et al. Reference Mallet, Schekochihin and Chandran2017b). Thus, each scale $\lambda _n$ will in fact be smeared over some range and, as the successive intervals $(\lambda _n,\lambda _{n+1})$ become narrower, this smear can easily exceed their width. Pending a detailed theory of intermittency in the tearing-mediated range, perhaps the best way to think of the spectrum and other scalings in this range is, therefore, in a ‘coarse-grained’ sense, focusing on the characteristic dependence of all interesting quantities on $\lambda _n$, treated as a continuous variable.

7.2.2. Restoration of Kolmogorov cutoff

To calculate the dissipative cutoff for the tearing-mediated cascade, one balances the larger of the viscous and resistive diffusion rates with the nonlinear cascade rate (7.13) (this is the longest of the nonlinear time scales involved): using (7.14), one gets

(7.17)\begin{align} \frac{\nu+\eta}{\lambda_{\eta}^2} \sim \gamma\quad\Rightarrow\quad \lambda_{\eta} \sim \varepsilon^{{-}3/4}\eta^{3/4}(1+\mathrm{Pm})^{3/2} \sim \lambda_{\mathrm{CB}}\mathrm{Rm}^{{-}3/4}(1+\mathrm{Pm})^{3/2} \equiv\lambda_{\eta}^\mathrm{tearing}. \end{align}

For $\mathrm {Pm}\lesssim 1$, this is the same scale as (7.16), where the $-11/5$ and $-5/3$ scalings meet; for $\mathrm {Pm}\gg 1$, (7.17) is reached before (7.16). The condition for the range $[\lambda _{\mathrm {D}},\lambda _{\eta }^\mathrm {K}]$ to be non-empty is

(7.18)\begin{equation} \frac{\lambda_{\mathrm{D}}}{\lambda_{\eta}^\mathrm{tearing}}\sim\left[\frac{\mathrm{Rm}}{(1+\mathrm{Pm})^{10}}\right]^{5/28}\gg1. \end{equation}

This is less stringent than (7.3), so will always be satisfied if the disruption occurs in the first place.

The good news (or, at any rate, the news) is that Kolmogorov's scaling of the dissipative cutoff is rehabilitated for $\mathrm {Pm}\lesssim 1$. Notably, it is not quite rehabilitated for $\mathrm {Pm}\gg 1$, but that is likely an illusion. Indeed, while viscosity destroys velocities, aligned magnetic structures survive unscathed (cf. § 11), so tearing can keep going in the viscously dominated regime. If it breaks up the aligned structures at the scale $\lambda _{\eta }^\mathrm {tearing}$ into flux ropes of the same width and lesser alignment, the turnover time of these structures will be shorter than the viscous-diffusion time – in just the same way as the debris of disruption had shorter turnover times than their mother sheets in § 7.2.1 – and the tearing-mediate cascade can continue. A rough estimate of how small the debris have to get to be killed completely by diffusion is to set their Lundquist number (7.11) obtained on the assumption of complete lack of alignment to $S_{\lambda _n}\sim 1+\mathrm {Pm}$ (in other words, the cutoff occurs when either $\mathrm {Rm}$ or $\mathrm {Re}$ associated with the $\lambda _n$-scale structures is order unity). This gives

(7.19)\begin{equation} \lambda_{\eta}^\mathrm{K} \sim \lambda_{\mathrm{CB}} \left(\frac{\mathrm{Rm}}{1+\mathrm{Pm}}\right)^{{-}3/4} = \lambda_{\mathrm{CB}}\widetilde{\mathrm{Re}}^{{-}3/4} \sim \frac{(\nu+\eta)^{3/4}}{\varepsilon^{1/4}}, \end{equation}

where $\widetilde {\mathrm {Re}}$, defined in (6.24), is $\mathrm {Rm}$ when $\mathrm {Pm}\lesssim 1$ and $\mathrm {Re}$ when $\mathrm {Pm}\gg 1$. This is the proper, classic Kolmogorov scale.

It is interesting to recall that it is the Kolmogorov scaling at (and of) the dissipation scale that was the strongest claim made by Beresnyak (Reference Beresnyak2011, Reference Beresnyak2012a, Reference Beresnyak2014b, Reference Beresnyak2019) on the basis of a convergence study of his numerical spectra (see § 6.3 and figure 10b). While he inferred from that an interpretation of these spectra as showing a $-5/3$ scaling in the inertial range, it is their convergence at the dissipative end of the resolved range that appeared to be the least negotiable feature of his work. He may well have been right. Indeed, his largest simulations (see figure 10) fall somewhere in between the condition (7.18) ($\lambda _{\mathrm {D}}/\lambda _{\eta }^\mathrm {K} \gtrsim 3$ would require perhaps $\mathrm {Rm}\gtrsim 10^3$ at $\mathrm {Pm}\sim 1$) and the more stringent condition (7.3) needed to stop Boldyrev's cutoff (6.24) from taking over ($\lambda _{\mathrm {D}}/\lambda _{\eta } \gtrsim 3$ if $\mathrm {Rm}\gtrsim 10^5$). Thus, $\lambda _{\mathrm {D}}$ in Beresnyak's simulations could not have been more than a factor of order unity larger than the dissipation scale. An optimist might argue that this could have been just about enough to pick up the Kolmogorov scaling of the latter.

7.2.3. Alignment in the tearing-mediated range

The structures corresponding to the lower (GS95) envelope (7.10) are unaligned, whereas the alignment corresponding to the upper envelope (7.14) is the tightest alignment sustainable in the tearing-mediated range and achieved by each aligned ‘mini-cascade’ just before it is disrupted by tearing at the scale $\lambda _n$. This is (cf. Boldyrev & Loureiro Reference Boldyrev and Loureiro2017)

(7.20)\begin{equation} \sin\theta_{\lambda_n^+} \sim \frac{\lambda_n/\delta Z_{\lambda_n^+}}{\tau_{\mathrm{nl}}} \sim \left(\frac{\lambda_{\mathrm{D}}}{\lambda_{\mathrm{CB}}}\right)^{1/4}\left(\frac{\lambda_n}{\lambda_{\mathrm{D}}}\right)^{{-}4/5}. \end{equation}

Equivalently, the fluctuation-direction coherence scale is

(7.21)\begin{equation} \xi_n \sim \frac{\lambda_n}{\sin\theta_{\lambda_n^+}}\sim \lambda_{\mathrm{CB}} \left(\frac{\lambda_{\mathrm{D}}}{\lambda_{\mathrm{CB}}}\right)^{3/4}\left(\frac{\lambda_n}{\lambda_{\mathrm{D}}}\right)^{9/5}. \end{equation}

The corresponding spectral exponent is again $-5/3$, which is automatically the case given the definitions of $\tau _{\mathrm {nl}}$, $\theta _\lambda$ and $\xi$ (see (6.4) and § 6.5).

Thus, the smallest possible alignment angle, having reached its minimum at $\lambda _{\mathrm {D}}$, gets larger through the tearing-mediated range, according to (7.20). It becomes order unity at the scale (7.16), which is the same as the Kolmogorov cutoff (7.19) for $\mathrm {Pm}\lesssim 1$ and a bit smaller than it for $\mathrm {Pm}\gg 1$ (but in fact there is no more alignment below the Kolmogorov cutoff).

To the (doubtful) extent that existing numerical evidence can be considered to be probing this regime, perhaps we can take heart from the numerical papers by both Beresnyak and by Boldyrev's group cited in § 6.3 all reporting that alignment fades away at the small-scale end of the inertial range – although this may also be just a banal effect of the numerical resolution cutoff.

7.2.4. Parallel cascade in the tearing-mediated range

As ever, CB should be an enduring feature of our turbulence. This means that the parallel spectrum (5.1) will not notice the disruption scale and blithely extend all the way through the tearing-mediated range. Since the unaligned (or less aligned) flux ropes produced in the wake of the disruption of aligned structures have a shorter decorrelation time than their aligned progenitors, they should break up in the parallel direction (cf. Zhou, Loureiro & Uzdensky Reference Zhou, Loureiro and Uzdensky2020 and § 12.6). The resulting parallel coherence scale, the same as the scale (5.5) in the GS95 theory, is the lower bound on $l_{\parallel }$ at each $\lambda _n$. The upper bound can be inferred by equating the nonlinear time (7.13) at $\lambda _n$ to the Alfvén time $l_{\parallel }/v_{\mathrm {A}}$. The result is

(7.22)\begin{equation} v_{\mathrm{A}}\varepsilon^{{-}1/3}\lambda_n^{2/3}\lesssim l_{{\parallel}}\lesssim v_{\mathrm{A}}\varepsilon^{{-}1/5}\eta^{{-}2/5}(1+\mathrm{Pm})^{1/5}\lambda_n^{6/5} \sim L_{{\parallel}}\left(\frac{\lambda_{\mathrm{D}}}{\lambda_{\mathrm{CB}}}\right)^{1/2}\left(\frac{\lambda_n}{\lambda_{\mathrm{D}}}\right)^{6/5}. \end{equation}

Thus, the upper bound on the parallel anisotropy $l_{\parallel }/\lambda$ decreases with scale in this range (turbulence becomes less anisotropic).

7.3. Is this a falsifiable theory?

With considerable difficulty.

Numerically, anything like a definite confirmation of the tearing disruption (§ 7.1) and the existence of a tearing-mediated cascade (§ 7.2) requires formidably large simulations: the condition (7.3) demands $\mathrm {Rm}\sim 10^5$ at least (estimated via the frivolous but basically sound principle that the smallest large number is 3) – and probably quite a bit larger if one is to see the scaling of the tearing-mediated-range spectrum (7.15). However, as was mooted above, an optimist might find cause for optimism in the evidence of the MHD turbulence cutoff seen in simulations appearing to obey the Kolmogorov scaling (7.19) (see discussion in § 7.2.2) or in the numerically measured alignment petering out at the small-scale end of the inertial range, as required by § 7.2.3. While the trouble to which I have gone to keep track of the $\mathrm {Pm}$ dependence of the tearing-mediated-range quantities did not yield anything qualitatively spectacular, there is perhaps an opportunity here for numerical tests: e.g., can one obtain Boldyrev's scaling (6.24) of the dissipation cutoff in the limit of moderate $\mathrm {Rm}$ and large $\mathrm {Pm}$? – which, in view of (7.3), is unlikely to need to be very large to take over and shut down tearing.

One way to circumvent the need for getting into a hyper-asymptotic regime is to simulate directly the dynamics of structures that resemble Alfvénic sheets deep in the inertial range. Such a study by Walker, Boldyrev & Loureiro (Reference Walker, Boldyrev and Loureiro2018), in 2D, of the decay of an Alfvénic ‘eddy’ highly anisotropic in the perpendicular plane, has shown it breaking up promisingly into plasmoids and giving rise to a steeper spectrum than exhibited by a larger-$\mathrm {Rm}$ case where tearing was too slow. Dong et al. (Reference Dong, Wang, Huang, Comisso and Bhattacharjee2018) went further and actually demonstrated a spectral break at the disruption scale and a $k_{\perp }^{-11/5}$ spectrum below it, with sheets in a turbulent system very vividly breaking up into plasmoids (figure 15) – but still in 2D. So far so good, but the real prize is always for 3D. The present review has spent so long going through revisions that 3D results, of which I originally wrote in remote future tense, are in fact now on the brink of arriving, again due to Dong et al. (Reference Dong, Wang, Comisso, Huang and Bhattacharjee2021), who have performed the mother of all MHD simulations (with $\mathrm {Rm}$ comfortably over $10^5$) and are claiming to see a ‘subinertial range’ with a $-11/5$ spectral slope, as per (7.15); the spectral break where this scaling starts appears to be in the right place, $\lambda _{\mathrm {D}}\propto \mathrm {Rm}^{-4/7}$ as per (7.2). If I read their results correctly (based on C. Dong's presentation at the 2021 APS DPP Meeting), hallelujah.

Figure 15. A snapshot of current density ($j_z$) from a 2D, $\mathrm {Rm}=10^6$ (64,000$^2$) MHD simulation by Dong et al. (Reference Dong, Wang, Huang, Comisso and Bhattacharjee2018) (reprinted with permission from Dong et al. Reference Dong, Wang, Huang, Comisso and Bhattacharjee2018, copyright (2018) by the American Physical Society; I am grateful to C. Dong for letting me have the original figure file). Zoomed areas show sheets broken up into plasmoids. The 3D versions of these visualisations reported by Dong et al. (Reference Dong, Wang, Comisso, Huang and Bhattacharjee2021) look generally messier, with islands less clearly delineated (they look a bit like the flux ropes in figure 50), but not, at first glance, qualitatively different in a paradigm-changing way.

Observationally, our best bet for fine measurements of turbulence is the solar wind and the terrestrial magnetosphere (e.g., the magnetosheath). However, these are collisionless environments, so, before any triumphs of observational confirmation can be celebrated, all the resistive reconnection physics on which the tearing-mediated-range cascade depends needs to be amended for the cornucopia of kinetic effects that await at the small-scale end of the cascade (see § 14.1). Once the tearing disruption in MHD was proposed, such generalisations were the ripe, low-hanging fruit quite a lot of which was immediately picked (Mallet, Schekochihin & Chandran Reference Mallet, Schekochihin and Chandran2017a; Loureiro & Boldyrev Reference Loureiro and Boldyrev2017a; Boldyrev & Loureiro Reference Boldyrev and Loureiro2019; see also Loureiro & Boldyrev Reference Loureiro and Boldyrev2018, where these ideas were ported to pair plasmas).

Pending all this validation and verification, the tearing-mediated cascade remains a beautiful fantasy – but one must be grateful that after half a century of scrutiny, MHD turbulence still has such gifts to offer.

7.4. Tearing disruption, plasmoid chains, fast reconnection, and reconnection-driven turbulence

This section deals with what many readers might feel are fairly esoteric details. They are right – and so skipping straight to § 8 will not subtract much from their experience.

In what until recently was a separate strand of research, much interest in the reconnection community has focused on stochastic plasmoid chains that arise in current sheets susceptible to the plasmoid instability (a subspecies of tearing), where a lively population of islands (plasmoids) are born, grow, travel along the sheet with Alfvénic outflows, occasionally eat each other (coalesce),Footnote ³⁷ and, as shown by Uzdensky, Loureiro & Schekochihin (Reference Uzdensky, Loureiro and Schekochihin2010), cause reconnection in the sheet that they inhabit to be fast,Footnote ³⁸ meaning independent of $\eta$ as $\eta \to +0$ (a derivation of the plasmoid instability, a long list of references on stochastic plasmoid chains, and, in figure 46, an example of one, can be found in Appendix D.4.2; the Uzdensky et al. Reference Uzdensky, Loureiro and Schekochihin2010 argument is reproduced in Appendix D.6). A stochastic chain can be viewed as a kind of ‘1D turbulence’, and has some distinctive statistical properties (see Appendices D.6.1 and D.6.2). Should one imagine the disrupted aligned structures spawning multiple instances of such a turbulence, and does the simple theory presented in § 7.2 describe this situation or does it need to be revised to represent a superposition of many fast-reconnecting, plasmoid-infested sheets (as attempted in three different ways by Loureiro & Boldyrev Reference Loureiro and Boldyrev2017b, Reference Loureiro and Boldyrev2020 and Tenerani & Velli Reference Tenerani and Velli2020b)?

7.4.1. Nature of tearing disruption

In considering this problem, I first want to return to the question of what the ‘disruption’ of the aligned structures actually consists of. There are two lines of thinking on this, articulated most explicitly in the papers by Mallet et al. (Reference Mallet, Schekochihin and Chandran2017b) and Boldyrev & Loureiro (Reference Boldyrev and Loureiro2017), of which I have so far stuck with the former. Namely, at the end of § 7.1, I followed Uzdensky & Loureiro (Reference Uzdensky and Loureiro2016), Mallet et al. (Reference Mallet, Schekochihin and Chandran2017b) and Loureiro & Boldyrev (Reference Loureiro and Boldyrev2017b) in invoking the collapse of the $X$-points separating the tearing-mode islands as a means of consummating the disruption of the aligned structure – the theory of the tearing-mediated cascade in § 7.2 was then presented as a corollary of this view.

The collapse of inter-island $X$-points is actually the first step in the formation of a stochastic chain. A nuance that I have previously elided is that comparing the $X$-point collapse rate (which is $\sim$ the tearing growth rate $\gamma$) with the growth rate of a secondary tearing instability of the same $X$-point shows that, at asymptotically large Lundquist numbers, the latter is always greater than the former. Therefore, the collapse may itself be disrupted by tearing, producing more islands and more $X$-points, followed by the collapse of those, also disrupted, and so on. My take on this recursive tearing is presented in Appendix D.5.2 (alongside a review of prior literature, starting with Shibata & Tanuma Reference Shibata and Tanuma2001). I argue there that the smaller-scale islands that are produced in this process are not energetically relevant and so we need not worry about including them to amend the ‘one-level’ scenario of tearing disruption presented in § 7.1.

The recursive tearing proceeds until inter-island current sheets are short enough to be stable, at which point the true nonlinear plasmoid chain can form, involving not just multiple tearings, but also nonlinear plasmoid growth by reconnection, their coalescence, and ejection from the sheet. While the multiscale statistics of such a chain may be different from that of a tearing-mediated cascade that I described in § 7.2, I assumed there, implicitly, that the chain could not survive for a long time, if at all: indeed, the characteristic time scale of the process of fully forming the sheet out of an aligned structure is the tearing time ($\sim \gamma ^{-1}$), and the time to break apart that aligned structure by ideal-MHD dynamics is of the same order ($\sim \tau _{\mathrm {nl}}$). So it seems that the ‘mother sheet’ (collapsed aligned structure) should break apart entirely shortly after (or even before) fully forming and release its plasmoids (flux ropes) into the general turbulent wilderness, where they are free to interact with each other or with anything else that comes along, and are thus no different from turbulent fluctuations of a particular size generically splashing around in a large nonlinear system. This gives rise to the ‘mini-cascades’ in § 7.2, with the overall $k_{\perp }^{-11/5}$ spectral envelope (7.15).

Boldyrev & Loureiro (Reference Boldyrev and Loureiro2017) also derive the $k_{\perp }^{-11/5}$ spectrum by using (7.13) as the operational prescription for the cascade time ($\tau _{\mathrm {nl}}\sim \gamma ^{-1}$ at each scale in the tearing-mediated range). However, they have a different narrative about what happens dynamically: they do not believe that inter-island $X$-points ever collapse, but that, rather, the tearing mode upsets alignment by order unity, changing the effective nonlinear cascade rate to the tearing rate. This is based on the (correct) observation that the alignment angle at the disruption scale (7.2)

(7.23)\begin{equation} \sin\theta_{\lambda_{\mathrm{D}}}\sim S_{\lambda_{\mathrm{D}}}^{{-}1/2}(1+\mathrm{Pm})^{{-}1/4} \end{equation}

is the same (at least for $\mathrm {Pm}\lesssim 1$) as the angular distortion of the field line caused by the tearing perturbation at the onset of the nonlinear regime: indeed, using (7.7) and (7.4) at $\lambda =\lambda _{\mathrm {D}}$,

(7.24)\begin{equation} \theta_\mathrm{tearing}\sim wk_* \sim (k_*\lambda_{\mathrm{D}})^2 \sim S_{\lambda_{\mathrm{D}}}^{{-}1/2}(1+\mathrm{Pm})^{1/4}. \end{equation}

They think that this is enough to make the aligned structure ‘cascade’, in some unspecified manner, without much reconnection, production of flux ropes, etc. In Loureiro & Boldyrev (Reference Loureiro and Boldyrev2020), they revise their view a little and allow that, since the collapse time, the tearing time and, therefore, the cascade time are comparable to each other, some aligned structures might, in fact, collapse into proper reconnecting sheets.Footnote ³⁹

It is hard to say whether the two pictures outlined above represent a disagreement in substance or merely in the style of presentation. While in the Mallet et al. (Reference Mallet, Schekochihin and Chandran2017b) interpretation, the collapse of the inter-island $X$-points is the way in which the distortion of alignment caused by tearing leads to faster nonlinear break-up of the aligned structures, Loureiro & Boldyrev (Reference Loureiro and Boldyrev2020) think this is not necessary but does happen with some finite probability. This might not be sufficiently quantifiable a difference to be testable.

7.4.2. Onset of fast reconnection

What may be consequential physically, however, is the onset of fast reconnection in those aligned, tearing structures that do manage to collapse into proper sheets. Loureiro & Boldyrev (Reference Loureiro and Boldyrev2020) conjecture that when that happens, one should start worrying about the resulting reconnecting sheets making a difference to the nature of the tearing-mediated (now reconnection-mediated) turbulence.Footnote ⁴⁰ In order for such a transition to be realisable, $\mathrm {Rm}$ must be large enough for the characteristic time of fast reconnection to be shorter than the cascade time:

(7.25)\begin{equation} \tau_{\mathrm{rec}} \sim \epsilon_{\mathrm{rec}}^{{-}1}\frac{\lambda}{\delta Z_\lambda} \lesssim \tau_{\mathrm{nl}} \sim \frac{\lambda}{\delta Z_\lambda} \left(\sin\theta_\lambda\right)^{{-}1} \quad\Leftrightarrow\quad \sin\theta_\lambda \lesssim \epsilon_{\mathrm{rec}}, \end{equation}

where the reconnection time has been estimated as the time for all of the flux in a magnetic structure of scale $\lambda$ to be reconnected, and $\epsilon _{\mathrm {rec}}\sim 10^{-2}(1+\mathrm {Pm})^{-1/2}$ is the dimensionless fast-reconnection rate in a Uzdensky et al. (Reference Uzdensky, Loureiro and Schekochihin2010) plasmoid chain (see Appendix D.6), or some modified version of it appropriate for a turbulent environment. The estimate (7.25) simply says that if the alignment angle (inverse aspect ratio) of a structure manages to become smaller than the reconnection rate, such a structure might be capable of becoming a fast-reconnecting sheet.

Let us recall that the alignment angle decreases with $\lambda$ according to (6.22) in an aligned MHD cascade, reaches its minimum at the disruption scale $\lambda _{\mathrm {D}}$ (given by (7.2)), and then increases with decreasing $\lambda$ according to (7.20) in a tearing-mediated cascade. The condition (7.25) is realisable if

(7.26)\begin{equation} \sin\theta_{\lambda_{\mathrm{D}}}\sim \left(\frac{\lambda_{\mathrm{D}}}{\lambda_{\mathrm{CB}}}\right)^{1/4} \lesssim \epsilon_{\mathrm{rec}} \quad\Leftrightarrow\quad \mathrm{Rm} \gtrsim \epsilon_{\mathrm{rec}}^{{-}7}(1+\mathrm{Pm})^{{-}1/2} \sim 10^{14}(1+\mathrm{Pm})^3. \end{equation}

Obviously, this is never going to be numerically (or indeed experimentally) achievable for resistive MHD – unless $\epsilon _{\mathrm {rec}}$ is substantially enhanced in a turbulent setting (Loureiro & Boldyrev Reference Loureiro and Boldyrev2020 show that this ‘nonlinear-reconnecting’ regime might also be more easily accessible in certain kinetic settings).

In the asymptotic world where the condition (7.25) can be realised, it will be realised for $\lambda \in [\lambda _{\mathrm {rec}}^<,\lambda _{\mathrm {rec}}^>]$, with the two scales lying below and above $\lambda _{\mathrm {D}}$: using (6.22) and (7.20) for the alignment angle in the aligned and tearing-mediated cascades, respectively, one gets

(7.27)\begin{equation} \lambda_{\mathrm{rec}}^> \sim \epsilon_{\mathrm{rec}}^4\lambda_{\mathrm{CB}},\qquad \lambda_{\mathrm{rec}}^< \sim \epsilon_{\mathrm{rec}}^{{-}5/4}\lambda_{\mathrm{D}}^{21/16}\lambda_{\mathrm{CB}}^{{-}5/16}. \end{equation}

One might argue that $\lambda _{\mathrm {rec}}^>$ is irrelevant because at that scale ideal nonlinear interactions are faster than tearing, so the instability that assists in the formation of a plasmoid chain is too slow to get it going before the energy cascades to small scales (if this is not true, then reconnection-mediated turbulence sets in at $\lambda _{\mathrm {rec}}^>$, an $\mathrm {Rm}$-independent scale, putting a hard lower bound on the allowed alignment angles, $\sin \theta _\lambda \gtrsim \epsilon _{\mathrm {rec}}$, and vindicating Beresnyak Reference Beresnyak2019). In contrast, the interval $[\lambda _{\mathrm {D}},\lambda _{\mathrm {rec}}^<]$ is tearing-dominated, so plasmoid chains could be on the cards.

What might turbulence in such a situation look like? To get some idea of that, it is perhaps wise to start with numerical evidence about turbulence in individual, ‘stand-alone’ plasmoid chains – a subject on which copious literature exists, quoted in Appendix D.4.2.

8.1. Is this the end of the road?

It never quite is (see §§ 7.4.2–7.4.3 and the second part of this review, starting from § 9), but the basic story looks roughly complete for the first time in years, at least as far as forced, balanced RMHD turbulence is concerned. The principle of critical balance gave us a sound ideology for dealing with anisotropic turbulence in a system that supports propagation of waves (§ 5) and a firm prescription for the parallel spectrum (§ 5.2). The aligned cascade (§ 6) produced a plausible prediction for the perpendicular spectrum but used to have an air of unfinished business about it, both in the sense that it gave rise to a state that appeared unsustainable at asymptotically small scales and in view of the objections, physical and numerical, raised by Beresnyak (Reference Beresnyak2011, Reference Beresnyak2012a, Reference Beresnyak2014b, Reference Beresnyak2019). With the revised interpretation of alignment as an intermittency effect (§ 6.4) and with the tearing-mediated cascade (§ 7) connecting the inertial-range, aligned cascade to the Kolmogorov cutoff (7.19), these issues appear to be satisfactorily resolved. In what is also an aesthetically pleasing development, the tearing-mediated cascade has emerged as an ingenious way in which MHD turbulence contrives to thermalise its energy while shedding the excessive alignment that ideal-MHD dynamics could not help producing in the inertial range. This development joins together in a most definite way the physics of turbulence and reconnection – arguably, this was always inevitable, but it is good that we now appear to have some grip on what happens specifically.

Before I move on to the miscellany of Part II, let me make a few comments about the robustness of the general picture presented above and its connection to other schools of thought on MHD turbulence.

8.2. What can go wrong?

It is only fair to spell out explicitly what is settled and what can go catastrophically wrong with this entire picture.

The principle of critical balance and, therefore, the theory of the parallel cascade (§ 5) are, in my view, quite safe. They are straightforward physically and have been quite convincingly verified both observationally and numerically (pace the ‘waves vs. structures’ confusion: see § 8.3.1). The CB approach also appears to offer an attractive and credible strategy for dealing with turbulence in wave-carrying systems other than MHD, in both plasmas and hydrodynamics (e.g., Cho & Lazarian Reference Cho and Lazarian2004; Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009, Reference Schekochihin, Parker, Highcock, Dellar, Dorland and Hammett2016; Schekochihin, Kawazura & Barnes Reference Schekochihin, Kawazura and Barnes2019; Nazarenko & Schekochihin Reference Nazarenko and Schekochihin2011; Barnes, Parra & Schekochihin Reference Barnes, Parra and Schekochihin2011; Boldyrev et al. Reference Boldyrev, Horaites, Xia and Perez2013; Chen & Boldyrev Reference Chen and Boldyrev2017; Passot, Sulem & Tassi Reference Passot, Sulem and Tassi2017; Loureiro & Boldyrev Reference Loureiro and Boldyrev2018; Avsarkisov Reference Avsarkisov2020; Adkins et al. Reference Adkins, Schekochihin, Ivanov and Roach2022; Skoutnev Reference Skoutnev2022), further bolstering its claim to being a universal physical principle.

In contrast, the aligned perpendicular cascade (§ 6), its tearing disruption (§ 7.1) and replacement by a tearing-mediated cascade (§ 7.2), which have formed the bulk of my story so far, are hotly debated concepts. Arguably, it is still subject to verification (requiring currently inaccessible resolutions) that alignment is not a transient, large-scale feature, as Beresnyak (Reference Beresnyak2019) would have it. It seems to me that we do know, however, that if we stir up unaligned turbulence, it will get aligned at smaller scales (see numerical studies cited in § 6.3), so its possible transient nature can only be due to some secondary instability of the aligned structures. The picture presented above relied on this being the tearing instability – but it is not entirely impossible that it is, in fact, an ideal MHD instability, e.g., some version of Kelvin–Helmholtz (KH) instability. The difference is that tearing required resistivity and so the disruption scale $\lambda _{\mathrm {D}}\propto \eta ^{4/7}$ was asymptotically separated from the outer scale $\lambda _{\mathrm {CB}}$ (see (7.2)), whereas the KH instability would kick in at some $\lambda \sim$ a finite fraction of $\lambda _{\mathrm {CB}}$. The usual expectation is that the KH instability is quenched by the magnetic field (and indeed hence perhaps the statistical preponderance of current sheets over shear layers; see § 10.4), but this can in principle turn out not to be enough. If it does, ideal MHD will take care of limiting alignment, without alignment there will be no need for, or dynamical tendency to, the tearing-mediated remedy to it at small scales, and, presumably, we will be back to GS95, in which case I apologise to my readers for having wasted their time (this scenario is explained in another language at end of § 8.3.3 and at the beginning of Appendix D.7.4).

8.3. What is lost in translation?

This section is devoted to misunderstandings and arguments at cross purposes. As any vibrant, fast-evolving field, MHD turbulence is a subject talked about in many languages, and differences in vocabulary are sometimes mistaken for fundamental disagreements.

9.1. Imbalance global and local

Since both incompressible MHD and RMHD conserve two invariants – the total energy and cross-helicity – each of the two Elsasser fields $\boldsymbol {Z}_{\perp }^\pm$ has its own conserved energy (see (3.3)). The energy fluxes $\varepsilon ^\pm$ of these fields are, therefore, independent parameters of MHD turbulence. Setting them equal to each other makes arguments simpler, but does not, in general, correspond to physical reality, for a number of reasons.

First, everyone's favourite case of directly measurable MHD turbulence is the solar wind, where the Alfvénic perturbations propagating away from the Sun are launched from the Sun (Roberts et al. Reference Roberts, Goldstein, Klein and Matthaeus1987), while the counterpropagating ones have to be supplied by some mechanism that is still under discussion and probably involves Alfvén-wave reflection as plasma density decreases outwards from the Sun (see Chandran & Perez Reference Chandran and Perez2019 and references therein). The counterpropagating component is usually energetically smaller, especially in the fast wind (Bruno & Carbone Reference Bruno and Carbone2013; Chen et al. Reference Chen, Bale, Bonnell, Borovikov, Bowen, Burgess, Case, Chandran, Dudok de Wit and Goetz2020).

Secondly – and, for a theoretical physicist interested in universality, more importantly – it is an intrinsic property of MHD turbulence to develop regions of local imbalance. This can be understood dynamically as a desire to evolve towards an Elsasser state, $\boldsymbol {Z}_{\perp }^+=0$ or $\boldsymbol {Z}_{\perp }^-=0$, which is an exact solution of RMHD equations (confirmed in simulations of decaying RMHD turbulence; see § 12.7), or statistically as a tendency for the local dissipation rates $\varepsilon ^\pm$ to fluctuate in space – a mainstay of intermittency theories since Landau's famous objection (see Frisch Reference Frisch1995) to Kolmogorov (Reference Kolmogorov1941b) and the latter's response in the form of the refined similarity hypothesis, accepting a fluctuating $\varepsilon$ (Kolmogorov Reference Kolmogorov1962) (the theories of intermittency for balanced MHD turbulence proposed by Chandran et al. Reference Chandran, Schekochihin and Mallet2015 and Mallet & Schekochihin Reference Mallet and Schekochihin2017, skimmed through in § 6.4.2, were based on the same premise). In this context, a complete intermittency theory for MHD turbulence must incorporate whatever local modification (if any) of the MHD cascade is caused by $\varepsilon ^+\neq \varepsilon ^-$, something that no existing theory has as yet accomplished or attempted. Another influential school of thought on the root causes of local imbalance connects it to a tendency for the cross-helicity to be less vigorously cascaded than energy, leading to local enhancements of $\boldsymbol {u}_{\perp }\cdot \boldsymbol {b}_{\perp }$ (see footnotes 15 and 16 and references therein).

That an intimate connection must exist between any verifiable theory of MHD turbulence and local imbalance is well illustrated (in figure 16) by the following piece of observational analysis, rather noteworthy, in my (not impartial) view. Wicks et al. (Reference Wicks, Roberts, Mallet, Schekochihin, Horbury and Chen2013b) took a series of measurements by Wind spacecraft of magnetic and velocity perturbations in fast solar wind and sorted them according to the amount of imbalance, both Elsasser and Alfvénic (§ 10), at each scale. They then computed structure functions conditional on these imbalances. While the majority of perturbations were imbalanced one way or the other (or both), there was a subpopulation with $\delta Z^+_\lambda \sim \delta Z^-_\lambda \sim \delta b_\lambda \sim \delta u_\lambda$. Interestingly, the structure function restricted to this subpopulation had what seemed to be a robust GS95 scaling (corresponding to a $k_{\perp }^{-5/3}$ spectrum), even though the structure functions of the imbalanced perturbations were consistent with Boldyrev's $k_{\perp }^{-3/2}$ aligned-cascade scaling and indeed exhibited some alignment, unlike the GS95 population (although not necessarily an alignment with the theoretically desirable scale dependence: see footnote 21 and Wicks et al. Reference Wicks, Mallet, Horbury, Chen, Schekochihin and Mitchell2013a; Podesta & Borovsky Reference Podesta and Borovsky2010 reported analogous results, conditioning on the presence of cross-helicity only). It is important to recognise that imbalance and alignment of Elsasser fields do not automatically imply each other, so balanced fluctuations are not absolutely required to be unaligned, or aligned fluctuations to be imbalanced (see Appendix B.1). However, as I argued in § 6.4.2, dynamical alignment is an intermittency effect and so there may be a correlation between the emergence of imbalanced patches at ever smaller scales and Elsasser fields shearing each other into alignment (cf. Chandran et al. Reference Chandran, Schekochihin and Mallet2015).

Figure 16. (a) Distribution of normalised cross-helicity ($\sigma _{\mathrm {c}}$) and residual energy ($\sigma _{\mathrm {r}}$) (defined in (B4)) in an interval of fast-solar-wind data taken by Wind spacecraft and analysed by Wicks et al. (Reference Wicks, Roberts, Mallet, Schekochihin, Horbury and Chen2013b), from whose paper both plots in this figure are taken (© AAS, reproduced with permission). (b) Structure functions corresponding to the total energy (sum of kinetic and magnetic) conditioned on values of $\sigma _{\mathrm {c}}$ and $\sigma _{\mathrm {r}}$ and corresponding to Regions 1 (balanced), 2 (Elsasser-imbalanced), and 3 (Alfvénically imbalanced towards magnetic perturbations), indicated in panel (a). The $f^{-2/3}$ slope corresponds to a $k_{\perp }^{-5/3}$ spectrum, the $f^{-1/2}$ slope to a $k_{\perp }^{-3/2}$ one.

Intuitively then, since patches of imbalance are locally ubiquitous even in globally balanced turbulence (Perez & Boldyrev Reference Perez and Boldyrev2009; see figure 17) and since the theory of balanced turbulence described in § 6.4 incorporates intermittency effects in the form of alignment, we might expect that this allows for local imbalance – and, therefore, that mildly imbalanced turbulence might look largely similar to the balanced one. Indeed, how would perturbations in the middle of inertial range ‘know’ that the local imbalance they ‘see’ is local rather than global? Obviously, on average, there will not be an imbalance and so the results for $\delta Z_\lambda$ that one derives for balanced turbulence (§§ 6 and 7) are effectively averaged over the statistics of the stronger and weaker Elsasser fields – which of $\delta Z^+_\lambda$ and $\delta Z^-_\lambda$ is which, depends on time and space.

Figure 17. Cosine of the angle between increments $\delta \boldsymbol {u}_{\boldsymbol {\lambda }}$ and $\delta \boldsymbol {b}_{\boldsymbol {\lambda }}$ in the $(x,y)$ plane, for $\lambda = L_{\perp }/6$ (left) and $\lambda = L_{\perp }/12$ (right, corresponding to the region demarcated by the white square within the left panel) in a balanced RMHD simulation by Perez & Boldyrev (Reference Perez and Boldyrev2009) (reprinted with permission from Perez & Boldyrev Reference Perez and Boldyrev2009, copyright (2009) by the American Physical Society). Since $\delta \boldsymbol {u}_{\boldsymbol {\lambda }}\cdot \delta \boldsymbol {b}_{\boldsymbol {\lambda }} =\left (|\delta \boldsymbol {Z}^+_{\boldsymbol {\lambda }}|^2 - |\delta \boldsymbol {Z}^-_{\boldsymbol {\lambda }}|^2\right )/4$, this is an illustration of patchy local imbalance, as well as of local alignment between the velocity and magnetic field.

If we now allow $\varepsilon ^+>\varepsilon ^-$ on average, it becomes reasonable to expect $\delta Z^+_\lambda > \delta Z^-_\lambda$ nearly everywhere or, at least, typically – unless $\varepsilon ^+/\varepsilon ^-$ is close enough to unity that fluctuations of local imbalance overwhelm the overall global one. In the latter case, presumably the global imbalance does not matter – at any rate, in the balanced considerations of §§ 6 and 7, we only ever required $\varepsilon ^+\sim \varepsilon ^-$, rather than $\varepsilon ^+=\varepsilon ^-$ exactly. What I am driving at here, perhaps with too much faffing about, is the rather obvious point that it is only the limit of strong imbalance, $\varepsilon ^+\gg \varepsilon ^-$, that can be expected to be physically distinct, in a qualitative manner, from the balanced regime.

9.2. Numerical and observational evidence

As usual, it is this most interesting limit that is also the hardest to resolve numerically and so we have little definitive information as to what happens in the strongly imbalanced regime. As in the case of the spectra of balanced turbulence, the debate about the numerical evidence regarding the imbalanced cascade and its correct theoretical interpretation has been dominated by the antagonistic symbiosis of Beresnyak and Boldyrev, so it is from their papers (Beresnyak & Lazarian Reference Beresnyak and Lazarian2008, Reference Beresnyak and Lazarian2009b, Reference Beresnyak and Lazarian2010; Beresnyak Reference Beresnyak2019; Perez & Boldyrev Reference Perez and Boldyrev2009, Reference Perez and Boldyrev2010a,Reference Perez and Boldyrevb) that I derive much of the information reviewed below. Perez & Boldyrev (Reference Perez and Boldyrev2010a,Reference Perez and Boldyrevb) argue that large imbalances are unresolvable and refuse to simulate them. Beresnyak & Lazarian (Reference Beresnyak and Lazarian2009b, Reference Beresnyak and Lazarian2010) do not necessarily disagree with this, but believe that useful things can still be learned from strongly imbalanced simulations, even if imperfectly resolved. Based (mostly) on both groups’ simulations, imbalanced MHD turbulence appears to exhibit the following distinctive features (which I recount with a degree of confidence as they have been reproduced in two sets of independent, unpublished RMHD simulations by Mallet & Schekochihin Reference Mallet and Schekochihin2011 and by Meyrand & Squire Reference Meyrand and Squire2020).

(i) The stronger field has a steeper spectrum than the weaker one, with the former steeper and the latter shallower than the standard balanced-case spectra (figure 18a). However, it is fairly certain that these spectra are not converged with resolution: as resolution is increased, the tendency appears to be for the spectral slopes to get closer to each other, both when the imbalance is weak (Perez & Boldyrev Reference Perez and Boldyrev2010a) and when it is strong (Mallet & Schekochihin Reference Mallet and Schekochihin2011). This led Perez & Boldyrev (Reference Perez and Boldyrev2010a) to argue that numerical evidence was consistent with the two fields having the same spectral slope in the asymptotic limit of infinite Reynolds numbers. There is no agreement as to whether the two fields’ spectra might be ‘pinned’ (i.e., equal) to each other at the dissipation scale: ‘yes’, it seems, in weakly imbalanced simulations of Perez & Boldyrev (Reference Perez and Boldyrev2010a), ‘no’ in the strongly imbalanced ones of Beresnyak & Lazarian (Reference Beresnyak and Lazarian2009b) and Meyrand & Squire (Reference Meyrand and Squire2020).Footnote ⁴² In fact, in the latter studies, the dissipation scales of the two fields do not appear to be the same: larger for the stronger field, smaller for the weaker field (see figure 18a and figure 35c).

Figure 18. A typical MHD simulation with large imbalance: (a) spectra, (b) anisotropy, $l_{\parallel }^\pm$ vs. $\lambda$. These plots are adapted from Beresnyak & Lazarian (Reference Beresnyak and Lazarian2009b) (© AAS, reproduced with permission). Mallet & Schekochihin (Reference Mallet and Schekochihin2011) and Meyrand & Squire (Reference Meyrand and Squire2020) have qualitatively similar results (although the difference in slopes between the weaker and the stronger fields’ spectra is much smaller in the higher-resolution simulations of Meyrand & Squire Reference Meyrand and Squire2020 – see an example of that in figure 35c, taken from Meyrand et al. Reference Meyrand, Squire, Schekochihin and Dorland2021).

(ii) The ratio of stronger to weaker field's energies, a crude outer-scale quantity that Beresnyak & Lazarian (Reference Beresnyak and Lazarian2008, Reference Beresnyak and Lazarian2009b, Reference Beresnyak and Lazarian2010) argue (reasonably, in my view) to be more likely to be numerically converged than inertial-range scalings, scales very strongly with $\varepsilon ^+/\varepsilon ^-$: it increases at least as fast as

(9.1)\begin{equation} \frac{\langle|\boldsymbol{Z}_{{\perp}}^+|^2\rangle}{\langle|\boldsymbol{Z}_{{\perp}}^-|^2\rangle} \sim \left(\frac{\varepsilon^+}{\varepsilon^-}\right)^2 \end{equation}

and possibly faster (which is inconsistent with the theory of Perez & Boldyrev Reference Perez and Boldyrev2009, another casus belli for the two groups; see § 9.4). Mallet & Schekochihin (Reference Mallet and Schekochihin2011, see figure 19) and Meyrand & Squire (Reference Meyrand and Squire2020) found the same scaling in their simulations for values of $\varepsilon ^+/\varepsilon ^-$ up to 10 (simulations with much higher imbalance are numerically suspect). In fact, early numerical evidence for (9.1) appears already in the (decaying) simulations of Verma et al. (Reference Verma, Roberts, Goldstein, Ghosh and Stribling1996) (see § 12.7).

Figure 19. Scalings found in (unpublished) $512^3$ RMHD numerical simulations by Mallet & Schekochihin (Reference Mallet and Schekochihin2011): perpendicular (parallel) spectral indices $\mu _\perp$ ($\mu _\parallel$) (inferred from structure functions calculated as explained in § 5.3) for both fields, denoted by the $\pm$ superscripts. In terms of the scaling exponents $\gamma ^\pm _{\perp,\parallel }$ of the field increments ($\delta Z^\pm _\lambda \propto \lambda ^{\gamma ^\pm _\perp }$, $\delta Z^\pm _{l_{\parallel }}\propto l _{\parallel }^{\gamma ^\pm _\parallel }$), these are $\mu ^\pm _{\perp,\parallel } = -2\gamma ^\pm _{\perp,\parallel } - 1$. The last column shows the overall Elsasser ratio $R_{\mathrm {E}}=\langle |\boldsymbol {Z}_{\perp }^+|^2\rangle /\langle |\boldsymbol {Z}_{\perp }^-|^2\rangle$. The parallel scalings of the weaker field were converged with resolution, while the perpendicular scalings of the stronger (weaker) field at $\varepsilon ^+/\varepsilon ^- = 10$ became shallower (steeper) as resolution was increased from $256^3$ to $512^3$ to $1024^2\times 512$; the parallel scaling of the stronger field also appeared to become shallower. Simulations with $\varepsilon ^+/\varepsilon ^- = 100, 1000$ should be viewed as numerically suspect.

(iii) According to Beresnyak & Lazarian (Reference Beresnyak and Lazarian2008, Reference Beresnyak and Lazarian2009b) and Meyrand & Squire (Reference Meyrand and Squire2020), the stronger field is less anisotropic than the weaker one, in the sense that $l_{\parallel }^+< l_{\parallel }^-$ and $l_{\parallel }^+$ drops faster with $\lambda$ than $l_{\parallel }^-$ (figure 18b). Beresnyak (Reference Beresnyak2019) notes that this is true in his simulations even though he forces the two fields with the same parallel scale, i.e., given an opportunity to keep $l_{\parallel }^+=l_{\parallel }^-$, the system refuses to do so.

(iv) Mallet & Schekochihin (Reference Mallet and Schekochihin2011, see figure 19) found that the parallel spectrum of the weaker field (measured via its local-field-parallel structure function, as in § 5.3) was very robustly $k_{\parallel }^{-2}$ – to be precise, the exponent varied between $-1.9$ and $-2.1$, but in a manner that evinced no systematic dependence on $\varepsilon ^+/\varepsilon ^-$. For the stronger field, they found a gradual steepening of the parallel spectrum with higher imbalance.

(v) Beresnyak & Lazarian (Reference Beresnyak and Lazarian2009b) found that the alignment angle between the Elsasser fields, defined as $\sin \theta$ in (B1), with numerator and denominator averaged separately, decreased with scale roughly as $\lambda ^{0.1}$, independently of the degree of imbalance. Mallet & Schekochihin (Reference Mallet and Schekochihin2011) measured the same exponent, quite robustly for a wide range of imbalances, but noticed also that the scaling exponent depended on the definition of the ‘alignment angle’: e.g., if root-mean-square numerator and denominator were used, the scaling was $\lambda ^{0.2\dots 0.25}$, closer to the familiar theoretical result (6.20). This is not special to the imbalanced cascades – the same is true in balanced turbulence (Mallet et al. Reference Mallet, Schekochihin, Chandran, Chen, Horbury, Wicks and Greenan2016).

(vi) The observational picture is only just emerging. A steeper scaling for the stronger field noted in item (i) appears to be consistent with the structure functions measured in the fast solar wind by, e.g., Wicks et al. (Reference Wicks, Horbury, Chen and Schekochihin2011), although, besides this, they also exhibit low Alfvén ratio (see § 10), which simulations do not, and a rather-hard-to-interpret (or, possibly, to trust) scale dependence of the anisotropy. In contrast, Podesta & Borovsky (Reference Podesta and Borovsky2010) report a scale-independent Elsasser ratio and $k^{-3/2}$ spectra for both fields in a number of reasonably imbalanced cases of solar-wind turbulence at 1 AU. The same result has been reported by Chen et al. (Reference Chen, Bale, Bonnell, Borovikov, Bowen, Burgess, Case, Chandran, Dudok de Wit and Goetz2020) from the very recent measurements by the Parker Solar Probe made closer to the Sun, where the imbalance gets larger ($\langle |\boldsymbol {Z}^+|^2\rangle /\langle |\boldsymbol {Z}^-|^2\rangle \approx 15$) – this may be damning for any theory or simulation where the two fields’ spectra scale differently, at least insomuch as these theories or simulations aspire to apply to the solar wind.

(vii) As the solar wind offers practically the only chance of observational testing of theory – a chance greatly enhanced by the launch of the Parker Solar Probe – there is a growing industry of direct numerical modelling of the generation of inward-propagating ($\boldsymbol {Z}^-$) perturbations by reflection of the outward-propagating ones ($\boldsymbol {Z}^+$), which is what is supposed to happen in the expanding solar wind. The latest and most sophisticated study of this kind is Chandran & Perez (Reference Chandran and Perez2019) (who also provide an excellent overview of previous work). Their results appear to be quite different from the idealised periodic-box, artificially-forced-turbulence studies discussed above: the stronger field's spectrum is actually shallower than the weaker one's (sometimes as shallow as $k^{-1}$), but both asymptote towards $k^{-3/2}$ with increasing heliocentric distance – good news for modelling, in view of what Chen et al. (Reference Chen, Bale, Bonnell, Borovikov, Bowen, Burgess, Case, Chandran, Dudok de Wit and Goetz2020) have found. Chandran & Perez (Reference Chandran and Perez2019) acknowledge, however, that they can break these results by fiddling with how their turbulence is forced in the photosphere. Thus, the nature of large-scale energy injection appears to matter,Footnote ⁴³ at least at finite resolutions, perhaps reinforcing the doubts expressed above about the convergence of even the more idealised simulations.

For a short while still, the field appears set to remain open to enterprising theoreticians.

9.3. Lithwick et al. (Reference Lithwick, Goldreich and Sridhar2007)

In what follows, in view of the discussion in § 6.1 and in Appendix B.2, I shall stick with my use of the Elsasser-field alignment angle $\theta$ in the expression (6.4) for $\tau _{\mathrm {nl}}^{\pm }$. This angle is obviously the same for both fields, so

(9.2)\begin{equation} \tau_{\mathrm{nl}}^\pm \sim \frac{\lambda}{\delta Z^\mp_\lambda\sin\theta_\lambda}\quad\Rightarrow\quad \frac{\tau_{\mathrm{nl}}^+}{\tau_{\mathrm{nl}}^-} \sim \frac{\delta Z^+_\lambda}{\delta Z^-_\lambda} > 1, \end{equation}

i.e., the cascade of the stronger field is slower (because it is advected by the weaker field). Assuming nevertheless that both cascades are strong, we infer immediately

(9.3)\begin{equation} \frac{(\delta Z^\pm_\lambda)^2}{\tau_{\mathrm{nl}}^\pm} \sim \varepsilon^\pm\quad\Rightarrow\quad \frac{\delta Z^+_\lambda}{\delta Z^-_\lambda} \sim \frac{\varepsilon^+}{\varepsilon^-}. \end{equation}

Thus, the two fields’ increments have the same scaling with $\lambda$ (the same $k_{\perp }$ spectra) and the ratio of their energies is $\sim (\varepsilon ^+/\varepsilon ^-)^2$, in agreement with (9.1). This is the conclusion at which Lithwick et al. (Reference Lithwick, Goldreich and Sridhar2007, henceforth LGS07) arrived – they considered unaligned GS95-style turbulence ($\sin \theta \sim 1$), but that does not affect (9.3) (note that this result already appeared in (5.4)).

Things are, however, not as straightforward as they might appear. LGS07 point out that it is, in fact, counterintuitive that the weaker $\delta Z_\lambda ^-$ perturbation, which is distorted by $\delta Z_\lambda ^+$ on a shorter time scale $\tau _{\mathrm {nl}}^-$, can nevertheless coherently distort $\delta Z_\lambda ^+$ for a longer time $\tau _{\mathrm {nl}}^+$. Their solution to this is to argue that, while the weaker field is strongly distorted in space by the stronger one, it remains correlated in time for as long as the stronger field does (coherence time of the long-correlated advector inherited by the advectee). In other words, during its (long) correlation time $\tau _{\mathrm {nl}}^+$, the stronger field (in its reference frame travelling at $v_{\mathrm {A}}$) sees a weak field that has been rendered multiscale by the spatial variation of the stronger field, but remains approximately constant for a time $\tau _{\mathrm {nl}}^+$ and so can keep distorting the stronger field in a time-coherent way.

This argument can only work, it seems, if the long-term coherence of the weaker field is not upset by the way in which it is forced, so one must assume that it is forced at the outer scale with the same (long!) correlation time as the cascade time of the stronger field, in the Alfvénic frame of the latter. Chandran & Perez (Reference Chandran and Perez2019) dub this the ‘coherence assumption’. While hard to justify in general,Footnote ⁴⁴ it is great for them as, in their model, the weaker field is generated by the reflection of the stronger one as the latter propagates outwards in an expanding solar wind, and so one should indeed expect the two fields to be tightly correlated at the outer scale. Their endorsement of LGS07 – perhaps with an amendment that $\theta _\lambda$ should have some scaling with $\lambda$ determined by alignment/intermittency (§ 6.4.2) – isbacked up by their numerical results, where both fields’ spectra approach $k_{\perp }^{-3/2}$ (and their alignment increases) with increasing heliocentric distance.

9.4. Perez & Boldyrev (Reference Perez and Boldyrev2009)

Perez & Boldyrev (Reference Perez and Boldyrev2009) disagree with the entire approach leading to (9.3): they think that the two Elsasser fields should have two different alignment angles $\theta ^\pm _\lambda$, both small, and posit that those ought to be the angles that they make with the velocity field.Footnote ⁴⁵ Why that should be the case they do not explain, but if one takes their word for it, then (as is obvious from the geometry in figure 38)

(9.4)\begin{equation} \delta Z^+_\lambda\sin\theta^+_\lambda \sim \delta Z^-_\lambda\sin\theta^-_\lambda \quad\Rightarrow\quad \tau_{\mathrm{nl}}^+ \sim \tau_{\mathrm{nl}}^- \sim \frac{\lambda}{\delta Z^\pm_\lambda\sin\theta^\pm_\lambda} \quad\Rightarrow\quad \frac{\delta Z^+_\lambda}{\delta Z^-_\lambda} \sim \sqrt{\frac{\varepsilon^+}{\varepsilon^-}}. \end{equation}

The last result follows from the first relation in (9.3) with $\tau _{\mathrm {nl}}^+\sim \tau _{\mathrm {nl}}^-$. The equality of cascade times also conveniently spares them having to deal with the issue, discussed above, of long-time correlatedness, or otherwise, of the weaker field (or with $l_{\parallel }^+\neq l _{\parallel }^-$; see § 9.5).

Perez & Boldyrev (Reference Perez and Boldyrev2009, Reference Perez and Boldyrev2010a,Reference Perez and Boldyrevb) are not forthcoming with any detailed tests of this scheme (viz., either of the details of alignment or of the energy-ratio scaling), while Beresnyak & Lazarian (Reference Beresnyak and Lazarian2010) present numerical results that contradict very strongly the expectation of the energy ratio scaling as $\varepsilon ^+/\varepsilon ^-$ (as implied by (9.4)) and possibly support $(\varepsilon ^+/\varepsilon ^-)^2$ (i.e., (9.3)).Footnote ⁴⁶ Perez & Boldyrev (Reference Perez and Boldyrev2010b) reply that (9.4) should only be expected to hold for local fluctuating values of the amplitudes and of $\varepsilon ^\pm$ and not for their box averages. It is not impossible that this could make a difference for cases of weak imbalance ($\varepsilon ^+/\varepsilon ^-\sim 1$), with local fluctuations of energy fluxes superseding the overall imbalance, although it seems to me that if it does, we are basically dealing with balanced turbulence anyway: I do not see any fundamental physical difference between $\varepsilon ^+=\varepsilon ^-$ and $\varepsilon ^+\sim \varepsilon ^-$ on the level of ‘twiddle’ arguments by which everything is done in these theories. At strong imbalance, (9.3) seems to work better (Beresnyak & Lazarian Reference Beresnyak and Lazarian2009b, Reference Beresnyak and Lazarian2010) for the overall energy ratio, but not for spectra, which do not have the same slope (figure 18a). Perez & Boldyrev (Reference Perez and Boldyrev2010a,Reference Perez and Boldyrevb) argue that such cases in fact cannot be properly resolved, the limiting factor being the weaker field providing too slow a nonlinearity to compete with dissipation and produce a healthy inertial range. If so, the interesting case is inaccessible and the accessible case is uninteresting, we know nothing.

9.5. Parallel scales and two flavours of CB

By the CB conjecture (§ 5.1), the parallel coherence lengths of the two fields are, in the ‘naïve’ theory leading to (9.3),

(9.5)\begin{equation} l_{{\parallel}}^\pm \sim v_{\mathrm{A}}\tau_{\mathrm{nl}}^\pm \quad\Rightarrow\quad \frac{l_{{\parallel}}^+}{l_{{\parallel}}^-} \sim \frac{\varepsilon^+}{\varepsilon^-} > 1, \end{equation}

whereas in the Perez & Boldyrev (Reference Perez and Boldyrev2009) theory (9.4), the equality of cascade times implies $l_{\parallel }^+\sim l _{\parallel }^-$, end of story. LGS07 argue that, in fact, also (9.5) should be replaced by

(9.6)\begin{equation} l_{{\parallel}}^+ \sim l_{{\parallel}}^- \sim v_{\mathrm{A}}\tau_{\mathrm{nl}}^- \end{equation}

because $\boldsymbol {Z}_{\perp }^+$ perturbations separated by distance $l_{\parallel }^-$ in the parallel direction are advected by completely spatially decorrelated $\boldsymbol {Z}_{\perp }^-$ perturbations, which would then imprint their parallel coherence length on their stronger cousins (the parallel coherence length of the short-correlated advector imprinted on the advectee).

Furthermore, if one accepts the LGS07 argument that the correlation time of the $\boldsymbol {Z}_{\perp }^-$ field is $\tau _{\mathrm {nl}}^+$, not $\tau _{\mathrm {nl}}^-$ (see § 9.3), then $l_{\parallel }^-\sim v _{\mathrm {A}}\tau _{\mathrm {nl}}^-$ must be justified not by temporal (causal) decorrelation but by the weaker field being spatially distorted beyond recognition on the scale $l_{\parallel }^-$, even if remaining temporally coherent. This is more or less what Beresnyak & Lazarian (Reference Beresnyak and Lazarian2008) call ‘propagation CB’ (the other CB being ‘causality CB’). They note that the typical uncertainty in the parallel gradient of any fluctuating field at scale $\lambda$ is

(9.7)\begin{equation} \delta k_{{\parallel}}\sim\frac{\boldsymbol{b}_{{\perp}}\cdot\boldsymbol{\nabla}_{{\perp}}}{v_{\mathrm{A}}} \sim \frac{\delta b_\lambda}{\xi_\lambda v_{\mathrm{A}}}. \end{equation}

In balanced turbulence, $\delta k _{\parallel }^{-1}\sim v _{\mathrm {A}}\tau _{\mathrm {nl}} \sim l_{\parallel }$ (cf. (6.32)), so this is just a consistency check. In imbalanced turbulence,

(9.8)\begin{equation} \delta b_\lambda\sim\delta Z^+_\lambda \quad\Rightarrow\quad \delta k_{{\parallel}}^{{-}1} \sim \frac{\xi_\lambda v_{\mathrm{A}}}{\delta Z^+_\lambda} \sim v_{\mathrm{A}}\tau_{\mathrm{nl}}^-, \end{equation}

where $\tau _{\mathrm {nl}}^-$, given by (9.2), is the spatial-distortion time of $\delta Z^-_\lambda$, not necessarily its correlation time. The parallel scale of any field will be the shorter of $\delta k _{\parallel }^{-1}$ and whatever is implied by the causality CB. In the LGS07 theory, the latter is $v_{\mathrm {A}}\tau _{\mathrm {nl}}^+$ for both fields. Since $\tau _{\mathrm {nl}}^+\gg \tau _{\mathrm {nl}}^-$, we must set $l_{\parallel }^+ \sim l_{\parallel }^- \sim \delta k _{\parallel }^{-1}$, which is the same as (9.6).

Thus, we end up with both Elsasser fields having $\tau _{\mathrm {A}} \sim l_{\parallel }^-/v_{\mathrm {A}}$ that is smaller than their correlation time $\tau _{\mathrm {nl}}^+$ (even though the weaker field has a shorter spatial distortion time $\tau _{\mathrm {nl}}^-\sim \tau _{\mathrm {A}}$), but their cascades are nevertheless strong. Whatever you think of the merits of the above arguments, neither (9.5) nor (9.6) appear to be consistent with any of the cases reported by Beresnyak & Lazarian (Reference Beresnyak and Lazarian2009b), weakly or strongly imbalanced, which all have $l_{\parallel }^+< l_{\parallel }^-$ (see, e.g., figure 18b). No other numerical evidence on the parallel scales in imbalanced turbulence is, as far as I know, available in print.

9.6. Towards a new theory of imbalanced MHD turbulence

The Beresnyak & Lazarian (Reference Beresnyak and Lazarian2008) argument was, in fact, more complicated than presented in § 9.5, because they did not agree with LGS07 about the long correlation time of the weaker field, assumed the stronger field to be weakly, rather than strongly, turbulent, and were keen to accommodate $l_{\parallel }^+ < l_{\parallel }^-$. Their key innovation was to allow interactions to be non-local. I will not review their theory here, because it depends on a number of ad hoc choices that I do not know how to justify, and does not, as far as I can tell, lead to a fully satisfactory set of predictions, but I would like to seize on their idea of non-locality of interactions, although in a way that is somewhat different from theirs. The resulting scheme captures most of the properties of imbalanced turbulence observed in numerical simulations (§ 9.2) and reduces to the already established theory for the balanced case when $\varepsilon ^+/\varepsilon ^-\sim 1$, so perhaps it deserves at least some benefit of the doubt.

9.6.1. Two semi-local cascades

Let me assume a priori that, as suggested by numerics (Beresnyak & Lazarian Reference Beresnyak and Lazarian2009b; Beresnyak Reference Beresnyak2019), $l_{\parallel }^+\lll _{\parallel }^-$ in the inertial range, viz.,

(9.9)\begin{equation} \frac{l_{{\parallel}\lambda}^+}{l_{{\parallel}\lambda}^-} \sim \left(\frac{\lambda}{L_{{\perp}}}\right)^\alpha, \end{equation}

where $\alpha > 0$ and $L_{\perp }$ is the perpendicular outer scale (so the two Elsasser fields are assumed to have the same parallel correlation length, $L_{\parallel }$, at the outer scale – e.g., by being forced that way, as was done by Beresnyak Reference Beresnyak2019).

This implies that, at the same $\lambda$, the stronger field $\delta Z_\lambda ^+$ oscillates much faster than the weaker field $\delta Z_\lambda ^-$. I shall assume therefore that the interaction between the two fields local to the scale $\lambda$ is not efficient: even though $\delta Z_\lambda ^-$ is buffeted quite vigorously by the stronger field $\delta Z_\lambda ^+$, most of this cancels out. Rather than attempting to pick up a contribution arising for the resulting weak interaction, let me instead posit that the dominant, strong nonlinear distortion of $\delta Z_\lambda ^-$ will be due to the stronger field $\delta Z_{\lambda '}^+$ at a scale $\lambda ' > \lambda$ such that

(9.10)\begin{equation} l_{{\parallel}\lambda'}^+ \sim l_{{\parallel}\lambda}^-. \end{equation}

In other words, the interaction is non-local in $\lambda$ but local in $l_{\parallel }$.Footnote ⁴⁷ The constancy of the flux of the weaker field then requires

(9.11)\begin{equation} \frac{\left(\delta Z_\lambda^-\right)^2\delta Z_{\lambda'}^+}{\xi_{\lambda'}} \sim \varepsilon^-, \end{equation}

where $\xi _{\lambda '}$ has been introduced to account for a possible depletion of the nonlinearity due to alignment:

(9.12)\begin{equation} \frac{\xi_{\lambda'}}{L_{{\perp}}} \sim \left(\frac{\lambda'}{L_{{\perp}}}\right)^\beta. \end{equation}

In the absence of alignment, $\beta = 1$. For aligned, balanced, locally cascading ($\lambda ' \sim \lambda$) turbulence, $\beta = 3/4$ (see (6.31)). By the usual CB argument, the parallel coherence scale of the weaker field is

(9.13)\begin{equation} l_{{\parallel}\lambda}^- \sim \frac{v_{\mathrm{A}}\xi_{\lambda'}}{\delta Z_{\lambda'}^+}. \end{equation}

Note that, in the terminology of § 9.5, this is both the causality CB and the propagation CB, because $\delta k _{\parallel }$ for $\delta Z_\lambda ^-$ is determined by the propagation of the latter along the ‘local mean field’ $\delta b_{\lambda '}$ (see (9.7)).

Now consider the cascading of the stronger field by the weaker one. Since $l_{\parallel \lambda }^+ \ll l_{\parallel \lambda }^-$, the $\delta Z_\lambda ^-$ fluctuations are, from the point of view of the $\delta Z_\lambda ^+$ ones, slow and quasi-2D, and so the weaker field can cascade the strong one locally, in the same way as it does in any of the theories described above:

(9.14)\begin{equation} \frac{(\delta Z_\lambda^+)^2\delta Z_\lambda^-}{\xi_\lambda} \sim \varepsilon^+. \end{equation}

Causality CB would imply $l_{\parallel \lambda }^+\sim v_{\mathrm {A}}\xi _\lambda /\delta Z_\lambda ^-$, but that is long compared to $\delta k _{\parallel }^{-1}$ given by (9.7), so I shall use propagation CB instead, just like LGS07 and Beresnyak & Lazarian (Reference Beresnyak and Lazarian2008) did:

(9.15)\begin{equation} l_{{\parallel}\lambda}^+ \sim \frac{v_{\mathrm{A}}\xi_\lambda}{\delta Z_\lambda^+}. \end{equation}

Reassuringly, this choice immediately clicks into consistency with the requirement of parallel locality (9.10) if $l_{\parallel \lambda }^-$ is given by (9.13).

There are two nuances here. First, in order for the $\delta Z_\lambda ^-$ field to be able to distort $\delta Z_\lambda ^+$ according to (9.14), it needs to remain coherent for a time $\sim \xi _\lambda /\delta Z_\lambda ^-$. To make it do so, let me invoke the LGS07 argument already rehearsed in § 9.3: according to (9.11), $\delta Z_\lambda ^-$ stays coherent as long as $\delta Z_{\lambda '}^+$ does, which, according to (9.14) with $\lambda =\lambda '$, is $\xi _{\lambda '}/\delta Z_{\lambda '}^-$ – long enough!

Secondly, in (9.14), I used the same fluctuation-direction scale $\xi _\lambda$ as in (9.11), except at $\lambda$, rather than at $\lambda '$. This may be a somewhat simplistic treatment of alignment in local vs. non-local interactions, but I do not know how to do better, and the scalings that I get this way will have all the right properties. A reader who finds this unconvincing may assume $\xi _\lambda \sim \lambda$ and treat what follows as a GS95-style theory that ignores alignment altogether.

To summarise, I am considering here an imbalanced turbulence that consists of two ‘semi-local’ cascades: that of the stronger field, local in $\lambda$ but not in $l_{\parallel }$, and that of the weaker one, local in $l_{\parallel }$ but not in $\lambda$ (figure 20).

Figure 20. Cartoon of the spectra of imbalanced turbulence. The interactions are shown by arrows: red (advection of the weaker field by the stronger field, non-local in $\lambda$ but local in $l_{\parallel }$) and blue (advection of the stronger field by the weaker field, local in $\lambda$ but non-localin $l_{\parallel }$). The inset shows the parallel scales $l_{\parallel }^\pm$ vs. the perpendicular scale $\lambda$.

9.6.2. Perpendicular spectra

In view of (9.13) and (9.15), (9.11) can be rewritten as follows:

(9.16)\begin{equation} \frac{\left(\delta Z_\lambda^-\right)^2\delta Z_\lambda^+}{\xi_\lambda} \sim \varepsilon^-\frac{\xi_{\lambda'}\delta Z_\lambda^+}{\xi_\lambda\delta Z_{\lambda'}^+} \sim \varepsilon^-\frac{l_{{\parallel}\lambda}^-}{l_{{\parallel}\lambda}^+} \sim \varepsilon^-\left(\frac{\lambda}{L_{{\perp}}}\right)^{-\alpha}, \end{equation}

the last step being a recapitulation of the assumption (9.9). Dividing (9.14) by (9.16), one gets

(9.17)\begin{equation} \frac{\delta Z_\lambda^+}{\delta Z_\lambda^-} \sim \frac{\varepsilon^+}{\varepsilon^-}\left(\frac{\lambda}{L_{{\perp}}}\right)^\alpha. \end{equation}

Thus, the ratio of the energies at the outer scale ($\lambda =L_{\perp }$) is $(\varepsilon ^+/\varepsilon ^-)^2$, likely the correct scaling (see § 9.2, item (ii) and § 9.4), and the spectrum of the stronger field is steeper than that of the weaker field, also in agreement with numerics (§ 9.2, item (i)).

Now, by using (9.14), (9.17) and the alignment assumption (9.12), it becomes possible to determine the scalings of both fields:

(9.18)\begin{equation} \delta Z_\lambda^+ \sim \left[\frac{(\varepsilon^+)^2}{\varepsilon^-}L_{{\perp}}\right]^{1/3} \left(\frac{\lambda}{L_{{\perp}}}\right)^{(\beta + \alpha)/3},\qquad \delta Z_\lambda^- \sim \left[\frac{(\varepsilon^-)^2}{\varepsilon^+}L_{{\perp}}\right]^{1/3} \left(\frac{\lambda}{L_{{\perp}}}\right)^{(\beta - 2\alpha)/3}. \end{equation}

Comparing the first of these with (9.11), one can also work out how non-local the interactions are:

(9.19)\begin{equation} \frac{\lambda}{\lambda'} \sim \left(\frac{\lambda}{L_{{\perp}}}\right)^{3\alpha/(2\beta-\alpha)}. \end{equation}

With $\alpha =0$ and $\beta =1$ in (9.18), we are back to GS95 (§ 5.3), whereas with $\alpha = 0$ and $\beta =3/4$, we recover the aligned theory of § 6.4.

9.6.3. Parallel spectra

Now, from (9.15), (9.12), (9.18), and (9.9), the parallel scales are

(9.20)\begin{equation} \frac{l_{{\parallel}\lambda}^+}{L_{{\parallel}}} \sim \left(\frac{\lambda}{L_{{\perp}}}\right)^{(2\beta-\alpha)/3},\qquad \frac{l_{{\parallel}\lambda}^-}{L_{{\parallel}}} \sim \left(\frac{\lambda}{L_{{\perp}}}\right)^{2(\beta-2\alpha)/3}, \end{equation}

where the parallel outer scale is (cf. (6.16))

(9.21)\begin{equation} L_{{\parallel}} = v_{\mathrm{A}} L_{{\perp}}^{2/3}\left[\frac{(\varepsilon^+)^2}{\varepsilon^-}\right]^{{-}1/3}. \end{equation}

Combining (9.20) with (9.18) gives us the parallel scalings of the field increments:

(9.22)\begin{equation} \delta Z_{l_{{\parallel}}}^+ \sim \left[\frac{(\varepsilon^+)^2L_{{\parallel}}}{\varepsilon^-v_{\mathrm{A}}}\right]^{1/2} \left(\frac{l_{{\parallel}}}{L_{{\parallel}}}\right)^{(\beta + \alpha)/(2\beta - \alpha)},\qquad \delta Z_{l_{{\parallel}}}^- \sim \left(\frac{\varepsilon^-l_{{\parallel}}}{v_{\mathrm{A}}}\right)^{1/2}. \end{equation}

Whereas the stronger field's scaling is (for small $\alpha$, slightly) steeper than $l_{\parallel }^{1/2}$, the weaker one's is exactly that, corresponding to a $k_{\parallel }^{-2}$ spectrum, as is indeed seen in numerical simulations (§ 9.2, item (iv) and figure 19). This makes sense because the weaker field was assumed to have a local parallel cascade with the usual CB conjecture, so the standard arguments for its parallel spectrum given in § 5.2 remain valid.

9.6.4. Pinning

It turns out that it is possible to determine $\alpha$ by considering what happens at the dissipation scale(s). The dissipation cutoffs $\lambda _{\eta }^\pm$ for the two Elsasser fields can be worked out by balancing their fluxes with their dissipation rates:

(9.23)\begin{equation} \varepsilon^\pm \sim \frac{\nu + \eta}{\left(\lambda_{\eta}^\pm\right)^2}\left(\delta Z^\pm_{\lambda_{\eta}^\pm}\right)^2. \end{equation}

Using (9.18) to work out the field amplitudes at $\lambda _{\eta }^\pm$, one gets

(9.24)\begin{equation} \frac{\lambda_{\eta}^+}{L_{{\perp}}} \sim \left(\frac{\varepsilon^-}{\varepsilon^+}\widetilde{\mathrm{Re}}\right)^{{-}1.5/(3-\beta-\alpha)}, \qquad \frac{\lambda_{\eta}^-}{L_{{\perp}}} \sim \widetilde{\mathrm{Re}}^{{-}1.5/(3-\beta+2\alpha)}, \qquad \widetilde{\mathrm{Re}} = \frac{\delta Z^+_{L_{{\perp}}}L_{{\perp}}}{\nu + \eta}, \end{equation}

where, as before, $\widetilde {\mathrm {Re}}$ is the smaller of $\mathrm {Re}$ and $\mathrm {Rm}$. There are two possibilities: either $\lambda _{\eta }^+ < \lambda _{\eta }^-$ or $\lambda _{\eta }^+ > \lambda _{\eta }^-$. The first of these is, in fact, impossible: if the weaker field is cut off at $\lambda _{\eta }^-$, there is nothing to cascade the stronger field at $\lambda < \lambda _{\eta }^-$ (locally in $\lambda$, as I assumed in § 9.6.1). The second possibility is $\lambda _{\eta }^+ > \lambda _{\eta }^-$. Since the weaker field is cascaded by the stronger one non-locally, a self-consistent situation would be one in which $\lambda _{\eta }^+$ were the scale $\lambda '$ corresponding to $\lambda =\lambda _{\eta }^-$. Using (9.19), we must therefore ‘pin’ the dissipation scales together in the following way:

(9.25)\begin{equation} \frac{\lambda_{\eta}^+}{L_{{\perp}}}\sim\left(\frac{\lambda_{\eta}^-}{L_{{\perp}}}\right)^{2(\beta-2\alpha)/(2\beta-\alpha)} \quad\Rightarrow\quad 1 -\frac{2(\beta-2\alpha)(3-\beta-\alpha)}{(2\beta-\alpha)(3-\beta+2\alpha)} = \frac{\ln(\varepsilon^+{/}\varepsilon^-)}{\ln\widetilde{\mathrm{Re}}}. \end{equation}

Assuming $\alpha \ll 1$, we get

(9.26)\begin{equation} \alpha \approx\frac{2(3-\beta)\beta}{3(3+\beta)}\frac{\ln(\varepsilon^+{/}\varepsilon^-)}{\ln\widetilde{\mathrm{Re}}}. \end{equation}

Thus, indeed, $\alpha \to 0$ as $\widetilde {\mathrm {Re}}\to \infty$, but very slowly, with very large $\widetilde {\mathrm {Re}}$ needed to achieve a modicum of asymptoticity at larger imbalances. For the record, this implies, from (9.24),

(9.27)\begin{equation} \frac{\lambda_{\eta}^+}{L_{{\perp}}} \sim \widetilde{\mathrm{Re}}^{{-}1.5/(3-\beta)}\left(\frac{\varepsilon^+}{\varepsilon^-}\right)^{0.5(9+\beta)/(9-\beta^2)}, \qquad \frac{\lambda_{\eta}^-}{L_{{\perp}}} \sim \widetilde{\mathrm{Re}}^{{-}1.5/(3-\beta)}\left(\frac{\varepsilon^+}{\varepsilon^-}\right)^{2\beta/(9-\beta^2)}. \end{equation}

It actually does appear to be the case that $\lambda _{\eta }^+>\lambda _{\eta }^-$ in figure 18(a) (Beresnyak & Lazarian Reference Beresnyak and Lazarian2009b) and figure 35(c) (Meyrand et al. Reference Meyrand, Squire, Schekochihin and Dorland2021) – perhaps the strongest evidence that we have of the non-locality of imbalanced cascades.

While what I have proposed above is a kind of ‘pinning’, it is not the conventional ‘pinning’ that means equating the amplitudes of the two fields to each other at the dissipation scale – one of the tenets of the theory of weak imbalanced turbulence (§ 4.3). Indeed, (9.23) implies that the ratio of the Elsasser amplitudes at their respective dissipation scales is

(9.28)\begin{equation} \frac{\delta Z^+_{\lambda_{\eta}^+}}{\delta Z^-_{\lambda_{\eta}^-}} \sim \sqrt{\frac{\varepsilon^+}{\varepsilon^-}}\frac{\lambda_{\eta}^+}{\lambda_{\eta}^-} \sim \left(\frac{\varepsilon^+}{\varepsilon^-}\right)^{(3+\beta/2)/(3+\beta)}, \end{equation}

where I have used (9.27), valid in the limit $\widetilde {\mathrm {Re}}\to \infty$. Notably, the amplitude ratio is independent of $\widetilde {\mathrm {Re}}$ in this limit, but is not equal to $(\varepsilon ^+/\varepsilon ^-)^{1/2}$, as one might have concluded from (9.23) for $\lambda _{\eta }^+=\lambda _{\eta }^-$ (as Beresnyak & Lazarian Reference Beresnyak and Lazarian2008 did).

Arguably, this is a rather attractive theory: asymptotically, the spectra are parallel, interactions are local, etc., but in any finite-width inertial range, there are finite-$\widetilde {\mathrm {Re}}$ logarithmic corrections to scalings, locality, etc., accounting for all of the distinctive features of imbalanced turbulence seen in non-asymptotic simulations (§ 9.2).

10.1. Observational and numerical evidence

Going back to figure 16(a), we see that real MHD turbulence observed in the solar wind is distributed between cases with a local Elsasser imbalance (cross-helicity) and those with an Alfvénic one – the latter in favour of the magnetic field. Thus, the imbalanced cascades are only half of the story. According to the second relation in (B6), in imbalanced turbulence ($|\sigma _{\mathrm {c}}|\approx 1$), it is a geometric inevitability that $|\sigma _{\mathrm {r}}|\ll 1$, as illustrated by figure 16(a) and confirmed directly in the statistical study of solar-wind data by Bowen et al. (Reference Bowen, Mallet, Bonnell and Bale2018). In contrast, when the cross-helicity is not large (i.e., when $\sigma _{\mathrm {c}}$ is not close to $\pm 1$), there is flexibility for the perturbations to have finite residual energy: in the event, $\sigma _{\mathrm {r}}<0$. The definitive observational paper on this is Chen et al. (Reference Chen, Bale, Salem and Maruca2013), confirming negative $\sigma _{\mathrm {r}}$ over a large data set obtained in the solar wind. They also report that residual energy has a spectrum consistent with $k_{\perp }^{-2}$ or perhaps a little shallower, but certainly steeper than either the kinetic- or magnetic-energy spectra: the scalings of all three are reproduced in figure 21. This seems to be in agreement with earlier observational and numerical evidence (Müller & Grappin Reference Müller and Grappin2005; Boldyrev et al. Reference Boldyrev, Perez, Borovsky and Podesta2011, figure 22).

Figure 21. Spectra of magnetic (red), kinetic (blue), total (black) and residual (green) energies measured by Chen et al. (Reference Chen, Bale, Salem and Maruca2013) (figures taken from Chen Reference Chen2016): (a) typical spectra; (b) average spectral indices vs. normalised cross-helicity $\sigma _{\mathrm {c}}$ (defined in (B4)).

Figure 22. Spectra of total, magnetic, kinetic ((a), solid, dashed and dot-dashed lines, respectively, compensated by $k_{\perp }^{3/2}$) and residual energy ((b), compensated by $k_{\perp }^{1.9}$) in an RMHD simulation by Boldyrev et al. (Reference Boldyrev, Perez, Borovsky and Podesta2011) (© AAS, reproduced with permission). Beresnyak (Reference Beresnyak2014b) does a convergence study by rescaling to the dissipative cutoff and finds instead a $k_{\perp }^{-1.7}$ scaling to be the best fit to the residual-energy spectrum (see his figure 3), so it is possible that what is displayed here is a large-scale transient.

There are two ways in which a shade of legitimate doubt extends over both numerical and observational evidence quoted above.

(i) Beresnyak (Reference Beresnyak2014b), analysing his (largest-ever) simulations, reports that the residual energy at the small-scale end of its spectrum scales approximately as $k_{\perp }^{-1.7}$, same as his kinetic- and magnetic-energy spectra. His conclusion is that residual energy is merely a scale-independent finite fraction $\approx 0.15$ of the total energy.

(ii) Solar wind's expansion with heliocentric distance and the resulting reflection of the outward-propagating Elsasser field leads to an increase in the negative residual energy, which has nothing to do with nonlinear interactions in the locally homogeneous turbulence (see, e.g., Perez & Chandran Reference Perez and Chandran2013). The steepening of the magnetic field's spectrum compared to velocity's due to this effect is captured quite clearly in, e.g., the expanding-box simulations of Squire, Chandran & Meyrand (Reference Squire, Chandran and Meyrand2020), to whose paper I also refer the reader requiring further space-physics references on this subject. In assessing solar-wind measurements for evidence of residual energy in MHD turbulence, it may be non-trivial to separate this expansion effect from the organic local tendency of the nonlinear interactions to favour negative residual energy.

Alas, nothing is ever clear-cut in this subject, but let me press on. Obviously, it cannot be true at asymptotically small scales that, as the solar-wind data suggests, the magnetic- and kinetic-energy spectra scale as $k_{\perp }^{-5/3}$ and $k_{\perp }^{-3/2}$, respectively, while their difference scales as $k_{\perp }^{-2}$ – the $\boldsymbol {b}$ and $\boldsymbol {u}$ spectra must meet somewhere, as they indeed do in figure 21(a). The residual energy appears to peter out at the same scale (although that is also where the noise effects kick in) – but it would not be asymptotically impossible for it to retain the $k_{\perp }^{-2}$ scaling as a subdominant correction to approximately equipartitioned $\boldsymbol {b}$ and $\boldsymbol {u}$ spectra (as suggested by Boldyrev et al. Reference Boldyrev, Perez, Borovsky and Podesta2011) – or perhaps to parallel, finitely offset ones (à la Beresnyak Reference Beresnyak2014b). I will discuss a plausible origin for such a correction in § 10.4, but first some history.

10.2. Old theories

The first awakening of the MHD turbulence community to the turbulence's tendency for residual-energy generation dates back to the dawn of time (Pouquet, Frisch & Leorat Reference Pouquet, Frisch and Leorat1976; Grappin et al. Reference Grappin, Frisch, Pouquet and Leorat1982, Reference Grappin, Leorat and Pouquet1983), when theories and simulations based on isotropic EDQNMFootnote ⁴⁸ closure models of MHD turbulence predicted a negative residual energy (i.e., an excess of magnetic energy) scaling as a $k^{-2}$ correction to the dominant $k^{-3/2}$ IK spectrum (see § 2.2). While the isotropic IK theory certainly cannot be relevant to MHD turbulence with a strong mean field (see § 2), the modern evidence (§ 10.1) looks very much like those old results, with $k$ replaced by $k_{\perp }$. This led Müller & Grappin (Reference Müller and Grappin2005) to claim a degree of vindication for the EDQNM-based theory. This vindication cannot, however, be any stronger than the vindication of IK provided by Boldyrev's theory (§ 6.2) and its variants (§ 6.4): same scaling, different physics.

Below the turgid layers of EDQNM formalism, the basic physical idea (best summarised by Grappin, Müller & Verdini Reference Grappin, Müller and Verdini2016) is that residual energy is generated from the total energy by nonlinear interactions that favour magnetic-field production (the ‘dynamo effect’)Footnote ⁴⁹ and removed by the ‘Alfvén effect’, which tends to equalise $\boldsymbol {u}_{\perp }$ and $\boldsymbol {b}_{\perp }$ perturbations. A balance of these two effects leads to a prediction for the residual-energy spectrum in the form

(10.1)\begin{equation} E_{\mathrm{res}} \sim \frac{\tau_{\mathrm{A}}}{\tau_b} E \sim \left(\frac{\tau_{\mathrm{A}}}{\tau_{\mathrm{nl}}}\right)^\alpha E, \end{equation}

where $E$ is the total-energy spectrum, $\tau _b$ is the characteristic time scale of the generation of excess magnetic energy at a given scale, $\tau _{\mathrm {A}}$ and $\tau _{\mathrm {nl}}$ are our old friends Alfvén and nonlinear times, and the exponent $\alpha$ depends on one's theory of how $\tau _b$ is related to these two basic times. For example, in the IK theory, $\tau _b\sim \tau _{\mathrm {nl}}^2/\tau _{\mathrm {A}}$ (because IK turbulence is weak; cf. (4.5) and footnote 5), so $\alpha =2$. Using the IK scalings (2.5) and $\tau _{\mathrm {A}}/\tau _{\mathrm {nl}} \sim \delta u_\lambda /v_{\mathrm {A}}$, one then gets from (10.1)

(10.2)\begin{equation} E_{\mathrm{res}}(k) \sim\frac{\varepsilon}{v_{\mathrm{A}}} k^{{-}2}. \end{equation}

I know of no unique or obvious way of adjusting this promising (but necessarily wrong because IK-based) result to fit a critically balanced cascade: indeed, the CB requires $\tau _{\mathrm {A}}\sim \tau _{\mathrm {nl}}$, implying $E_{\mathrm {res}}\sim E$, i.e., a scale-independent ratio between the residual and total energy (this was also the conclusion of Gogoberidze et al. Reference Gogoberidze, Chapman and Hnat2012, who undertook the heroic but thankless task of constructing an EDQNM theory of anisotropic, critically balanced MHD turbulence). Neither solar wind nor MHD simulations appear to agree with this (§ 10.1).

Obviously, once we enter the realm of intermittent scalings of the kind described in § 6.4.1, i.e., allow the outer scale to matter, there is a whole family of possibilities admitted by the RMHD symmetry and dimensional analysis: by exactly the same argument as the one that led to (6.14) (and noting that spectrum $\sim$ amplitude$^2/k_{\perp }$), we must have

(10.3)\begin{equation} E_{\mathrm{res}}(k_{{\perp}}) \sim\varepsilon^{2(1+\delta)/3}\left(\frac{L_{{\parallel}}}{v_{\mathrm{A}}}\right)^{2\delta}k_{{\perp}}^{-(5-4\delta)/3}, \end{equation}

where $\delta$ is some new exponent. In order to determine it, one must input some physical or mathematical insight.