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

Abstract

We continue our study of the transition of ideal, homogeneous, incompressible, magnetohydrodynamic (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 field and represent the five cases of ideal, homogeneous, incompressible, and MHD turbulence: Case I (Run 1), with no rotation or mean field; Case II (Runs 2a and 2b), where only rotation is imposed; Case III (Run 3), which has only a mean magnetic field; Case IV (Run 4), where rotation vector and mean magnetic field direction are aligned; and Case V (Run 5), which has non-aligned rotation vector and mean field directions. Statistical mechanics predicts that dynamic Fourier coefficients are zero-mean random variables, but largest-scale coherent magnetic structures emerge and manifest themselves as Fourier coefficients 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’.

1. Introduction

The ‘dynamo problem’ is the problem of understanding how magnetofluids within planets and stars produce a quasi-stationary, energetic dipole magnetic field. Over a century ago, Larmor [1] hypothesized that the solution lay in a ‘self-excited dynamo’ within the Earth or the Sun. The conjectured dynamo was not mechanical, like Faraday’s dynamo, but was thought to arise from ‘convective circulation’ and ‘electric currents.’ Four decades later, Elsässer [2] saw that for ‘the dynamo problem, that is … the problem of generating and maintaining magnetic fields which draw their energy from the mechanical energy of the fluid, the nonlinear character of the equations is altogether essential’, as it produces ‘turbulence, the most conspicuous of the nonlinear phenomena of fluid dynamics.’ Thus, he realized that statistically stationary hydromagnetic turbulence was fundamental to the solution of the dynamo problem, rather than ‘rigorously stationary flow’ (i.e., a kinematic dynamo). He added that there were ‘qualitative conditions, three in number, requisite for the operation of … dynamo models’; these conditions were (1) large linear dimensions; (2) rotation; and (3) convection. The first condition implies that Reynolds numbers are sufficiently large; the second that rotation axis and magnetic dipole vector appear on average to be aligned in planets and stars; and the third, that convection is the dynamo’s source of energy. MHD turbulence, due to conditions (1) and (3), is expected to occur in planetary liquid cores and stellar interiors when global magnetism is observed, and must be an integral part of any dynamo model. Condition (2), rotation, appears to cause a close alignment between the dipole moment vector and the rotation axis; however, despite being prevalent in planets and stars, is not essential for dynamo action, whereas MHD turbulence is. Our conclusion, as discussed herein, is that MHD turbulence is the dynamo.

Earth is the primary example we have of a planet with a quasi-steady, mostly dipole magnetic field; current knowledge about the external geomagnetic field is contained in the International Geomagnetic Reference Field [3]. Reynolds numbers are predicted to be very large within the Earth’s outer core, so that convective forcing is able to produce MHD turbulence, which must be involved in creating and maintaining the geodynamo [4]. The realistic numerical simulations of the geodynamo [5,6,7], by including dissipation and forcing, have established that MHD processes within the Earth are capable of creating a magnetic field similar to the actual geomagnetic field, including the reversals of the dominant dipole component. There have also been laboratory experiments involving magnetofluids [8], some of which have shown the growth of a self-generated magnetic field, i.e., a dynamo effect [9,10,11]. Thus, the ‘dynamo problem’ falls within the purview of MHD theory, experiment and computation. However, in spite of successful numerical simulations and laboratory experiments, the fundamental origin of the geomagnetic dipole field still appeared to be a theoretical mystery [12]. Here, we present results aimed at resolving this mystery. Our focus will be on ideal MHD theory and computation as this allows us to examine large-scale dynamics without the complications of dissipation and forcing, or the associated artefacts that arise in sub-grid scale modeling.

In that regard, our new results come from extending into equilibrium the numerical simulations reported in [13], where relevant background and references are given. In brief, we use a Fourier spectral transform method as a surrogate for a planetary interior that contains a turbulent magnetofluid, such as the Earth’s outer core. The connection between the periodic box of a Fourier model and the spherical shell of an Earth-like model is that they have essentially the same ideal MHD statistical mechanics [14]. However, no spectral transform methods for studying MHD turbulence in a spherical shell exist yet, and although non-transform methods [15,16,17] have been used, they will always have insufficient resolution. Thus, we use Fourier spectral transform methods in a periodic box because they are efficient and, again, serve as a surrogate for transform methods more explicitly designed for spherical shells, until such methods are available.

In our previous computer runs [13], the transition of ideal MHD turbulence from non-equilibrium initial conditions to near equilibrium was studied and the results showed that a quasi-steady, energetic, coherent structure (equivalent to a dipole magnetic field) arose in all the cases defined in

Table 1. Cases and associated runs with different integral invariants for ideal MHD turbulence. When ${\Omega }_{\mathrm{o}}=\sigma {\mathbf{B}}_{\mathrm{o}}$, the ‘parallel helicity’ of Case IV is $H_P=H_C-\sigma H_M$.

Case	Mean Field	Rotation	Invariants	${128}^3$ Runs
I	${\mathbf{B}}_{\mathrm{o}}=0$	${\Omega }_{\mathrm{o}}=0$	E, $H_C$, $H_M$	1
II	${\mathbf{B}}_{\mathrm{o}}=0$	${\Omega }_{\mathrm{o}}\ne 0$	E, $H_M$	2a,b
III	${\mathbf{B}}_{\mathrm{o}}\ne 0$	${\Omega }_{\mathrm{o}}=0$	E, $H_C$	3
IV	${\mathbf{B}}_{\mathrm{o}}\ne 0$	${\Omega }_{\mathrm{o}}=\sigma {\mathbf{B}}_{\mathrm{o}}$	E, $H_P$	4
V	${\mathbf{B}}_{\mathrm{o}}\ne 0$	${\Omega }_{\mathrm{o}}\times {\mathbf{B}}_{\mathrm{o}}\ne 0$	E	5

. The appearance of coherent structure is an example of ‘broken ergodicity’ [18]. As explained in the Appendix to [13], ergodicity is dynamically broken when some of the largest-scale magnetic Fourier coefficients become very large, as statistically predicted, and then become essentially constant as they are slightly affected by the myriad of very small random perturbations due to all the other Fourier modes.

However, these coherent structures were not statistically expected to arise in Cases III and V of Table 1, which had a mean magnetic field (i.e., constant in space and time) embedded in the magnetofluid, with no rotation in Case III and with an unaligned rotation vector in Case V. In order to see whether these persisted, we extended these runs (and the others) by a sufficient amount so as to fully enter equilibrium. As we will show here, those unexpected coherent structures were, in fact, transitory, while the coherent structures seen, as expected, in those cases without a mean magnetic field, i.e., Cases I (non-rotating) and II (rotating) of Table 1, remained energetic and quasi-steady. In Case IV, which has a mean magnetic field aligned with the rotation axis, there was also a quasi-steady coherent structure, as expected, but with energy which was an order of magnitude less than those in Cases I and II.

Cases I and II have no mean magnetic field and are thus most analogous to a planetary interior. In those cases, the statistical mechanics of ideal MHD turbulence, as reviewed in [13], explains dipole alignment and predicts that ‘dipole energy’ $E_d$ and magnetic helicity $H_M$ are related by

$$E_d=k_{min}\left|H_M\right|.$$

(1)

This was indeed seen to hold in Runs 1, 2a and 2b, and, approximately, in Run 4; it was also seen to apply in forced, dissipative Fourier method simulations of MHD turbulence on a ${64}^3$ grid [19,20]. (regarding (1), in the Fourier case, $k_{min}=1$, while in a spherical shell model of the Earth’s outer core [14], $k_{min}\cong 1.8638$, both in dimensionless units).

Since the necessary details of the model system and the statistical theory appear in [13], we avoid repeating this information as much as possible here, in order to focus on our new computational results and how they verify the statistical theory of ideal MHD turbulence; again, dissipative and forced simulations have shown that this theory is applicable to real MHD turbulence [19,20]. In essence, what is presented here is a solution to the ‘dynamo problem’.

2. Global Quantities

There are various important global quantities that can be expressed as averages over either $\mathbf{x}$ space or, equivalently, $\mathbf{k}$ space. We define the volume average of a quantity $\Phi (\mathbf{x},t)$ multiplied by a quantity $\Psi (\mathbf{x},t)$ as $\left\{\Phi \Psi \right\}$, which is an integral over the periodic box of side length $2\pi$:

$$\begin{array}{ccc}\hfill \left\{\Phi \Psi \right\}& \equiv & {\left(2\pi \right)}^{-3}\int \Phi (\mathbf{x},t)\Psi (\mathbf{x},t)d^{3}x=\frac{1}{N^3}\sum _{\mathbf{k}}{\tilde{\Phi }}^{*}(\mathbf{k},t)\tilde{\Psi }(\mathbf{k},t).\hfill \end{array}$$

(2)

${\tilde{\Phi }}^{*}(\mathbf{k},t)$ and $\tilde{\Psi }(\mathbf{k},t)$ are Fourier transforms of $\Phi (\mathbf{x},t)$ and $\Psi (\mathbf{x},t)$, as discussed in [13].

Using (2), the volume-averaged energy E, enstrophy $\Omega$, mean-squared current J, cross helicity $H_C$, magnetic helicity $H_M$ and mean-squared vector potential A (the last two defined in terms of the magnetic vector potential $\mathbf{a}$, where $\nabla \times \mathbf{a}=\mathbf{b}$, $\nabla ·\mathbf{a}=0$) are

$$\begin{array}{c}\hfill E_K=\frac{1}{2}\left\{u^2\right\},E_M=\frac{1}{2}\left\{b^2\right\},\Omega =\frac{1}{2}\left\{{\omega }^2\right\},J=\frac{1}{2}\left\{j^2\right\},\end{array}$$

(3)

$$\begin{array}{c}\hfill E=E_K+E_M,H_C=\frac{1}{2}\left\{\mathbf{u}·\mathbf{b}\right\},H_M=\frac{1}{2}\left\{\mathbf{a}·\mathbf{b}\right\},A=\frac{1}{2}\left\{a^2\right\}.\end{array}$$

(4)

Since all functions are periodic in $\mathbf{x}$ space, we can use the (incompressible) MHD equations, along with integration by parts to derive the following relations [13]:

$$\begin{array}{ccc}\hfill \frac{dE}{dt}& =& -2\left(\nu \Omega +\eta J\right),\hfill \\ \hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill \frac{d{H}_C}{dt}& =& {\Omega }_{\mathrm{o}}·\left\{\mathbf{b}\times \mathbf{u}\right\}-\frac{1}{2}(\nu +\eta )\left\{\mathbf{j}·\mathbf{\omega }\right\},\hfill \\ \hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill \frac{d{H}_M}{dt}& =& {\mathbf{B}}_{\mathrm{o}}·\left\{\mathbf{b}\times \mathbf{u}\right\}-\eta \left\{\mathbf{j}·\mathbf{b}\right\}.\hfill \end{array}$$

(7)

When $\nu =\eta =0$ and ${\Omega }_{\mathrm{o}}=0$, the quantities E, $H_C$ and $H_M$ are the traditional ideal integral invariants of MHD turbulence [2,21,22]. If ${\Omega }_{\mathrm{o}}\ne 0$, then (6) indicates that $H_C$ is no longer an ideal invariant. If an external mean magnetic field ${\mathbf{B}}_{\mathrm{o}}$ is imposed, then $H_M$ would also no longer be an ideal invariant [23].

The helicity $H_P$ arises when ${\Omega }_{\mathrm{o}}=\sigma {\mathbf{B}}_{\mathrm{o}}$; if (6) is added to $-\sigma$ times (7), we obtain

$$\begin{array}{ccc}\hfill \frac{d{H}_P}{dt}& =& -\frac{1}{2}(\nu +\eta )\left\{\mathbf{j}·\mathbf{\omega }\right\}+\sigma \eta \left\{\mathbf{j}·\mathbf{b}\right\},H_P=H_C-\sigma H_M.\hfill \end{array}$$

(8)

$H_P$ has been called the ‘parallel helicity’ [23] and is an invariant when ${\Omega }_{\mathrm{o}}=\sigma {\mathbf{B}}_{\mathrm{o}}$ and $\nu =\eta =0$. Ref. [24] calls $H_P$ the ‘hybrid’ helicity and tries to apply this case to the geodynamo, where ${\mathbf{B}}_{\mathrm{o}}$ is identified with the Earth’s dipole field. However, the geomagnetic dipole field is dynamic and not external and the application of the results of [24] seems more appropriate to a tokamak [25].

Although the kinetic helicity $H_K$ is not an ideal MHD invariant, it is an ideal invariant of fluid turbulence [26]. As discussed in [13], we find that, for a periodic box,

$$\frac{d{H}_K}{dt}=\left\{\mathbf{\omega }·\mathbf{j}\times \mathbf{b}\right\}-\nu \left\{{\partial }_i\mathbf{u}·{\partial }_i\mathbf{\omega }\right\},H_K=\frac{1}{2}\left\{\mathbf{u}·\mathbf{\omega }\right\}.$$

(9)

Thus, if $\mathbf{b}(\mathbf{x},t)\equiv 0$ and $\nu =0$, then $H_K$ is an ideal invariant, but only for hydrodynamic turbulence and not MHD.

The dynamic variables of interest are the turbulent velocity $\tilde{\mathbf{u}}(\mathbf{k},t)$ and magnetic field $\tilde{\mathbf{b}}(\mathbf{k},t)$, but these can be written as

$$\begin{array}{ccc}\hfill \tilde{\mathbf{u}}(\mathbf{k},t)& =& {\tilde{u}}_{+}(\mathbf{k},t){\widehat{\mathbf{e}}}_{+}\left(\mathbf{k}\right)+{\tilde{u}}_{-}(\mathbf{k},t){\widehat{\mathbf{e}}}_{-}\left(\mathbf{k}\right),\hfill \\ \hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill \tilde{\mathbf{b}}(\mathbf{k},t)& =& {\tilde{b}}_{+}(\mathbf{k},t){\widehat{\mathbf{e}}}_{+}\left(\mathbf{k}\right)+{\tilde{b}}_{-}(\mathbf{k},t){\widehat{\mathbf{e}}}_{-}\left(\mathbf{k}\right).\hfill \end{array}$$

(11)

The ${\tilde{u}}_{+}(\mathbf{k},t)$ and ${\tilde{u}}_{-}(\mathbf{k},t)$ are coefficients with positive and negative kinetic helicity, respectively; similarly, the ${\tilde{b}}_{+}(\mathbf{k},t)$ and ${\tilde{b}}_{-}(\mathbf{k},t)$ are coefficients with positive and negative magnetic helicity, respectively. The ${\widehat{\mathbf{e}}}_{\pm }\left(\mathbf{k}\right)$ are modal unit vectors associated with positive and negative helicity (see [13] for more details).

The magnetic field $\mathbf{b}$ is the curl of the magnetic vector potential $\mathbf{a}$, where $\mathbf{b}=\nabla \times \mathbf{a}$ and $\nabla ·\mathbf{a}=0$; thus, $\tilde{\mathbf{b}}(\mathbf{k},t)=i\mathbf{k}\times \tilde{\mathbf{a}}(\mathbf{k},t)$; again, see [13]. Auxiliary variables, such as $\mathbf{a}$, also have a helical representation in $\mathbf{k}$ space:

$$\tilde{\mathbf{a}}(\mathbf{k},t)=\frac{i}{k^2}\mathbf{k}\times \tilde{\mathbf{b}}(\mathbf{k},t)=\frac{1}{k}\left[{\tilde{b}}_{+}(\mathbf{k},t){\widehat{\mathbf{e}}}_{+}\left(\mathbf{k}\right)-{\tilde{b}}_{-}(\mathbf{k},t){\widehat{\mathbf{e}}}_{-}\left(\mathbf{k}\right)\right].$$

(12)

It proves very useful to use helical coefficients in the definition of various quadratic forms.

Explicitly, in $\mathbf{k}$ space, the total energy E, magnetic energy $E_M$, kinetic energy $E_K$, kinetic helicity $H_K$, mean-squared vector potential A, cross helicity $H_C$, magnetic helicity $H_M$, parallel helicity $H_P$, enstrophy $\Omega$ and mean-squared current J can be represented as

$$\begin{array}{ccc}\hfill E& =& E_M+E_K,\hfill \\ \hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill E_M& =& \frac{1}{2N^3}\sum _{\mathbf{k}}\left[{\left|{\tilde{b}}_{+}(\mathbf{k},t)\right|}^2+{\left|{\tilde{b}}_{-}(\mathbf{k},t)\right|}^2\right],\hfill \\ \hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill E_K& =& \frac{1}{2N^3}\sum _{\mathbf{k}}\left[|{\tilde{u}}_{+}{(\mathbf{k},t)|}^2+{\left|{\tilde{u}}_{-}(\mathbf{k},t)\right|}^2\right],\hfill \\ \hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill H_K& =& \frac{1}{2N^3}\sum _{\mathbf{k}}k\left[|{\tilde{u}}_{+}{(\mathbf{k},t)|}^2-{\left|{\tilde{u}}_{-}(\mathbf{k},t)\right|}^2\right],\hfill \\ \hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill A& =& \frac{1}{2N^3}\sum _{\mathbf{k}}\frac{1}{k^2}\left[{\left|{\tilde{b}}_{+}(\mathbf{k},t)\right|}^2+{\left|{\tilde{b}}_{-}(\mathbf{k},t)\right|}^2\right],\hfill \\ \hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill H_C& =& \frac{1}{2N^3}\sum _{\mathbf{k}}\left[{\tilde{u}}_{+}(\mathbf{k},t){\tilde{b}}_{+}^{*}(\mathbf{k},t)+{\tilde{u}}_{-}(\mathbf{k},t){\tilde{b}}_{-}^{*}(\mathbf{k},t)\right],\hfill \\ \hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill H_M& =& \frac{1}{2N^3}\sum _{\mathbf{k}}\frac{1}{k}\left[{\left|{\tilde{b}}_{+}(\mathbf{k},t)\right|}^2-{\left|{\tilde{b}}_{-}(\mathbf{k},t)\right|}^2\right],\hfill \\ \hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill H_P& =& H_C-\sigma H_M,\hfill \\ \hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill \Omega & =& \frac{1}{2N^3}\sum _{\mathbf{k}}{k}^2\left[{\left|{\tilde{u}}_{+}(\mathbf{k},t)\right|}^2+{\left|{\tilde{u}}_{-}(\mathbf{k},t)\right|}^2\right],\hfill \\ \hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill J& =& \frac{1}{2N^3}\sum _{\mathbf{k}}{k}^2\left[{\left|{\tilde{b}}_{+}(\mathbf{k},t)\right|}^2+{\left|{\tilde{b}}_{-}(\mathbf{k},t)\right|}^2\right].\hfill \end{array}$$

(22)

Again, ${\tilde{u}}_{\pm }(\mathbf{k},t)$ and ${\tilde{b}}_{\pm }(\mathbf{k},t)$ are coefficients with positive (+) and negative (−) helicities. In the ideal MHD case with homogeneous b.c.s, E, $H_M$ (when ${\mathbf{B}}_{\mathrm{o}}=0$), $H_C$ (when ${\Omega }_{\mathrm{o}}=0$) and $H_P$ (when ${\Omega }_{\mathrm{o}}=\sigma {\mathbf{B}}_{\mathrm{o}}$) are ideal invariants, as shown in Table 1. The other quadratic forms involving ${\tilde{u}}_{\pm }(\mathbf{k},t)$ and ${\tilde{b}}_{\pm }(\mathbf{k},t)$ given above will generally be time-dependent, particularly in numerical simulations during the transition from initial conditions to an equilibrium state, at which point they may become quasi-stationary.

3. Numerical Procedure

A Fourier fully de-aliased spectral transform method [27] on an $N^3$ grid with $N=128$ is used; the minimum wave number is $k=\left|\mathbf{k}\right|=1$ and the maximum wave number is $K=\sqrt{3640}\cong 60.34$. Time integration is performed with a third-order Adams–Bashforth–Adams–Moulton method [28] with a time-step of $\Delta t=0.0005$. The cpu-time required for one time-step in any run is $\Delta T\cong 7.34$ s, so the total cpu-time for all the runs listed in

Table 2. Time averages (avg) and standard deviations (std) for global quantities over the last quarter of each run appear below for six new ideal MHD turbulence long-time ${128}^3$ Runs 1, 2a, 2b, 3, 4 and 5. The 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 $\Omega$ and mean-squared current J; these are defined in (2)–(4) and more explicitly in (13)–(22). The angle ${\theta }_D$ at $t=t_{end}$, is defined in (30) and is analogous to the geomagnetic dipole angle (if all dipole moment components were equal, then ${\theta }_D=54.7^{\circ }$).

Run:	1	2a	2b	3	4	5
$t_{end}$:	1475	1526	1426	1132	1276	846
${\Omega }_{\mathrm{o}}$	0	$10\widehat{\mathbf{z}}$	$10\widehat{\mathbf{z}}$	0	$2\widehat{\mathbf{z}}$	$1\widehat{\mathbf{z}}$
${\mathbf{B}}_{\mathrm{o}}$	0	0	0	$1\widehat{\mathbf{z}}$	$1\widehat{\mathbf{z}}$	$\frac{1}{2}\widehat{\mathbf{y}}$
${\theta }_D$	$53.7^{\circ }$	$14.6^{\circ }$	$10.9^{\circ }$	$61.2^{\circ }$	$67.1^{\circ }$	$63.1^{\circ }$
$E^{avg}$	1.01100	1.01246	1.01226	1.02980	1.03506	1.00699
$E^{std}$	1.2609 × 10${}^{-3}$	1.1801 × 10${}^{-3}$	1.2758 × 10${}^{-3}$	2.1086 × 10${}^{-3}$	1.7358 × 10${}^{-3}$	7.7732 × 10${}^{-4}$
$E_K^{avg}$	0.44822	0.47564	0.44539	0.51489	0.51648	0.50349
$E_K^{std}$	6.8596 × 10${}^{-4}$	8.3465 × 10${}^{-4}$	9.4303 × 10${}^{-4}$	1.1379 × 10${}^{-3}$	9.1951 × 10${}^{-4}$	5.0655 × 10${}^{-4}$
$〈E_M〉$	0.56277	0.53665	0.56648	0.51490	0.51961	0.50350
$E_M^{avg}$	0.56278	0.53682	0.56687	0.51491	0.51857	0.50350
$E_M^{std}$	7.0116 × 10${}^{-4}$	5.4670 × 10${}^{-4}$	5.1064 × 10${}^{-4}$	1.0846 × 10${}^{-3}$	9.4653 × 10${}^{-4}$	5.0999 × 10${}^{-4}$
$A^{avg}$	1.1570$\times {10}^{-1}$	6.0735 × 10${}^{-2}$	1.2060$\times {10}^{-1}$	4.5031 × 10${}^{-4}$	2.1962 × 10${}^{-3}$	4.1079 × 10${}^{-4}$
$A^{std}$	6.6740$\times {10}^{-6}$	1.2467 × 10${}^{-4}$	1.9983 × 10${}^{-4}$	5.9076 × 10${}^{-5}$	3.5112 × 10${}^{-5}$	5.2188$\times {10}^{-6}$
$H_K^{avg}$	6.0009 × 10${}^{-3}$	−4.1736 × 10${}^{-4}$	−6.2487 × 10${}^{-5}$	3.5946 × 10${}^{-4}$	3.2657 × 10${}^{-2}$	−5.2689 × 10${}^{-4}$
$H_K^{std}$	2.1744 × 10${}^{-2}$	2.3103 × 10${}^{-2}$	2.1481 × 10${}^{-2}$	2.5246 × 10${}^{-2}$	2.4646 × 10${}^{-2}$	2.3873 × 10${}^{-2}$
$H_M^{avg}$	1.1570$\times {10}^{-1}$	6.0840 × 10${}^{-2}$	1.2070$\times {10}^{-1}$	−1.5986 × 10${}^{-5}$	2.2020 × 10${}^{-3}$	−1.6685 × 10${}^{-5}$
$H_M^{std}$	2.4011$\times {10}^{-15}$	2.3803$\times {10}^{-15}$	5.7460$\times {10}^{-15}$	1.5901 × 10${}^{-5}$	3.5315 × 10${}^{-5}$	1.7279 × 10${}^{-5}$
$H_C^{avg}$	5.6359 × 10${}^{-2}$	−9.9448 × 10${}^{-5}$	4.5145 × 10${}^{-4}$	5.2491 × 10${}^{-2}$	−1.2703$\times {10}^{-1}$	−1.0241 × 10${}^{-4}$
$H_C^{std}$	6.5119 × 10${}^{-5}$	2.5384 × 10${}^{-4}$	4.5670 × 10${}^{-4}$	2.1550 × 10${}^{-5}$	6.9637 × 10${}^{-5}$	1.7423 × 10${}^{-4}$
$H_P^{avg}$	⋯	⋯	⋯	⋯	−1.3144$\times {10}^{-1}$	⋯
$H_P^{std}$	⋯	⋯	⋯	⋯	5.2471 × 10${}^{-5}$	⋯
${\Omega }^{avg}$	9.7631 × 10${}^2$	1.0387 × 10${}^3$	9.7270 × 10${}^2$	1.1245 × 10${}^3$	1.1277 × 10${}^3$	1.0996 × 10${}^3$
${\Omega }^{std}$	1.5758	1.9354	2.0917	2.5683	2.0880	1.2396
$J^{avg}$	9.7686 × 10${}^2$	1.0393 × 10${}^3$	9.7328 × 10${}^2$	1.1244 × 10${}^3$	1.1282 × 10${}^3$	1.0996 × 10${}^3$
$J^{std}$	1.5981	1.9306	2.1060	2.5085	2.0831	1.2747

is 3.6 cpu-years; thus, the average cpu-time per run presented here is 0.6 cpu-years (each run has its own cpu). However, we have been able to run to an average simulation time of ${\overline{t}}_{end}=1280$, and this has allowed us to see important phenomena in the transition of ideal MHD turbulence from initial conditions to an equilibrium state. In comparison, there is a very interesting run on a ${512}^3$ grid performed by [29] studying energy fluxes in the inertial range of forced, dissipative MHD turbulence using a method similar to ours. This ${512}^3$ run ends at a simulation time of $t_{end}=15$; if we ran our code on a ${512}^3$ grid for this length of simulation time, it would take approximately 14 years, and to have one run going to $t_{end}=1280$ would take about 150 years. However, in the work of [29], the run time was greatly shortened by massively parallel processing [30]; for us, however, this is a future project. Thus, with our current single-cpu approach, there is a trade-off: use a large grid-size such as ${256}^3$ and run for a relatively short simulation time, or a smaller one such as ${128}^3$ and run for much longer simulation times. Each has its advantages and we choose the latter in order to have very-long-duration runs.

Six long-duration Fourier spectral transform method runs on ${128}^3$ grids are considered in this paper. They are the continuations of runs that we reported on in [13], where their transition from non-equilibrium initial conditions was studied. Parameters for the continued runs are given here in Table 2, which is similar in form to Table 2 of [13], except that some new rows were added, including one indicating the statistical expectation value $\langle E_M\rangle$. This has been added in order to compare it to the temporal mean value seen in Table 2; the expectation values and the temporal means match to fractions of a percentage point, providing a validation of the statistical theory.

Initial, non-equilibrium magnetic and kinetic modal energy values (spectra) are ${\overline{\mathcal{E}}}_M\left(k^2\right)\sim {\overline{\mathcal{E}}}_K\left(k^2\right)\sim k^{4}exp(-k^2/k_{\mathrm{o}}^2)$, where $k_{\mathrm{o}}=6$. Viscosity and magnetic diffusivity are set to zero so that the flow is ideal. Again, 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$/\Delta t$; thus, a single run of 2 × 10${}^6\Delta t$s requires about 36 weeks of cpu time.

As mentioned previously, six computer simulations covering the five cases in Table 1 were 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 ${\tilde{u}}_{\pm }(\mathbf{k},t)$ and ${\tilde{b}}_{\pm }(\mathbf{k},t)$, with $0<k\le K$, as defined in (13)–(22). The partial differential equations for MHD in $\mathbf{x}$ space are transformed into a set of ordinary differential equations in $\mathbf{k}$ space. The $\mathbf{k}$ space equations, explicitly shown in [13], are numerically integrated to advance the $\tilde{\mathbf{u}}(\mathbf{k},t)$ and $\tilde{\mathbf{b}}(\mathbf{k},t)$, with nonlinear terms being evaluated by a de-aliased transform method.

4. Computational Results

As seen in Table 1, the integral invariants of ideal MHD turbulence are the volume-averaged energy E and magnetic helicity $H_M$ when ${\mathbf{B}}_{\mathrm{o}}=0$, as well as the cross-helicity $H_C$ when ${\Omega }_{\mathrm{o}}=0$ and $H_P$ when ${\Omega }_{\mathrm{o}}=\sigma {\mathbf{B}}_{\mathrm{o}}\ne 0$. In our numerical simulations, 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 shown in Table 2. We stated that these runs came to ‘near equilibrium.’ To obtain an idea as to what this means, consider Figure 1 where we show, at different times for Runs 1 and 4, the average magnetic and kinetic energy spectra, ${\overline{\mathcal{E}}}_M\left(k^2\right)$ and ${\overline{\mathcal{E}}}_K\left(k^2\right)$, for modes $\mathbf{k}$ having the same value of $k^2={\left|\mathbf{k}\right|}^2$. There are 3036 different values of $k^2$, where $1\le k^2\le 3640$, for $N=128$; the number of independent $\mathbf{k}$ (i.e., $\mathbf{k}$ but not $-\mathbf{k}$) is $\mathsf{n}\left(k^2\right)$.

The definitions of the averaged spectra ${\overline{\mathcal{E}}}_M\left(k^2\right)$ and ${\overline{\mathcal{E}}}_K\left(k^2\right)$ are

$$\begin{array}{ccc}\hfill {\overline{\mathcal{E}}}_M\left(k^2\right)& =& \frac{1}{\mathsf{n}\left(k^2\right)}\sum _{\mathbf{k}}^{{\left|\mathbf{k}\right|}^2=k^2}\frac{{\left|\tilde{b}\left(\mathbf{k}\right)\right|}^2}{N^3},\hfill \\ \hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill {\overline{\mathcal{E}}}_K\left(k^2\right)& =& \frac{1}{\mathsf{n}\left(k^2\right)}\sum _{\mathbf{k}}^{{\left|\mathbf{k}\right|}^2=k^2}\frac{{\left|\tilde{u}\right(\mathbf{k}\left)\right|}^2}{N^3}.\hfill \end{array}$$

(24)

Here, $\mathsf{n}\left(k^2\right)$ is the number of independent $\mathbf{k}$ that satisfies ${\left|\mathbf{k}\right|}^2=k^2$. The number $\mathsf{n}\left(k^2\right)$ jumps around as $k^2$ increases; for example,

$$\begin{array}{cccc}\mathsf{n}\left(1\right)=3,\hfill & \mathsf{n}\left(2\right)=6,\hfill & \mathsf{n}\left(3\right)=4,\hfill & \mathsf{n}\left(4\right)=3,\hfill \\ \mathsf{n}\left(5\right)=12,\hfill & \mathsf{n}\left(6\right)=12,\hfill & \mathsf{n}\left(7\right)=0\hfill & \cdots \hfill \\ \mathsf{n}\left(3628\right)=36,\hfill & \mathsf{n}\left(3629\right)=552,\hfill & \mathsf{n}\left(3630\right)=264,\hfill & \mathsf{n}\left(3631\right)=0,\hfill \\ \mathsf{n}\left(3632\right)=60,\hfill & \mathsf{n}\left(3633\right)=144,\hfill & \mathsf{n}\left(3634\right)=144\hfill & \cdots \hfill \end{array}$$

(25)

We have $\mathsf{n}\left(k^2\right)=0$ whenever $k^2=4^a(8m+7),a,m=0,1,2,\cdots$; see, for example, [31]. The full magnetic and kinetic energy spectra, $\mathsf{n}\left(k^2\right){\overline{\mathcal{E}}}_K\left(k^2\right)$ and $\mathsf{n}\left(k^2\right){\overline{\mathcal{E}}}_M\left(k^2\right)$, jump wildly as $k^2$ increases, which is why we prefer to look at the average energy spectra (23) and (24), as in Figure 1. However, when averaged over nearest neighbors, $\overline{\mathsf{n}}\left(k^2\right)$ is well approximated by $\overline{\mathsf{n}}\left(k^2\right)=\pi k$ and this can be used to give smoothed energy spectra, if desired.

Referring to [23], we find that the ideal expectation values of ${\overline{\mathcal{E}}}_M\left(k^2\right)$ and ${\overline{\mathcal{E}}}_K\left(k^2\right)$ are

$$\begin{array}{ccc}\hfill 〈{\overline{\mathcal{E}}}_M\left(k^2\right)〉& =& \frac{2}{N^3}\frac{\widehat{\alpha }{\widehat{\delta }}^2}{{\widehat{\delta }}^4-{\widehat{\alpha }}^2{\widehat{\gamma }}^2/k^2},\hfill \\ \hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill 〈{\overline{\mathcal{E}}}_K\left(k^2\right)〉& =& \frac{2}{N^3}\frac{\widehat{\alpha }({\widehat{\delta }}^2-{\widehat{\gamma }}^2/k^2)}{{\widehat{\delta }}^4-{\widehat{\alpha }}^2{\widehat{\gamma }}^2/k^2}.\hfill \end{array}$$

(27)

Here, ${\widehat{\delta }}^2={\widehat{\alpha }}^2-{\widehat{\beta }}^2/4$ and $\widehat{\alpha }$, $\widehat{\beta }$ and $\widehat{\gamma }$ are the normalized inverse temperatures related to inverse temperatures appearing in the phase space probability density; please see [13] for details. For each of the cases, $\widehat{\alpha }$, $\widehat{\beta }$ and $\widehat{\gamma }$ are determined by numerically finding the minimum of the entropy functional with proviso that, for Case II (Runs 2a and 2b), $\widehat{\beta }=0$; for Case III (Run 3), $\widehat{\gamma }=0$; and for Case V (Run 5), $\widehat{\beta }=\widehat{\gamma }=0$. The values of $\widehat{\alpha }$, $\widehat{\beta }$ and $\widehat{\gamma }$ for the different runs are given in

Table 3. Values of the inverse temperatures $\widehat{\alpha }$, $\widehat{\beta }$ and $\widehat{\gamma }$ for the Runs in Table 1 with proviso that, for Runs 2a and 2b, $\widehat{\beta }\equiv 0$; for Run 3, $\widehat{\gamma }\equiv 0$; and for Run 5, $\widehat{\beta }\equiv 0$ and $\widehat{\gamma }\equiv 0$. When needed, $\mathcal{E}$, ${\mathcal{H}}_C$ and ${\mathcal{H}}_M$ took their values from $E^{avg}$, $H_C^{avg}$ and $H_M^{avg}$, respectively, in Table 2 (the need for precision here is due to the possible smallness of the denominators in (26) and (27)).

Run	$\widehat{\mathit{\alpha }}$	$\widehat{\mathit{\beta }}$	$\widehat{\mathit{\gamma }}$
1	0.99101483544	$-0.19849208643$	$-0.98106334761$
2a	0.92181618056	0	$-0.92179251089$
2b	0.98391439550	0	$-0.98390250359$
3	0.83658011283	$-0.17565170681$	0
4	0.90153576833	0.42548700154	$-0.85097400308$
5	0.87722015381	0	0

.

The spectra are shown in Figure 1 for various times. The spectra are approaching their expectation values and are near equilibrium; when they match the ideal expectation values (within small fluctuations), then they are in true equilibrium. The spectra ${\overline{\mathcal{E}}}_M\left(k^2\right)$ and ${\overline{\mathcal{E}}}_K\left(k^2\right)$ for Runs 3 and 5 will both be similar to those in Figure 1d, i.e., becoming flat, while the spectra for Runs 2a and 2b will be similar to Figure 1a,b, although not as highly peaked at $k^2=1$.

Now, consider Figure 1, where average magnetic and kinetic energy spectra, ${\overline{\mathcal{E}}}_M\left(k^2\right)$ and ${\overline{\mathcal{E}}}_K\left(k^2\right)$, for Runs 1 and 4 are given. These are averages over modes $\mathbf{k}$ having the same value of $k^2={\left|\mathbf{k}\right|}^2$; there are 3036 different values of $k^2$, where $1\le k^2\le 3640$, for $N=128$. In (a) and (b), the spectra are shown for times $t=0$, 1095 and 1475 for Run 1 ($t=1095$ was the end of Run 1 in [13] and $t=1475$ is the end of the continuation); and in (c) and (d) the spectra are shown for times $t=0$, 895 and 1276 for Run 4 ($t=895$ was the end of Run 4 in [13] and $t=1276$ is the end of the continuation). The spectra in Figure 1 appear to be reaching expected values, although when MHD turbulence is in equilibrium, spectra will be seen to fluctuate about these ideal expectation values; moreover, for Run 4, an extra term may be needed in the perturbation expansion we used to more closely align computed and expected values at $k^2=1$; however, this effort will be deferred. The spectra for Runs 2a, 2b, 3 and 5 behave similarly. The correlation between computed and predicted spectra is another validation of the statistical theory.

In Figure 2, 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) $\Omega$ and mean-squared electric current J; for Run 4, parallel helicity $H_P$ is also shown; again, please see (13)–(22) for precise definitions of these quantities. In Figure 2, the integral invariants for each run appear to be horizontal lines, i.e., have constant values. Time averages and standard deviations for these quantities over the course of the runs is given in Table 2. For each run in Table 2, the quantities that are supposed to be canonical invariants are again seen to be well conserved. 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, where it is part of the invariant parallel helicity $H_P$. Also note that statistics for the enstrophy $\Omega$ and mean-squared current J are essentially equal in each run, indicating that smaller length-scales effectively have equipartition in magnetic and kinetic energy.

In Figure 2, near the ends of some runs, notice that one or both of the helicities $H_M$ and $H_C$ appear to become constants for a period of time. These helicities, as well as mean-average vector potential A for (a) Run 2a, (b) Run 2b, (c) Run 3 and (d) Run 5 are shown again in Figure 3. Here, we more clearly see that $H_M$ or $H_C$ or both become constant in these runs for a significant period of time: (a) $\Delta t=267.4$ (from $t=1140.1$ to 1407.5); (b) $\Delta t=267.1$ (from $t=1041.1$ to 1308.2); (c) $\Delta t=268.4$ (from $t=718.5$ to 986.9) and (d) $\Delta t=266.6$ (from $t=465.4$ to 732.0)—also notice that, during these periods in (c) and (d), A makes a transition to a lower value. The essential equivalence of the $\Delta t$s is intriguing; however, to date, it remains unexplained and a puzzle for the future.

As mentioned earlier, continuation has resolved what seems to be an anomaly occurring in the transition phase, i.e., the appearance of mildly energetic coherent magnetic structures in Runs 3 and 5, as reported in [13]. These quasi-stationary structures are not predicted and, as we will see, eventually disappear, signaling the end of transition for those runs. Another novel feature is the behavior of the helicities $H_M$ and $H_C$ in those runs where one or both were not integral invariants. Specifically, this feature is that of the helicity or helicities that were not invariant became quasi-invariant for a significant interval of time and then resumed their fluctuating behavior. This happened in Runs 2a, 2b, 3 and 5, and signals the end of transition from initial conditions to equilibrium; for Runs 3 and 5, the quasi-invariance of the helicities coincided with the collapse of their transitory coherent structures, as seen in Figure 4.

In the previous work [13], there were many figures showing the trajectories of various complex Fourier coefficients as 2D phase plots (real vs. imaginary parts), but we will not look at these explicitly here. These 2D phase plots are quite interesting in that they are projections onto different complex planes of the random walk through a phase space of dimension 3,679,328 which the dynamical system takes as it evolves in time. These phase trajectories would not really look that different for the continued runs discussed here and the reader is referred back to [13] if there is a desire to view them. Here, we show behavior of the $k=1$ modes through 3D plots of a ‘dipole moment vector’, as described below.

For all the runs in Table 1, the values of all the $\tilde{\mathbf{u}}(\mathbf{k},t)$ and $\tilde{\mathbf{b}}(\mathbf{k},t)$ with $k^2\le 3$ were saved every 0.1 units of simulation time t (i.e., every 200 $\Delta t$s). From these, we can calculate a time history of modal energies

$$\begin{array}{ccc}\hfill E_K(\mathbf{k},t)& =& N^{-3}{\left|\tilde{\mathbf{u}}(\mathbf{k},t)\right|}^2,E_K^{\left(k\right)}\left(t\right)=\sum _{k=\left|\mathbf{k}\right|}{E}_K(\mathbf{k},t),\hfill \\ \hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill E_M(\mathbf{k},t)& =& N^{-3}{\left|\tilde{\mathbf{b}}(\mathbf{k},t)\right|}^2,E_M^{\left(k\right)}\left(t\right)=\sum _{k=\left|\mathbf{k}\right|}{E}_M(\mathbf{k},t).\hfill \end{array}$$

(29)

We now define a ‘dipole moment vector’ $\mathbf{D}$ as one that has components $D_x$, $D_y$, $D_z$ and an angle ${\theta }_D$ with respect to the z axis:

$$D_x=\sqrt{E_M\left(\widehat{\mathbf{x}}\right)},D_y=\sqrt{E_M\left(\widehat{\mathbf{y}}\right)},D_z=\sqrt{E_M\left(\widehat{\mathbf{z}}\right)},{\theta }_D={tan}^{-1}\sqrt{\frac{D_x^2+D_y^2}{D_z^2}}.$$

(30)

The evolution of $\mathbf{D}$ is shown in Figure 4 for each of the runs in Table 1. In Figure 4, the initial growth of $\mathbf{D}$ from $t=0$ is shown in red, while the last 10% of time for each run is shown in black; the periods shown in Figure 3 when $H_M$ or $H_C$ or both become constant appear in Figure 4 in gray for (b) Run 2a, (c) Run 2b, (d) Run 3 and (f) Run 5. Table 1 also has ${\theta }_D$ at the end of each run and there we see that Run 2a (${\theta }_D=14.6^{\circ }$) and 2b (${\theta }_D=10.9^{\circ }$) have vectors $\mathbf{D}$ that are relatively well aligned with the rotation axis. We discuss how this happens in [13].

In Figure 5, the temporal evolution of the magnetic energies defined in (29) at $k=1$, i.e., $E_M\left(\widehat{\mathbf{k}}\right)$, $\widehat{\mathbf{k}}=\widehat{\mathbf{x}},\widehat{\mathbf{y}},\widehat{\mathbf{z}}$, are presented for (a) Run 1, (b) Run 2a, (c) Run 2b and (d) Run 4. In Figure 4, the square roots of the $E_M\left(\widehat{\mathbf{k}}\right)$ defined a ‘dipole moment vector’ which was seen there to become quasi-stationary for these runs. In Figure 5, energies are divided by the corresponding mean values of $|H_M|$ over Runs 1, 2a and 2b, and the mean value of $|H_M|$ over the last 5% of Run 4. Here, we see a verification of the theoretical result that the energy of the dipole is $E_d=k_{min}\left|H_M\right|$, where, for the Fourier case, $k_{min}=1$. This result is predicted for Runs 1, 2a and 2b, where $H_M$ is an ideal invariant, but also seems to apply approximately to Run 4, as seen in (d).

In Figure 6, the temporal evolution of the magnetic energies at $k=1$, i.e., $E_M\left(\widehat{\mathbf{k}}\right)$, $\widehat{\mathbf{k}}=\widehat{\mathbf{x}},\widehat{\mathbf{y}},\widehat{\mathbf{z}}$, are presented for (a) Run 3 and (b) Run 5. The theoretical prediction when $H_M$ is not an invariant is that energy will be equipartioned among all Fourier modes and all $E_M\left(\mathbf{k}\right)$ and $E_K\left(\mathbf{k}\right)$ will have the same value, indicated above by $\langle E_M\left(\widehat{\mathbf{k}}\right)\rangle =\langle E_K\left(\widehat{\mathbf{k}}\right)\rangle =E/2\mathcal{M}$. Here, we see that, during the initial phase of growth starting from the initial spectra shown in Figure 1, the temporary, relatively energetic, coherent structures arise but collapse in the final transition to equilibrium. This collapse occurs during the gray parts of the trajectories shown in Figure 4d,f. In Figure 6, we see that the $E_M\left(\widehat{\mathbf{k}}\right)$ fluctuate about their expectation values $\langle E_M\left(\widehat{\mathbf{k}}\right)\rangle$ and that the anomalously large values that occurred during transition were, in fact, transitory.

Note in Figure 6a,b that the $E_M\left(\widehat{\mathbf{k}}\right)$ with $\widehat{\mathbf{k}}$ parallel to ${\mathbf{B}}_{\mathrm{o}}$ is very small during the transition phase while the other $E_M\left(\widehat{\mathbf{k}}\right)$ are much more energetic. This is an example of anisotropy induced by the mean magnetic field ${\mathbf{B}}_{\mathrm{o}}$ [32,33,34,35], an anisotropy at the largest-scale that is indicated by the ${\theta }_D$ in Table 2.

5. Discussion

The results presented herein were found by continuing to run the simulations first presented in [13]. The continued runs have entered a state of apparent equilibrium and resolved an important open question from that work, specifically, we have seen that the coherent structures that appeared during transition from initial conditions in Runs 3 and 5—which had mean non-aligned magnetic fields—disappeared. Run 4, in which the mean magnetic field was aligned with the rotation axis, maintained a coherent structure, but one that was an order of magnitude less energetic that those in Runs 1, 2a and 2b, which had no mean magnetic field. Runs 1, 2a and 2b are the pertinent ones for modeling those planets and stars that contain a turbulent magnetofluid. The magnetic energy at $k=1$ (the ‘dipole’ energy $E_d$) matched the theoretical prediction that $E_d=k_{min}\left|H_M\right|$ and that the effective ‘dipole moment vector’ $\mathbf{D}$, with dipole angle ${\theta }_D$ and defined in (30), was quasi-stationary, and, Runs 2a and 2b—which most closely modeled a rotating planetary core—showed alignment with the rotation axis, as seen in Table 2 and graphically illustrated in Figure 4. This quasi-stationarity was initially unsuspected [36] since the ensemble prediction is that all variables have zero-mean, so what we have seen, again, is a dynamically broken ergodicity.

Furthermore, the statistical theory of ideal MHD turbulence developed for a periodic box has the same form as that developed for a spherical shell [14]. 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 these dipole magnetic fields arise and is thus a solution to ‘the dynamo problem.’ The statistical theory of ideal MHD turbulence began when [36] recognized that an inverse cascade of energy may move from intermediate scales to the largest scale in the system, as confirmed by [37]. After this initial discovery, it was eventually noticed in 1982 (as discussed in [13]) that this cascade created a coherent structure; the history of the subsequent development built on it is given in the list of related references appearing in [13]. Once it was realized that MHD turbulence, per se, was a dynamo, it was clear, at least to the author, that a deeper understanding of the statistical nature of ideal MHD turbulence, along with numerical simulations of real MHD turbulence [19,20], held the key to solving the ‘dynamo problem.’ That has been our path over many decades and it has led, we believe, to a fundamental understanding of self-generated planetary magnetic fields.

The exact nature of the turbulence within a planet, of course, is hidden from view in its deep interior, but observations of the global magnetic fields of the Earth and other planets indicate that they have long periods of quasi-stationarity in which a magnetic dipole field exists and is closely aligned with a rotation axis. The only requirements that we need to relate our results to these objects is for them to contain a sufficiently spherical shell of relatively incompressible, turbulent magnetofluid that is in a state of quasi-equilibrium, and that this state has an effectively constant (i.e., quasi-stationary) magnetic helicity; how this state is maintained is an important question, but extraneous to our simple (but not too simple) model. Our unique prediction that $E_d=k_{min}\left|H_M\right|$ tells us that if a planet has a magnetic dipole field and we measure its value, then this also tells us that we know its magnetic helicity.

6. Conclusions

Although recent numerical simulations [13] of MHD turbulence corresponding to each of the five cases of ideal MHD turbulence in Table 1 had shown the quasi-stationary largest-scale coherent structure to arise during transition, once these runs were allowed to continue on and enter equilibrium, the transitory structures in Runs 3 and 5, where they were not predicted to exist, were seen to collapse. On the other hand, Runs 1, 2a and 2b had very energetic largest-scale coherent structures, both in transition and equilibrium, that were robust and behaved as predicted. Run 4, where the mean field ${\mathbf{B}}_{\mathrm{o}}$ and the rotation vector ${\Omega }_{\mathrm{o}}$ were parallel, also showed a coherent structure, although an order of magnitude less energetic. The largest-scale coherent structures seen in Runs 2a and 2b, in particular, represent the dipole magnetic fields that are seen to occur around various planets and stars.

The origin of these largest-scale coherent structures in ideal MHD turbulence has been understood through the application of classical statistical mechanics, with extensions, to the interacting set of Fourier coefficients; again, this is detailed in [13]. The difference between the statistical mechanics of an ideal gas and that of MHD turbulence, is that the former has only one invariant, the total energy, while the latter can have one or two invariant helicities. In the case of ideal MHD turbulence, an entropy functional can be defined and the minimization of this functional with respect to magnetic energy leads, through a series of logical steps, to an expression linking the energy of the dipole field to the magnetic helicity contained in the turbulent magnetofluid (as well as to a specific value for the phase entropy). Thus, in dynamo theory, as we have presented it, the statistical mechanics of ideal MHD turbulence plays an essential role and is a key to theoretical understanding of the geodynamo, as well as other planetary and stellar dynamos. Numerical simulation is also another key factor as it offers a validation of the statistical theory.

In this paper, we presented new computational data, discussed new results and resolved an open question left over from previous work. The importance of our theoretical and numerical results, both now and in the past [13], is that they show that MHD turbulence, per se, is a dynamo that produces an energetic, quasi-steady, magnetic dipole that is closely aligned with a rotation axis when one is present. The energy in the ‘dipole’ field and magnetic helicity were predicted and seen to satisfy Equation (1), the fundamental relation $E_d=k_{min}\left|H_M\right|$. This relation is also a prediction for model systems in which a turbulent magnetofluid is confined within a spherical shell, since ideal MHD turbulence in this geometry has the same statistical mechanics as in the Fourier case. In summary, broken ergodicity in MHD turbulence appears to be a viable solution to the ‘dynamo problem.’

A future direction that our investigations can take is to add dissipation and forcing to our ${128}^3$ simulations as was previously performed on ${64}^3$ grids [19,20]; we also wish to move to the larger grid-sizes enabled by massively parallel processing. 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 [19,20], we have always seen behavior analogous to the ideal MHD turbulence at the largest-length scale, specifically the emergence of a largest-scale coherent magnetic structure whose energy is determined by the amount of magnetic helicity contained by the magnetofluid turbulence in those cases (I and II) that are most pertinent to planets and stars. We would expect to see this in ${128}^3$ and greater grid-size simulations that were dissipative and forced. Perhaps some novel effects would also occur, but this, and much else, remains worked that is yet to be performed.