Systematics of the magnetic-Prandtl-number dependence of homogeneous, isotropic magnetohydrodynamic turbulence

We present the results of our detailed pseudospectral direct numerical simulation (DNS) studies, with up to 10243 collocation points, of incompressible, magnetohydrodynamic (MHD) turbulence in three dimensions, without a mean magnetic field. Our study concentrates on the dependence of various statistical properties of both decaying and statistically steady MHD turbulence on the magnetic Prandtl number PrM over a large range, namely 0.01⩽PrM⩽10. We obtain data for a wide variety of statistical measures, such as probability distribution functions (PDFs) of the moduli of the vorticity and current density, the energy dissipation rates, and velocity and magnetic-field increments, energy and other spectra, velocity and magnetic-field structure functions, which we use to characterize intermittency, isosurfaces of quantities, such as the moduli of the vorticity and current density, and joint PDFs, such as those of fluid and magnetic dissipation rates. Our systematic study uncovers interesting results that have not been noted hitherto. In particular, we find a crossover from a larger intermittency in the magnetic field than in the velocity field, at large PrM, to a smaller intermittency in the magnetic field than in the velocity field, at low PrM. Furthermore, a comparison of our results for decaying MHD turbulence and its forced, statistically steady analogue suggests that we have strong universality in the sense that, for a fixed value of PrM, multiscaling exponent ratios agree, at least within our error bars, for both decaying and statistically steady homogeneous, isotropic MHD turbulence.


Introduction
The hydrodynamics of conducting fluids is of great importance in many terrestrial and astrophysical phenomena. Examples include the generation of magnetic fields via dynamo action in the interstellar medium, stars and planets [1]- [11] and in liquid-metal systems [12]- [18] that are studied in laboratories. The flows in such settings, which can be described at the simplest level by the equations of magnetohydrodynamics (MHD), are often turbulent [5]. The larger the kinetic and magnetic Reynolds numbers, Re = U L/ν and Re M = U L/η, respectively, the more turbulent is the motion of the conducting fluid; here L and U are typical length and velocity scales in the flow, ν is the kinematic viscosity and η is the magnetic diffusivity. The statistical characterization of turbulent MHD flows, which continues to pose challenges for experiments [19], direct numerical simulations (DNS) [20] and theory [21], is even harder than its analogue in fluid turbulence, because (i) we must control both Re and Re M , and (ii) we must obtain the statistical properties of both the velocity and the magnetic fields.
The kinematic viscosity ν and the magnetic diffusivity η can differ by several orders of magnitude, so the magnetic Prandtl number Pr M ≡ Re M /Re = ν/η can vary over a large range. For example, Pr M 10 −5 in the liquid-sodium system [15,16], Pr M 10 −2 at the base of the Sun's convection zone [22] and Pr M 10 14 in the interstellar medium [8,20]. Furthermore, two dissipative scales play an important role in MHD; they are the Kolmogorov scale d (∼ν 3/4 at the level of Kolmogorov 1941 (K41) phenomenology [23,24]) and the magnetic-resistive scale M d (∼η 3/4 in K41). A thorough study of the statistical properties of MHD turbulence must resolve both of these dissipative scales. Given current computational resources, this is a daunting task at large Re, especially if Pr M is significantly different from unity. Thus, most 4

Magnetohydrodynamic (MHD) equations
The hydrodynamics of a conducting fluid is governed by the MHD equations [1]- [5], [7], in which the Navier-Stokes equation for a fluid is coupled to the induction equation for the magnetic field, ∂u ∂t Here, u, b, ω = ∇ × u and j = ∇ × b are, respectively, the velocity field, the magnetic field, the vorticity and the current density at the point x and time t; ν and η are the kinematic viscosity and the magnetic diffusivity, respectively, and the effective pressure isp = p + (b 2 /8π), where p is the pressure; f u and f b are the external forces; while studying decaying MHD turbulence, we set f u = f b = 0. The MHD equations can also be written in terms of the Elsässer variables z ± = u ± b [7,25]. We restrict ourselves to low-Mach-number flows, so we use the incompressibility condition ∇ · u(x, t) = 0; and we must, of course, impose ∇ · b(x, t) = 0. By using the incompressibility condition, we can eliminate the effective pressure and obtain the velocity and magnetic fields via a pseudospectral method that we describe in section 2.1. The effective pressure then follows from the solution of the Poisson equation,

Direct numerical simulation
Our goal is to study the statistical properties of homogeneous and isotropic MHD turbulence, so we use periodic boundary conditions and a standard pseudospectral method [44] with N 3 collocation points in a cubical simulation domain with sides of length L = 2π ; thus, we evaluate spatial derivatives in Fourier space and local products of fields in real space. We use the 2/3 dealiasing method [44] to remove aliasing errors; after this dealiasing, k max is the magnitude of the largest-magnitude wave vectors resolved in our DNS studies. We have carried out extensive simulations with N = 512 and N = 1024; the parameters that we use for different runs are given in table 1 for decaying and statistically steady turbulence.
We use a second-order, slaved, Adams-Bashforth scheme, with a time step δt, for the time evolution of the velocity and magnetic fields; this time step is chosen such that the Courant-Friedrichs-Lewy (CFL) condition is satisfied [45].
In our decaying-MHD-turbulence studies, we have taken the initial (superscript 0) energy spectra E 0 u (k) and E 0 b (k), for velocity and magnetic fields, respectively, to be the same; specifically, we have chosen where E 0 , the initial amplitude, is chosen in such a way that we resolve both fluid and magnetic dissipation scales η u d and η b d , respectively: in all, except for a few, of our runs, k max η u d 1 and k max η b d 1. The initial phases of the Fourier components of the velocity and magnetic fields are taken to be different and chosen such that they are distributed randomly and uniformly between 0 and 2π. In such studies, it is convenient to pick a reference time at which various statistical Table 1. List of parameters for our 16 DNS runs R1-R5, R3B-R5B, R1C-R4C and R1D-R4D: N 3 is the number of collocation points in our simulation, ν is the kinematic viscosity, Pr M is the magnetic Prandtl number, δt is the time step; and u rms , I , λ and Re λ are the root-mean-sqare velocity, the integral scale, the Taylor microscale and the Taylor-microscale Reynolds number, respectively. These are obtained at t c for our decaying-MHD turbulence runs R1-R5, R3B-R5B and R1C-R4C; and for statistically steady MHD turbulence (runs R1D-R4D), these are averaged over the statistically steady state; here, t c (iteration steps multiplied by δt) is the time at which the cascades for both the fluid and the magnetic fields are completed (see text); η u d and η b d are, respectively, the Kolmogorov dissipation length scales for the fluid and magnetic fields. k max is the magnitude of the largest-magnitude wave vectors resolved in our DNS studies which use the 2/3 dealiasing rule; k max 170. 67  properties can be compared. One such reference time is the peak that occurs in a plot of the energy dissipation versus time; this reference time has been used in studies of decaying fluid turbulence [46,47], decaying fluid turbulence with polymer additives [48,49] and decaying MHD turbulence [25,26,50]. Such peaks are associated with the completion of the energy cascade from large length scales, at which energy is injected into the system, to small length scales, at which viscous losses are significant. In the MHD case, these peaks occur at slightly different times, t u and t b , respectively, in plots of the kinetic ( u ) and magnetic ( b ) energydissipation rates. In our decaying-MHD-turbulence studies, we store velocity and magnetic fields at time t c ; if t u > t b , t c = t u ; and t c = t b otherwise; from these fields we calculate the statistical properties that we present in the next section.
In the simulations in which we force the MHD equations to obtain a nonequilibrium statistically steady state (NESS), we use a generalization of the constant-energy-injection method described in [51]. We do not force the magnetic field directly, so we choose f b = 0.
6 The force f u (x, t) is specified most simply in terms off u (k, t), its spatial Fourier components, as follows,f where (k f − k) is 1 if k k f and 0 otherwise, P is the power input, and E u (k f , t) = k k f E u (k, t); in our DNS we use k f = 2. This yields a statistically steady state in which the mean value of the total energy dissipation rate per unit volume balances the power input, i.e. = P; (6) once this state has been established, we save 50 representative velocity-and magneticfield configurations over 36.08t I , 29.29t I , 32.61t I and 30.95t I , for R1D, R2D, R3D and R4D, respectively, where t I = I /u rms is the integral-scale eddy-turnover time. We use these configurations to obtain the statistical properties we describe below. For decaying MHD turbulence, we have carried out eight simulations with 512 3 collocation points and four simulations with 1024 3 collocation points. The parameters used in these simulations, which we have organized into three sets, are given in table 1.
In the first set of runs, R1-R5, we set the magnetic diffusivity η = 10 −3 and use five values of ν, namely 10 −4 , 5.0 × 10 −4 , 10 −3 , 5.0 × 10 −3 and 10 −2 , which yield Pr M = 0.1, 0.5, 1, 5 and 10. These runs have been designed to study the effects, on decaying MHD turbulence, of an increase in Pr M , with the initial energy held fixed: in particular, we use E 0 (4) for runs R1-R5. Given that this initial energy and η are both fixed, an increase in Pr M leads to a decrease in Re and thus an increase in k max η u d and k max η b d , as we discuss in detail later.
In our second set of decaying-MHD-turbulence runs, R3B, R4B and R5B, we increase E 0 in equation (4) as we increase ν and thereby Pr M , so that k max η u d 1 and k max η b d 1. Thus, in these runs, the inertial ranges in energy spectra extend over comparable ranges of the wavevector magnitude k.
Our third set of decaying-MHD-turbulence runs, R1C, R2C, R3C and R4C, uses 1024 3 collocation points and Pr M = 0.01, 0.1, 1 and 10, respectively. By comparing the results of these runs with those of R1-R5, R3B, R4B and R5B, we can check whether our qualitative results depend significantly on the number of collocation points that we use.
We have carried out another set of four runs, R1D, R2D, R3D and R4D, in which we force the MHD equations, as described above, until we obtain a NESS. These runs help us to compare the statistical properties of decaying and statistically steady turbulence. In these runs, we use 512 3 collocation points and ν and η such that Pr M = 0.01, 0.1, 1 and 10, respectively.
We calculate the kinetic, magnetic and total energy spectra E u (k) = k |k|=k |ũ(k)| 2 , E b (k) = k |k|=k |b(k)| 2 and E(k) = E u (k) + E b (k), respectively, the kinetic, magnetic and total energies E u = k E u (k)/2, E b = k E b (k)/2 and E = E u + E b and the ratio E b /E u . 7 Spectra for the Elsässer variables, energy dissipation rates and the effective pressure are, Our MHD simulations are characterized by the Taylor-microscale Reynolds number Re λ = u rms λ/ν, the magnetic Taylor-microscale Reynolds number Rm λ = u rms λ/η and the magnetic Prandtl number Pr M = Rm λ /Re λ = ν/η, where the root-mean-square velocity u rms = √ 2E u /3 and the Taylor microscale λ = [ k k 2 E(k)/E] −1/2 . We also calculate the integral length scale I = [ k E(k)/k]/E, the mean kinetic energy dissipation rate per unit mass, , the mean magnetic energy dissipation rate per unit , the mean energy dissipation rate per unit mass = u + b and the dissipation length scales for velocity and magnetic fields η u We calculate the eigenvalues u n and the associated eigenvectorsê u n , with n = 1, 2 or 3, of the rate-of-strain tensor S whose components are 3 denote the eigenvalues of the tensile magnetic stress tensor T, which has components T i j = −b i b j ; the corresponding eigenvectors are, respectively,ê b 1 ,ê b 2 andê b 3 . For incompressible flows n u n = 0, so at least one of the eigenvalues u n must be positive and another negative; we label them in such a way that u 3 is positive, u 1 is negative and u 2 lies in between them; note that u 2 can be positive or negative. We obtain PDFs of these eigenvalues; furthermore, we obtain PDFs of the cosines of the angles that the associated eigenvectors make with vectors such as u, ω, etc. These PDFs and those of quantities such as the local cross helicity H C = u · b help us to quantify the degree of alignment of pairs of vectors such as u and b [53]. We also compare PDFs of magnitudes of local vorticity ω, the current density j and local energy dissipation rates u and b to obtain information about intermittency in velocity and magnetic fields.
We also obtain several interesting joint PDFs; to the best of our knowledge, these have not been obtained earlier for MHD turbulence. We first obtain the velocity-derivative tensor A, also known as the rate-of-deformation tensor, with components A i j = ∂ i u j , and then the invariants Q = − 1 2 tr(A 2 ) and R = − 1 3 tr(A 3 ), which have been used frequently to characterize fluid turbulence [58]- [60]. The zero-discriminant line D ≡ 27 4 R 2 + Q 3 = 0 and the Q and R axes divide the Q R plane into qualitatively different regimes. In particular, regions in a turbulent flow can be classified as follows: when Q is large and negative, local strains are high and vortex formation is not favoured; furthermore, if R > 0, fluid elements experience axial strain, whereas if R < 0, they feel biaxial strain. In contrast, when Q is large and positive, vorticity dominates the flow; if, in addition, R < 0, vortices are compressed, whereas if R > 0, they are stretched. Thus, some properties of a turbulent flow can be highlighted by making contour plots of the joint PDF of Q and R; these Q R plots show a characteristic, tear-drop shape. We explore the forms of these and other joint PDFs, such as joint PDFs of u and b , in MHD turbulence.
To characterize intermittency in MHD turbulence, we calculate the order-p, equal-time, longitudinal structure functions S a p (l) ≡ |δa (x, l)| p , where the longitudinal component of the field a is given by δa (x, l) ≡ [a(x + l, t) − a(x, t)] · l l , where a can be u, b or one of the Elsässer variables. From these structure functions, we also obtain the hyperflatness F a 6 (l) = S a 6 (l)/[S a 2 (l)] 3 . For separations l in the inertial range, i.e. η u d , η b d l L, we expect S a p (l) ∼ l ζ p a , where ζ a p are the inertial-range multiscaling exponents for the field a; the Kolmogorov phenomenology of 1941 [23]- [25], henceforth referred to as K41, yields the simple scaling

Results
To set the stage for the types of studies we carry out for MHD turbulence, we begin with a very brief summary of similar and well-known results from studies of homogeneous, isotropic Navier-Stokes turbulence, which can be found, e.g., in [24,46,47], [61]- [67].

Overview of fluid turbulence
For ready reference, we show here illustrative plots from a DNS study that we have carried out for the three-dimensional Navier-Stokes equation by using a pseudospectral method, with 512 3 collocation points and the 2/3 rule for removing aliasing errors; here, ν = 0.001, Re λ 340 and k max η u d 0.3. In decaying fluid turbulence, energy is injected at large spatial scales, as described in the previous section for the MHD case. This energy cascades down till it reaches the dissipative scale at which viscous losses are significant. We study various statistical properties; these are given in points (i)-(vi) below.
(i) Plots of the energy E and the mean energy dissipation rate versus time show, respectively, a gentle decay and a maximum, as shown, e.g., by the full red and dotted blue curves in figure 1(a). This maximum in is associated with the completion of the energy cascade at a time t c ; the remaining properties (ii)-(vi) are obtained at t c . (ii) If Re λ is sufficiently large and we have a well-resolved DNS (i.e. k max η u d > 1), then at t c , the spectrum E(k) shows a well-developed inertial range, where at the K41 level E(k) ∼ k −5/3 , and a dissipation range, in which the behavior of the energy spectrum is consistent with E(k) ∼ k α exp(−βk), where α and β are non-universal, positive constants [62,65] and 5k d < k < 10k d , with k d = 1/η u d . An illustrative energy spectrum is shown by the dashed red line in figure 1(b); the blue dotted curve shows the compensated spectrum k 5/3 E(k); the associated dissipation or enstrophy spectrum (k) is shown in figure 1(c) and the inertial-range pressure spectrum [68], P(k) ∼ k −7/3 at the K41 level, is shown in figure 1(d). (Note that our DNS for the Navier-Stokes equation, which suffices for our purposes of illustration, does not have a well-resolved dissipation range because k max η u d 0.3 < 1; this is also reflected in the lack of a well-developed maximum in the enstrophy spectrum of figure 1(c).) (iii) Illustrative PDFs of the eigenvalues u n of the rateof-strain tensor S are given for n = 1, 2 and 3, respectively, by the full red, dashed green and dotted blue curves in figure 2(a); PDFs of the cosines of the angles that the vorticity ω and the velocity u make with the associated eigenvectorsê u n are given, respectively, in figures 2(b) and (c) via full red (n = 1), dashed green (n = 2) and dotted blue (n = 3) curves; these show that both ω and u have a tendency to be preferentially aligned parallel or antiparallel toê u 2 [60]; the PDF of the cosine of the angle between u and ω also indicates preferential alignment or antialignment of these two vectors, but with a greater tendency towards alignment, as found in experiments with a small amount of helicity [69] and as illustrated in figure 2(d). Finally, we give representative PDFs of the pressure p, the modulus of vorticity ω = |ω| and the local energy dissipation in figures 3(a)-(c), respectively; note that the PDF of the pressure is negatively skewed. (iv) Inertial-range structure functions S u p (l) ∼ l ζ p u show significant deviations [24] from the K41 result ζ u K 41 p = p/3, especially for p > 3. From these structure functions, we can obtain the hyperflatness F u 6 (l); this increases as the length scale l decreases, a clear signature of intermittency, as shown, e.g., in [49,66]. This intermittency also leads to non-Gaussian tails, especially for small l, in PDFs of velocity increments (see e.g. [66,70,71]), such as δu (l).
(v) Small-scale structures in turbulent flows can be visualized via isosurfaces [72] of, say, ω, and p, illustrative plots of which are given in figures 4(a)-(c); these show that regions of large ω are organized into slender tubes, whereas isosurfaces of look like shredded sheets; pressure isosurfaces also show tubes [37,47] but some studies have suggested the term 'cloud-like' for  them [61]. (vi) Joint PDFs also provide useful information about turbulent flows; in particular, contour plots of the joint PDF of Q and R, as in the representative figure 5, show a characteristic tear-drop structure.
The properties of statistically steady, homogeneous, isotropic fluid turbulence are similar to those described in points (ii)-(vi) in the preceding paragraph for the case of decaying fluid turbulence at cascade completion at t c . In particular, the strong-universality [42] hypothesis suggests that the multiscaling exponents ζ u p have the same values in decaying and statistically steady turbulence.
The remaining part of this section is devoted to our detailed study of the MHD-turbulence analogues of the properties (i)-(vi) summarized above; these are discussed, respectively, in the six sections 3.2-3.7.

Temporal evolution
We examine the time evolution of the energy, the energy-dissipation rates and related quantities, first for decaying and then for statistically steady MHD turbulence.   Figure 6 also shows similar plots for the mean energy dissipation rate (red full line), the mean kinetic-energy dissipation rate u (green dashed line) and the mean magnetic-energy dissipation rate b (blue dotted line) versus time t for runs R1-R5 (figures 6(a.2)-(e.2)) and runs R3B-R5B (figures 6(f.2)-(h.2)). In addition, figure 6 depicts the time evolution of the ratio E b /E u for runs R1-R5 (figures 6(a.3)-(e.3)) and runs R3B-R5B (figures 6(f.3)-(h.3)). We see from these figures that, for all the values of Pr M we have used, the energies E and E u decay gently with t but E b rises initially such that the ratio E b /E u rises, nearly monotonically, with t over the times we have considered; this is an intriguing trend that does not seem to have been noted earlier. The times over which we have carried out our DNS are comparable to the cascade-completion time t c that can be obtained from the peaks in the plots of , u and b versus t (figures 6(a.2)-(h.2)); by comparing these plots we see that, as we move from Pr M = 0.1 to Pr M = 10, with fixed η, we find that ( u − b ) and (t b − t u ) grow from negative values to positive ones because u increases with Pr M , where t b and t u are the positions of the cascade-completion maxima in b and u , respectively. We do not pursue the time evolution of our system well beyond t u and t b because the integral scale begins to grow thereafter and, eventually, can become comparable to the linear size of the simulation domain [46].

Spectra
We now discuss the behaviors of the energy, kinetic-energy, magnetic-energy, Elsässer variable, dissipation-rate and effective-pressure spectra, first for decaying and then for statistically steady MHD turbulence. In the former case, spectra are obtained at the cascade-completion time t c ; in the latter, they are averaged over the statistically steady state that we obtain.
We present compensated spectra of the total energy E c (k) = −2/3 k 5/3 E(k) (red full line), the kinetic energy E u c (k) = −2/3 k 5/3 E u (k) (green dashed line) and the total magnetic energy   d close to 1, so the dissipation ranges of these spectra span comparable ranges of k; however, as Pr M increases, more and more of the energy is concentrated in the magnetic field. These trends are not affected (a) if we increase the number of collocation points, as can be seen from the compensated spectra in figures 8(a.2)-(d.2) for runs R1C-R4C, which use 1024 3 collocation points, or (b) if we study energy spectra for statistically steady MHD turbulence, as can be seen from the compensated spectra in figures  have Pr M = 1; so they provide a convenient way of comparing the Re λ dependence of these spectra with a fixed value of Pr M = 1. All of the spectra in the subfigures of figure 8 have been compensated for by the 5/3 power of k and, to the extent that they show small, flat parts, their inertial-range, energy-spectra scalings are consistent with k −5/3 behaviors; other powers, such as −3/2, can also give small, flat parts in compensated spectra. A detailed error analysis is required to decide which power is most consistent with our data; we defer such an error analysis to section 3.5, where we carry it out for structure functions.
Compensated spectra of the Elsässer variables, namely show these spectra for statistically steady MHD turbulence in runs R1D-R4D, respectively. Note that the dissipation ranges of E + c (k) and E − c (k) overlap nearly on the scales of these figures. Differences between these are most pronounced at small k, where, typically, E − c (k) lies below E + c (k); these differences decrease with increasing Pr M if we hold the initial energy fixed as in figures 9(a.1)-(e.1) for runs R1-R5.
Next we come to the energy-dissipation (or enstrophy) spectra 1, respectively, are well suited for uncovering the functional forms of E u (k) and E b (k) in their dissipation ranges. In figures 11(a) and (b), we show, respectively, the kinetic-and magnetic-energy spectra E u (k) and E b (k) deep in their dissipation ranges for runs R5 and R1, respectively; our data for these spectra can be fitted to the form ∼k α exp(−βk) for k deep in the dissipation range and α and β non-universal numbers that depend on the parameters of the simulation; similar results have been obtained for fluid turbulence [62,65]. In particular, our data (figures 11(a) and (b)) for runs R5 and R1C are consistent with E u (k) ∼ k 2.68 exp(−0.235k), for 5k u We now turn to the spectra for the effective pressure P(k) (red full lines) and their compensated versions k 7/3 P(k) (blue dashed lines) that are shown at t c for runs R1-R5 (figures 12(a.1)-(e.1)) and R3B-R5B (figures 12(f.1)-12(h.1)) for decaying MHD turbulence; and for statistically steady MHD turbulence they are shown in figures 12(a.2)-(d.2) for runs R1D-R4D. Pressure spectra have been studied for fluid turbulence as, e.g., in [47,68]; to the best of our knowledge they have not been obtained for MHD turbulence. The compensated spectra here show that, for all of our runs, the inertial-range behaviors of these effective-pressure spectra are consistent with the power law k −7/3 ; this is consistent with the k −5/3 behaviors of the energy spectra discussed above. Furthermore, as Pr M increases from 0.1 to 10 in runs R1-R5, P(k) falls more and more rapidly, as can be seen from the vertical scales in figures 12(a.1)-(e.1).

Probability distribution functions
We calculate several PDFs to characterize the statistical properties of decaying and statistically steady MHD turbulence. In the former case, PDFs are obtained at the cascade-completion time t c ; in the latter, they are averaged over the statistically steady state that we obtain. The PDFs we consider are of two types: the first type are PDFs of the cosines of angles between various  vectors, such as u and ω; these help us to quantify the degrees of alignment between such vectors; the second type are PDFs of quantities such as u , b and the eigenvalues of the rate-ofstrain tensor.
In figure 13, we show plots of the PDFs of cosines of the angles between the vorticity ω and the eigenvectors of the fluid rate-of-strain tensor S, namelyê 1 u (red full line),ê 2 u (green dashed lines) andê 3 u (blue dotted lines) for runs R1-R5 and R3B-R5B at the cascade-completion time t c for the case of decaying MHD turbulence. In figure 14, we show similar plots of the PDFs of cosines of the angles between the current density j and the eigenvectors of the fluid rate-ofstrain tensor S. The most important features of these figures are sharp peaks in the green dashed lines; these show that there is a marked tendency for the alignment or antialignment of ω andê 2 u , as in fluid turbulence, and of a similar tendency for the alignment or antialignment of j andê 2 u ; these features do not depend very sensitively on Pr M . Furthermore, the PDFs of cosines of the angles between ω andê 1 u (blue dotted lines) and ω andê 3 u (red full lines) show peaks near zero in figure 13; in contrast, analogous PDFs for the cosines of the angles between j andê 1 u (red full lines) and ω andê 3 u (blue dotted lines) show nearly flat plateaux in the middle with very gentle maxima near −0.5 and 0.5 ( figure 14). Runs R1C-R4C and R1D-R4D yield similar PDFs, for the cosines of these angles, so we do not give them here.
Plots of the PDFs of cosines of the angles between the velocity u and the eigenvectors of the fluid rate-of-strain tensor S are given in figure 15; their analogues for b are given in figure 16. Again, the most prominent features of these figures are sharp peaks in the green dashed lines; these show that there is a marked tendency for the alignment or antialignment of u andê 2 u and of a similar tendency for the alignment or antialignment of j andê 2 u ; these features do not depend very sensitively on Pr M . The PDFs of cosines of the angles between u andê 1 u (red full line) and u andê 3 u (blue dotted lines) show gentle, broad peaks that imply a weak preference for angles close to 45 • or 135 • ; these peaks are suppressed as we increase Pr M (figures 15(a.1)-(e.1) for runs R1-R5) with fixed initial energy, but they reappear if we compensate for the increase in Pr M by increasing the initial energy (figures 15(f.1)-(h.1)). Similar, but sharper, peaks appear in the PDFs of cosines of the angles between u andê 1 u (red full lines) and u andê 3 u (blue dotted lines); these show a weak preference for angles close to 47 • or 133 • (figure 16). Some simulations of compressible MHD turbulence have noted the presence of such peaks [35] for Pr M = 1. The PDFs of figures 13-16 have a marginal dependence on Pr M . Furthermore, they look quite similar to those obtained earlier in convection-driven dynamos (see figure 15 of [36]).
Only one of the eigenvalues b 1 of the tensile magnetic stress tensor T is non-zero; and the corresponding eigenvectorê 1 b is identically aligned with b. Thus PDFs of cosines of angles between u, ω, j and b and the eigenvectors of T are simpler than their counterparts for S and are not presented here. Figure 17 shows     However, these PDFs are quite broad and distinctly non-Gaussian; this can be seen easily from the values of the mean µ H C , standard deviation σ H C , skewness γ 3,H C and kurtosis γ 4,H C given in table 2. Thus fluctuations of H C away from the mean are very significant. Table 2 also gives the value of the mean energy E and the ratio E/µ H C , which does not appear to be universal; for runs R1-R5 and R3B-R5B it lies in the range 0.23-0.26, for R1C-R2C in the range −0.04-0.04 and for R1D-R4D in the range 0.05-0.2. For all of our runs, with the exception of R2C, the mean µ H C and the skewness γ 3,H C are positive. Even if the PDF of H C had been a Gaussian, its mean value would have been within one standard deviation of 0; the actual PDF is much broader than a Gaussian. On symmetry grounds, there is no reason for the system to display a non-zero value for µ H C unless there is some bias in the forcing or in the initial condition (the latter for the case of decaying turbulence). In any given run, if there is some residual H C , it is reflected in a slight asymmetry in alignment (or antialignment) of u and b, which we have studied above via the PDF of the cosine of the angle between u and b. When we consider the ratio µ H C /E, it seems to be substantial in some runs but, given the arguments above, we expect it to vanish in runs with a very large number of collocation points; indeed, it is very small in runs R1C-R4C.  PDFs of 1 u and 3 u have long tails on the right-and left-hand sides, respectively. These tails shrink as we increase Pr M (figures 19(a.1)-(e.1) for runs R1-R5, respectively), by increasing ν while holding the initial energy fixed; thus, there is a substantial decrease in regions of large strain. However, if we compensate for the increase in ν by increasing the energy in the initial condition such that k max η u d and k max η b d are both 1, we see that these tails stretch out, i.e. regions of large strain reappear.
We show PDFs of the kinetic-energy dissipation rate u (blue dashed lines) and the magnetic-energy dissipation rate b (red full lines) that are obtained at t c for runs R1-R5 Table 2. The mean µ H C , standard deviation σ H C , skewness µ 3,H C and kurtosis µ 4,H C of the PDF of the cross helicity H C for our runs R1-R5 and R3B-R5B for decaying MHD turbulence at cascade completion; columns 6 and 7 give, respectively, the mean energy E and ratio of the means of the cross helicity and the energy, i.e. µ H C /E.   3) for runs R1, R1C and R1D, respectively). This indicates that large values of b are more likely to appear than large values of u and, given the long tails of these PDFs, suggests that, except at the smallest values of Pr M we have used, we might obtain more marked intermittency for the magnetic field than for the velocity field. Table 3. The mean µ u , standard deviation σ u , skewness γ 3, u and kurtosis γ 4, u of the PDFs of the local fluid energy dissipation u , and their analogues for b , for runs R1-R5, R3B-R5B, R1C-R4C and R1D-R4D. Furthermore, as we expect, the tail of the PDF of u is drawn in towards small values of u as we increase Pr M (figures 20(a.1)-(e.1) for runs R1-R5, respectively) while holding η and the initial energy fixed. However, if we compensate for the increase in ν by increasing the initial energy so that k max η u d and k max η b d are both 1, we see that the tails of the PDFs of b and u get elongated as we increase Pr M , e.g. in figures 20(f.1)-(h.1) for runs R3B-R5B, respectively. The values of the mean µ u , standard deviation σ u , skewness γ 3, u and kurtosis γ 4, u of the PDFs of the local fluid energy dissipation u are given for all our runs, and their counterparts for b are given in table 3. From these values, we see that the right tails of these distributions fall much more slowly than the tail of a Gaussian distribution.
Similar trends emerge if we examine the PDFs of the moduli of the vorticity and the current density, ω (blue dashed lines) and j (red full lines), respectively: these are presented at t c for runs R1- 3) for runs R1, R1C and R1D, respectively), so large values of j are more likely than large values of ω. Thus, given that these PDFs have long tails, it is reasonable to expect that, except at the smallest values of Pr M we have used, intermittency for the magnetic field might be larger than that for the velocity field. Moreover, the tail of the PDF of ω is drawn in towards small values of ω as we increase Pr M (figures 21(a.1)-(e.1) for runs R1-R5, respectively) while holding η and the initial energy fixed; but if, while increasing ν, we also increase the initial energy so that k max η u d and k max η b d are 1, we see that the tails of the PDFs of j and ω get stretched out as we increase Pr M , e.g. in figures 21(f.1)-(h.1) for runs R3B-R5B, respectively. The values of the mean µ ω , standard deviation σ ω , skewness γ 3,ω and kurtosis γ 4,ω of the PDFs of the modulus of the local vorticity ω for all of our runs and their counterparts for j are given in table 4. From these values, we see that the right tails of these distributions fall much more slowly than the tail of a Gaussian distribution.
We move now to PDFs of the local effective pressure (green full lines), which are shown at t c for runs R1-  Table 4. The mean µ ω , standard deviation σ ω , skewness γ 3,ω and kurtosis γ 4,ω of the PDFs of the modulus of the local vorticity ω, and their analogues for j, for runs R1-R5, R3B-R5B, R1C-R4C and R1D-R4D. skewness γ 3, p and kurtosis γ 4, p of the PDFs of the local effective pressure p are given for these runs in table 5. These have mean µ p = 0 but are distinctly non-Gaussian as can be seen from the values of γ 3, p and γ 4, p . Pressure PDFs are negatively skewed in pure fluid turbulence, as we have mentioned above; however, for MHD turbulence, we find that the PDFs of the effective pressurep can be positively skewed, as in runs R1-R5, R3B-R5B and run R4D, or negatively skewed, as in runs R1D-R3D; negative skewness seems to arise at low values of Pr M . The scale dependence of PDFs of velocity increments provides important clues to the nature of intermittency in fluid turbulence. To explore similar intermittency in MHD turbulence [73], we present data for the scale dependence of PDFs of velocity and magneticfield increments. As mentioned above, these increments are of the form δa (x, l) ≡ a(x + l, t) − a(x, t)] · l l , with a being either u or b, l = |l| the length scale and x an origin over which we can average to determine the dependence of the PDFs of δa on the scale l; for notational convenience, such velocity and magnetic-field increments are denoted by δu(l) and δb(l) in our plots. These PDFs are obtained at t c for runs R1-  of the red and black dashed lines in these plots indicates that the PDFs of the magnetic-field increments are broader than their velocity counterparts in most of our runs; this suggests, as we had surmised from the PDFs of energy-dissipation rates given above, that the magnetic field displays stronger intermittency than the velocity field at all but the smallest values of Pr M (figures 23(a.1), (a.2) and (a.3) for runs R1, R1C and R1D); the general trend that we note from these figures is that the magnetic-field intermittency is stronger than that of the velocity field at large magnetic Prandtl numbers but the difference between these intermittencies decreases as Pr M is lowered. We will try to quantify this when we present structure functions in section 3.5.

Structure functions
We continue our elucidation of intermittency in MHD turbulence by studying the scale dependence of order-p equal-time, velocity and magnetic-field longitudinal structure functions  and S b p (l) ∼ l ζ p b , where ζ u p and ζ b p are inertial-range multiscaling exponents for velocity and magnetic fields, respectively; if these fields show multiscaling, we expect significant deviations from the K41 result ζ u K 41 p = ζ bK 41 p = p/3. (Note that we do not expect any Iroshnikov-Kraichnan [74] scaling because we have no mean magnetic field in our simulations.) Given large inertial ranges, the multiscaling exponents can be extracted from slopes of log-log plots of structure functions versus l. However, in practical calculations, inertial ranges are limited, so we use the ESS procedure [40,41] in which we determine the multiscaling exponent ratios ζ u p /ζ u 3 and ζ b p /ζ b 3 , respectively, from slopes of log-log plots of (a) S u p versus S u 3 and (b) S b p versus S b 3 ; we refer to these as ESS plots. Our data for structure functions are averaged over 51 and 400 origins, respectively, for simulations with 512 3 and 1024 3 collocation points.
We begin with data from our decaying-MHD-turbulence runs R1C-R4C, which use 1024 3  ; the local slopes of these ESS curves are shown in the insets of these figures. Flat portions in these plots of local slopes help us to identify the inertial ranges. The regions that we have chosen for our fits are indicated by black horizontal lines with vertical ticks at their ends. In such a region, the mean value and the standard deviation of the local slope of the ESS plot for S u p (r ) (or S b p (r )) yield, respectively, our estimates for the exponent ratio ζ u p /ζ u 3 (or ζ b p /ζ b 3 ) and its error bars. Figures 24(a.3)-(d.3) show plots of these exponent ratios versus p for the velocity field (blue dotted line with thick error bars) and the magnetic field (red dashed line with thin error bars); the black solid line shows the K41 result for comparison. Although earlier studies [25,39] have obtained such exponents from DNS studies, they have done so, to the best of our knowledge, only for Pr M = 1; furthermore, they have not reported error bars. Although our (conservative) error bars are large, our plots of exponent ratios suggest the following: (a) deviations from the K41 result are significant, especially for p > 3, as in fluid turbulence; ; the local slopes of these ESS curves are shown in the insets of these figures. We obtain estimates for the exponent ratio ζ u p /ζ u 3 and ζ b p /ζ b 3 and their error bars, as in figure 24. Figures 25(a.3)-(d.3) show plots of these exponent ratios versus p for the velocity field (blue dotted line with thick error bars) and the magnetic field (red dashed line with thin error bars); the black solid line shows the K41 result for comparison. Plots versus l of the hyperflatnesses F u 6 (l) (red line) and F b 6 (l) (blue dashed line) are presented Table 6. Multiscaling exponent ratios ζ u bars. Thus, at least at this level of resolution and accuracy, we have strong universality of these exponent ratios, for a given value of Pr M , in as much as the ratios from decaying-MHD turbulence agree with those from the statistically steady case. The dependence on Pr M will be examined in section 4. We have tried to fit our data for the multiscaling exponent ratios to the generalized She-Leveque formula used in [25,75] to extract the codimensions of dissipative structures; these fits are not very good; they yield values for these codimensions that are close to or slightly lower than 1 and are closer to the estimates of [25] than those of [75]. (See [75] for a discussion of possible reasons for the discrepancies between these estimates.) Given the large error bars in multiscaling exponent ratios, we suggest that such fits are fraught with considerable uncertainty.

Isosurfaces
As we have mentioned in our discussion of fluid turbulence, isosurface plots of quantities such as ω, the modulus of the vorticity, give us a visual appreciation of small-scale structures in a turbulent flow; in fluid turbulence, iso-ω surfaces are slender tubes if ω is chosen to be well above its mean value [66,72]. For the case of MHD turbulence, it is natural to consider isosurface plots [76] of ω, the modulus j of the current density, energy dissipation rates and the effective pressure. 3) for runs R1D-R4D; these isosurfaces go through points at which the value of ω is two standard deviations above its mean value (for any given plot). For Pr M = 1, it has been noted in several DNS studies that such isosurfaces are sheets [5,25,76,77] and that there is a general tendency for such sheet formation in MHD turbulence; our results show that this tendency persists even when Pr M = 1. The number of high-intensity isosurfaces of ω shrinks as we increase Pr M (figures 26(a.1)-(e.1) for runs R1-R5, respectively), by increasing ν while holding the initial energy fixed. However, if we compensate for the increase in ν by increasing the energy in the initial condition such that k max η u d and k max η b d are both 1, we see that high-ω sheets reappear (figures 26(f.1)-(h.1) for runs R3B-R5B, respectively). These trends are also visible in our high-resolution, decaying-MHD-turbulence runs R1C-R4C (figures 26(a.2)-(d.2)) and the statistically steady ones, namely R1D-R4D (figures 26(a.3)-(d.3)). One interesting point that has not been noted before is that some tube-type structures appear along with the sheets at small values of Pr M , as can be seen by enlarging figure 26(a.3)  3) for runs R1D-R4D; these isosurfaces go through points at which the value of j is two standard deviations above its mean value (for any given plot). Again, the dominant features in these isosurface plots are sheets; their number goes down as Pr M increases with ν while the initial energy is held constant; but if this energy is increased, the number of high-intensity sheets increases.
Isosurfaces   3) for runs R1D-R4D; the isosurfaces go through points at which the value of b is two standard deviations above its mean value (for any given plot).
Here too, the isosurfaces are sheets; they lie close to, but are not coincident with, isosurfaces of ω and j; changes in Pr M affect these isosurfaces much as they affect isosurfaces of ω and j.  go through points at which the value ofp is two standard deviations above its mean value (for any given plot). The general form of these isosurfaces is cloud-type, to borrow the term that has been used for isosurfaces of the pressure in fluid turbulence [61]. Here also changes in Pr M affect these isosurfaces much as they affect isosurfaces of ω and j, in as much as high-intensity isosurfaces are suppressed as Pr M increases via an increase in ν, unless this is compensated for by an increase in the initial energy (in the case of decaying MHD turbulence) or Re λ .

Joint probability distribution functions (PDFs)
In this subsection, we present three sets of joint PDFs that have, to the best of our knowledge, not been used to characterize MHD turbulence previously. The first of these is a Q R plot that is often used in studies of fluid turbulence, as we have discussed in sections 2.2 and 3.1; the next is a joint PDF of ω and j; and the last is a joint PDF of u and b .
We show Q R plots, i.e. joint PDFs of Q and R, via filled contour plots; these are obtained at t c for runs R1- 3) for runs R1D-R4D; the black curve in these plots is the zero-discriminant line D ≡ 27 4 R 2 + Q 3 = 0. These Q R plots retain overall, aside from some distortions, the characteristic tear-drop structure familiar from fluid turbulence (see section 3.1 and figure 5). If we recall our discussion of Q R plots in section 2.2 and we note that, as we increase Pr M (figures 31(a.1)-(e.1) for runs R1-R5, respectively) while holding η and the initial energy fixed, there is a general decrease in the probability of having large values of Q and R, i.e. regions of large strain or vorticity are suppressed; this corroborates what we have found from the PDFs and isosurfaces discussed above. However, if we compensate for the increase in ν by increasing the initial energy, or Re λ , so that k max η u d and k max η b d are both 1, we see that Q and R can increase again. Note that when Pr M is very small, as in run 3) for runs R1D-R4D. All of these joint PDFs have long tails; as we move away from Pr M = 1, they become more and more asymmetrical. Furthermore, as we expect, the tails of these PDFs are drawn in towards small values of ω and j as we increase Pr M (figures 32(a.1)-(e.1) for runs R1-R5, respectively) while holding η and the initial energy fixed. However, if we compensate for the increase in ν by increasing the initial energy or Re λ , so that k max η u d and k max η b d are both 1, we see that the tails of the PDFs become elongated again. In the end, we consider joint PDFs of u and b that are obtained at t c for runs R1-R5 in figures 33(a.1)-(e.1), runs R3B-R5B in figures 33(f.1)-(h.1) and runs R1C-R4C in figures 33(a.2)-(d.2) for decaying MHD turbulence; and for statistically steady MHD turbulence they are shown in figures 33(a.3)-(d.3) for runs R1D-R4D. The trends here are similar to the ones discussed in the previous paragraph. In particular, these joint PDFs have long tails; as we move away from Pr M = 1, they become more and more asymmetrical; and the tails of these PDFs are drawn in towards small values of u and b as we increase Pr M (figures 33(a.1)-(e.1) for runs R1-R5, respectively) while holding η and the initial energy fixed. But if we make up for the increase in ν by increasing the initial energy or Re λ so that k max η u d and k max η b d are both 1, we see that the tails of the PDFs become elongated again.

Discussions and conclusion
We have carried out an extensive study of the statistical properties of both decaying and statistically steady homogeneous, isotropic MHD turbulence. Our study, which has been designed specifically to study the systematics of the dependence of these properties on the magnetic Prandtl number Pr M , uses a large number of statistical measures to characterize the statistical properties of both decaying and statistically steady MHD turbulence. Our study is restricted to incompressible MHD turbulence; we do not include a mean magnetic field as, e.g., in [33]; furthermore, we do not study Lagrangian properties considered, e.g., in [34]. In our studies, we obtain (a) various PDFs, such as those of the moduli of the vorticity and current density, the energy dissipation rates, of cosines of angles between various vectors and scaledependent velocity and magnetic-field increments, (b) spectra, e.g. those of the energy and the effective pressure, (c) velocity and magnetic-field structure functions that can be used to characterize intermittency, (d) isosurfaces of quantities such as the moduli of the vorticity and current, and (e) joint PDFs, such as Q R plots. The evolution of these properties with Pr M has been described in detail in the previous section.
To the best of our knowledge, such a comprehensive study of the Pr M dependence of incompressible, homogeneous, isotropic MHD turbulence, both decaying and statistically steady, has not been attempted before. Studies that draw their inspiration from astrophysics often consider anisotropic flows [78]- [85], flows that are compressible [35,86] or flows that include a mean magnetic field [33,83,87,88]. Yet other studies concentrate on the alignment between various vectors such as u and b as, e.g., in [28,35,53]; some of these include a few, but not all, of the PDFs we have studied; and, typically, these studies are restricted to the case Pr M = 1. Some of the spectra we study have been obtained in earlier DNS studies but, typically, only for the case Pr M = 1; a notable exception is [38], which examines the Pr M dependence of energy spectra but with a relatively low resolution. The papers [6,32,89] have also considered some Pr M dependence but not for low Pr M . Isosurfaces of the moduli of vorticity and current density have been obtained earlier [25,39,76] for the case Pr M = 1. The Pr M dependence of these and other isosurfaces is presented here for the first time. The joint PDFs we have shown above have also not been investigated previously.
Here we wish to highlight, and examine in detail, the implications of our study for intermittency. Some earlier DNS studies, such as [30], noted that, for the case Pr M = 1, the magnetic field is more intermittent than the velocity field. This is why we have concentrated on velocity and magnetic-field structure functions. Our study confirms this finding, for the case Pr M = 1. This can be clearly seen from the comparison of our exponent ratios, for Pr M = 1, with those of the recent DNS of decaying-MHD turbulence in [30] in table 8; the error bars that we quote for our exponent ratios have been calculated as described in the previous section; we have obtained exponent ratios of the paper [30] by digitizing 5 the data in their plot (figure 3 of [30]) of multiscaling exponents versus the order p (error bars are not given in their plot). Thus, at least given our error bars, there is agreement between our exponent ratios, both for decaying and statistically steady MHD turbulence, and those of [30] for Pr M = 1. We note in passing that the latter DNS is one of decaying MHD turbulence but with a very special initial condition, which allows an effective resolution greater than that we have obtained; however, the initial condition we use in our decaying-MHD-turbulence DNS is more generic than that of [30]. It is Table 8. A comparison of multiscaling exponent ratios ζ u our expectation that non-universal effects, associated with different initial conditions [26,50], might not affect multiscaling exponent ratios, except if we use non-generic, power-law initial conditions [26] in which E(k) grows with k (at least until some large-k cut-off).
Direct numerical simulations of decaying-MHD turbulence, e.g. those of [25,30], often average data obtained from field configurations at different times that are close to the time at which the peak appears in plots of the energy-dissipation rate. This is a reasonable procedure, for Pr M = 1, because the temporal evolution of the system is slow in the vicinity of this peak. We have not adopted this procedure here because, as we move away from Pr M = 1, the cascadecompletion peaks occur at different times in plots of u and b , as we have discussed in detail in earlier sections of this paper.
Let us now turn to the Pr M dependence of the multiscaling exponent ratios shown in tables 6 and 7 and in figures 24(a.3)-(d.3) and 25(a.3)-(d.3). Even though our error bars are large, given the conservative, local-slope error analysis we have described in the previous section, a trend emerges: at large values of Pr M , the magnetic field is clearly more intermittent than the velocity field, in as much as the deviations of ζ b p from the simple-scaling prediction are stronger than their counterparts for ζ u p . However, the velocity field becomes more intermittent than the magnetic field as we lower Pr M . Could this result, namely the dependence of our multiscaling exponent ratios on Pr M , be an artifact? We believe not. As we have discussed above, dissipation ranges in our spectra are adequately resolved; furthermore, we have determined exponent ratios from a rather stringent local-slope analysis, which is rarely presented in earlier DNS studies of MHD turbulence. Ultimately, of course, this Pr M dependence of multiscaling exponents in MHD turbulence must be tested in detail in very-highresolution DNS studies of MHD turbulence; such studies should become possible with the next generation of supercomputers. In particular, such supercomputers should allow us to achieve high Reynolds numbers along with higher values k max η u d and k max η b d than we have been able to obtain in runs R1D and R2D so that we have both well-resolved inertial and dissipation ranges.
For a fixed value of Pr M , we conjecture that, as in fluid turbulence, the energy-spectrum exponents and multiscaling exponent ratios do not depend on the Re λ if the latter is large enough to ensure that we have fully developed, homogeneous, isotropic MHD turbulence with a substantial inertial range. Our results for the spectral exponents (see e.g. figures 8(c.1), (g.1), (c.2) and (c.3) for runs R3 (Re λ = 121), R3B (Re λ = 210), R3C (Re λ = 172) and R3D (Re λ = 239), respectively) are consistent with this conjecture, as are the multiscaling exponent ratios given in table 8 for Pr M = 1.
It is useful to note at this stage that a recent experimental study of MHD turbulence in the solar wind [57] provides evidence for velocity fields that are more strongly intermittent than the magnetic field; this study does not give the value of Pr M . However, its data for multiscaling exponents are qualitatively similar to those we obtain at low values of Pr M . Furthermore, PDFs of H C have also been obtained from solar-wind data [56]; these are similar to the PDFs we obtain for H C . Of course, we must exercise caution in comparing results from DNS studies of homogeneous, isotropic, incompressible MHD turbulence with measurements on the solar wind, where anisotropy and compressibility can be significant; and, for the solar wind, we might also have to consider kinetic effects that are not captured by the MHD equations.
The last point we wish to address is the issue of strong universality of exponent ratios. In the fluid-turbulence context, such strong universality [42,43] implies the equality of exponents (and, therefore, their ratios) determined from decaying-turbulence studies (say at the cascadecompletion time) or from studies of statistically steady turbulence. Does such strong universality have an analogue in MHD turbulence? Our data, for any fixed value of Pr M in tables 6 and 7, are consistent with such strong universality of multiscaling exponent ratios in MHD turbulence; but, of course, our large error bars imply that a definitive confirmation of such strong universality in MHD turbulence must await DNS studies that might become possible in the next generation of high-performance computing facilities.