MHD Turbulence

We review the current status of research in MHD turbulence theory and numerical experiments and their applications to astrophysics and solar science. We introduce general tools for studying turbulence, basic turbulence models, MHD equations and their wave modes. Subsequently, we cover the theories and numerics of Alfvenic turbulence, imbalanced turbulence, small-scale dynamos and models and numerics for supersonic MHD turbulence.


Introduction
Turbulence is a time-dependent, stochastic flow commonly found in fluids with low viscosity. Shear viscosity is associated with microscopic phenomena, and the size of the turbulent system is much larger than the viscous scale. Turbulence develops from laminar flow due to instabilities and has many degrees of freedom. Despite its complexity, researchers investigate turbulence because of practical importance. The effect of turbulence is not only unpredictability of each realization of the flow, but often very important, quantifiable and predictable effects which are attractive to scientists and engineers. For example, ideal equations of motion, such as Euler's equation can be used, under certain conditions, to derive conservation laws. The conservation of energy and the conservation of the velocity circulation along the path frozen into the fluid (Kelvin's theorem) are notable examples of these "ideal invariants". However, physicists realized very early on that moving through the fluid involves drag and the loss of energy. Despite there is always a stationary ideal flow that produces zero drag (d'Alembert paradox), in practice such flows are not realizable due to instabilities and finite viscosity. Turbulence research elucidated this energy loss process and argued that it could happen for arbitrarily small viscosity due to the conserved quantity forming a "cascade" through scales finally dissipating on sufficiently small scales (Richardson-Kolmogorov picture). Likewise, Kelvin's circulation theorem is broken for flow around the wing, making possible lift force and the airplane flight.
Compared to turbulence on Earth, astrophysical turbulence is characterized by even larger scale separation between the problem size and the dissipative size, this makes turbulence in space almost unavoidable. Unlike the flows of non-conductive fluids on Earth, well-described by the Navier-Stokes equations, astrophysics deals with flows of ionized plasmas which, in most cases, can be considered perfectly conducting. These flows are described by equations with currents, magnetic fields and the Lorentz force, namely magnetohydrodynamic (MHD) equations.
In most astrophysical environments magnetic fields are observed and often are dynamically important. In our Galaxy, similar to other spiral galaxies, the magnetic field has a regular as well as random components. The value of the magnetic field, around 5 µG, suggests equipartition between magnetic and kinetic energies. In the galaxy clusters, the magnetic field is of order 1-3 µG, which is around 1/20th of the equipartition. Another example is the convective cell on the Sun. Its magnetic field is also somewhat close to the equipartition with the motions. Our Universe would have been very boring if it had both electric and magnetic charges, so that both electric and magnetic fields are screened out on large scales. Fortunately, this is not the case, the magnetic field and large-scale motions result in the acceleration of particles and the Universe if filled with non-thermal radiation in all wavebands.
Considering the space is filled with ionizing radiation it is not surprising that astrophysical plasmas are well-conductive. However, do they always have to be well-magnetized as well? The process of generation or amplification of the field is known as a dynamo, and this process seems to work sufficiently fast to do its job everywhere. If we start with zero magnetic field in the MHD equations, this produces precisely zero field in the future in an apparent contradiction with the ubiquity of magnetic fields. Do we always have to rely on primordial magnetic fields or the effects beyond simple MHD equations? In this review, among other things, we will emphasize that the growth of magnetic energy can be described in a framework somewhat similar to the loss of kinetic energy in the nearly ideal hydrodynamic flows. In other words, fast dynamo is an inherent property of turbulence.
Magnetic turbulence is also the primary cause of accretion onto gravitating objects, in particular accretion onto black holes is estimated to be the most potent source of energy in the Universe, exceeding thermonuclear burning in stars. Thin stationary accretion disks in a Keplerian potential are hydrodynamically stable, so in order to generate accretion one has to rely on the excitation of the the magnetic degree of freedom, the problem known as magnetorotational instability (MRI). Related to MRI-unstable disks are astrophysical jets, highly collimated flows perpendicular to the accretion disks in which magnetic field is essential in the process of launching and collimation of the flow.
Alfven theorem of perfectly conducting magnetohydrodynamics states that magnetic field lines are perfectly frozen into the conductive fluid, which places a severe restriction on the process of the so-called magnetic reconnection -the change of topology of the magnetic configuration by magnetic field lines crossing and moving through magnetic null. By way of restricting such change, Alfven theorem also precludes fast release of magnetic energy in highly conductive environments. This is in gross contradiction with high-energy phenomena above the solar surface known as X-ray flares. Again, turbulence comes to the rescue and allows for the radical breaking of the Alfven theorem, even in near-perfectly conducting fluids, very much like breaking of Kelvin's theorem of hydrodynamics.
In recent years the two traditional pillars of physics -the theory and the experiment has been complemented by a new method, numerical simulations. Numerics is valuable because it covers the gap between the real world, experimental data, and the idealized theory, a gap often being too great and precluding discovery. Numerics solves idealized equations directly, in this aspect it is similar to theory. On the other hand, numerics may be referred to as "numerical experiment", measurement of physical quantities without Figure 1: Faraday rotation measure maps of radio sources within a galaxy cluster. From [1,2]. The cluster electrons act as a foreground for the radio source. These maps indicate random magnetic fields of several µG in the cluster, changing on scales of 10-40 kpc.
invoking much of the assumptions or prerequisites. Compared to the real-life experiment, in numerics, it is easier to study idealized cases such as statistically homogeneous or statistically stationary turbulence for which the theory has something to say At the same time numerics reduces almost infinite space of theoretical ideas by weeding out theories which are incompatible with numerical measurements. In the studies of turbulence, the strength of numerics is manifested in the high statistical accuracy of the results, especially on small scales. Compared to the experimental measurements or observations which have high statistical and systematic uncertainties this helps to discriminate between theories and make quicker progress.
One type of numerics, direct numerical simulations (DNS) will be highlighted in this review. DNS refers to "fully resolved" numerical experiment, where numerics is very accurate and faithfully reproduce solutions of the original equations on all scales. On the other end, there are Implicit Large Eddy Simulations (ILES), calculations aiming to get the large-scale features of the flow correct without caring about the details of the dissipation in small-scale turbulence or shocks. These are very common is astrophysics, allowing to simulate large objects which are indeed out of reach of DNS, however as we Figure 2: Fluctuation of density in the interstellar medium (ISM), "big power law in the sky" from [3,4]. This 3D spectral density with slope of −11/3 correspond to Kolmogorov's 1D spectrum with slope of −5/3. Whether this power law, obtained from a variety of observations of different ISM components, is part of a single turbulent cascade is still an open question.
will show in the dynamo section this should be used with caution.
This, primarily theoretical, review mostly deals with homogeneous (although usually anisotropic) turbulence. Homogeneous models benefit from the opportunity to average quantities over volume, instead of averaging over ensemble (see more in Section 3). Few examples of important physical problems involving inhomogeneous turbulence: 1) large-scale dynamo, where inhomogeneity is required to break the statistical symmetry of turbulence, typically mirror symmetry, to produce large-scale magnetic fields, see also Section 9; 2) generation of imbalanced turbulence with a localized source of perturbations, see also Section 8; 3) large-scale dynamics of expanding solar wind; 4) magnetic shear as a driver of turbulence, see Section 11; 5) MRI, mentioned above. Very often inhomogeneous problems are treated with the scale-separation technique, where turbulence is described is some sort of "local box" approximation, within the box it is assumed homogeneous, but have overall driving, for example, shear boundary condition, as in the case of MRI.

MHD turbulence in astrophysics
Turbulence results from instabilities of large-scale fluid motions experiencing low friction forces. Dimensionless Reynolds number characterizes the relative importance of viscosity where L is the characteristic scale of the flow, often called "outer scale," e.g., the diameter of a jet, V is its velocity, and ν is fluid kinematic viscosity (in units of [L] 2 /[T ]). Likewise, one can introduce similar magnetic Reynolds number where η = c 2 /4πσ is magnetic diffusivity, c is a speed of light, and σ is a conductivity and Lundquist number where is Alfven speed, in units of velocity. Re, Re m and S are typically very large in astrophysics, meaning that viscous and resistive effects should be very small, the numbers of Figure 4: Density, shown in color, in a simulated accretion disk around a black hole, subject to MRI. The length unit is two gravitational radii of the black hole. From [6].
order 10 10 or larger are common. A notable caveat of this simple picture is that astrophysical plasmas are very often collisionless and the rigorous derivation of simple diffusive transport coefficients, such as Chapman-Enskog expansion, simply fails. So, our transport coefficients refer to some "effective" diffusivities, the physical meaning of which can be understood as follows. In molecular physics, the kinematic viscosity can be estimated as ul, a product of thermal speed and the mean free path. It is then clear that the Reynolds number is the ratio of the product of velocity and scale corresponding to macro-and microscales. A suitably chosen "effective" mean free path will allow estimating Re and the scale at which fluid motion transitions into the dissipative or dispersive regime, for example, the Kolmogorov scale that we introduce in Section 3. Such trickery works extremely well for fluid flows with turbulence or shocks and also the logic behind ILES. Throughout this review, the reader will see many examples of the so-called scale locality of turbulence in action, in particular, large-scale properties of the flow are insensitive to the diffusivities. One way to understand astrophysical turbulence is to understand the source of energy and how it is converted to turbulent motions. The biggest source of energy in the Universe is gravity. Turbulence can be driven by cosmological flows when gravity amplifies initially small density perturbations and cause structure formation. This is a very slow process, however, and typically dynamical times of voids are less then unity, in units of the age of the Universe, the dynamical times of filaments (superclusters) are of order unity, while dynamical times in the intra-cluster medium (ICM) of the galaxy clusters are of order 20, meaning they are expected to be turbulent. The size of a typical large galaxy cluster is of the order of several megaparsecs (Mpc), and the main source of turbulence is the infalls of large chunks of matter into it, so-called major mergers. While the direct evidence of kinetic turbulence motions in clusters is still scarce, the evidence of magnetic fields produced by dynamo action is available, see, e.g., Fig. 1. The source of turbulence in the plasma of ordinary galaxies, called interstellar medium (ISM), is likely multiple. At the present time, collisions with other galaxies are fairly rare. Several other mechanisms of driving turbulence can be identified, however: a) galactic disk is conductive and as such is subject to MRI, b) supernova explosions expand to the scale of several parsecs colliding with density inhomogeneities of the ISM, producing large-scale irregular motions, c) ISM turbulence is subject to Parker's instability when volumes of the ISM which are more magnetized and filled with more cosmic rays (CRs) are more buoyant compared to volumes poorly magnetized and scarce in CRs, so convective instability against the gravity of the disk ensues, d) CRs and starlight heat ISM at the same time cools itself by atomic and molecular emission on lower frequencies, and this produces thermal instability in the gas, e) jets and winds from young stars collide with ISM inhomogeneities producing irregular motions, f) at high redshifts also accretion and merger. The complexity of the ISM turbulence is rather overwhelming, and we refer to [8,9] for further reading. The evidence of turbulence present in the ISM was compiled from different sources by [3] and sometimes is referred to as a "Big power law in the sky", see Fig. 2. One amusing property of ISM turbulence is the large-scale dynamo, which produces the magnetic field on the scales of the disk, tens of kpc, while the outer scale of turbulence is only 10-100 pc. The evidence of large-scale magnetic fields in other galaxies is abundant, see, e.g. Fig. 3, while the evidence of fluctuating component of the magnetic field is mostly limited to our own Galaxy, due to the limited resolution of the observations. Another source of kinetic energy derived from gravity is jets from disks around black holes. This process is actually more efficient than thermonuclear burning in stars, in- Figure 6: Left: A sample of fast solar wind at a distance of 0.9 AU measured by the Helios 2 spacecraft. Right: Power density spectra of magnetic field fluctuations observed by Helios 2 between 0.3 and 1 AU. From [10].
cluding explosive burning in supernovae. The disk is conductive and unstable to MRI, turbulence in the disk creates rather large-scale structures, see, e.g., Fig. 4 and magnetic fields. The jet may be driven from the rotating black hole directly or may be driven centrifugally at larger distances by gas escaping along open magnetic field lines.
In the solar system, we can make in-situ measurements in the solar wind, the flow of tenuous magnetized plasma emitted from the Sun at speeds 400-800 km/s and propagating outwards to the boundaries of the solar system. Such direct measurements of the solar wind parameters and fluctuations in different regions from 0.3 to 5 AU distance to the Sun are especially valuable because they convey much more precise information about turbulent fluctuations compared to astrophysical observations of ISM and ICM mired by limited resolution and projection effects. Ion and electron counters on the satellite provide information about the flow, while magnetometers measure magnetic fields. Solar wind properties widely vary depending on the flow angle with respect to the ecliptic, see Fig. 5. The measurements by a single spacecraft represent time-sequence, demonstrating fluctuations on timescales from days to seconds, which can be Fourier-analyzed. The velocity of the solar wind is much larger than the local Alfvén speed of around 30 km/s so that the measurement can be interpreted as the spacial spectrum, see, e.g., Fig. 6 for spectra obtained from measurement by Helios 2 spacecraft. The f −1 part of the spectrum corresponds to the shot-noise statistics of features emitted by the Sun, while the f −5/3 part is the evidence of well dynamically evolved turbulence, the characteristic timescales on these scales are indeed shorter than the time of flight from the Sun.
The Sun's outer envelope transports energy to the surface by convection, also generating magnetic fields in the process. The magnetic field is distributed extremely unevenly on the surface reaching several kilogauss in sunspots. Sunspots are connected by magnetic arcs visualized by structure because hot plasma has high thermal conductivity along the field and low conductivity perpendicular to it. Interaction of strong magnetic flux tubes above the solar surface leads to magnetic reconnection which results in two spectacular phenomena: X-ray flares (see Fig. 7) and coronal mass ejections (CME). It is conjectured that reconnection and the release of magnetic energy are due to the thin current sheet at the intersection of flux tubes becomes unstable and generate turbulence, see Section 11.

Statistical description of turbulence
In this section, we will briefly introduce a statistical description. For a more in-depth review of this subject, we highly recommend the monograph by Monin and Yaglom [11]. While single realizations of turbulent flow are chaotic and unpredictable, there's some order. This order can mostly be described statistically, at the same time engineers and scientists are not interested in individual realizations, but rather in averaged quantities, such as an averaged lift of an airfoil. Turbulence is a volume-filling and persistent process, its realizations filling configuration space densely so that the statistical ensemble measure- Figure 8: Simulated MHD turbulence visualized by magnetic field magnitude shown in grayscale. This statistically homogeneous, isotropic turbulence with zero net magnetic flux was driven by volumetric force to statistically stationary state, one snapshot of which is shown on the picture. ments sometimes can be replaced with time-and volume-averaging (ergodic hypothesis.) The theory relies typically on ensemble averaging, but numerical experiments mostly use volume and time averaging, see Fig. 8. We will designate averaging <> without specifying whether its statistical, time or volume averaging.
The spatial variability of a physical variable, e.g., v(r), over some scale l can be described as some function of the difference of v between points separated by a distance l. Second-order statistics can be related to energy, for example second order structure function (SF) of velocity, In the limit of large l this equals to four times kinetic energy, while for smaller l it is four times "characteristic energy" on this scale and all smaller scales. Fourier-transformed SF can be related to "energy spectrum" E(k) (see below). The energy spectrum is the energy distributed in wavenumber space, with dE = E(k)dk being the energy at a particular wavenumber and E(k)dk being total energy. When turbulence is statistically self-similar we expect a power-law scaling of statistical quantities, e.g., E(k). The SF above represents the sum of the longitudinal and transverse components of the velocity with respect to direction perpendicular and parallel to l. Naturally, longitudinal and transverse functions can be calculated separately. The longitudinal SF is historically important in experimental research of hydrodynamic turbulence due to being the primary quantity measured by the heated wire technique. The large-scale flow v 0 around the wire carries smaller fluctuations which cause fluctuations of the absolute value of v, which is what is measured by the changing resistance of the wire, however, fluctuations in |v| are mostly due to fluctuations parallel to the average v 0 . The Taylor hypothesis assumes that time variations correspond to variations in space, i.e. measurements separated by t correspond to l = v 0 δt. Thus we measure only the component parallel to l. In the solar wind measurements, all three vector components are recovered so that the transverse, longitudinal and full structure functions can be calculated.
In the case of isotropic turbulence SF(l) is only a function of l, MHD turbulence is not isotropic, however, so there is a wider variety of structure functions that we can measure. However, as we show below, in the reduced MHD limit there is a particular structure function which plays the similar role as the isotropic SF in hydrodynamics, the perpendicular SF where n is a vector perpendicular to the magnetic field. The turbulent quantity u(r) can be Fourier-transformed: with the square of the transform called power spectrum: This function can be integrated over the sphere in k-space, e.g. if F (k) depends only on the magnitude of k we have E(k) = 4πk 2 F (k), the resulting quantity we will call three-dimensional spectrum. Similar procedure is possible when sampling the field along the line, i.e. in one dimension, this quantity will be called a one-dimensional spectrum E 1 (k) = 2F 1 (k). Note that 2 comes from F(k) being defined for positive and negative wavenumbers. E(k), E 1 (k) can be related in isotropic case by while the above mentioned parallel one-dimensional spectrum, which we designate as E (k) where only parallel component of velocity is used, for a solenoidal isotropic field: In practice spectras are never exact power laws so the shape of these spectra are different Fig. 9 shows three types of spectra from a simulation of MHD turbulence. Spectra and structure functions have one-to-one correspondence by Fourier transforms: If the spectrum has a power-law dependence k α , then by substitution k = x/r we obtain provided that remaining dimensionless integral converge. This relation is satisfied for α between −3 and −1.
From statistical viewpoint turbulence self-similarity, i.e. the assumption that turbulence have a single-fractal structure would mean that for structure functions of arbitrary orders n and m one can write: Some exact relations for structure functions in turbulence are known for hydrodynamics and MHD, which helps to test numerics. In the subsequent Section, we explain in more detail the concept of the inertial range -a range of scales where energy is being overall conserved and is being transferred from one scale to another. From the dynamical viewpoint, these are scales at which dissipation term can be ignored, and the energy is only injected from large-scale motions but not from an external force.
The Kolmogorov −4/5 law relates a parallel signed structure function for velocity in the inertial range with the turbulent dissipation rate: Another exact relation, similar to the Yaglom's -4/3 law for incompressible hydro exists for axially symmetric MHD turbulence: where l is taken perpendicular to the axis of statistical symmetry -the direction of the mean magnetic field B ( [12,13]). The testing of numerics involves measuring SFs and comparing them with theoretical predictions which help to establish which part of the spectrum is the inertial range, and which scales are dissipative and driving scales. The inertial range in a simulation is often defined as a range of scales where −SF 3 /l is closest to its theoretical value, i.e., where the influence of energy injection from driving and energy dissipation from the viscous term is minimized. Fig. 10 shows several structure functions, compensated by various powers of l. The ratio of different structure functions can test turbulence self-similarity. If this ratio is dimensionless, it is supposed to be constant through scales. For the test of self-similarity on Fig. 10 we show the ratio of parallel third order structure function and full third order SF, SF 3 = |v(r − l) − v(r)| 3 . Fig. 10 shows that hydrodynamic turbulence is rather self-similar at the same time the scaling of the second-order structure function in the inertial range is around l 0.7 , i.e., close to the Kolmogorov scaling (see next section).

Kolmogorov cascade model
In hydrodynamic turbulence, a useful starting point is the Kolmogorov model [14] for incompressible turbulence. The incompressible case has constant density so that the energy dissipation can be defined per unit mass and assumed statistically homogeneous as well. This quantity, ǫ have units of cm 2 /s 3 and plays a crucial role in many situations and will be used plenty through this review. The Kolmogorov model assumes that the statistical properties of turbulence are uniquely determined by the amount of energy available in this stationary homogeneous system, i.e., by the ǫ alone. Furthermore, it is argued that the energy self-similarly cascades through the series of scales known as the inertial range. Cascade means that the energy is being transferred from one scale to another without dissipation.
The dimensional derivation of Kolmogorov scaling involves noting that the spectrum defined in the previous section have units of cm 3 /s 2 and the wavenumber have units of cm −1 , so that where C K is a dimensionless Kolmogorov constant. From the previous Section we also know that the 3D spectrum The hand-waving derivation involves introducing "characteristic velocity on scale l" u l and imagining that the energy rate is constant for all scales: where t casc is the "cascading timescale" the time it takes for nonlinearity to remove energy from scale l and transfer it to smaller scales. It is further assumed that in the hydrodynamic cascade t casc is a dynamical time on each particular scale, i.e. t casc ≈ l/u l , which results in From the definition of the spectrum and its relation to the SF, we can argue that E(k)k ∼ u 2 l so that the two formula for Kolmogorov scaling agree. A compilation of experimental results for hydrodynamic turbulence [15] suggests that a Kolmogorov constant C K is universal for a wide variety of flows. High-resolution numerical simulations of isotropic incompressible hydrodynamic turbulence, see Fig. 11 and [16] suggest the value around 1.6. See Fig. 12 for a spectrum from 1024 3 simulation from [17].
The expression Eq. 19 can be written for the largest scale in the system, the outer scale: This can be regarded as scaling with the outer scale velocity v, and/or the scale of the system L. Three different things can be done to study this law, known as the zeroth law of turbulence, empirically. One can scale an experimental apparatus from L to L ′ , one can increase or decrease velocity, and one can change the fluid to vary viscosity. From the symmetries of the hydrodynamic equations, we know that the only real change would be a change in Reynolds number. The same type of turbulent flow results in approximately the same dimensionless coefficient of C ′ K . In the systems that generate turbulence easily, e.g., flow past the grid the above expression is reasonably precise for Re > 200. Note that in the statistically stationary case ǫ is also the energy dissipation rate, which happens on small scales due to viscosity. This fact illustrates that the outer scale L the dissipative scale only know each other through ǫ.
Outer scale is sometimes formally defined through integral over the spectrum, e.g.
. Usually, this is around a scale where energy is injected into the system. Inverse cascade of energy in two-dimensional hydrodynamic turbulence is one counter-example of this..
The energy "cascades" down to smaller scales until it hits the so-called Kolmogorov scale, where dissipative processes overcome nonlinear transfer of energy. The Kolmogorov scale can be expressed as a combination of viscosity/diffusivity and energy dissipation rate, which gives a unit of length.
where n is the order of the viscosity, e.g. n=2 for classic molecular viscosity, ν n is the value of the diffusivity, so that we obtain Navier-Stokes equation by adding to the RHS of the Euler equation the dissipation operator −ν n (−∇ 2 ) n/2 . Dimensionless ratio L/η could serve as a "length of the inertial range", although in practice spectrum is around an order of magnitude shorter.
Criticism of the Kolmogorov model points to the fact that the assumption of selfsimilarity is quite arbitrary and points to the examples of turbulence which are notably not self-similar. In the three-dimensional hydrodynamic turbulence deviations from selfsimilarity in the second-order measurements, such as energy spectrum, are fairly small, however. The more precise formula can be obtained by multiplying RHS of Eq. 17 by the "intermittency correction" (kL) α , where α ≈ 0.035. For more details see [18].

Lagrangian spectrum
Lagrangian measurements are performed by following a fluid element. Lagrangian viewpoint offers a simpler conceptual picture as the model above is conceptually simpler in Lagrangian formulation. The Euler's equation is a third Newton's law for the fluid element: Here D/Dt is the advective (Lagrangian) derivative corresponding to changes of a fluid element's properties over time, The work per unit mass, done upon a fluid element by pressure of surrounding fluid elements will be expressed, therefore as v · dv/dt. The Kolmogorov theory would therefore assume that, given a characteristic time interval τ , the work done per unit mass upon a fluid element during this interval, δv τ · δv τ /τ , will be constant when τ corresponds to inertial-range timescales and equal to the turbulence energy cascade rate per unit mass ǫ. Formally, in stationary turbulence the second-order Lagrangian structure function of velocity should satisfy: in the inertial range, where v(t) is a velocity as a function of time for a given fluid element. This time structure function will correspond to the frequency spectrum of see Eq. (13). This first appeared in the texbook [19] and also in [20,21]. The scaling ω −2 and the fact that the energy spectrum is proportional to energy injection rate ǫ appear to be conceptually simpler than ǫ 2/3 scaling of the standard Eulerian Kolmogorov scaling. This spectrum has a dissipation timescale associated with the lifetime of critically damped eddies, also called the Kolmogorov timescale (for n=2): this is the location of the dissipative cutoff in Lagrangian spectrum. The direct measurement of the Lagrangian frequency spectrum is fairly challenging, however, as the probe has to be embedded in the flow. Temporal measurement of spectra from a wind tunnel or channel flow (e.g. [22]) does not correspond to the Lagrangian spectrum but can be connected, by Taylor hypothesis, to spatial spectrum (see Section 3). Below, in Section 7.3 we explain how a parallel spectrum in MHD turbulence can act as a surrogate of the Lagrangian spectrum.

More general Kolmogorov phenomenology
More general phenomenology is possible assuming cascading time scale relating to the dynamical timescale as in which case using Eq. 18 we get so that, assuming self-similarity, the spectral slope will be −1 − 2(1 − α)/3 = −5/3 + 2α/3 (Eq. 13). This will also result in a different Kolmogorov scale which we get by equating cascading time above and the viscous time l n /ν n : This reduces to Eq. 22 for α = 0. Alternatively, one can also assume that interaction is reduced by u l /c s : where c s is the sound speed. This model is called acoustic/wave turbulence and gives the −3/2 spectral slope.

Scaling convergence in turbulence: numerics and experiments
works similarly in the compared experiments. The discretization error and other numerical inaccuracies of statistically averaged quantities should depend only on the ratio of Kolmogorov scale to the grid scale provided that the timestep is also determined by the grid scale. So we need to keep the grid scale as the fixed fraction of the Kolmogorov scale. Keep in mind, that Kolmogorov scale itself is determined based on particular phenomenology of the cascade (see Section 4.2) and may not be known apriori. In this case, rigorous scaling convergence would require going through available hypotheses and checking each in turn.
We express the spectra of several simulations in dimensionless units corresponding to the expected scaling, for example, a E(k)k 5/3 ǫ −2/3 for the Kolmogorov model and plot it versus dimensionless wavenumber kη, where dissipation scale η again, corresponds to the same phenomenology. On the plot, the two spectra should collapse onto the same curve on the viscous scales, as long as the model works. The method has been used extensively in hydrodynamics [23,16,24] with great success. In numerics, it is especially efficient since, while experimental data may suffer from systematic uncertainties, numerics does not, and it collects tremendously large statistics on small scales, driving statistical error virtually to zero. Let us understand why this is the case. If we refer to the Kolmogorov cascade picture, described above, the energy cascade is local in scale and the only information that is being transferred from large scales to small scales is the local cascade rate ǫ. Now, assuming that at each scale, each eddy is independently created, its energy content on this scale should only depend on ǫ as ǫ 2/3 . So if we normalize the measurement by ǫ −2/3 , each eddy will represent, presumably, independent estimate of such normalized energy content at each scale. Given a characteristic eddy scale l, the number of eddies in a datacube goes as l −3 , while the number of correlation timescales for strong turbulence goes as l −2/3 , so the statistical error due to volume and time-averaging should decrease as l −11/6 . The plotted normalized spectrum I(x) = I(kη) = E(k)ǫ −2/3 k 5/3 should be "pinned" on the dissipation scale, because it should satisfy The precision of the convergence method was demonstrated in [24], where 4096 3 simulations allowed to capture the intermittency correction, which is a correction of −0.04 to the −5/3 spectral slope.

MHD equations, modes
Below we write ideal MHD equations that describe perfectly conducting, inviscid fluid. It should be kept in mind that solving ideal equations often require including special treatment of shocks and turbulence. Figure 12: Spectrum of hydrodynamic turbulence compensated by Kolmogorov scaling to give approximate constant function vs. wavenumber. Statistics from smaller datacube should largely repeat statistics from larger datacubes, as we see on collapsing spectra on the left and visually represented on the right.
with current j = ∇ × B and vorticity ω = ∇ × v, P (ρ, s) is an equation of state.
In the ideal case, specific entropy s is decoupled from the rest of the equations and can be described as a passive scalar. Here we use Heaviside units, redefining electric charge with a factor of 1/4π and getting rid of 4π factors in Maxwell's equations.
Introducing sound speed c 2 s = ∂P/∂ρ, linearized MHD equations reveal four perturbation modes: 1) Alfvén mode -transverse waves with v and B perturbations along k × B and dis- while its group velocity ∂ω/∂k = ±v A , hence the term Alfvén velocity. 2,3) Fast and slow modes -compressible waves with perturbations in the k, B plane propagating correspondingly faster and slower than v A , with the dispersion relation 4) Entropy mode -non-propagating passive scalar perturbations of specific entropy [25,12].
In this section, we will skip eigenvectors for the three modes for brevity and write them out in Section 10.

Numerical methods to simulate MHD turbulence
Several tools are available to simulate turbulence numerically:

Pseudospectral codes
The pseudospectral code solves MHD equations as a set of ordinary differential equations in time for each spacial Fourier harmonic. The coupling of harmonics is through the nonlinear term which is calculated in real space (hence "pseudo") and then converted back to Fourier space. Since the pseudospectral method approximates derivatives non-locally, using all data points, it does not suffer from dispersion error. Also, if some care is taken with timestep integration, for example, a symplectic integrator is used, it also preserves energy, i.e., does not suffer from dissipation error. The explicit dissipation, e.g. viscosity or resistivity are done with simple algebraic operations in Fourier space and can be made unconditionally stable, irrespective of the time step with no numerical expense. A typical pseudospectral code have symplectic integrator, corrects for aliasing error and has explicit dissipation in the form a t+∆t = a t exp(−νk 2 ∆t). Aliasing error comes from frequencies above 2/3 or below −2/3 of the Nyquist frequency if the nonlinear term is second order. E.g., if the Nyquist frequency is π, and keeping frequencies within [−2/3π, 2/3π], the sum or difference will still be within the interval modulo 2π (e.g. 2/3π + 2/3π − 2π = −2/3π). Another advantage of pseudospectral code for the incompressible case is that divergencefree condition for velocity and magnetic field can be done with simple algebraic operations in Fourier space. The spacial reconstruction uses all Fourier harmonics. As a result, the method's precision increases exponentially with the number of points in one dimension. This makes it practical to do "fully resolved" simulations, i.e., when the viscosity and magnetic diffusivities are explicit, and all scales of interests are represented with reasonable precision. The usual rule of thumb for well-resolved simulation in a periodic box with size 2π and number of points N in one direction, when the wavenumbers are represented by integers [−N/2 + 1, ..., 0, 1, ..., N/2] is k max η 1, where k max = N/3 for a 2/3 dealiased code. The main disadvantage is difficulty in introducing arbitrary boundary conditions. This method's periodic box comes naturally when we try to simulate homogeneous isotropic turbulence, however. One example of a publicly available pseudospectral code is Snoopy: http://ipag.osug.fr/~lesurg/snoopy.html.

Finite difference codes
Finite difference codes estimate derivatives by finite differencing. The precision of the code increases, typically, as a power law with the number of points, the index of the power law is the "order" of the code. The main advantage is simplicity and numerical speed. Disadvantages include special treatment of shocks with "shock viscosity". High order finite difference codes with explicit diffusivities can be rather precise in simulating turbulence. The divergence-free condition can be kept with "divergence cleaning" or with equations formulated in terms of magnetic potential. Pencil code is popular publicly avalable high order finite difference code: https://github.com/pencil-code

Finite volume codes with Riemann solvers
Also known as Godunov codes. Finite volume codes keep values of cell averages then reconstruct (interpolate) values on the interface of the cells both from the right and from the left. Thus the interface value is discontinuous and may be evolved for a short time as a "Riemann problem" -initial value problem with a single discontinuity. The time-average fluxes of conserved quantities through the interface are then computed from, typically, the approximate solution of the Riemann problem by the "Riemann solver". Finally, fluxes are used to advance cell averages in time. The inherent ability to describe discontinuities makes Godunov codes very robust and a code of choice to simulate supersonic turbulence. See also documentation of the publicly available code Athena++: https://github.com/PrincetonUniversity/athena-public-version

Lagrangian codes
Lagrangian codes refer to codes which use grid or material elements that are moving with the fluid. These include N-body codes (e.g., collisionless particles moving under the action of gravity), smooth particle hydrodynamics (SPH) codes -particle codes simulating hydrodynamics, moving mesh codes, among which purely Lagrangian (mesh moving with the fluid) or arbitrary Lagrangian-Eulerian (ALE, mesh moving in an arbitrary way). These codes are often used to simulate collapse under gravitational forces, but less common to simulate turbulence. Publicly avalable code Gadget2: https://wwwmpa.mpa-garching.mpg.de/gadget/

Theory of Alfvenic Turbulence
In Section 4 we introduced the standard Kolmogorov description of the inertial range of incompressible hydrodynamic turbulence. It is clear, however, that this picture is not applicable to MHD. Turbulence spectra are typically steeper than k −1 meaning that RMS fields are dominated by large scales. In hydrodynamics, however, large-scale velocity can be nullified by an appropriate choice of the reference frame. In MHD large-scale magnetic field cannot be nullified and will be dynamically important on all scales, including very small scales. This combination of sizable RMS large-scale field and small-scale fluctuations of the fields is the main difference from hydrodynamics. Also known as a "strong field limit", it was pointed out by Iroshnikov and Kraichnan [26,27] and it was suggested that inertial-range MHD turbulence is weak turbulence. Here weak turbulence refers to the picture of wave turbulence where wave packets propagate almost freely, and collision between waves leads to the small perturbation in their structure so that the perturbation theory is applicable [28]. The interaction of wave packets in MHD, however, is very different from the collision of sound waves. Introducing wavevector components parallel and perpendicular to the mean field, k and k ⊥ we see that the wave frequency ω = k v A depends only on k . This anisotropic dispersion relation results in anisotropic turbulence. The subsequent analytic work demonstrated that MHD turbulence tends to become stronger and not weaker during the cascade [29], as we will show below. We will also show that the Alfvenic part of MHD perturbations governs this highly anisotropic turbulence, hence the terms "Alfvenic Turbulence".
The rationale of working with simplified incompressible equations is similar to hydrodynamics. Assuming that a) turbulence have no shocks, b) no sizable energy is carried by sound waves (in MHD case, fast MHD mode), c) the Mach number M s = V L /c s is small, we can argue that a scale-wise Mach number M s = δv/c s should also be small and decrease with scale. The fluid compressibility will, therefore, be small in the inertial range.
Incompressible MHD equations consist of two dynamical equations and two constraints: Here we normalized magnetic field to velocity units in the same manner as in previous The dynamical equations are known as the momentum equation and the induction equation. The induction equation honors the divergence-free constraint for the magnetic field, this effectively results in no new constraint, if the initial condition is chosen divergence-free. The divergence-free constraint for velocity is satisfied by the appropriate choice of the scalar function P ′ = P + ρv 2 /2. The pressure, therefore, is a dummy variable.
Introducing solenoidal projectionŜ = (1 − ∇∆ −1 ∇) we can rewrite the equations without explicit constraints: Very useful change of variables to the Elsässer variables w ± = v ± b makes these equations even more compact: We introduce total energy, 1/2 (v 2 + b 2 ) dr and cross-helicity v · b dr which are conserved in this incompressible formulation. By taking the sum and difference of these quantities, we obtain conservation of each of Elsasser energies 1/2 (w ± ) 2 dr.

From weak to strong turbulence
Keeping in mind the above argument of a sizable mean field let us explicitly write it down as the constant field v A in a given volume, and perturbations as δw ± = w ± v A : Let us denote and ⊥ as directions parallel and perpendicular to v A and the subscript to the vector means projection to v A or the perpendicular plane, respectively.
In the limit of small δw's they represent perturbations, propagating along B or in the opposite direction, with the nonlinear term describing their interaction. Note that "self-interaction" of δw + or δw − is absent, both being an exact solution in the absence of another. The dominant nonlinear interaction is a three-wave process, so writing the dispersion relation and the conservation laws for energy and momentum, we see that one of the ω n must be zero. Let us choose ω 3 = 0, this means |k 1 | = |k 2 |, but there's no restrictions on k ⊥1,2 . The cascade preserves frequencies and goes forward by increasing only k ⊥ .
In wave turbulence theory the interaction strength ξ is the ratio of the nonlinear shear rate k ⊥ δw to the wave frequency k v A , it describes a fractional perturbation during one wave period: It is also the estimate of the ratio of the nonlinear term to the mean-field term in Eq. 46. In MHD turbulence the dynamical timescale τ dyn = 1/k ⊥ δw does not have to be proportional to the cascade timescale as in hydrodynamic turbulence. Instead, τ casc is increased by a factor of 1/ξ. This can also be understood in terms of perturbations of a wave packet being a random walk. Each individual perturbation is ξ strong, so it takes (1/ξ) 2 steps to destroy the wavepacket completely: The energy cascade rate is the energy on each scale divided by the cascade time on this scale. This rate is expected to be constant through scales and we designate it ǫ: Note k here is constant, so the phenomenological cascade spectrum is determined by . This argument can be followed rigorously by perturbation collision integral approach, used in wave turbulence [28] and solved exactly by Zakharov transformation, which was accomplished in [29,30].
One consequence of this solution is that turbulence grows anisotropic, with k ⊥ /k ∼ k ⊥ . Interestingly, it becomes stronger and not weaker on smaller scales, in other words, ξ is an increasing function of k ⊥ . Indeed, if we maintain k constant, this will result in Our two conclusions from this simple perturbation theory is that: a) the resonance condition results in a "perpendicular cascade", making MHD turbulence anisotropic, b) turbulence becomes stronger along the cascade until ξ ∼ 1. One can wonder if weak MHD turbulence is ever realized in nature. We can hypothesize that this is the case in astrophysical objects where the strong magnetic field is anchored in a heavy object, i.e., a star and is extended into the magnetosphere where perturbations of the field are much smaller than this anchored field. The empirical evidence for this case and specifically for the k −2 perpendicular spectrum is lacking, however. One can argue that large-scale dynamo (which we consider in subsequent Sections) can generate a mean field which is much stronger than perturbations, but empirically we know from the ISM observations that they are of the same order. This results in MHD turbulence being strong on the outer scale. Equation 46 can be further simplified assuming anisotropy k ⊥ ≫ k and the fact that δw ≪ v A . This allows to neglect parallel gradients in the nonlinear term, indeed, the mean field term with the parallel gradient (v A ∇ )δw ± is always much larger than similar contribution from the nonlinear term, (δw ∓ ∇ )δw ± and the latter could be ignored. So the three vector components of equation 46 are split into interdependent equations for the scalar δw ± and vector δw ± ⊥ :

Reduced MHD approximation
Note that Equation 55 depends on Eq. 56, but not vice-versa. Since Equation 55 represent passive dynamics and does not have essential nonlinearity, the nonlinear cascade is completely governed by Eq. 56. This latter equation is known as reduced MHD. In this anisotropic limit, the δw ± ⊥ is purely the Alfvén mode, and δw ± is the amplitude of the slow mode. Turbulence in Eq. 56 is called Alfvénic turbulence.
Slow mode for δw + is a passive scalar to δw − ⊥ and vice versa. If Alfvén and slow modes will be injected similarly from large scales, they will have the same statistics. In practice, the slow mode content can be determined from numerics.
It turns out, reduced MHD is more general than incompressible MHD and can be used beyond collisional fluid description. Alfvénic perturbations are transverse and rely only on the tension of the magnetic field line as a restoring force, the charged particles tied to this magnetic field line provide inertia. The [E × B] drift waves with wavelengths much smaller than the ion skin depth are indeed just Alfvén waves, and they exist regardless of the collisionality of the plasma [31], which is useful for the description of the collisionless solar wind. The anisotropy of MHD turbulence has been known empirically for a while, and RMHD had been formulated for perturbations in plasma in strongly magnetized case some time ago [32,33]. Since RMHD motions do not require plasma pressure we assume the results that we find for Alfvenic turbulence in this section do not depend on the ratio of plasma pressure to magnetic pressure "β", despite we started the derivation assuming infinite β.
Introducing parallel length Λ = 2π/k and perpendicular length λ = 2π/k ⊥ we see that reduced MHD has a two-parametric symmetry: A and B are arbitrary parameters of the transformation. This is the same symmetry as in hydrodynamics, except for the parallel scale Λ transforms similar to time, not to length. Λ being similar to time is very important and leads to analogies between dynamics in time and parallel structure in space as we show below. MHD equations do not have such symmetry, so Kolmogorov self-similarity arguments, technically, can not be applied to the MHD case. In practice, this regime for MHD can be achieved within the inertial range where δw ≪ v A condition and anisotropy condition are satisfied. In numerics, it is challenging to reach these universal dynamics directly from isotropic scales with δw ∼ v A . Instead, one can directly solve RMHD equations. As a practical comment, the statistics from the full MHD with δw ± ∼ 0.1v A is very close of that one of RMHD, see [34]. Another symmetry is evident in RMHD, related to the value of v A , The equations are unchanged under transformation v A → v A A, Λ → ΛA. The parallel scale and the Alfvén speed can be rescaled simultaneously without changing the dynamics.

Strong MHD turbulence
As we demonstrated above in Section 7.1, the perpendicular cascade will result in the growth of ξ and will naturally lead to strong turbulence, with ξ ∼ 1. Goldreich and Sridhar [35] proposed that the growth of ξ will be limited by the uncertainty relation between the cascading timescale and the wave-packet frequency, namely that the cascade time cannot be shorter than the wave period: τ casc ω ≥ 1. Using Eq. 52 we get ξ ≤ 1. This will make ξ to be stuck around unity, which was termed as "critical balance" by Goldreich and Sridhar. As far as ξ ∼ 1 we have τ dyn ∼ τ casc ∼ 1/ω, we can regard turbulence as "strong" and apply Kolmogorov phenomenology. For the cascade of the two Elsasser energies: These are two independent cascades, but in a theory with ǫ + = ǫ − this becomes standard Kolmogorov phenomenology in the k ⊥ direction.
We will return to the more general "imbalanced" case with ǫ + = ǫ − in the next section, but briefly note that such theory is non-trivial since it is impossible to maintain critical balance for resonant waves with different amplitude.
The assumption of critical balance ξ ∼ 1 allow us to estimate perturbation anisotropy directly. The "wavevector anisotropy" relates two wavevectors at which the one-dimensional spectrum along the field and perpendicular to the field have the same power. A similar relation can be obtained between parallel and perpendicular scales Λ and λ vis SFs. Using ξ = 1, and the −5/3 scaling δw ∼ λ 1/3 we obtain k ∼ k 2/3 ⊥ , which is known as GS95 anisotropy.
There is a different argument, however, that is sufficient to obtain this anisotropy. This argument is based on RMHD symetry Λ ∼ v A we discussed in Section 7.2 and dimensional grounds. Indeed, if we have Λ ∼ v A , the rest of the expression for Λ must have units of time, which is uniquely obtained from λ and ǫ as λ 2/3 ǫ −1/3 : where we introduced a dimensionless "anisotropy constant" C A . The perpendicular SF which correspond to k −5/3 ⊥ spectrum will have the scaling SF ⊥ ∼ λ 2/3 (Eq. 13), while inserting Λ ∼ λ 2/3 , we get parallel structure function as SF ∼ λ 2/3 ∼ Λ.
The parallel spectrum, which corresponds to such SF is E(k ) ∼ k −2 and from dimensional arguments we recover the prefactor as where we introduced dimensionless constant C . Equations 59 and 60 (or, alternatively 61) describe the spectrum and anisotropy of MHD turbulence, which may still be corrected for intermittency.
A modification which leads to a shallower and not steeper spectrum was proposed in [36,37], henceforth B06 suggesting that GS95 scalings are modified by a scale-dependent factor that decreases the strength of the interaction, effectively, the theory described in Sec. 4.2 with α = 1/4. Different arguments to the same effect were proposed in [38]. In this case the spectrum will be expressed as E(k) = C K2 ǫ 2/3 k −3/2 L 1/6 , see Eq. 28, Figure 14: Second order SF scaling for M1-3. Left: SF plotted vs dimensional distance r compensated by r −0.58 . Right: scaling convergence study for the r 2/3 (Kolmogorov) scaling, described in Sec. 4.3, both axes are dimensionless. We see convergence, i.e., the overall scaling is r 2/3 . the factor ξ is modified by (l/L) 1/4 , so that anisotropy follows modified critical balance with k ∼ k ⊥ 1 /2. The Kolmogorov scale of B06 model is obtained from Eq. 29: It turns out the anisotropy can be argued from the Lagrangian frequency spectrum without postulating critical balance or involving uncertainty relations. In the incompressible MHD all modes propagate with the same speed, the Elsässer components propagate either along or against the local magnetic direction, i.e., along with the magnetic field line. This propagation will be described by the functional form f (s ∓ v A t), where s is a distance along the field line. The nonlinear interaction will contribute to the slower time evolution of f and the trajectory s = ±v A t will be analogous to following hydrodynamic fluid element in the Lagrangian formulation. Let us record w + and w − along the field line in a fixed time. The positive direction s will be equivalent to following the evolution of w + backward in time and w − forward in time. In measuring the frequency spectrum, the sign of time will be unimportant. So the measurement of power spectrum along the field line will be analogous to Lagrangian frequency spectrum with frequency ω replaced by the wavenumber k = ω/v A [39]: This is the same expression as obtained in Section 7.3 from phenomenological considerations. The parallel structure function SF (l) ∼ ǫlv −1 A The dimensional argument involving Alfvén symmetry of reduced MHD arrive at the same result [40]. This symmetry allows E(k )dk to depend only on k v A , which will require that E(k ) ∼ v −1 A . The rest of the expression can be obtained from units.  Table 1 as M1-3 and M1-3H. These have a strong mean field we denote B 0 , RMS fields v rms ≈ B rms ≈ 1, perpendicular box size of 2π and parallel box size of 2πB 0 . The driving was anisotropic with anisotropy B 0 /B rms so that turbulence starts being strong from the outer scale. Technically, B 0 is arbitrary. However, the RMHD limit is only applicable to very large B 0 as we showed above. In simulations, we see a rapid decrease of parallel correlation length right after the driving scale, which indicates the efficiency of nonlinear interaction and the regime of strong turbulence. The correlation timescale for v and B was around τ ≈ 0.97, so the box contained around 6.5 parallel correlation lengths. Each simulation was started from long-evolved low-resolution simulation and was subsequently evolved for ∆t = 13.5 in code units in high resolution, and we used the last 7 dynamical times for averaging. In earlier work [17,41] it was found that averaging over ∼ 7 correlation timescales gives a reasonably good statistic on the outer scale and very good statistics on smaller scales.
Numerically, we used k max η > 1 resolution criterion, with η being classic Kolmogorov scale. Additionally, we checked the precision of the spectra by performing a resolution study on lower resolutions. In particular, we saw spectral error lower than 8 × 10 −3 , up to kη = 0.5 when increasing resolution from 576 3 to 960 3 and the spectral error lower than 3 × 10 −3 when we increased parallel resolution in a 1152 3 simulation by a factor of two. We presume this error is a mostly systematic error, associated with grid effects because the statistical error is likely to be vanishingly small, see the end of Section 4.3. We did not use any data above kη = 0.5 for fitting as the spectrum sharply declines after this point and contains negligible energy.
On Fig. 13 we plotted first order perpendicular structure functions of velocity and We see that the scaling is not obvious, with higher-resolution SF having the shallower slope of the flat part, the 4096 3 seemingly having r 0.58 scaling not expected from theory. On the right of Fig. 14 we use rigorous scaling convergence study (Sec. 4.3) to show that the overall scaling is r 2/3 (Kolmogorov). Fig. 15 presents a convergence test of the perpendicular 3D spectrum for the −5/3 model, and the convergence is reasonable, while the best convergence is reached at the -1.7 scaling. Fig. 16 shows a convergence study of the residual energy spectrum (magnetic energy minus kinetic energy). The best convergence is, again, near -1.7 slope. In all cases the convergence is consistent across two simulation groups with different dissipation prescriptions, M1-3 and M1-3H.
The flat part of the normalized spectrum can be used to obtain a Kolmogorov constant of C KA = 3.3 ± 0.1, which was first reported in [17]. The total Kolmogorov constant for both Alfvén and slow mode in the above paper was estimated as C K = 4.2 ± 0.2 for isotropically driven turbulence with zero mean field. This was obtained using an empirical energy ratio between slow and Alfvén mode, C s which is between 1 and 1.3. This larger Figure 16: Residual energy convergence. Best convergence is k −1.70 scaling for M1-3 and k −1.69 scaling for M1-3H. From [39] value of Kolmogorov constant, C K = C KA (1 + C s ) 1/3 is due to slow mode being passively advected and not contributing to nonlinearity.
We also from these simulations that the residual energy, E B − E v have the same spectral slope as the total energy, i.e., there is a constant fraction of residual energy in the inertial range. The results in Fig. 16 show that residual energy scaling is the same as for total energy so that residual energy is a constant fraction of the total energy. Our best estimate for this fraction is σ r = 0.15 ± 0.03. More commonly used in the solar wind community, Alfven ratio r A = E v /E B = (1 − σ r )/(1 + σ r ) ≈ 0.74. Residual energy and its scale-dependence has been discussed in the past and has recently been associated with the so-called called alignment measures in simulations [43] and in the solar wind measurements [44,42]. Explaining previously reported −2 scaling [45] for the residual energy is challenging from the theoretical standpoint. Assuming particular residual energy on the outer scale, and the −2 scaling, its value in the inertial range will depend on the scale separation. This would mean a nonlocal character of residual energy. Our simulations, showing that the residual energy is just a fraction of the total energy in the inertial range, resolve this conceptual difficulty and make theories suggesting different scalings for magnetic and kinetic energies obsolete.
The solar wind spectra often feature different kinetic and magnetic scalings, see Fig. 17. The amount of residual energy changes from measurement to measurement and is different for the fast and the slow solar wind [44,42]. These deviations are not observed in numerics Figure 17: Power spectra of magnetic field, velocity and residual energy measured in the solar wind. Alfven ratio was strongly fluctuating, and the average was around 0.71. From [42] and will be the subject of future study. We optimistically believe that RMHD is valid for large-scale solar wind fluctuations. However, the disagreement between simulations and the measurements could also be due to the solar wind being inhomogeneous, expanding and accelerating [46], anisotropic with respect to the sunward direction [47], and having the large number of discontinuities [48].

Numerics: parallel spectrum
We plotted the parallel spectrum E(k ) vs dimensionless wavenumber kv A τ η , compensated by k 2 ǫ −1 v A to see how the scaling is consistent with (61). This measurement is presented on Fig. 18. For the RMHD case the spectra collapsed, meaning the overall scaling of k −2 .
Reduced MHD can be performed with different the mean field strength, which in practice requires a particular choice of ǫ to generate strong turbulence from the outer scale. The Alfven symmetry of numerical RMHD formulation ensures that E(k ) scale precisely linearly ǫ. However, statistically isotropic MHD simulations with B 0 = 0 MHD1-2 do not have this symmetry, and the inertial range scaling (61) cannot be rigorously argued based on units. Our test of Eq. (61) is the test not only of the Lagrangian spectrum idea but also the Kraichnan hypothesis of dominant local v A . We substituted the RMS field instead of v A in Eq. (61). Fig. 18 demonstrates that there's convergence to ǫk −2 in this zero mean field case as well. Another spectral measurement is along the direction of the global mean field in M1-3, M1-3H. We expect these scalings to be the same as the perpendicular scalings, i.e., Kolmogorov because while Alfvén waves propagate along the local field direction which deviates by an angle of δB L /B 0 from B 0 , the angular anisotropy in this frame is δB l /B 0 , with inertial range values of δB l much smaller than the outer scale value of δB L . It follows that the anisotropy will be washed out. Fig. 19 presents a measurement of the spectrum along the global mean field direction, which is consistent with −5/3.
Worth noting that the application of the critical balance fails in the imbalanced turbulence (more on this in Section 8). A more rigorous Lagrangian argument does not have this problem. We can imagine that the energy cascade is manifested both in space and time domains, with the parallel direction is equivalent to the time domain, the anisotropy relation k ∼ k 2/3 ⊥ being the correspondence between space domain (Eulerian) and frequency domain (Lagrangian) spectra. Observational data from the solar wind points to the k −2 parallel spectrum, e.g. [50]. Numerical studies overwhelmingly support k −2 , as long as the measurements was along the local field direction see, e.g., [51,52,34,43,41], while the measurements in the global frame usually demonstrate scale-independent anisotropy, see, e.g. [53]. Figure 19: The spectra along the global mean field in M1-3, M1-3H. The M1-3H spectra have been shifted by a factor of two. The energy spectrum scales as k −5/3 , i.e., in the same way as the perpendicular scaling. From [49].

Numerics: anisotropy
We alluded above that the anisotropy should be universal in the inertial range due to the relation between Lagrangian and Eulerian spectra. We expect the relation between parallel and perpendicular scales to follow Eq. 60. Both Alfvénic and slow modes are expected to have the same anisotropy. Similar relation is expected to hold between parallel and perpendicular wavenumber: On Fig. 20 we plotted wavevector anisotropy. We determine anisotropy relation by solving the equation for k , given a range of k ⊥ : A similar procedure is done with parallel and perpendicular second order structure functions to obtain the relation between Λ and λ of Fig. 21. We also did a convergence study in the same spirit that was done for spectra in previous sections, which are on the bottom parts of each two figures above. We see that the theoretical scalings are followed fairly well. To summarize the above three sections, we showed that the spectrum and anisotropy of Alfvenic turbulence follows the two relations below: In this Section, we argued that the properties of Alfvén and slow components of incompressible MHD turbulence in the inertial range would be determined only by the Alfvén speed v A , dissipation rate ǫ and the scale of interest λ. The energy spectrum and anisotropy of Alfvén mode will be expressed as Also, we found numerically that the ratio of kinetic to magnetic energies in the inertial range is constant, ⊥ , which is expected from theory and represented by the solid line. The dashed line corresponds to the theory with 3/2 spectral scaling. Top: x-axis is a dimensional wavenumber, bottom: x-axis is a dimensionless k ⊥ η, so this plot corresponds to scaling study for anisotropy.

Dynamic alignment models
MHD has more degrees of freedom than hydro, which results, in first-order measures, in two independent dimensionless quantities (four degrees of freedom of w ± ⊥ minus rotational freedom minus normalization). One example of this is the fraction of residual energy, introduced earlier. These dimensionless quantities may, in principle, have a non-trivial scaling in the inertial range.
After spectral scaling of k −3/2 has been found in [45], a number of models have been proposed suggesting that strong turbulence phenomenology have to be modified along the lines of Sec. 4.2 to become consistent with this scaling. Among these models are [36] and [38]. A sizable confusion ensued, however to which alignment measure represent the scale-dependent weakening of interaction more accurately. The original [36] idea was analyzed in [54] and no significant alignment was found for the averaged angle between w + and w − , AA = |δw + λ × δw − λ |/|δw + λ ||δw − λ | , but when this angle was weighted with the amplitude PI = |δw + λ × δw − λ | / |δw + λ ||δw − λ | , some scale-dependent alignment was found. Later [37] proposed the alignment between v and b and subsequently [55] suggested a particular amplitude-weighted measure, DA = |δv λ × δb λ | / |δv λ ||δb λ | , that was claimed to depend of perpendicular scale as λ 1/4 . In a sense, DA is very similar to PI but uses B and v instead of w ± . The measure for local imbalance was introduced in [43] as While reevaluating the logic in [37] it becomes clear that the choice of DA as a measure exclusively responsible for the interaction weakening is arbitrary, at the same time the argument that DA scales as l 1/4 due to field line wandering, is invalid keeping in mind the Alfven symmetry of RMHD equations. The argument in [37] suggests that DA will tend to increase, but will be bounded by field wandering, i.e., the alignment on each scale will be created independently of other scales and will be proportional to the relative perturbation amplitude δB/B 0 . However, this violates Alfvén symmetry of RMHD equations (see Section 7.2), which requires that B 0 can be factored out of the dynamics and appear only in combination with k . The most apparent contradiction we find while following [37] is that a perfectly aligned state, Figure 22: Scaling study of alignment measures DA = |δv×δb| / |δvδb| and IM = |δ(w + ) 2 − δ(w − ) 2 | / δ(w + ) 2 + δ(w − ) 2 from M1-3H (top) and M1-3 (bottom). The alignment slopes converge to relatively small values, e.g., 0.06 for DA which is smaller than 1/4, predicted in alignment theories. From [39], see also [43,17,41]. e.g., with δw − = 0 is a precise solution of MHD equations and it is not destroyed by its field wandering. Empirically we know that alignment measures showed very little or no dependence on δB L /B 0 (see, e.g., [43]). Fig. 22 studies scale-dependency of DA and IM by the method of scaling study. The asymptotic scale-dependency slope for DA for our data is found around 0.06, which is way below its value of 1/4 suggested in [37]. From this figure it is evident that the scale-dependency of DA seems is tied to the outer scale, i.e., it is non-universal.

Imbalanced MHD turbulence
While hydrodynamic turbulence have only one energy cascade, the incompressible MHD turbulence has two, due to the exact conservation of the Elsässer (oppositely going wave packets') "energies". The situation of zero total cross-helicity, which we considered in previous sections has been called "balanced" turbulence as the amount of oppositely moving wavepackets balance each other, the alternative being "imbalanced" turbulence. Imbalanced turbulence, however, is very common, it is enough to have a mean magnetic field (like in the ISM) and a localized source of perturbations. Another example is the solar wind, where the dominant perturbation component propagate away from the Sun, see Fig. 23. Likewise, we expect similar phenomena in active galactic nuclei (AGN), where the jet has a mean magnetic field component and the perturbations will propagate primarily away from the central engine.

Theoretical considerations
A conceptual difficulty in the imbalanced case arises from the application of the critical balance idea. If δw + l critical balance depends on δw − l amplitude, their parallel scales and frequencies are mismatched, and the cascade cannot proceed in a normal manner. Below we describe three models that tried to deal with this issue, references [57,58,59 In short, only BL08 model is consistent with all numerical evidence, taking into account cascading rates, spectra and anisotropy. Below we shortly describe these theories.
PB09 employs dynamic alignment which depends on the scale as l 1/4 , this alignment, however, is acting differently on w + and w − so that it effectively results in the same nonlinear timescales for both components. It could be rephrased that PB09 predicts turbulent viscosity on each scale which is equal for both components. this results in an expression for the ratios of energies It is not clear how this is consistent with the limit of large imbalances, where the weak component will not be able to produce any sizable interaction.
LGS07 the authors proposed that the parallel scale for both components is determined by the shear rate of the stronger component, despite the cascading timescale is different Figure 24: Upper: a w + wavepacket, produced by cascading by w − wavepacket is aligned with respect to w − wavepacket, but misaligned with respect to the local mean field on scale λ 1 , by the angle θ. Lower: the longitudinal scale Λ of the wavepackets, as a function of their transverse scale, λ; Λ + , Λ − , λ 1 , λ 2 are the notations used in this review. From [58].
for both components. The energy cascade is still a strong cascade: This results in the prediction Both PB09 and LGS07 predict the same anisotropy for both components. In BL08 the authors proposed a new formulation of critical balance for the stronger component. This model is described in greater detail below. BL08 relaxes the assumption of local cascading for the strong component w + , while saying the w − classic critical balance and local cascading. In BL08 picture the waves have different anisotropies (see Fig. 24) and the w + wave have smaller anisotropy than w − , which is opposite to what a naive application of critical balance would predict. The anisotropies of the waves are determined by  Table 2) and A7 and A5 from [34]. The solid line is the LGS07 prediction, and the dashed line is a PB09 prediction, also a prediction for purely viscous/non-turbulent dissipation of eddies.
where λ * = √ λ 1 λ 2 , and the energy cascading is determined by weak cascading of the dominant wave and strong cascading of the subdominant wave: BL08 model, unlike LGS07, does not produce self-similar (power-law) solutions when turbulence is driven with the same anisotropy for w + and w − on the outer scale. BL08, however, claims that, on sufficiently small scales, the initial non-power-law solution will transit into asymptotic power law solution that has Λ − 0 /Λ + 0 = ǫ + /ǫ − and λ 2 /λ 1 = (ǫ + /ǫ − ) 3/2 . The larger imbalance will require larger transition to this asymptotic regime. Table 2 summarizes RMHD simulations with imbalanced driving. In these simulations, we kept the energy injection rate constant. All experiments were evolved into the stationary state. The imbalanced runs have to be were evolved for a longer time to achieve stationary state due to longer cascading timescales for the stronger component.
from [34]. I1 and I3 are simulations with the normal viscosity similar to I2 and I4. They show slightly less energy imbalances than I2 and I4. We see that most data points are above the prediction of LGS07, which is consistent with BL08. In other words, numerics strongly suggest that Although there is a tentative correspondence between LGS07 and the data for small degrees of imbalance, the deviations for large imbalances are significant. As to PB09 prediction, it is inconsistent with data for all degrees of imbalance including those with small imbalances. The approximate equality in Eq. 74 for very small imbalances, however, is an excellent test of the expression in Eq. 58 that we assumed in the balanced case. So in some sense, empirical study of the imbalanced case validated the theory of the balanced case as well. Fig. 26 shows spectra from low-imbalance simulation I2, compensated by the predictions of PB09 and LGS07. We see that the collapse of two curves for w + and w − is much better for the LGS07 model. Fig. 27 shows spectra from all I1-6 simulations, compensated by the prediction of LGS07. For lower imbalances, the collapse is reasonably good and become progressively worse for larger imbalances. This deviation, however, does not entirely follow the prediction of the asymptotic power-law solutions from BL08, which will predict that the solid curve will go above C KA and the dashed curve -below it. This is possibly explained by the fact that asymptotic power-law solutions were not reached in these limited resolution experiments, this is also observed for anisotropies.
We measured parallel and perpendicular structure functions in simulations I1-I6, in order to quantify the anisotropies of eddies. Fig. 28 shows anisotropies for I1-6 simulations. All simulations were driven by the same anisotropies on the outer scale, which is unfavorable for obtaining the asymptotic power law solutions of BL08. It is, however, favorable to the LGS07 model, which predicts the same w + and w − anisotropies for all scales. Therefore, these simulations are a sensitive test between LGS07 and BL08 models, both of which are roughly consistent in terms of energy ratios and spectra for small imbalances. If the LGS07 model is correct, we would observe the same anisotropy on all scales, but this is not what is observed in Fig. 28, where anisotropies start to diverge on smaller scales. The ratio of anisotropies is roughly consistent with the BL08 prediction of ǫ + /ǫ − for small imbalances and somewhat smaller for larger imbalances.

MHD dynamo
One of the main questions of MHD dynamics is how conductive fluid generates its magnetic field, a process known broadly as "dynamo". Turbulent dynamo is known as "largescale/mean-field dynamo" and "small-scale/fluctuation dynamo", in the first case magnetic fields are amplified on scales larger than the outer scale of turbulence in the seconds on smaller scales. An example of the flow generating no dynamo is an axisymmetric situation [60], the natural turbulence, however, possesses no exact symmetries and is expected to amplify magnetic field by stretching, due to the particle separation in a turbulent flow. For the large-scale dynamo, a "twist-stretch-fold" mechanism was introduced in [61].
If the turbulent flow possess statistical isotropy, it can not generate a large-scale field, i.e., the field with scales larger than the outer scale of turbulence. To generate the observed large-scale fields, such as fields in the disk galaxies, large-scale asymmetries of the system must break the statistical symmetry of turbulence. Further analysis of the induction equation shows that the rotation (described by the pseudovector of angular velocity) is insufficient to provide large-scale dynamo, stratification or shear should provide a real vector, so that the pseudoscalar alpha-effect result in a large-scale turbulent EMF [62,63].
One approach to large-scale dynamo is mean field theory, see, e.g. [64], where the magnetic and velocity fields are decomposed into the mean and fluctuating part. The equations for the mean field are closed using statistical or volume averaging over the fluctuating part. The traditional theories of mean field dynamo, however, often fail due to issues related to magnetic helicity [62,65].
The studied of the large-scale dynamos is big and complex science due to the variety of ways, large-scale symmetries are broken in different astrophysical objects. In this review, we propose the reader follow other reviews focused on large-scale dynamos, e.g., [66,67] and instead concentrate on the small-scale or fluctuating dynamo as more universal, generic and crucial to understand the overall level of magnetization in astrophysical environments. We also emphasize that proper understanding of turbulent magnetic fields is crucial as it is subsequently slowly ordered and made large-scale by the large-scale dynamo process.

Kinematic dynamo
If the magnetic energy is less than the kinetic energy of turbulent motions, the turbulence may generate the magnetic field, which is referred to as "turbulent dynamo". There are two distinct regimes: 1) the magnetic energy is much less than the kinetic energy of driving eddies at all scales down to the dissipation scale and 2) the kinetic and magnetic energies come to the equipartition at some scale. The first regime is called "kinematic" or "linear" dynamo, referring to induction equation (35) being linear with respect to the magnetic field. If the magnetic field is so weak that it provides no back-reaction to velocity, the problem reduces to studying the solutions of the induction equation with a prescribed velocity field. This kinematic regime was, historically, the most studied (see Kazantsev model outlined in [68,69,70,71]). For the Kolmogorov-type turbulence, the fastest magnetic field amplification comes from the fastest turbulent eddies, i.e., the eddies at the dissipation scales. It would follow that the growth rate γ = 1/τ η (see Eq. 26). The spectrum of the magnetic field will be E(k) ∼ k 3/2 rising sharply to the magnetic dissipation scale. Kazantzev picture was not without drawbacks since it had relied on an artificial delta-correlated velocity field instead of a realistic turbulent velocity field. This resulted in an overestimated γ compared to reality. Precise numerical experiments with Pr m = Re m /Re = 1 have found a prefactor γτ η = 0.0326 [72,73,40] which is much smaller than unity. This small number, as had been suggested in [40] is due to turbulence being time-asymmetric.
From kinematic models, it is not clear whether magnetic energy will continue to grow after the end of the kinematic regime, however, keeping in mind γ ∼ Re 1/2 and very large astrophysical Re, the kinematic growth is incredibly fast and occupies a tiny fraction of the dynamical time of the system. For example, while the galaxy clusters form on timescales of 15 billion years, the characteristic growth time of kinematic dynamo is less than a million years [74], so we do not expect kinematic dynamo to operate in present time. The magnetic spectrum of the kinematic dynamo, with its positive spectral index, is incompatible with observations in galaxy clusters [75], see Fig. 1. These observations indicate steep spectrum with negative power index at small scales.
Summarizing, the kinematic dynamo is inapplicable in most astrophysical environments, since the observed magnetization corresponds to Alfvén speed which is many orders of magnitude higher than the Kolmogorov velocity ǫ 1/4 ν 1/4 .

Nonlinear dynamo
Before numerical simulations were commonplace, it was proposed that after saturation of kinematic dynamo the magnetic energy will stop growing. If we agree to this proposition and assume that the magnetic energy indeed saturates as soon as the dynamo becomes nonlinear, then the saturation level, in this case, will be of order ρv 2 η /2, where v η is a Kolmogorov velocity scale. This is a factor of Re −1/2 smaller than the kinetic energy density Figure 29: A cartoon of kinetic and magnetic spectra in small-scale dynamo, at a particular moment of time when equipartition wavenumber is k * . and will be completely dynamically unimportant in high-Re astrophysical environments. In fact, observations indicate the opposite -a sizable energy density of the magnetic field in large-scale systems, see Fig. 1 An alternative had been proposed in the early work by Schlüter and Bierman [76], who suggested that dynamo will not stop and will continue to grow, saturating on each subsequent scale after a dynamical time. Recently small-scale dynamo underwent revival due to the availability of direct numerical simulations. Simulations of the dynamo saturated state produced steep spectra and significant outer-scale fields. The saturated state was only weakly dependent on Re and Pr m as long as Re was large, see, e.g., [72].
Apart from the saturated state, the growth stage was suggested to have growth of magnetic energy which is linear in time [77,78,79,80,40]. In [40] the locality of the dynamo, which is necessary for the linear growth picture was argued analytically. Let us imagine that the magnetic and kinetic spectra at a particular moment of time are similar to what is presented in Fig. 29. Magnetic and kinetic spectra cross at some "equipartition" scale 1/k * , below which both spectra are steep, typically k −5/3 due to the MHD cascade (Sec. 7). A number of arguments suggest this assumption. Firstly, we expect kinematic dynamo [70] to proceed until the moment when magnetic spectrum intersect the kinetic spectrum at the viscous scales (assuming Pr m = 1), which will correspond to the beginning of the nonlinear regime. Secondly, this is supported by the wealth of numerics, e.g. [43,78]. Thirdly, this is also somewhat supported by observations of magnetic fields in clusters [75].
At large scales magnetic spectrum is shallow, k α , α > 0, while kinetic spectrum is steep due to the hydrodynamic cascade. Most of the magnetic energy is contained at the scale of 1/k * . We designate C K and C M as Kolmogorov constants of hydro and MHD respectively. The hydrodynamic cascade rate is ǫ and the MHD cascade rate as ǫ 2 . Due to the conservation of energy in the inertial range, magnetic energy will grow at a rate of ǫ − ǫ 2 . We will designate C E = (ǫ − ǫ 2 )/ǫ as an "efficiency of the small-scale dynamo" and will argue that this is a true constant, since a) turbulent dynamics is local in scale in the inertial range; b) ideal MHD or Euler equations do not contain any scale explicitly. Magnetic energy will grow linearly with time if ǫ = const: 1 8π The equipartition scale L B = 1/k * will grow with time as t 3/2 [80]: here we introduced dimensionless constant c l . Alternatively one can say that small-scale dynamo saturates at several dynamical times at scale 1/k * and proceeds to a twice larger scale [76,77]. If magnetic energy grows approximately till equipartition [72,78], the whole process will take around several outer timescales of the system, or more quantitatively, (C It was demonstrated analytically in [40] that as long as the kinetic spectrum is steep (spectral slope between -1 and -3), the magnetic spectrum is steep below the equipartition scale and magnetic spectrum is shallow (slope higher than 0) above the equipartition scale, the dynamo is indeed local and the picture described above can indeed be rigorously argued. So, besides the fact that local interactions dominate the kinetic cascade at large scales and the MHD cascade at small scales, we also know that dynamo is governed by local interactions. Assuming Kolmogorov phenomenology, it is also possible to estimate the upper limit on the wavevector interval in which nonlinear dynamo operates, namely the interval k * [C , which not very restrictive as a practical upper limit for numerical simulations, but undoubtedly essential for astrophysical situations where the inertial range is many orders of magnitude in scale. Given the sufficiently large inertial range, C E a universal dimensionless constant of MHD dynamics, much like Kolmogorov constant and since it relates energy fluxes, not energies, it is also unaffected by intermittency. Interestingly, magnetic energy dynamics at k ≪ k * is likely dominated by nonlocal interactions with k * but this part of the spectrum contains negligible magnetic energy and the universality claim is unaffected by this nonlocality.

Numerical results
Numerical simulations of statistically homogeneous isotropic small-scale dynamo in [40] were performed by solving MHD equations with stochastic non-helical driving and explicit dissipation with P r m = 1. The results in Fig. 30 is the statistical average over three different simulations. Growth is initially exponential and smoothly transition into the linear stage. Note that scatter is initially small, but grows with time, which is consistent with the picture of the magnetic field that grows at the progressively larger scales, and has progressively less independent realizations in a single datacube. The value of the dynamo efficiency that we measure C E = 0.05 is much smaller than unity. One would expect this quantity of order unity because this is a universal number, determined only by strong interaction on the equipartition scale. If we refer to the ideal incompressible MHD equations, written in terms of Elsässer variables, ∂ t w ± +Ŝ(w ∓ · ∇)w ± = 0, dynamo could be understood as decorrelation of w ± which are originally equal to each other in the hydrodynamic cascade. In our case, this decorrelation is happening at the equipartition scale 1/k * . Being time-dependent, it propagates upscale, while, ordinarily, energy cascade goes downscale. The small value of C E might be due to this. As opposed to picture with multiple reversals and dissipation due to microscopic diffusivity, typical for the kinematic case, in our picture we appeal to turbulent diffusion which helps to create the large-scale field. Both stretching and diffusion depend on turbulence at the same designated scale 1/k * , so, in the asymptotic regime of large Re, one of these processes must dominate. As C E is small, we conclude that stretching and diffusion are close to canceling each other. In [40], the interplay of stretching and mixing was studied by simulations of kinematic dynamo forward and backward in time. Basically, forward in time stretching is less efficient, while mixing is more efficient. This also tells a cautionary tale that using artificial statistics of velocity, such as delta-correlated statistics, may be grossly misleading when dealing with real turbulence.
The arguments of the previous section can be applied even if the energy injection rate is not stationary. On Fig. 31 we presented simulations with intermittent driving checking relations (75,76).  [40]), while the lower plot correspond to intermittent driving with the period of 8τ c . In this relation we used the dissipation rate ǫ averaged over 2τ c . Outer scale was determined by the peak wavenumber of the spectrum. Best fit c l ≈ 0.18 in Eq. (76). From [74] 9.4 Caveats in simulating astrophysical dynamo As we alluded above in Section 9.1 the timescales of initial growth of the magnetic field are extremely short compared to the dynamical timescale of the system -these two are separated by the factor of Re 1/2 , which is very large in astrophysical environments. Not so in numerical simulations, which are limited in terms of Re. This situation is further exacerbated by the small prefactor 0.0326 in the growth rate of γ = 0.0326/τ η . For example, simulation with Re = 1600 will result in characteristic growth time of 1/(0.0326Re 1/2 ) ≈ 0.77 of the outer timescale of the system, while in astrophysical reality this is a negligible fraction of the outer timescale 1 . The lesson from this is that if we take a very small initial magnetic field while simulating rather dynamically young object, e.g., collapsing cloud or a forming galaxy cluster, this may artificially delay the onset of the nonlinear dynamo and grossly underestimate the magnetic field at the end of the evolution [74]. This delay effect can also lead to the artificial dependence on the initial field value or direction, which should not normally appear in nonlinear dynamo, which, as any turbulence, erases any traces of the initial condition. The solution to this issue could be injecting initial field based on the amount of energy in the cascade, along the lines of Section 9.2.
Another issue is that popular ILES codes lack any knowledge or prescription of the microscales. While very often this is not an issue, because ILES code would simply absorb cascade energy into the thermal energy on the grid scales, but there is a qualitative difference between nature and an ILES simulation with zero initial field, which would produce zero magnetic field for all subsequent times. In nature, the small-scale nonideal contributions to the induction equation, such as Biermann battery term, will always jump-start the dynamo and result in non-zero fields, also large Re in nature will ensure that average magnetization will be close to the estimate given by nonlinear dynamo in Section 9.2, i.e. the magnetic energy will be a certain fraction of the dissipated cascade energy. From MHD point of view, taking diffusivities to zero allows us to have two types of solutions -completely unmagnetized and fairly strongly magnetized, however in nature only magnetized solution will be realized 2 . We conclude that in dynamo situations ILES code should sometimes use subgrid dynamo prescription.

Application to galaxy clusters
Galaxy clusters constitute an interesting case of the small-scale dynamo. All properties of the cluster continue to evolve during the cosmic time, its mass also determining virial velocity which is always around the inflow velocity. Clusters are heated by the major mergers, which is also the primary source of energy for the intracluster turbulence. On average, around 0.4 of heating comes from dissipation of solenoidal turbulence [81]. At the same time, as we saw in the previous section dynamo converts a fraction of 0.05 of this energy into magnetic energy. Thus we may conclude that magnetization β, the ratio of magnetic to thermal energy, should stay approximately constant through cosmic time, ∼ 40 for the past 10 Gyr. Similarly, the outer magnetic scale will stay a constant fraction of the cluster virial radius, ∼ 1/200.

Compressible and Supersonic Turbulence
The study of supersonic ISM turbulence is vital for understanding the structure of molecular clouds and subsequent star formation. In this respect, the studies of density scalings and thermal instability in DNS have become commonplace. One way to look at it is assuming that in addition to the Alfven mode, additional compressible modes are excited as well. How this approach is valid in the regime where Alfven waves couple strongly among themselves is still a matter of debate. Below we will show a mode decomposition techniques, as well as study density perturbations which are perturbed nonlinearly, by orders of magnitude in the supersonic regime.
Compressible turbulence is characterized by two dimensionless numbers, sonic Mach number M s = δv/c s and Alfvenic Mach number M A = δv/v A . Plasma beta, the ratio of gas pressure to magnetic pressure can be expressed as β = 8πP gas /B 2 = 2γM 2 A /M 2 s . In a turbulent environment, these numbers tend to evolve towards "preferred" values, depending on the physics. For example in the absence of cooling, supersonic turbulence will heat itself to trans-sonic state M s ∼ 1 in one dynamical time. The persistent supersonic regime, therefore, is limited only to regions which have very efficient atomic or molecular cooling, such as molecular clouds, where M s 10.
Likewise, the Alfvenic Mach number is limited by the dynamo process that we described in Section 9. Turbulence with a very weak initial field tends to generate equipartitionstrength fields, with M A ∼ 1, in around 20 dynamical times. We expect large (but not very large) M A only in dynamically young systems, such as in collapsing molecular clouds.
In galaxy clusters, which are continuously heated up by accretion and turbulence, its mass, and, therefore, virial velocity, also continue to increase. This leads to the situation when both δv and c s are determined by virial velocity and M s ∼ 1. Due to the continuous growth of the cluster mass, it never approaches equipartition, so M A ∼ 7.

Decomposition into modes
Here we will use expressions for the phase wave speeds of the Alfven and slow/fast modes (37,38) from Section 5 and the same notation for cos θ = (k ·B). Our coordinate system will be defined by the unit vector of the magnetic fieldB, the unit vector of the Alfven perturbationê A = k × B 0 /|k × B 0 |, and the third unit vector perpendicular to the three, in the (k, B) plane:ê ⊥ =ê A ×B 0 , see Fig. 32.
The perturbations of the velocity for the Alfven mode are alongê A , also from induction equation −ωδB = k × (δv × B 0 ), however the sign of δB depend on the sign of ω, i.e. whether the wave is propagating along the field direction, or in the opposite direction. The waves where δB ∼ −δv the waves propagate along the field. This also corresponds to δw + propagating opposite to the field and δw − along the field.
This two vectors andê A form an orthogonal coordinate system in wavenumber space. The perturbation δB for both modes is in the direction, however again the sign of the contribution for each mode will depend on the direction of the wave propagation. Introducing ∆v f,s as perturbations where oppositely propagating wave contribute with different signs, the contribution to δB for each mode will be proportional to |∆v f,s /u f,s | (see, e.g., [82]). Alternatively, one can decompose in the coordinate system rotated aroundê A by θ, which is made ofk,ê A andê ⊥2 =ê A ×k, in which case the velocity and magnetic field vectors along fast and slow modes can be expressed as It was suggested in [83,82] that the Alfven mode have an independent cascade (Sec. 7), slow mode is passive to Alfven cascade and has the same spectra and anisotropy, and fast mode has an independent isotropic acoustic/wave turbulence cascade (Sec. 4.2): Fast: This was broadly consistent with simulations in [83,82], however, steeper spectrum ∼ k −2 was reported in [84] for fast modes. Regarding the amplitude of each mode in realistic turbulence, it stronly depend on the way the turbulence is driven. For a particular incompressible driving of [83,82] they suggested a scaling relation Figure 33: Anisotropy of the Alfvén, slow and fast modes as evidenced by the contours of the second order structure function. Here we used the SF decomposition method. The Alfvén and slow mode exhibit scale-dependent anisotropy, while the fast mode is almost isotropic.
where δE fast and δE Alf are energy of fast and Alfvén modes, respectively.

Decomposition in real space
Another way to decompose into modes is by using structure functions. In this method, the separation vector l of the structure function plays the role of the wavenumber, because of the correspondence between one-dimensional structure function along the particular line and the power spectrum along the same line (Sec. 3). Fig. 33 shows the contours of the structure function corresponding to each mode obtained in datacubes from M s = 10 supersonic simulations used earlier in [85]. The anisotropies of each mode show the same behavior as in the global decomposition method discussed above. The advantages in using the SF decomposition method is that it is a local measurement, so we can apply it even when the global average magnetic field is zero, e.g., it has been applied to the decomposition of MHD turbulence obtained in the high-resolution cosmological simulation of a galaxy cluster [86].

Density scalings
The properties of density in supersonic turbulence is interesting to astronomers due to its relation to star formation. The density is primarily perturbed by the slow mode; however, in supersonic case, these perturbations are not small. Instead, density varies by several orders of magnitude within the box (Fig. 34). At the same time, the statistics of density is very different from the statistics of the slow mode velocity. It is worth noting that supersonic isothermal hydrodynamics predicts the log-normal distribution of density [87] due to the gauge symmetry of log-density in this case. This symmetry is broken in MHD. However the PDF still resemble log-normal law (Fig. 35). The paper [85] pointed out that in the low beta supersonic regime the perturbation of the slow mode velocity is almost along B 0 , so the dynamics of the slow mode is quasi-one-dimensional. This results in the generation of many slow shocks, which is indeed observed in simulations.  From [85] The perturbations of density, e.g., its log-normal PDF, are created by these random slow shocks, similar to hydrodynamics. The Alfén mode only mixes these perturbations by shearing motions, without affecting PDF. On the other hand, the structure of density (SFs) is almost entirely determined by this shearing and is expected to have an anisotropic structure like velocity and magnetic field. In other words, two distinct physical processes act simultaneously and affect different statistical measures of the turbulent density field. The random multiplication of density induced by shocks, affect the PDFs, while the other, Alfvénic shearing, affects anisotropy and scaling of the structure function of the density. Revealing the structure created by shearing requires overcoming the effect of high-density clumps which will dominate SF 2 (ρ). So, instead, we should use statistics of log(ρ), which have approximate Gaussian PDFs. This exercise is shown in Fig. 35. Indeed, the statistics of log(ρ) shows anisotropy characteristic of Alfvén shearing.

Turbulence driven by magnetic field
Turbulence in reconnection can appear as a result of instabilities, for example, resistive tearing [88]. The paper [89] demonstrated that the instability becomes faster and not slower with decreasing resistivity above a critical Lundquist number around 10 4 . Observing effects of the feedback of the release of magnetic energy in numerics is challenging because currently available 3D MHD numerics are limited by the Lundquist numbers of several of 10 4 . Physically, periodic box simulations, like numerics in [90] correspond to early times in the current sheet disruption when the outflow did not develop. Importantly enough, it did demonstrate fast (resistively-independent) reconnection rates, defined as mixing rates of the fluid. Simulations with open boundaries in [91] have been performed for a sufficiently long time to allow for the establishment of the stationary state. They correspond to later times when the stationary inflow/outflow appears.
One of the simplest setups to study the development of turbulence in the thin current layer is a periodic setup with the mean field B z0 threading the box, reconnecting field ±B y0 changing the sign in the x direction, see Fig. 36. We consider the incompressible case, in which case the only dimensionless parameters of the problem are the Lundquist number S and the ratio B y0 /B z0 . We use the planar sheet in an attempt to simulate a zoomed-in portion of a very large and unstable Sweet-Parker current layer. Lundquist number is defined with the box size, as S = v Ay L/η. We imagine that this box is a part of a bigger system with larger system size L global . We aim to simulate early times, Figure 37: Left: The y-z power spectra of velocity and magnetic perturbations of turbulence in the current layer. Right: Anisotropy from the ratio of parallel to perpendicular scales obtained from equating 2-order SFs (see Section 3.12). We used simulations with 2nd as well as 4th order diffusivities (hyper diffusivities) to evaluate the effect of the dissipation on the statistics of turbulence in the layer. From [90].
t < L global /V A , when the global outflow did yet develop. We also assume that the global Lundquist number, determined by the larger system, is asymptotically large so that we can ignore large-scale gradients. The simulation end is determined by the development of structures with the size comparable to the box size, at which point our artificial periodic boundary starts influencing the result.
The free energy in the system is the energy density of the opposing fields B 2 y0 /8π, which is declining in the turbulent current layer due to dissipation. After t ≈ 0.3L/v Az the fraction of the dissipated energy w d becomes approximately constant, around w d ≈ 0.4. We calculate the reconnection rate as the speed of growth of the turbulent current layer width ∆, i.e. we define V r = d∆/dt. The evolution of d and the inferred reconnection rate are shown in Fig. 36. V r was around 0.015v Ay for high Lundquist numbers and is rather insensitive to the imposed mean field B z0 (Fig. 36). The dissipation rate per unit area of the current sheet can be calculated from w d and v r as Note that we arrived at the expression not only for "fast reconnection" (independent on resistivity and viscosity) but also for "fast dissipation". This expression, modulo numerical coefficient, can be obtained by dimensional analysis using only ρ and v Ay .
The field in the current layer can be analyzed statistically. We show the spectrum for one time slice on Fig. 37. The peak of the spectrum moves towards smaller wavenumbers, i.e., the outer scale of this turbulence is growing in time. This is unlike driven turbulence (Sec 7) where this scale was determined by forcing and fixed. Another difference with driven turbulence is that magnetic spectrum is above kinetic on all scales, but closer to equipartition on smaller scales. This is similar to decaying MHD turbulence described, e.g., in [12]. Qualitatively reconnection turbulence is very similar to decaying turbulence created by the initial random magnetic field.
Scale-locality is an important component of turbulent reconnection. Our spontaneous reconnection numerics corroborate scale-locality, because the spectral slope of perturbations is between -1 and -3. In the real world, we expect the reconnection rate to be independent of system size as long as ion Larmor radius r L and ion skip depth d i are both much smaller than the layer width ∆. On the right-hand side, Fig. 37 shows anisotropy expressed as a ratio of parallel to perpendicular scale λ /λ ⊥ , obtained by a method we explain in Sec 7. We can also estimate the interaction strength parameter ξ = δvλ /v A λ ⊥ and see that for this case it is around unity, i.e., we are dealing with critically balanced strong turbulence. Note that the anisotropy of our turbulence, being around k /k ⊥ ∼ 1/20 is very different from the tangent of the fastest growing oblique tearing mode, k /k ⊥ = B z /B y = 1. So turbulence forgets the properties of the oblique tearing that started it. From simulations with higher B z one also confirm Alfvén symmetry: increasing B z only increases parallel lengthscale, while keeping dynamics essentially unchanged (see [90]).

Conclusion
This review covers a set of topics that were chosen because they are either: a) basic topics that are essential for the understanding of subsequent material or b) have seen rapid progress recently, which is otherwise not covered in books or reviews. Due to this choice, many things, especially astrophysical and space applications of MHD turbulence, has been omitted or mentioned only in passing. This document, however, is a living review and will be evolving, below we overview several topics that we expect to add or expand. For the impatient, we mention older books and reviews that can be used to expand and deepen the reader's knowledge of the topic. Most of these have been already mentioned in the course of the review.
Mathematical tools to work with the statistical ensemble of turbulent realization can be found in the monograph by Monin and Yaglom [11]. Overview of turbulence as a nonlinear dynamical process can be found in the book of Frisch [18], Falkovich (2011) [92], Davidson (2015) [93]. Comprehensive, although older, book dedicated to MHD turbulence is by Biskamp (2003) [12], a few topics in MHD turbulence are also covered in Davidson (2013) [94]. Older book on mean-field dynamo is Krause and Raedler (1980) [64], a more modern approach to the same topic, primarily for solar dynamo applications is a living review by Charbonneau (2010) [67]. A more broad review on dynamo is by Brandenburg and Subramanian (2005) [66]. In future editions, we plan to cover large-scale dynamo as well. For an in-depth review of solar wind turbulence, see living review by Bruno and Carbone [10]. We plan to expand the section related to the solar wind and cover energy flux [95] as well as magnetic helicity measurements [96]. An interesting case of energy cascade with applications to cosmological-scale magnetic fields and its dynamical evolution is a freely decaying homogeneous MHD turbulence, see, e.g. [97], which we also plan to cover in the future.
Several topics connecting astrophysical plasmas and MHD turbulence can be found in a book by Kulsrud (2005) [98]. In the ISM, MHD turbulence coexists with cosmic rays.
Cosmic ray interaction with MHD turbulence is a fairly large topic, for an introduction to cosmic rays as well as quasilinear scattering theory one can start with Schlickeiser (2002) [99]. One particularly important application of mutual interaction between cosmic rays and MHD turbulence which will be covered in future editions of this review is the acceleration of cosmic rays in supernova remnants. In front of strong shocks, MHD turbulence is self-generated by fast particles. This is supported by estimates of diffusion coefficient D of cosmic rays in supernova remnants. D in front of the shock is estimated to be many orders of magnitude smaller than D in the ambient ISM, i.e. cosmic rays create their own MHD turbulence and dynamo and scatter themselves.
Supersonic turbulence with applications to ISM and star formations is covered in Mac Low and Klessen 2004 [8] and McKee and Ostriker (2007) [9]. The physics of turbulent energy cascade in the supersonic case has been an open question for some time, but recently we saw progress in deriving exact analytic relations in supersonic case [100,101,102], as well as empirical findings and numerical verification [103,104]. The earlier phenomenological approach of replacing statistics of velocity u with the statistics of ρ 1/3 u in the compressible case [105,106] have found a firmer foundation [107]. This has also been used to explain observed statistical correlations, such as Larson's laws [108] in starforming clouds [104].
Intermittency is the deviation from self-similarity of turbulence and is an important property that reminds us of the richness of the field of nonlinear fluid dynamics. While intermittency in hydrodynamics has been long studied as-is, the intermittency in different variables in MHD turbulence may give us something to think about its dynamical origin. On the other hand, extreme intermittency had been suggested to explain heating and molecular synthesis in the ISM. This will be covered in more detail in the future editions of this review. The numerical section will be expanded with mention of Lagrangian-Eulerian (moving mesh) codes and recent progress in this area. In future editions, we also will pay more attention to the connection between theory and observations. Big progress has been achieved in the area of cosmological structure formation by massive ab-initio simulations including ΛCDM initial conditions with Λ, dark matter and ordinary matter using grid refinement down to the scales of galaxies (e.g. Illustris project). Some of these results are relevant to understand magnetization in filaments and, possibly, clusters and will be added later.