Transition to Equilibrium and Coherent Structure in Ideal MHD Turbulence, Part 2

: We continue our study of the transition of ideal, homogeneous, incompressible, magne-tohydrodynamic (MHD) turbulence from non-equilibrium initial conditions to equilibrium using long-time numerical simulations on a 128 3 periodic grid. A Fourier spectral transform method is used to numerically integrate the dynamical equations forward in time. The six runs that previously went to near equilibrium are here extended into equilibrium. As before, we neglect dissipation as we are primarily concerned with behavior at the largest scale where this behavior has been shown to be essentially the same for ideal and real (forced and dissipative) MHD turbulence. These six runs have various combinations of imposed rotation and mean magnetic ﬁeld and represent the ﬁve cases of ideal, homogeneous, incompressible, and MHD turbulence: Case I (Run 1), with no rotation or mean ﬁeld; Case II (Runs 2a and 2b), where only rotation is imposed; Case III (Run 3), which has only a mean magnetic ﬁeld; Case IV (Run 4), where rotation vector and mean magnetic ﬁeld direction are aligned; and Case V (Run 5), which has non-aligned rotation vector and mean ﬁeld directions. Statistical mechanics predicts that dynamic Fourier coefﬁcients are zero-mean random variables, but largest-scale coherent magnetic structures emerge and manifest themselves as Fourier coefﬁcients with very large, quasi-steady, mean values compared to their standard deviations, i.e., there is ‘broken ergodicity.’ These magnetic coherent structures appeared in all cases during transition to near equilibrium. Here, we report that, as the runs were continued, these coherent structures remained quasi-steady and energetic only in Cases I and II, while Case IV maintained its coherent structure but at comparatively low energy. The coherent structures that appeared in transition in Cases III and V were seen to collapse as their associated runs extended into equilibrium. The creation of largest-scale, coherent magnetic structure appears to be a dynamo process inherent in ideal MHD turbulence, particularly in Cases I and II, i.e., those cases most pertinent to planets and stars. Furthermore, the statistical theory of ideal MHD turbulence has proven to apply at the largest scale, even when dissipation and forcing are included. This, along with the discovery and explanation of dynamically broken ergodicity, is essentially a solution to the ‘dynamo problem’.


Introduction
"One cannot deal with the physics of cosmic fluids without encountering at almost every turn the problem of turbulence."[1] In particular, one encounters the problem of magnetohydrodynamic (MHD) turbulence, since 'cosmic fluids', such as those found in planetary interiors and in stars, are electrically conducting and stirred into turbulence by buoyant forces.A universal feature of these objects is that they are rotating and often possess a quasi-stationary, mostly dipole magnetic field, i.e., a magnetic coherent structure.Just over a hundred years ago, ref. [2] hypothesized that internal magnetic fields coupled to fluid motions within the Sun, and by extension, the Earth, were responsible for the creation and maintenance of global magnetic dipole fields.Determining exactly how this was done came to be called 'the dynamo problem'.
In regard to Earth, numerical simulations of the geodynamo [3][4][5] established that MHD processes within the Earth are capable of creating a magnetic field similar to the actual geomagnetic field, including reversals of the dominant dipole component.In addition, there have been many laboratory experiments [6], some of which have shown the growth of a self-generated magnetic field, i.e., a dynamo [7][8][9].Although numerical simulations have been successful in proving the MHD nature of the geodynamo and experiments have shown a dynamo effect, the fundamental MHD origin of the dominant, quasi-steady geomagnetic dipole field still appear to be a theoretical mystery [10].
Our approach to resolving this mystery is to study three-dimensional (3-D) homogeneous, incompressible MHD turbulence, i.e., to examine a turbulent magnetofluid confined to a 3-D periodic box (a three-torus).This allows for the use of Fourier expansions to represent turbulent velocity and magnetic fields in order to study the nonlinear dynamics of a magnetofluid without the complicating factor of boundary conditions; in fact, introducing boundaries with no-slip conditions also requires introducing compressibility [11] and further complicates the problem.Fourier methods transform the problem from the partial differential equations of MHD (presented in Section 2.1) in x-space (physical space) to a dynamical system of nonlinear, coupled ordinary differential equations in k-space (Fourier space), which describe the evolution of the Fourier modes associated with velocity and magnetic field; this is discussed in Section 2.3.In x-space, the continuum magnetofluid is represented on an N 3 grid of velocity u(x, t) and magnetic field b(x, t) values, while in k-space, it is represented by the velocity and magnetic vector Fourier coefficients ũ(k, t) and b(k, t), where the wavevectors k define length-scales λ = 2π/k, k = |k| ≥ 1, in the turbulent magnetofluid and explicitly separate the largest length-scale modes from all of the smaller length-scale modes.This is particularly useful as the largest length-scale modes are analogous to the magnetic dipole components in a planet or star.
Viscous and ohmic dissipation occur in a real magnetofluid but 'cosmic fluids' generally have very large Reynolds numbers so that one may remove dissipation from a model system and study ideal homogeneous MHD turbulence for the following reason.We neglect dissipation because we are primarily concerned with behavior of a turbulent magnetofluid at the largest scale; large-scale dynamics have been shown to be essentially the same for ideal and real (forced and dissipative) MHD turbulence [12,13].The real MHD equations appear in (1) and (2) of Section 2.1 and become the ideal MHD equations when we set the viscosity ν and magnetic diffusivity η equal to zero.There are differences between solutions of the real and ideal MHD equations as the highest derivative terms (the dissipative terms) in the real MHD equations methods are removed to produce the ideal MHD equations; the connection between real and ideal solutions is examined by singular perturbation theory, e.g., [14].These differences are critical in real fluid and MHD flows since there may be thin regions, such as boundary layers, where the second-order dissipative terms become dominant, or in regions far from boundaries where there is a direct cascade of energy to the dissipation length-scale and where, again, dissipative terms are dominant.However, ideal results are often very useful.For example, in studying fluid turbulence, while it is usually necessary to retain the dissipative term in the Navier-Stokes equation because it determines the wavenumber spectrum in the medium-scale inertial and small-scale dissipative ranges, it is nevertheless possible to use ideal results to draw conclusions on the nature of real, helical fluid turbulence [15].
The fundamental difference between fluid and MHD turbulence is induced by the presence of a dynamic magnetic field, in that fluid turbulence only has a direct cascade to higher wavenumbers [16], while MHD turbulence has an inverse cascade as well as a direct one [17].This difference between fluid and MHD turbulence, due to the MHD inverse cascade, can create an energetically dominant large-scale magnetic field, depending on the particular case of MHD turbulence under consideration; these cases are given in Table 1 and the various helicities listed there are defined by equations described in the next two paragraphs and discussed in more detail in Appendix A. Ideal MHD statistical theory Appendix B predicts that in Cases I and II of Table 1, the large-scale magnetic field is quasi-stationary and has energy that is essentially equal to magnetic helicity, a result that also seems to apply to Case IV.These ideal MHD spectra are strongly peaked at the largest length scale, i.e., smallest wavenumber, where dissipation is weakest in the real case.That this survives the addition of forcing and dissipation to the magnetofluid has been shown in studies of helically forced, dissipative MHD turbulence using Fourier method simulations [12,13].Thus, if we wish to study the large-scale dynamics and structure of MHD turbulence, either ideal or real MHD turbulence simulations can be employed.Here, we choose ideal MHD.

Case
Mean Field Rotation Invariants 128 3 Runs Ideal MHD turbulence, as represented by a finite set of Fourier modes, is a conservative system with energy and up to two different helicities being integral invariants; again, the different cases are shown in Table 1.The number and form of the helical invariants depend on whether or not rotation or a mean magnetic field or both are imposed on the system (a mean magnetic field B o is one that is spatially and temporally constant and frame rotation is given by a vector Ω o ).Energy E, cross helicity H C , magnetic helicity H M and parallel helicity H P , as well as other quantities, are fully defined in Appendix A; here, we will, for easy reference and to serve as a 'nomenclature section', briefly discuss how E, H C , H M and H P , as well as other integrals, are determined.
First, we define the volume average of a quantity Φ(x, t) multiplied by a quantity Ψ(x, t) as {ΦΨ}, which is an integral over the volume of interest, here, a periodic box of side length 2π: The k-space functions Φ * (k, t) and Ψ(k, t) are Fourier transforms of the x-functions Φ(x, t) and Ψ(x, t) (Fourier transforms are defined in Section 2.2).Energy and the various helicities listed in Table 1 are then (H M is defined in terms of the magnetic vector potential a, where ∇ × a = b, ∇ • a = 0) Again, velocity u and magnetic fields b can be represented in x-space or k-space; their basic equations are given in (1) and (2) for x-space and ( 23) and (24) for k-space, and either set can be used to prove invariance.In terms of the mean magnetic field and rotation vector magnitudes Other quadratic forms, in addition to those given above, are also useful measures of the turbulent state: Above, we have E M (magnetic energy), E K (kinetic energy), H K (kinetic helicity, where ω = ∇ × u is the vorticity), A (mean squared vector potential), Ω ('enstrophy', i.e., mean squared vorticity) and J (mean squared electric current, where j = ∇ × b); please note that enstrophy Ω is different from rotation vector magnitude Ω o .For the computer runs discussed here, time-averages of these appear in Table 2.
Ω std 1.7385e+00 3.6272e+00 2.7515e+00 4.4614e+00 5.4717e+00 1.3995e+00 The ideal invariants given in Table 1 can be used to create a multi-invariant canonical probability density function and a statistical mechanics applicable to ideal MHD turbulence (Appendix B).This statistical theory was initiated by [18], applied to fluid turbulence by [16] and to MHD turbulence by [17,19]; these initial results were extended later, e.g., [20][21][22][23][24].This extension is relevant to the planets and stars as it provides an, at least qualitative, explanation for the emergence of a dominant, quasi-stationary, dipole magnetic field through broken ergodicity, which is defined by [25] as occurring if 'In a system that is nonergodic on physical timescales the phase point is effectively confined in one subregion or component of phase space'.The phase space itself is defined by the real and imaginary parts of the coefficients of the independent Fourier modes in the model system.
Numerical simulations have played an essential part in MHD turbulence research, starting with [26] using a 64 3 grid (which was large for that time), while advances in computer hardware have allowed more recent simulations to have much larger grid sizes, e.g., 1024 3 [27].However, these simulations were only run for short simulation times, with t max = 3 in the former and t max = 10 in the latter, as the emphasis was on large grid sizes rather than on longer simulation time.In contrast, we have focused on long-time runs with t max ∼10 3 to gather statistics and this dictated smaller grid sizes of 64 3 [22][23][24]; there is always this trade-off between grid size and simulation time.To illustrate the trade-off between grid size and maximum simulated time, ref. [28] presented runs from 64 3 up to 512 3 , but as the grid size went up, the time-step size decreased and number of time-steps increased, so that simulation time (number of 'eddy turnover times') decreased and only for the 64 3 run was quasi-equilibrium attained (see their figure 5).
Improved computational resources have allowed us to increase our grid size to 128 3 while still keeping t max ∼ 10 3 .We have recently run six of these 128 3 simulations, each for over a million time-steps, in a reasonable amount of calendar time (over half a year each) and they allow us to examine the transition of ideal MHD turbulence for each of the cases in Table 1, where the runs associated with each case are indicated.All of our 128 3 runs have moved through or are near the end of their transition periods and appear either to achieve equilibrium or to approach it asymptotically.Here, we focus on this initial period of transition to near-equilibrium and plan to examine the full equilibrium phase and present results in the future when these become available.
From our earliest numerical simulations [20] to those presented here, we have consistently seen the lowest-wavenumber (largest length-scale) modes grow from very small initial values into energetic, quasi-stationary, coherent structures.To understand this, it is necessary to develop the statistical theory of ideal MHD turbulence beyond that introduced by [17].In particular, we analyzed the entropy of the model dynamical system to discover the origin of this energetic structure; we have presented this analysis before [24,29] and will place it here in Appendix B. Furthermore, the statistical mechanics of an ideal, turbulent magnetofluid in a periodic box is essentially the same as that of one contained in a spherical shell [30] and thus, the smallest wavenumber Fourier modes are surrogates for geomagnetic dipole components.
Key discoveries in our past work are: (i) That the energy of a quasi-stationary magnetic dipole in Cases I and II is directly proportional to the absolute value of the magnetic helicity of the equilibrium turbulent magnetofluid [29,31]; this also appears to apply to Case IV, as will be shown here; and (ii) that the assumed ergodicity of the statistical ensemble is actually broken at the largest scale [20,24].Again, these ideal results have been seen to apply at the largest length-scale to real (i.e., forced and dissipative) MHD turbulence [12,13]; thus, modeling and theory provide a solution to the 'dynamo problem,' at least for the model systems considered here.(Again, ideal MHD statistical theory is reviewed in Appendix B of this paper for ease of referral).
Next, following a brief discussion of the mathematical model and numerical procedure, new computational results drawn from 128 3 simulations are presented.These are followed by a discussion of these results and their relevance to the dynamo problem.Lastly, we offer a conclusion to summarize our results and emphasize their importance.

Mathematical Model 2.1. Basic Equations
Magnetohydrodynamics (MHD) is the physics of electrically conducting fluids and principles and applications can be found in, e.g., [32,33], as well as in many other texts.The nondimensional form of the 3-D incompressible MHD equations in a rotating frame of reference with constant angular velocity Ω o and mean magnetic field B o (i.e., constant in space and time) can be written as These are described, for example, by [22]; these equations can be augmented to include thermal effects in the Boussinesq approximation [34], but we will not consider this here.Above, u(x, t) and b(x, t) are the turbulent velocity and magnetic fields, respectively.Velocity and magnetic fields are solenoidal: ∇ • u = ∇ • b = 0, as is appropriate for laboratory experiments and the Earth's outer core [3].The vorticity ω(x, t) and electric current density j(x, t) are defined by Nondimensional density does not appear in (1) because it equals unity.The symbols ν in (1) and η in (2) are shorthand for 1/R E and 1/R M , i.e., the inverses of the kinetic and magnetic Reynolds numbers, respectively.In the dimensional form of the equations, ν is the kinematic viscosity, while η is the magnetic diffusivity and ν = η = 0 for ideal MHD turbulence.Again, we avoid the complication of boundary conditions by placing the magnetofluid in a periodic box and expanding the various fields in terms of Fourier series.
In studying MHD turbulence, it is also possible to use other orthogonal function expansions, such as those appropriate to a cylinder [35], sphere [36] or spherical shell [30].However, the statistical mechanics of ideal MHD turbulence has the same form and predictions as in the Fourier case.In addition, the Fourier spectral transform methods used to simulate MHD turbulence in a periodic box are much more computationally efficacious in that they allow for larger grid sizes (more dynamical resolution) and longer runs (for time averaging).In fact, no spectral transform methods for these other geometries yet exist (the creation of which may be considered a 'grand computational challenge') and nontransform methods [36,37], in comparison, will always have insufficient resolution.In essence, Fourier transform methods in a periodic box are efficient and also serve as a surrogate for methods based on other orthogonal function expansions more suited to other geometries.

Fourier Representation
Discrete Fourier transforms for u and b, connecting x-space to k-space, are Here, N is the number of grid points in each x-space dimension, so we have a grid of N 3 points.The set of positions x and wave vectors (modes) k appearing in (4) and ( 5) are The components n j and k j , j = x, y, z, are integers.The n j satisfy 0 ≤ n j < N, while the integers k j lie in the range −N/2 + 1 and +N/2; thus, there are N 3 points in both spaces.The dot product k • x in (4) and ( 5) is The Fourier coefficients ũ(k, t) and b(k, t) are nonzero only for 1 ≤ k ≤ K < N/2, where k = |k|.The exact value of K is set by a de-aliasing requirement of [38]: K 2 is the largest integer such that K 2 ≤ 2N 2 /9, e.g., for N = 128, K 2 = 3640.
Note that the reality of u(x, t) and b(x, t) The vector fields u(x, t) and b(x, t) can be defined in a similar manner.
The periodic box of the model system is a cube of equal length edges.Although we will not do so here, we could have allowed for the possibility that it might be elongated or compressed in the x, y or z directions, which simulates a difference in oblateness that may occur in the Earth's inner and outer cores.We ignore this possibility because we have examined it and the dynamical effect is insignificant when compared to rotation [29].
Dynamically interacting coefficients fit within a sphere of radius K in k-space; all coefficients outside this k-sphere, as well as at k = 0, are initially zero and remain so during any numerical simulation.Since ũ(k, t) = ũ * (−k, t) and b(k, t) = b * (−k, t), only half of the k that satisfy 1 ≤ k 2 ≤ K 2 identify independent coefficients.Therefore, the number of independent modes k is M ∼ = 2πK 3 /3; for example, with N = 128 and K 2 = 3640, we have M ∼ = 459,950, while the actual numerical count is M = 459,916.
In k-space, the requirements ∇ Thus, ũ(k, t) and b(k, t) have two independent complex vector coefficients each, which can be defined as follows: First, determine a coordinate system for each k by starting with a unit vector ê3 (k) = k/k = k; then choosing a unit vector ê1 (k) orthogonal to ê3 (k), and then the remaining unit vector ê2 (k) is a vector product of the first two, forming a right-handed orthonormal basis for each k: (This is similar to a Craya-Herring decomposition [39]).Explicit choices [23] for any of the êj (k) will be mentioned as needed.
In terms of the êj (k) defined above, the Fourier vector coefficients (9) and ( 10) are However, an equivalent but very useful helical representation is: Here, the positive and negative helicity unit vectors and components are Note that ê * ± (k) = ê∓ (k).The orthonormality properties of the ê± (k An important property of the helical unit vectors concerns the curl operation: The vorticity and current in helical form are then Thus, the helical ± components of vorticity ω± (k, t) = ±k ũ± (k, t) and current j± (k, t) = ±k b± (k, t) are directly connected to velocity and magnetic field ± helical components.The helical representation is very useful for theoretical analysis and graphical presentation of numerical data, although numerical simulation is actually performed using the Cartesian components ( 9) and (10).

A Dynamical System
The Fourier-transformed 3-D MHD equations are found by placing expansions for ω(x, t) and b(x, t) of the form ( 4) into ( 1) and ( 2).The result is a set of coupled, nonlinear ordinary differential equations that represents a dynamical system in the sense of [40]: The nonlinear terms denoted by S are vector convolutions: The double summation in ( 25) is over all wavevectors p and q inside the truncation volume in k-space that satisfy p + q = k.In 3-D numerical simulations, these equations are integrated as described in Section 3.For N = 128, the phase space Γ of the dynamical system ( 23) and ( 24) has a dimension of N Γ = 3,679,328, as discussed in Appendix B. For ideal MHD, Table 1 indicates that the dynamical system conserves energy and various helicities depending on whether B o and Ω o are zero or not; these quantities and others are defined in Appendix A. The important Case II of Table 1, where Ω o = 0 and B o = 0, applies to the Earth and most other planets (and stars), so in these cases, energy and magnetic helicity are of primary interest.The Cases III, IV and V, where B o = 0, are appropriate for laboratory experiments and perhaps compact regions of solar wind turbulence where B o is the local interplanetary magnetic field.

Linear Modes
Using ( 11)-( 18), the MHD Equations ( 23) and ( 24) can be linearized to yield These can be put into a matrix form: It is a straightforward procedure to find linear eigenmodes of ( 28) by diagonalizing M ± ; there is one transform matrix for each k: The solution to (28) are the eigenmodes where the time dependence of the linear modes is in the matrix In the linear case, V ± (k, t) is constant, while in the nonlinear case, it becomes timedependent.In slightly more detail, V ± (k, t) has the form
If there is no nonlinearity, Ṽ± 1 (k, t) and Ṽ± 2 (k, t) of (34) will have constant, complex values determined by the initial conditions; they are the coefficients of traveling helical waves∼exp If there is nonlinearity, then these coefficients will generally change with time and it is these changes that contain information about MHD turbulence.

Numerical Procedure
A Fourier spectral transform method [41] on an N 3 grid with N = 128 is used; the minimum wave number is k = |k| = 1 and the maximum wave number is K = √ 3640 60.34.Time-integration is performed with a third-order Adams-Bashforth-Adams-Moulton method [42] with a time-step of ∆t = 0.0005; initial, non-equilibrium magnetic and kinetic modal energy values (spectra where k o = 6.Viscosity and magnetic diffusivity are set to zero so that the flow is ideal.A grid size of 128 3 was used so that the six single core Runs listed in Table 1 could be completed in a reasonable amount of time with the resources available, the Hopper Cluster at George Mason University, with each simulation running at ≈11 s/∆t; thus, a single run of 2 × 10 6 ∆ts requires about 36 weeks of CPU time. As mentioned, six computer simulations covering the five cases in Table 1 have been run and are identified in that table.The ideal invariants for each case are quadratic forms (global quantities) with terms that are scalar products of the vector Fourier coefficients ũ(k, t) and b(k, t), with 0 < k ≤ K, as defined in Appendix A. The partial differential equations for MHD in x-space are given by ( 1) and (2), while the transformed set of ordinary differential equations in k-space are given by ( 23) and (24).The set of equations in k-space is a finite dynamical system, as discussed in Section 2.3.The k-space Equations ( 23) and ( 24) are numerically integrated to advance the ũ(k, t) and b(k, t), as described in Section 3.
As seen in Table 1, the integral invariants of ideal MHD turbulence are the volumeaveraged energy E and magnetic helicity H M when B o = 0, as well as the cross helicity H C when Ω o = 0 and H P when Ω o = σB o = 0.In a numerical simulation, these ideal invariants typically have a standard deviation of less than 1% per million time-steps, while kinetic helicity H K , though an ideal invariant for hydrodynamic turbulence, falls to zero very quickly and then has small fluctuations about that value, as Table 2 shows.We mentioned that these runs came to 'near equilibrium'.To understand what this means, consider Figure 1 where we show, at different times for Runs 1 and 4, the average magnetic and kinetic energy spectra, ĒM (k 2 ) and ĒK (k 2 ), for modes k with the same value of k 2 = |k| 2 .There are 3036 different values of k 2 , where 1 ≤ k 2 ≤ 3640, for N = 128; the number of independent k (i.e., k but not −k) is n(k 2 ).
Average magnetic and kinetic energy spectra, ĒM (k 2 ) and ĒK (k 2 ), for Runs 1 and 2b.These are averages over modes k with the same value of k 2 = |k| 2 ; there are 3036 different values of k 2 , where 1 ≤ k 2 ≤ 3640, for N = 128.The spectra are shown for times t = 0, 0.01t end , 0.1t end 0.5t end and t end , where t end = 1095 for Run 1 and 1037 for Run2b.Spectra are approaching their expectation values and are near equilibrium; they are in true equilibrium when they fully match the ideal expectation values.
The definitions of the averaged spectra ĒM (k 2 ) and ĒK (k Here, n(k 2 ) is the number of independent k that satisfy |k| 2 = k 2 .The number n(k 2 ) jumps around as k 2 increases; for example, (We have n(k 2 ) = 0 whenever k 2 = 4 a (8m + 7), a, m = 0, 1, 2, . . .; see [43]).The full energy spectra are n(k 2 ) ĒM,K (k 2 ) at each value of k 2 and thus jumps wildly as k 2 increases, which is why we prefer to look at the average energy spectra (35) and (36), as in Figure 1.(However, its running average over near neighbors n(k 2 ) is well approximated by n(k 2 ) = πk.) If we use the results in Appendix C or in [24], we find that the ideal expectation values of ĒM (k 2 ) and ĒK (k Here, 2 = α 2 − β 2 /4 and α, β and γ are the normalized inverse temperatures related to inverse temperatures appearing in the phase space probability density (A20); please see Appendix B for details.For each of the cases, α, β and γ are determined by numerically finding the minimum of the entropy functional (A38) with the proviso that for Case II (Runs 2a and 2b), β = 0; for Case III (Run 3), γ = 0; and for Case V (Run 5), β = γ = 0.The values of α, β and γ for the different runs are given in Table 3.The spectra are shown in Figure 1 for times t = 0, 0.01t end , 0.1t end , 0.5t end and t end , where t end = 1095 for Run 1 and 1037 for Run2b.The spectra are approaching their expectation values and are near equilibrium; when they fully match the ideal expectation values, then they are in true equilibrium.The spectra ĒM (k 2 ) and ĒK (k 2 ) for Runs 3 and 5 will both be similar to those in Figure 1d, i.e., becoming flat, while the spectra for Run 4 will be similar to Figure 1a,b, though not as highly peaked at k 2 = 1.

Computational Results
In Figures 2 and 3, we show how the six runs in Table 1 change over time with respect to volume-averaged energy E, magnetic energy E M , kinetic helicity H K , mean squared magnetic vector potential A, cross helicity H C , magnetic helicity H M , enstrophy (mean squared vorticity) Ω and mean squared electric current J; for Run 4, parallel helicity H P is also shown; again, please see Appendix A for precise definitions of these quantities.
The time-averages and standard deviations for these quantities over the course of the runs are given in Table 2.For each run in Table 2, the quantities that are supposed to be conserved are, in fact, conserved as shown in Figures 2 and 3.The most conserved quantity, by far, is the magnetic helicity H M in Runs 1, 2a and 2b.Those quantities in Table 2 that have standard deviations larger in magnitude than their average values are basically just fluctuating around zero; in particular, this seems to be true for the kinetic helicity H K in every run except for Run 4. Note also that statistics for the enstrophy Ω and mean squared current J are essentially equal in each run, indicating that smaller length-scales have equipartition in magnetic and kinetic energy.
For all the runs in Table 1, the values of all the ũ(k, t) and b(k, t) with k 2 ≤ 3 were saved every 0.1 units of simulation time t (i.e., every 200 ∆ts).From these, we can calculate a time history of modal energies In Figures 4 and 5, we see how the E M ( k, t), k = x, ŷ, ẑ, for each run varies with time.In Figure 4, E M ( k, t) and E (1) M (t) are compared to H M (t) in each of Runs 1, 2a, 2b and 4; the result is compelling: M (t) is tending to become equivalent to H M (t), even for Run 4 where H M (t) is not an integral invariant of the dynamical equations (instead, we use its value at t = t end for Run 4).This is again verification of the theoretical result E (1) given by [29,31], which also seems to apply to Run 4 because H M has reached a steady constant value.Figure 5 shows the evolution of E M ( k, t) for Runs 3 and 5; in these runs, H M (t)∼0, as Table 2 shows.Furthermore, in Figure 5a (Case III) and Figure 5b (CaseV), we see E M ( k, t) ≈ 0 for k aligned with B o ; additionally, in Figure 5b, E M ( k, t) is largest for k aligned with Ω o , a 'dipole alignment' that has been seen before [29].The ensemble prediction for both (a) and (b) in Figure 5 is E (1) M = 3/(4M) 2 × 10 −7 , which is to be equally divided among the three E M ( k, t); thus, the steady values during transition here are several orders of magnitude too large, perhaps being supported by an increased magnetic energy inverse cascade during transition (and they do not correlate with |H M |).
Continuing these runs well into equilibrium (to be done) will determine whether these anomalous values persist.In all six runs, a running time-average of the vector coefficients ũ(k, t) and b(k, t) was kept: ũavg (k, t) and bavg (k, t).The expectation has been that these will tend to zero; however, this is not the case for k = 1.These time-averaged vector coefficients can be used to calculate the 'coherent total energy' E C , 'coherent magnetic energy' E C M and 'coherent kinetic energy' E C K .The coherent kinetic and magnetic energies for the six runs are shown in Figure 6; the coherent kinetic energy E C K levels off to a constant nonzero only for Run 1, as seen in Figure 6a, while the coherent magnetic energy E C M levels off to a constant nonzero for Runs 1, 2a, 2b and 4, as seen in Figure 6b and these values of E C M appear to approach their respective values of |H M |, which are shown in Figure 6c.With regard to Runs 3 and 5, Figure 6a,b indicates what seems to be exponential fall-off, while Figure 6c shows that H M for Runs 3 and 5 is fluctuating about zero; thus, the anomalously high values of E (1) M for Runs 3 and 5 seen in Figure 5 may only be transitory.
After numerical integration, components of the vectors ũ(k, t) and b(k, t) can be transformed into helical components ũ+ (k, t), ũ− (k, t), b+ (k, t) and b− (k, t), as discussed in Section 2.2.These are useful but can be further transformed into cyclic linear modes whose noncyclic factors are Ṽ+ 1 (k, t), Ṽ+ 2 (k, t), Ṽ− 1 (k, t) and Ṽ− 2 (k, t), as defined in Section 2.4.If there were no nonlinear interactions in the MHD equations, the Ṽ± 1,2 (k, t) would be complex constants; however, the MHD equations are nonlinear, so the Ṽ± 1,2 (k, t) may wander around with time.Since we have recorded data that can give us the time-histories of these for k 2 ≤ 3, we can plot their trajectories on a complex plane and clearly see the nature of ideal MHD turbulence, at least at the larger length-scales (i.e., smaller wave numbers k).These 'phase portraits' are projections of the dynamical trajectory in the high-dimensional phase space onto a 2-D plane which enables us to visualize the concepts of coherent structure, broken ergodicity and broken symmetry (please see Appendix E for a review of these concepts).For the six runs discussed here, some phase portraits are shown in Figures 7-13   1 remain constant in their associated runs.  1 remain constant in their associated runs.
helical ± components of vorticity ω± (k, t) = ±kũ ± (k, t) and current j± (k, t) = ±k b± (k nnected to velocity and magnetic field ± helical components.The helical representation theoretical analysis and graphical presentation of numerical data, although numerical sim done using the cartesian components (9) and (10).
ynamical System urier-transformed 3-D MHD equations are found by placing expansions for ω(x, t) and b( (4) into ( 1) and ( 2).The result is a set of coupled, nonlinear ordinary differential equatio a dynamical system in the sense of Birkhoff [1927]: Phase portraits from Run 5 are presented in Figure 7. Referring to Section 2.4, we see that (with Ω 1 ( ŷ) = 0.5 and Ω 2 ( ŷ) = −0.5) The transformed variables Ṽ± 1,2 ( ŷ, t) have no cyclic time variation and their evolution from initial values (which are all close to zero) shows the effect of the nonlinear terms in (23) and (24).In Figure 7, the Ṽ± 1,2 ( ŷ, t) appear to be zero-mean random variables whose fluctuations match statistical expectations as to mean and standard deviation; in fact, all Ṽ± 1,2 (k, t) for Case V of Table 1 are expected to be zero-mean random variables with the same standard deviation.(Please see Appendix B of a review of the statistical theory.)However, if we take a look at Figure 8 corresponding to Run 5, we see the unexpected nonlinear behavior of (a) Ṽ± 2 (x, t) and (c) Ṽ± 2 (ẑ, t); we compare this with the trajectory of Ṽ± 1,2 ( ŷ, t) that is shown in Figure 8b (which is the same as Figure 7c).In other words, what we see in Figure 8a,c are not the expected zero-mean random variables of small fluctuation, but very large nonzero mean random variables of small fluctuation; we see 'broken ergodicity' (as discussed in Appendix E).The quasi-steady end stage of the trajectories of Figure 8a,c represent a quasi-stationary coherent structure in the magnetofluid.This is, at first, surprising, but nonetheless, it occurs: ergodicity is dynamically broken for Case V, at least during transition to equilibrium; however, longer run-times (which are planned) are needed to examine behavior in equilibrium.At the largest length-scales, where k 2 = 1, these coherent structures seem very robust during transition in Run 5 and in other runs where they appear, while at smaller length-scales, k 2 > 1, 'random walk' trajectories appear to be converging towards zero-mean, both for Run 5 and for the other Runs in Table 1; again, longer run-times are planned to study equilibrium behavior.Phase portraits from Run 4 are shown in Figure 9. Parts (a) and (b) show phase trajectories of ũ± (ẑ, t) and b± (ẑ, t) before cyclic behavior is removed.In Figure 9c, an example is given of Ṽ± 1 (ẑ, t) trajectories where cyclicity is removed.In Figure 9d, quasi- stationarity indicates coherent structure.Again, statistical expectations are that all Case IV variables Ṽ± 1,2 ( k, t) have zero-mean; yet, they do not behave as expected and indicate that broken ergodicity has occurred.In fact, the behavior seen in Figure 9d is the same as that expected of Case II in Table 1, i.e., the three k 2 = 1 magnetic modes altogether have energy |H M | and thus, average magnitudes each of |H M |/3 = 0.0269.(In Figure 9d, H avg M from Table 2 is used for H M .)In Figure 9, these Run 4 Ṽ± 1,2 ( k, t) are, explicitly (with In Figure 9d, we also see 'broken symmetry' in that the magnitudes of the quasistationary values of Ṽ+ 2 (x, t), Ṽ+ 2 ( ŷ, t) and Ṽ+ 2 (ẑ, t) are unequal.Phase portraits from Run 3 are shown in Figure 10.The yellow circles around the origin show the expected average magnitudes; in (a) and (b), these circles have also been placed around the ends of the trajectories, indicating that the fluctuation levels are as predicted, but the values are far from the predicted zero-mean, compared with the expected behavior seen in Figure 10c.In Figure 10 The statistical ensemble expectation is that the Ṽ± 1,2 (k, t) for all k, 0 < k ≤ K, in addition to being zero-mean, will have the same standard deviation, (4M) −1/2 = 2.6 × 10 −4 .In Figure 10a,b, we see that this expectation is not met and that, again, we have broken ergodicity manifested during transition as large-scale coherent structure.Additionally, the magnitudes at the ends of the trajectories in Figure 10a,b are much larger than predicted; whether or not these persist will be determined by allowing Run 3 to continue running, as is planned for the next phase of our work.
Phase portraits from Runs (a) 1, (b) 2a and (c) 2b are presented in Figure 11.Light parts are the full trajectory, while dark ends are from the last half of the trajectory, i.e., from 0.5t end to t end .The black circles have radius √ H M /3, where H M = H avg M , is given for Runs 1, 2a and 2b in Table 2.In Figure 11, all Ṽ+ 2 ( k, t) = b+ ( k, t).These phase trajectories all clearly show broken ergodicity and broken symmetry, for although their rms magnitudes are predicted, their quasi-stationarity is not and is, instead, a dynamical phenomenon.Now, we turn to Figure 12 for close-ups of ends of Run 1 trajectories in Figure 11a; these also serve as examples of the level of quasi-stationarity in Runs 2a and 2b.
In Figure 12, we see close-up phase portraits of Ṽ+ 2 ( k, t), k = (a) x, (b) ŷ and (c) ẑ, from the ends of the Run 1 phase trajectories shown in Figure 11a.Light gray is part of the full trajectory, while dark gray is the last half of the trajectory, i.e., from 0.5t end to t end , and the black part is from 0.9t end to t end .Again, in these figures, Ṽ+ 2 ( k, t) = b+ ( k, t).The small extent of the parts from 0.9t end to t end indicate the high level of quasi-stationarity in the coherent structure that has been created as MHD turbulence transitioned from an initially small-scale distribution of energy at t = 0 to the states appearing in Figure 11; in this figure, Runs 2a and 2b are particularly relevant to understanding the origin of the MHD dynamo in rotating planets and stars: it is an inherent property of MHD turbulence.Again, Figure 12 presents a close-up example of the energetic, quasi-stationary, coherent structures that arise in MHD turbulence.
Figure 10.Phase portraits from Run 3; + marks the origin; the variables are normalized by N 3/2 , N = 128.The yellow circles around the origin are the expected average magnitudes; in (a-c), these circles have also been placed around the ends of the trajectories, indicating that the fluctuation levels are as predicted, but the mean values are far from zero.Here, Ṽ±   Close-up phase portraits of Ṽ+ 2 ( k, t), k = (a) x, (b) ŷ and (c) ẑ, from Run 1; these are the ends of the trajectories shown in Figure 11.(The variables are normalized by N 3/2 , N = 128.)Light gray is part of the full trajectory, while dark gray is the last half of the trajectory, i.e., from 0.5t end to t end , and the black part is from 0.9t end to t end .In these figures, Ṽ+ 2 ( k, t) = b+ ( k, t).
The energetic coherent structure generated by MHD turbulence is essentially magnetic.However, some coherent structure is exhibited by the velocity field, specifically in Runs 1 and 4. In Figure 13, you will find phase portraits of Ṽ+ 1 ( k, t), where k = x, ŷ and ẑ, for (a) Run 1 and (b) Run 4. These runs are the only ones with Ṽ+ 1 (k, t) variables (these are inertial waves) exhibiting coherent structure.These are all purely helical velocity variables Ṽ+ 1 ( k, t) = ũ+ ( k, t), except in (b), where the one Run 4 variable exhibiting coherence is mostly kinetic but with a small part that is magnetic: Once again, details of the variables Ṽ± 1,2 (k, t) and their relation to the primary variables ũ± (k, t) and b± (k, t) can be found in Section 2.4.thods used here, numerical simulation and statistical analysis, may also profitably

Discussion
The primary discovery here is that all of the five Cases of ideal MHD turbulence in Table 1 make a transition from small-scale initial conditions into a near-equilibrium state that contains coherent structure.Since the ensemble prediction is that all variables have zero-mean, what we have is broken ergodicity (and broken symmetry; please see Appendix E).The most surprising phenomena are that Run 3 (Case III) and Run 5 (Case V) developed coherent structure during transition.This was not expected, although coherent structure was expected in Run 1 (Case I), as well as Runs 2a and 2b (both Case II), since similar coherent structure has been seen in numerical simulations of both ideal [24,44] and real [12,13] MHD turbulence.
Coherent structure is identified by quasi-stationary Ṽ± 1,2 ( k, t), as Figures 2-13 show.This discovery was facilitated through Fourier spectral method numerical simulations on a 128 3 grid.It has been shown with similar 64 3 simulations that coherent structure also appears in dissipative MHD turbulence forced at an intermediate wavenumber and that these dissipative, forced results obey ideal MHD statistical theory at the largest wave lengths [12,13].Furthermore, the statistical theory of ideal MHD turbulence developed for a periodic box has the same form as that developed for a spherical shell [30].Therein lies the importance of ideal results to the real MHD turbulence contained in planetary liquid cores and in the convective layers of stars: it is a basic, qualitative explanation of how astrophysical dipole magnetic fields arise and thus, a solution to 'the dynamo problem'.
MHD turbulence in astrophysical objects is covered by Cases I and II in Table 1 and as has been shown here with 128 3 simulations, these are the cases with the most energetic coherent structure, i.e., quasi-stationary 'dipole magnetic fields.'The energy in the large-scale magnetic field is directly proportional to the total magnetic helicity in the turbulent magnetofluid; in the ideal case, this is set by initial conditions, while in the real case, the dipole field grows through inverse cascade and dynamic alignment as long as whatever forcing mechanism is present injects magnetic helicity of primarily one sign into the magnetofluid [12,13].The exact nature of the helical forcing is hidden from view deep within an astrophysical object, but our results imply that it must be there.
Cases III, IV and V of Table 1 all contain an external magnetic field and are thus related to terrestrial experiments with turbulent magnetofluids, e.g., tokamaks, or perhaps to concentrated regions of plasma within stellar winds, for example, rather than whole planets and stars.The statistical mechanics of ideal MHD turbulence (see Appendix B) tells us that the modal energies of Cases III and V are expected to be equipartitioned, while Case IV is expected to have much more energy in the k = 1 modes.Thus, the coherent structure seen in Runs 3 and 5 is a surprise, while in Run 4, the surprise is that it appears to be like Case I or II once H M reaches an essentially constant value.

Conclusions
In summary, the primary discovery here is that all of the five Cases of ideal MHD turbulence in Table 1 have developed quasi-stationary coherent structure.These coherent structures appear robust as the runs appear to have finished or are close to finishing transition from initial conditions to an equilibrium state.However, the appearance of coherent structure for Cases III and V during transition to equilibrium may be transitory and further long-term computation is needed to examine these cases once they fully reach equilibrium.The importance of the computational and theoretical results, particularly for Cases I and II which are most relevant for planets and stars, is that they show that MHD turbulence, per se, is a dynamo.We intend to continue these runs beyond the values of t end given in Table 2 in order to examine MHD turbulence while in full, long-term equilibrium; we will report on the results when we have them.
Another direction future investigations can take is to add dissipation and forcing to our 128 3 simulations as was done on 64 3 grids [12,13].The methods employed to force a model magnetofluid at intermediate wavenumbers are speculative because physical forcing mechanisms are hidden from observation within a planet or star.Nevertheless, having used a variety of mid-wavenumber, helical forcing methods on 64 3 grids, we have always seen behavior analogous to ideal MHD turbulence at the smallest wavenumber, which is the emergence of large-scale coherent magnetic structure whose energy is determined by the amount of magnetic helicity in the magnetofluid turbulence.We would expect to see this in 128 3 simulations, along with perhaps some novel effects, but this remains to be done.
The methods used here, numerical simulation and statistical analysis, may also profitably be used in examining other dynamical systems that have associated integral invariants.One example is a model of quantum electrodynamic turbulence [45] in which the coupled Dirac and Maxwell equations were numerically solved in a manner similar to that used here for a magnetofluid.This work was based on classical field theory but allowed for a nonperturbative approach.Another possible application is to study geometrodynamic turbulence, i.e., the evolution of space-time from arbitrary initial conditions; a preliminary numerical simulation could be based on the Einstein equations in a synchronous reference frame [46] before moving on to more general reference frames.These seem to be 'grand challenge problems' but should yield interesting results in terms of mathematics with perhaps some compelling physical relevance.Since all functions are periodic in x-space, we can use Equations ( 1) and ( 2), along with integration by parts to derive the following relations [22]: When ν = η = 0 and Ω o = 0, the quantities E, H C and H M are the traditional ideal integral invariants of MHD turbulence [1,18,47].If Ω o = 0, then (A5) indicates that H C is no longer an ideal invariant.If an external mean magnetic field B o is imposed, then H M would also no longer be an ideal invariant [22].We note that 'generalized helicities' G C and G M related to H C and H M can be defined and are ideal invariants even in the presence of nonzero B o and Ω o , but these are not useful for canonical ensemble theory [48].Thus, we will only need to refer to H M , H C and another possible invariant H P in this paper.
The helicity H P arises when Ω o = σB o ; if (A5) is added to −σ times (A6), we obtain H P has been called the 'parallel helicity' [22] and is an invariant when Ω o = σB o and ν = η = 0. [49] calls H P the 'hybrid' helicity and tries to apply this case to the geodynamo, where he identifies B o with the Earth's dipole field.This does not seem appropriate as the geomagnetic dipole field is dynamic and not external; application of his results seems more appropriate to a tokamak [50].
Although kinetic helicity H K is not an ideal MHD invariant, it is an ideal invariant of fluid turbulence [16].To show this, we use (1) to find that, for a periodic box, Thus, if b(x, t) ≡ 0 and ν = 0, then H K and E are the ideal invariants of hydrodynamic turbulence.
Now, b = ∇ × a, with ∇ • a = 0, so b(k, t) = ik × ã(k, t).Using ( 11), ( 16) and ( 18), we see that The total energy E, magnetic energy E M , kinetic energy E K , kinetic helicity H K , meansquared vector potential A, cross helicity H C , magnetic helicity H M , parallel helicity H P , enstrophy Ω and mean-squared current J are then represented in k-space by the quadratic sums: ) when apparent non-ergodicity was first noticed and reported, and confirmed later [21].
As already mentioned, this non-ergodicity is actually 'broken ergodicity' [25]; a review of broken ergodicity for ideal MHD turbulence is given by [24].In general, there is no reason to expect ergodicity in any dynamical system, as this can only be determined by experimentation or numerical simulation.This is because the statistical theory averages over all probable realizations while a single experiment or simulation only produces one dynamical realization.Remember that ergodicity is defined as the equivalence of statistical ensemble predictions with statistics drawn from a single dynamical realization; sometimes, this definition is ignored or unknown and erroneous conclusions result [54].Furthermore, one must use large enough grid sizes in simulations (see [23,24]) since turbulence is high-dimensional; otherwise, nonergodic behavior will be missed, as in the work of [55].
In the various cases of ideal MHD turbulence, expectation values of the global quantities in (A2) and (A3), as well as any other phase functions, can be determined with respect to the probability density function (A20).Expectation values are integrations over the phase space Γ, i.e., for all possible values of the coordinates in Γ, which are the real and imaginary parts of ũ± (k) and b± (k) (treated as phase coordinates and not dynamical variables).Given a quantity Q, the expectation value Q is defined by (A21) As an example, the velocity and magnetic field coefficients are expected to have zero mean values: ũ In the ideal case, the integral invariants E, H M , and possibly H C , should have timeindependent values E , H M , and H C that are equal to their expectation values: In fact, we require that (A23) be true and use this to determine the 'inverse temperatures' α, β and γ in (A20).Whereas (A22) is an 'ergodic hypothesis', (A23) is actually an a priori axiom on which the theory of ideal MHD turbulence is based, though justified by a posteriori numerical results.Ergodicity, or rather the lack of it, will be discussed in Appendix E.
Although the M k in (A26) can also be expressed as 8 × 8 real symmetric matrices and the ỹ(k) as 8 × 1 real arrays [17], finding eigenvalues and eigenvariables is facilitated by using the 4 × 4 matrices M k and 4 × 1 complex arrays ỹ, along with the properties of block matrices given by [56].
The real, symmetric matrices M k can be diagonalized (and more easily than the Hermitian matrices used previously [23,24,57]) to yield the modal PDFs dipole fields) that spontaneously arise within a turbulent magnetofluid such as is found in the Earth's outer core.
In Cases III and V in Table A1, ϕ o ≡ E/2, and modal spectra are expected to become equipartitioned in energy [24].Case IV in Table A1 can be evaluated by setting γ = −σβ where σ = Ω o /B o ; one can then use analysis similar to that discussed in [24] to show that, instead of the simpler looking set in Equation (A37), we obtain The exact value of ϕ = ϕ o must be determined by numerically finding the minimum of (A38).However, in (A49), a choice of ± must be made.The correct choice is + if H P > 0 and − if H P < 0; the choice − if H P > 0 and + if H P < 0 leads to a computed ϕ o outside the range (A50); and for H P ≡ 0, we have ϕ = E /2 and β = 0.The two choices for H P = 0 are due to the need for parity invariance in the statistical theory of ideal MHD turbulence: (A20) becomes D = Z −1 exp(−αE − βH P ) for Case IV and the product βH P must be invariant under a parity transformation (β and H P are both pseudoscalars, while βH P ≤ 0 is a scalar).λ (4) The eigenvariables have real (R) and imaginary (I) parts, i.e., ṽn (k) = ṽR n (k) + i ṽI n (k), n = 1, 2, 3, 4, each have the same eigenvalue λ As defined in (A51), the index n = 1 refers to negative and n = 2 to positive kinetic helicity coefficients; similarly, the index n = 3 refers to negative and n = 4 to positive magnetic helicity coefficients.The relations (A54) and (A55) tell us that the expected energies with respect to helicity are The sum of these over independent modes k is E plus a term of O(M −1 ), as it should be.Again, for cubical periodic boxes or symmetrical spherical shells, the three lowestwavenumber modes are expected to have the same energy.However, for the nonrotating case, and especially for the rotating case, there is always some dynamical symmetry breaking so that one of the lowest-wavenumber modes dominates dynamically, as will be discussed further shortly.
The statistical results given above are directly related to Case II of Table A1, but also apply approximately to Case I if H C is small compared to H M .In Case III, H M is not an ideal invariant but H C is and the prediction is an equipartitioned distribution of energy amongst all the kinetic variables and a different equipartition amongst the magnetic variables, while in Case V, all the variables are predicted to have the same energy.The statistical predictions for Case IV are discussed by [24] but not completely worked out; however, the spectrum is expected to peak strongly at k = 1.(Please note the numbering of Cases in [24] is different than Table A1 in that II and III are switched).or ẑ, does not become as large as predicted by (A60).These predictions are just average values over the ensemble and to see what is really happening, we must consider the sum of the expectation values b+ ( k, t).The smallest wavenumber k min = | k| = 1 occurs for the wavevectors k = x, ŷ or ẑ.The ensemble prediction (A61) tells us that the three complex vector coefficients b+ ( k, t) are very large and unique from those for k > 1; using these, we can define a six component vector in a 6-D real space or a three-component vector in a complex 3-D space; for compactness, we define a vector ṽd and dot product |ṽ d | 2 = ṽ † d ṽd in a 3-D complex space: The endpoint of ṽd is, for the reasons given in Appendix E.1, a quasi-stationary point on the surface of the hypersphere of radius |H M | in a 6-D subspace of the 8M-D phase space Γ.The coordinates of the 6-D space can undergo an arbitrary orthogonal rotation (or a unitary transformation in the complex 3-D space) that puts the point of the vector on any desired point on the surface.Thus, although (A62a) predicts that all b+ ( k, t) will have the same magnitude, this does not take into account that a suitable orthogonal transformation of initial conditions in the whole phase space will lead to the evolution of ṽd pointing in any direction its 6-D space that we choose, at least for the nonrotating case.Canonical ensemble predictions can give us the mean-squared expectation values (A61) but cannot predict the direction of ṽd ; this directional symmetry is broken dynamically.
The vector ṽd will exhibit broken ergodicity, i.e., its components are not zero-mean random variables; it will also exhibit broken symmetry because the quasi-stationary direction of ṽd is not predictable; instead, initial conditions impose a particular direction on it.The appearance of broken ergodicity has been noted many times before [20,21,24,44], as has the phenomenon of broken symmetry [44,58].Here, we see how these aspects of MHD turbulence are connected.
In equilibrium, the magnitudes of the b+ ( k, t), for k = x, ŷ and ẑ are often not that different from each other in nonrotating ideal MHD Case I but may vary appreciably when rotation is imposed, i.e., Case II, where the eigenfunction b+ ( k, t) with k parallel to Ω o has essentially all the dipole energy, as we have consistently seen numerically.In the real case of forced, dissipative MHD turbulence, this phenomenon is also usually observed numerically, although some forcing mechanisms may disrupt this process, as seen in Figure 7 of [13].
The generation of quasi-stationary, energetic dipole magnetic fields, along with dipole moment alignment with the rotation axis seems fairly ubiquitous in numerical simulations, as well as in planets and stars.Perhaps, the theory of ideal MHD turbulence and its relevance to real MHD turbulence at the largest scales will be accepted as a viable explanation of these planetary and stellar phenomena.Verbum sat sapienti est.

Figure 2 .
Figure 2. Global quantities (see Appendix A) versus time t for 128 3 homogeneous MHD Runs (a) 1, (b) 2a and (c) 2b.Notice how the integral invariants in Table1remain constant in their associated runs.

Figure 3 .
Figure 3. Global quantities (see Appendix A) versus time t for 128 3 homogeneous MHD Runs (a) 3, (b) 4 and (c) 5. Notice how the integral invariants in Table1remain constant in their associated runs.

Figure 4 .
Figure 4. Magnetic energies for the k = 1 modes, normalized by H M > 0, versus time t for Runs (a) 1, (b) 2a, (c) 2b and (d) 4. Notice how the total magnetic energy for k = 1, E

Figure 5 .Figure 6 .
Figure 5. Magnetic energies for the k = 1 modes versus time t for Runs (a) 3 and (b) 5.In these runs, B o = 0 so that H M is not an invariant, although H P is an invariant in Run 4 where Ω o = σB o = 0.In Run 5, B o is not aligned with Ω o , and the only invariant is E. The ensemble prediction for both (a,b) is E (1) M

Figure 9 .
Figure 9. Phase portraits from Run 4; the white + marks the origin; the variables are normalized by N 3/2 , N = 128.Parts (a,b) show phase trajectories of ũ± (ẑ, t) and b± (ẑ, t) before cyclic behaviors is removed.In (c) is an example of Ṽ± 1 ( k, t) trajectories when cyclicity is removed.In (d), stationarity indicates coherent structure.The yellow circle in (c) indicates predicted average magnitude of each

Figure 11 .
Figure 11.Phase portraits from Runs (a) 1, (b) 2a and (c) 2b; + marks the origin; the variables are normalized by N 3/2 , N = 128.Light parts are the full trajectory, while dark ends are from the last half of the trajectory, i.e., from 0.5t end to t end .The black circles have radius √ H M /3.Please see Figure 12 for close-ups of ends of trajectories in (a).In these figures, Ṽ+ 2 ( k, t) = b+ ( k, t).
Figure 12.Close-up phase portraits of Ṽ+ 2 ( k, t), k = (a) x, (b) ŷ and (c) ẑ, from Run 1; these are the ends of the trajectories shown in Figure 11.(The variables are normalized by N 3/2 , N = 128.)Light gray is part of the full trajectory, while dark gray is the last half of the trajectory, i.e., from 0.5t end to t end , and the black part is from 0.9t end to t end .In these figures, Ṽ+ 2 ( k, t) = b+ ( k, t).

Table 1 .
Cases and associated Runs with different integral invariants for ideal MHD turbulence.The 'parallel helicity' of Case IV is H P = H C − σH M and occurs when Ω o = σB o , i.e., Ω o is parallel to B o .

Table 2 .
Time averages and standard deviations (avg ± std) for various global quantities over the last half of each run are given below for six new ideal MHD turbulence long-time 128 3 runs A, B, C, D, E and F. These global quantities are: energy E, kinetic energy E K , magnetic energy E M , mean squared vector potential A, kinetic helicity H K , magnetic helicity H M , cross helicity H C , parallel helicity H P , enstrophy Ω and mean squared current J.

Table 3 .
Values of the inverse temperatures α, β and γ for the Runs in Table1with the proviso that for Runs 2a and 2b, β ≡ 0; for Run 3, γ ≡ 0; and for Run 5, β ≡ 0 and γ ≡ 0. When needed, E , H C , and H M took their values from E avg , H

Table 2 .
(38) need for precision here is due to the possible smallness of the denominators in(38)and (39))