Next Article in Journal
The Tailored CFD Package ‘containmentFOAM’ for Analysis of Containment Atmosphere Mixing, H2/CO Mitigation and Aerosol Transport
Next Article in Special Issue
Physical Background, Computations and Practical Issues of the Magnetohydrodynamic Pressure Drop in a Fusion Liquid Metal Blanket
Previous Article in Journal
Triple Decomposition of Velocity Gradient Tensor in Compressible Turbulence
Previous Article in Special Issue
An Optimized Method for 3D Magnetic Navigation of Nanoparticles inside Human Arteries
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Magnetic Helicity and the Geodynamo

Department of Physics and Astronomy, George Mason University, Fairfax, VR 22030, USA
Fluids 2021, 6(3), 99; https://doi.org/10.3390/fluids6030099
Submission received: 13 January 2021 / Revised: 10 February 2021 / Accepted: 19 February 2021 / Published: 2 March 2021
(This article belongs to the Special Issue Fluids in Magnetic/Electric Fields)

Abstract

:
We present theoretical and computational results in magnetohydrodynamic turbulence that we feel are essential to understanding the geodynamo. These results are based on a mathematical model that focuses on magnetohydrodynamic (MHD) turbulence, but ignores compressibility and thermal effects, as well as imposing model-dependent boundary conditions. A principal finding is that when a turbulent magnetofluid is in quasi-equilibrium, the magnetic energy in the internal dipole component is equal to the magnetic helicity multiplied by the dipole wavenumber. In the case of the Earth, measurement of the exterior magnetic field gives us, through boundary conditions, the internal poloidal magnetic field. The connection between magnetic helicity and dipole field in the liquid core then gives us the toroidal part of the internal dipole field and a model value of 3 mT for the average core dipole magnetic field. Here, we present the theoretical analysis and numerical simulations that lead to these conclusions. We also test an earlier assertion that differential oblateness may be related to dipole alignment, and while there is an effect, rotation appears to be far more important. In addition, the relationship between dipole quasi-stationarity, broken ergodicity and broken symmetry is clarified. Lastly, we discuss how inertial waves in a rotating magnetofluid can affect dipole alignment.

1. Introduction

Explaining the origin of planetary and stellar magnetic fields has been a scientific quest since Larmor [1] hypothesized that magnetohydrodynamic (MHD) motions within the Sun, and by extension, the Earth, were responsible for the creation and maintenance of global magnetic fields. Numerical simulations of the geodynamo [2,3,4] have established that MHD processes within the Earth are capable of creating magnetic fields similar to the actual geomagnetic field, including reversals of the dominant dipole component.
However successful these numerical simulations were in proving the MHD nature of the geodynamo, the fundamental origin of the dominant, quasi-steady geomagnetic dipole field still remained a mystery [5]. In previous work, we showed that extending the statistical mechanics of ideal MHD turbulence opened the door to understanding the fundamental origin of the geodynamo [6,7] and then demonstrated that these ideal results appeared to apply to real (i.e., forced and dissipative) MHD turbulence [8,9]. In Reference [9], we put forward the hypothesis that differences in the oblateness of the outer and inner boundaries of a spheroidal shell could cause dipole alignment. In the present work, we studied the analogue of this with a Fourier model by using stretching-shrinking factors q x , q y , and q z ( q m i n is the smallest and q m a x the largest of these factors). These factors change the relative lengths of the x-, y- and z-directions in both x-space and k-space, as will be described in Section 3.3.
We use a Fourier model as a surrogate for a spherical model because the statistical mechanics of ideal, rotating MHD turbulence are equivalent in the two cases [7]. Both models are based on the MHD equations transformed from physical, or x-space, to a dynamical system of interacting modes in k-space. In either x- or k-space, there is no imposed magnetic field; only overall rotation is imposed. Here, the magnetic field is dynamically generated by MHD turbulence.
We find, in both Fourier and spherical cases, that the largest-scale (smallest wavenumber k) magnetic field modes are very energetic and quasi-steady, i.e., random variables with large mean values subject to only very small fluctuations. The explanation of this is that there is a broken ergodicity: we expected all modes to be fluctuating zero-mean random variables and they mostly are—except for certain largest-scale modes, which will have a very large mean and comparatively very small fluctuations [6]. (Not all of the smallest-k modes need to have the same value—there is a dynamically broken symmetry.) The only requirement for this to occur is that the system have a non-zero magnetic helicity H M . In the case of ideal rotating MHD turbulence, H M is an integral invariant, along with the energy E. In the forced, dissipative, quasi-equilibrium case, driven so that it has a large, relatively constant H M , behavior is equivalent in that a largest-scale, quasi-steady ‘dipole’ component arises [8,9]. Although the Fourier and spherical models have equivalent ideal statistical mechanics and dynamical evolution in k-space, they appear, of course, very different and are not comparable in x-space, where the largest-scale mode in a periodic box has a sinusoidal shape, while in the spherical case it is a spherical harmonic dipole. (At times, we will refer to the largest-scale mode in the Fourier model as the ‘dipole’, although this is meant figuratively and not meant to imply that the x-space structures are similar, only that the corresponding modes in k-space are equivalent).
The smallest wavenumber k magnetic field quasi-steady mode (or modes) may be thought of as a ‘mean field’, but it is not the mean field that is assumed to exist in ‘mean field electrodynamics’ (MFE) [10]. The a priori MFE mean-field is introduced by altering the magnetic field evolution equation, thereby changing the fundamental nature of the MHD equations; the result is the MFE equations. Instead, basing our work on the MHD equations, we see that a quasi-steady, largest-scale component of the magnetic field emerges dynamically and has an explanation by a statistical mechanics with broken ergodicity and broken symmetry, as will be discussed at length in this paper. Thus, a ‘mean field’ emerges from MHD turbulence, per se, and there is no need to impose one a priori.
Let us add that the interesting question of geomagnetic reversals is not addressed here. This presumably depends on the equilibrium of an existing geomagnetic dipole orientation being destabilized by unseen processes in the Earth’s interior or by some natural bi-stability. Here, we are only concerned with a system that is in a state of equilibrium or quasi-equilibrium and important questions as to the cause of reversals are beyond the scope of this paper.
Our principal result is that without magnetic helicity there is no dipole field. Thus, if the Earth has a dipole field, we require that it have magnetic helicity. In fact, by measuring the strength of the Earth’s external dipole magnetic field and by mapping this mathematically onto the core-mantle boundary, we can estimate the value of the magnetic helicity in the core.
This principal result and all of the related results, we believe, further enhance the view that magnetic helicity and MHD turbulence play essential roles in understanding the fundamental cause of the geodynamo (and perhaps other astrophysical dynamos [11,12]). Next, we briefly review our numerical method, then summarize our new theoretical and computational results, followed by detailed explanations of how they were found.

2. Numerical Procedure

A Fourier spectral transform method [13] on an N 3 grid with N = 64 is used; the minimum wave number is k ̲ = q m i n and the maximum wave number is k ¯ = 30.2 q m a x . Time-integration is performed with a third-order Adams Bashforth–Adams Moulton method [14] with a time-step of Δ t = 0.001 ; initial magnetic and kinetic energy spectra are E M ( k ) E K ( k ) k 4 exp ( k 2 / k o 2 ) , where k o = 6 . Viscosity and magnetic diffusivity are set to zero so that the flow is ideal. (A grid size of 64 3 was used so that the ten runs of Table 1 could be completed in a reasonable amount of time with the resources available, an 8-core Linux computer, each core running at 8 s / Δ t ). The integral invariants of ideal MHD turbulence are the volume-averaged energy E and magnetic helicity H M , as well as the cross helicity H C when there is no rotation. The ideal invariants changed less than 0.1% during the runs of Table 1, while the kinetic helicity H K , though it is an ideal invariant for hydrodynamic turbulence, fell to zero very quickly and had small fluctuations about that value. Further details concerning mathematical modeling and numerical procedure may be found in [6,8,9].

3. Summary of Results

Here, we present a brief summary of the major results to be detailed later in this paper. In what follows, both a Fourier model (periodic box) [9] and a spherical harmonic model (spherical shell) are employed, with the Fourier model serving numerically as a surrogate for the spherical harmonic model. While the periodic box is a 3-torus with no boundaries, the spherical shell model [7] has homogeneous boundary conditions, i.e., the normal components of all fields are zero on the spherical surfaces bounding the internal, turbulent magnetofluid; these boundary conditions were introduced by [15]. Here, these boundaries are not coincident with the Earth’s outer or inner core boundaries (whose precise topography is not well known, anyhow), but are set conceptually inside of them a short distance away at the top of each laminar boundary layer; at the core-mantle boundary (CMB), this laminar layer is estimated to be ∼160 m thick [16] or about 0.07% of the thickness of the outer core. It is within this laminar boundary layer that the magnetic field gains a radial component, eventually becoming the external magnetic field. The external field is connected to the internal, turbulent magnetic field by assuming tangential continuity at the model’s upper boundary. This view of magnetic structure appears to be reflected in Figure 1 of Reference [2]. The efficacy of this model is discussed in Section 3.2, below.
Please note that what is meant by the ‘dipole part’ of a Fourier expansion is comprised of those coefficients that have the smallest wavenumbers. In addition, in regard to the magnetic field, the Fourier coefficients that have positive magnetic helicity are b ˜ + ( k ) and those with negative magnetic helicity are b ˜ ( k ) ; similarly, the Fourier coefficients that have positive kinetic helicity are u ˜ + ( k ) and those with negative kinetic helicity are u ˜ ( k ) . The set of coefficients with the same wavevector k will be called a ‘Fourier mode’ and the mode k and the mode k are not independent because, for example, b ˜ + ( k ) = b ˜ + * ( k ) ; coefficients can have the same wavenumber k = | k | and there are no coefficients with k = 0 . The number of independent modes for k ̲ k k ¯ is M ; in our simulations, N = 64 and M = 57,656. In what follows, for definiteness, we choose the global magnetic helicity H M to be positive. The integral invariants of ideal MHD turbulence are, again, the energy E and magnetic helicity H M , as well as the cross helicity H C when there is no rotation; these are all volume-averaged quantities as described in Section 5, as are the magnetic energy E M , kinetic energy E K and related quantities.

3.1. Magnetic Dipole Energy and Magnetic Helicity

The energy in the dipole part of the magnetic field generated by MHD turbulence is equal to the absolute value of the magnetic helicity multiplied by the smallest wavenumber. This result was implicit in our previous formulations [6,7], but further analysis, discussed in more detail in Section 9, reveals that
E M d = k ̲ | H M | .
Here, E M d is the expectation value E M d of the dipole magnetic energy and the expectation value of the magnetic helicity is H M = H M , where H M is a constant for ideal MHD turbulence; expectation values are defined in Section 6. As discussed in Section 4.2, the smallest possible wavenumber is k ̲ and is associated with the dipole, while the corresponding wavevector is k ̲ (in the symmetric case, there are three of them).
The result (1) has been tested on ten new ideal MHD turbulence Fourier method simulations of Table 1, as well as ten real (forced, dissipative) simulations, drawn from [9] that are given in Table 2. Let us denote the ratio of the right and left sides of (1) as R = E M d / k ̲ | H M | . Averaging R over the last 20% of each run gives a value R ¯ for each run, and then averaging the averages R ¯ gives us a mean R ¯ a v g and standard deviation R ¯ s t d for each of Table 1 and Table 2. Using the notation R = R ¯ a v g ± R ¯ s t d , for the ten runs in Table 1, R = 1.000 ± 0.036 and for the ten runs in Table 2, R = 0.998 ± 0.041 . The significance of Table 2 (real MHD turbulence) is that the values of R ¯ are commensurate with those found in Table 1 (ideal MHD turbulence), i.e., real and ideal MHD turbulence behave similarly with regard to the theoretical prediction (1).
The prediction (1) can be visually verified, as seen in Figure 1 and Figure 2, where the lowest-wavenumber positive magnetic helicity coefficients b ˜ + ( k ̲ ) are normalized by N 3 / 2 (since the expectation value of any coefficient’s energy is | b ˜ ± ( k ) | 2 / N 3 or | u ˜ ± ( k ) | 2 / N 3 , where N 3 is the grid size) and by H M in accordance with (1); if all dipole coefficients had equal energies, their doubly normalized magnitudes would be 1 / 3 and if one of these coefficients had all the energy predicted by (1), its magnitude would be unity. Thus, the theoretical prediction (1) appears to be well confirmed through the Fourier method (periodic box) numerical simulations of both ideal and real MHD turbulence.
In addition to the energy E M d contained in the dipole modes of the magnetic field, given by (1), the rest of the magnetic field has an energy E M h (‘h’ stands for higher-wavenumber, non-dipole positive magnetic helicity components, as well as all negative magnetic helicity ones), while the velocity field has kinetic energy E K ; the total energy is E = E M d + E M h + E K . In addition to (1), we find
E M h = E K = 1 2 ( E k ̲ | H M | ) .
The essential result (2), along with (1), is discussed in Section 9. Using (1) and (2), we clearly have E = E M d + E M h + E K .
Next, we discuss the application of the principal result (1) to the Earth’s outer core. We use the spherical shell model of [7] and assume that ideal results are qualitatively applicable to the low-wavenumber (large length-scale) behavior of a real turbulent magnetofluid. This assumption was tested previously numerically on a Fourier model and seen to be a fairly good one at low wavenumber [9]. Forced, dissipative MHD turbulence behaves very similarly to ideal MHD turbulence at the largest wavelengths and important measures of MHD turbulence are commensurate, as application of the ideal result (1) to forced, dissipative numerical data shows us in Table 2.

3.2. Knowledge of the Internal Toroidal Field

The Fourier model examines MHD turbulence in a periodic box, for which there is no ‘outside’ but only an ‘inside’. Thus, the results above are pertinent to the flow interior to the Earth’s fluid (outer) core. We have considered MHD turbulence in the outer core by developing a spherical shell model that uses spherical Bessel function-spherical harmonic expansions and have found that the statistical theory of ideal MHD turbulence developed for the Fourier model system also applies to this case [7]. In this spherical shell model, there are poloidal a l m n and toroidal b l m n magnetic field coefficients, and it is the poloidal coefficients a l m n that can be matched to the Gauss coefficients g l m and h l m at the outer boundary of the numerical model. The toroidal coefficients b l m n are internal to the outer core and it had long been assumed that these were unknowable. However, consider the following results.
A turbulent magnetofluid in the rotating spherical shell model of [7] has the ideal invariants energy E , where E = E M + E K , and magnetic helicity H M . In that model, the magnetic energy E M and helicity H M are given by
E M = 1 2 l , m , n | b l m n | 2 + k l n 2 | a l m n | 2 ,
H M = l , m , n a l m n * b l m n .
The indices l and m are the degree and order of a spherical harmonic Y l m ( θ , ϕ ) , and n denotes the nth zero k l n of g ^ l ( k l n r ) at r = 1 , where g ^ l ( k l n r ) is a linear combination of spherical Bessel and Neumann functions. Please see [7] for more details; in that paper, the smallest spherical shell dimensionless wavenumber for an Earth-like model is k ̲ = k 11 1.8638 (the next smallest wavenumber is k 21 2.1497 ). The spherical shell version of (1) is
E M d = k 11 | H M | .
The geomagnetic field is measured at the surface out to spherical harmonic degree l = 13 [17], where about 93% of the magnetic energy at the Earth’s surface is in the dipole, l = 1 . However, at the CMB, only about 38% of its magnetic energy is locked up in the dipole l = 1 components; this percentage is probably high, as there is presumably more energy that is unaccounted for in the l > 13 components. The dipole’s Gauss coefficients g 1 m and h 1 m given by [17] indicate that 97% of the dipole’s energy is held by g 1 0 , so the dipole moment is closely aligned with the rotation axis.
Theoretically, we expect that the components a 1 m 1 and b 1 m 1 that have the smallest wave number k ̲ = k 11 contribute most to H M in (4) [7]. Truncating the summations (3) and (4) to the l = 1 and n = 1 terms and setting | H M | = + H M , for H M > 0 , we have
m = 1 1 b 1 m 1 k 11 a 1 m 1 2 = 0 .
Therefore, if we know the major part of the poloidal dipole magnetic field, then we also know the major part of the toroidal dipole field:
b 1 m 1 = k 11 a 1 m 1 .
Note that, if we had set | H M | = H M , the result would be b 1 m 1 = k 11 a 1 m 1 .
The important point here is that because we can relate the coefficients a 1 m 1 to the Gauss coefficients g 1 m and h 1 m through the boundary conditions, we know both the poloidal coefficients a 1 m 1 and, through (7), the toroidal coefficients b 1 m 1 in the outer core. Using the results in References [7,17], we find, with H M > 0 ,
b 101 = k 11 a 101 = 2 C g 1 0 , b 111 = k 11 a 111 = C ( g 1 1 i h 1 1 ) .
Formulas in Reference [7] tell us that the conversion factor is C 19.5 for degree l = 1 , leading to an estimated mean value of the dipole magnetic field in the outer core of 3.02 mT, or 11.7 times that of average geomagnetic dipole field at the CMB, which is 0.259 mT, based on the 2015 International Geomagnetic Reference Field (IGRF) coefficients [17]. Although the value of 3.02 mT comes from a model with homogeneous boundary conditions [7], the model appears to have some robustness because 3.02 mT is compatible with the value of 2.5 mT estimated by [18] as the mean strength of the total magnetic field in the Earth’s outer core. In addition, the IGRF data tells us that E M d = 45.7 J/m 3 and H M = 24.5 J/m 2 .
In the case of the dipole, l = 1 , using k 11 is appropriate because it is the minimum wavenumber k ̲ of the spherical shell model. However, for estimating higher multipole ( l 2 ) field strengths, it is unclear which wavenumber k l n , or even how many, to use because the IGRF field components are indexed by degree l and order m, and, in the spherical shell model of [7], there are many wavenumbers k l n associated with each spherical harmonic degree l, but, at present, no obvious procedure for picking out which ones are important.

3.3. Differential Oblateness vs. Rotation in Dipole Alignment

In Reference [9], we put forward the idea that a greater oblateness of the inner core boundary (ICB) with respect to the CMB in the z-direction might influence dipole alignment. In the Fourier case, the analogy is that the periodic box would be stretched in the z-direction compared to the x- and y-directions in x-space; in k-space, this would entail a shrinking in the k z -direction compared to the k x and k y -directions (or some permutation of this). Thus, k j q j k j , j = x , y , z , while x x / q x , y y / q y and z z / q z . A number of runs were performed, using various values of q x , q y and q z , with q x q y q z = 1 , as Table 1 shows. In a numerical simulation with H M > 0 , the effect of making one of the q j smaller than the others is to give the positive magnetic helicity Fourier coefficient, say b ˜ + ( q x x ^ ) with smallest wave number k = q x = q m i n , more energy than it would have if all the q j = 1 ; the reason for this is discussed in Section 8.
Let us again compare the evolution of the b ˜ + ( k , t ) , q m i n k q m a x , for the four runs C0, S10, S20, and S30 of Table 1. This comparison is made using the trajectories shown in Figure 1, where we plot the real vs. imaginary parts of the doubly normalized b ˜ + ( k , t ) . In Figure 1, we see that, when the smallest wavevector is in the (b) x-, (c) y- or (d) z-directions, it induces an increase in the magnitude of the corresponding b ˜ + ( k , t ) compared to (a), where all q j = 1 . Although it was initially assumed that all coefficients would be zero-mean random variables [19,20], we clearly see ‘broken ergodicity,’ a concept defined by [21]. Broken ergodicity in MHD turbulence was first noted by [22]; in Figure 1, the phenomenon of ‘broken symmetry’ [23] also appears, in that nonlinear evolution can lead to a lack of equipartition even in the q x = q y = q z = 1 run C0.
Next, consider Figure 2 where (a) pertains to reference run S10 and Figure 2b–d, which pertain to three runs whose parameters in Table 1 are the same as run S10, except that they have an imposed rotation: in Figure 2b, S1X has Ω o = 10 x ^ , in (c), S1Y has Ω o = 10 y ^ and in (d), S1Z has Ω o = 10 z ^ . Figure 2 shows that the effect of rotation overrides the effect of differential oblateness. The low-wavenumber, positive magnetic helicity coefficients b ˜ + ( k ) with k parallel to Ω o develop much greater magnitude than those with k perpendicular to Ω o . We discuss a possible reason for this in Section 3.5 below.

3.4. Broken Ergodicity and Broken Symmetry

Figure 1 and Figure 2 show that the dipole components b ˜ + ( k ̲ ) exhibit quasi-stationarity, i.e., coherent structure due to broken ergodicity and broken symmetry. We discuss these phenomena in greater detail in Section 10. We believe we have a clearer understanding of their relationship and wish to summarize this here.
Assume q x = q y = q z = 1 and define a vector v ˜ d in a 3D complex space by
v ˜ d T = N 3 / 2 b ˜ + ( x ^ ) b ˜ + ( y ^ ) b ˜ + ( z ^ )
This has the dot product | v ˜ d | 2 = v ˜ d v ˜ d = | H M | . At least one of the components of v ˜ d will have a large mean value of magnitude | H M | that is much larger than the fluctuations it experiences, which are O M 1 / 2 . The ‘dipole vector’ v ˜ d will thus evolve to become quasi-stationary and, while the expectation value v ˜ d is zero, the time-average of v ˜ d is almost constant and has a magnitude M 1 / 2 greater than its standard deviation, as is seen in Figure 1. This is the origin of broken ergodicity.
The phase space of independent Fourier components can also be subject to an arbitrary unitary transformation. While the choice of a particular set of initial conditions leads to a specific direction for the quasi-stationary vector v ˜ d , a unitary transformation can be made so as to point it in any direction we chose: The canonical ensemble represents all directions, which average to zero, but the dynamical system can only evolve so that v ˜ d points in one direction. This is the origin of broken symmetry.
Rotation, however, induces dipole alignment, though the reason why is not yet completely clear.

3.5. Rotation and Dipole Alignment

As mentioned above and shown in Figure 1, differential oblateness can cause alignment, but, as shown in Figure 2, rotation appears to be the absolute controlling factor. The rotation vector Ω o enters explicitly into the basic equations only in the velocity Equation (10) and not in the magnetic induction Equation (11). However, it does enter into the magnetic induction equation implicitly through the electromotive field E = u × b because of the hydrodynamic inertial waves present in the velocity field u of a rotating system. While our discussion here will be brief, inertial waves, both Alfvén and hydrodynamic, are well known to be important in planetary cores, as discussed in detail by [24] and [25], respectively.
Here, we examine the inertial wave effects in the context of our Fourier model system. Since inertial waves may introduce a sinusoidal time variation into the electromotive field, the evolution of some of the Fourier magnetic dipole modes can be influenced: if Ω o = Ω o z ^ , then b ˜ + ( x ^ ) and b ˜ + ( y ^ ) are strongly affected, while b ˜ + ( z ^ ) is not affected at all, to leading order in M 1 / 2 . A qualitative appraisal is given more detail in Section 11.

4. Mathematical Model

4.1. Basic Equations

The non-dimensional form of the 3D incompressible MHD equations in a rotating frame of reference with constant angular velocity Ω o can be written as
ω t = × u × ω + 2 Ω o + j × b + ν 2 ω ,
b t = × u × b + η 2 b .
These are described, for example, by References [26,27]. Here, u ( x , t ) and b ( x , t ) are the turbulent velocity and magnetic fields, respectively. The velocity and magnetic fields are solenoidal: · u = · b = 0 , as are the vorticity ω ( x , t ) and electric current density j ( x , t ) , defined by
ω = × u , j = × b .
Non-dimensional density does not appear in (10) because it equals unity. The symbols ν in (10) and η in (11) 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; ν = η = 0 for ideal MHD turbulence.

4.2. Fourier Representation

Discrete Fourier transforms for u and b , connecting x-space (with vectors x ) to k-space (with vectors k ), are
u ( x , t ) b ( x , t ) = k u ˜ ( k , t ) b ˜ ( k , t ) exp i x · k N 3 / 2 ,
u ˜ ( k , t ) b ˜ ( k , t ) = x u ( x , t ) b ( x , t ) exp i x · k N 3 / 2 .
Here, N is the number of grid points in each x-space dimension. Note that the reality of u ( x , t ) and b ( x , t ) imply that u ˜ ( k , t ) = u ˜ * ( k , t ) and b ˜ ( k , t ) = b ˜ * ( k , t ) .
Although the periodic box of the model system is usually taken to be a cube, here we will allow for the possibility that it might be elongated or compressed differently in the x-, y-, and z-directions; this simulates the difference in oblateness that may occur in the Earth’s inner and outer cores and might have dynamical importance [9]. To accomplish this elongation or compression, we define x and k by
x = x x ^ + y y ^ + z z ^ = 2 π N n x q x x ^ + n y q y y ^ + n z q z z ^ ,
k = k x x ^ + k y y ^ + k z z ^ = q x m x x ^ + q y m y y ^ + q z m z z ^ .
The components n j and m j , j = x , y , z , are integers. The n j satisfy 0 n j < N , while the integers m j lie in the range N / 2 + 1 and + N / 2 ; thus, there are N 3 points in both spaces. No matter what the parameters q j > 0 are, the dot product k · x in Equations (13) and (14) is then always
k · x = 2 π N n x m x + n y m y + n z m z .
The Fourier transforms (13) and (14) are thus unaffected by the parameters q j , j = x , y , z . Here, we will choose q x q y q z = 1 so that the volume of the periodic box does not change even when any of the q j 1 .
For computational purposes, we define the following vectors with integer components defined above:
n = n x x ^ + n y y ^ + n z z ^ , m = m x x ^ + m y y ^ + m z z ^ .
The Fourier coefficients u ˜ ( k , t ) and b ˜ ( k , t ) are nonzero only for 1 | m | M < N / 2 . The exact value of M is set by a de-aliasing requirement of [28]: M 2 is the largest integer such that M 2 2 N 2 / 9 . (Note: Although we will use k as the magnitude of k , we will not use n as the magnitude of n , nor m as the magnitude of m because we will use n and m as general indices that appear in various mathematical terms).
Let q m i n equal the smallest of q x , q y , and q z , and let q m a x be the largest; then, the largest length-scale corresponds to smallest wavenumber k ̲ = q m i n , associated with k = k ̲ , and the smallest length-scale to then largest wavenumber k ¯ = q m a x M , associated with k = k ¯ . Dynamically interacting coefficients fit within a sphere in m -space (ellipsoid in k-space); all coefficients outside this m -sphere, as well as at m = 0 , are initially zero and remain so during any numerical simulation. Since u ˜ ( k , t ) = u ˜ * ( k , t ) and b ˜ ( k , t ) = b ˜ * ( k , t ) , only half of the m that satisfy 1 | m | M identify independent coefficients. Therefore, the number of independent modes k is M 2 π M 3 / 3 ; here, for N = 64 , M 57,515, while the actual numerical count is M = 57,656.
In k-space, the requirements · u = · b = 0 become i k · u ˜ ( k , t ) = i k · b ˜ ( k , t ) = 0 . Thus, u ˜ ( 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 e ^ 3 ( k ) = k / k = k ^ ; then choosing a unit vector e ^ 1 ( k ) orthogonal to e ^ 3 ( k ) , and then the remaining unit vector e ^ 2 ( k ) is a vector product of the first two, forming a right-handed orthonormal basis for each k :
e ^ 1 ( k ) · e ^ 3 ( k ) = 0 , e ^ 2 ( k ) = e ^ 3 ( k ) × e ^ 1 ( k ) , e ^ i ( k ) · e ^ j ( k ) = δ i j , e ^ 1 ( k ) · e ^ 2 ( k ) × e ^ 3 ( k ) = 1 .
(This is similar to a Craya–Herring decomposition [29]). Explicit choices [30] for any of the e ^ j ( k ) will be mentioned as needed.
In terms of the e ^ j ( k ) defined above, the Fourier vector coefficients u ˜ ( k ) and b ˜ ( k ) are
u ˜ ( k , t ) = u ˜ 1 ( k , t ) e ^ 1 ( k ) + u ˜ 2 ( k , t ) e ^ 2 ( k ) ,
b ˜ ( k , t ) = b ˜ 1 ( k , t ) e ^ 1 ( k ) + b ˜ 2 ( k , t ) e ^ 2 ( k ) .
However, an equivalent, but more efficacious helical representation can be used:
u ˜ ( k , t ) = u ˜ + ( k , t ) e ^ + ( k ) + u ˜ ( k , t ) e ^ ( k ) ,
b ˜ ( k , t ) = b ˜ + ( k , t ) e ^ + ( k ) + b ˜ ( k , t ) e ^ ( k ) .
Here, the positive and negative helicity unit vectors and components are
e ^ ± ( k ) = 1 2 e ^ 1 ( k ) ± i e ^ 2 ( k ) ,
u ˜ ± ( k , t ) = 1 2 u ˜ 1 ( k , t ) i u ˜ 2 ( k , t ) ,
b ˜ ± ( k , t ) = 1 2 b ˜ 1 ( k , t ) i b ˜ 2 ( k , t ) .
Note that e ^ ± * ( k ) = e ^ ( k ) . The orthonormality properties of the e ^ ± ( k ) are
e ^ ± ( k ) · e ^ ± * ( k ) = 1 , e ^ ± ( k ) · e ^ ± ( k ) = 0 = e ^ 3 ( k ) · e ^ ± ( k ) .
An important property of the helical unit vectors concerns the curl operation:
i k × e ^ ± ( k ) = ± k e ^ ± ( k ) ,
The vorticity and current in helical form are then
ω ˜ ( k , t ) = k [ u ˜ + ( k , t ) e ^ + ( k ) u ˜ ( k , t ) e ^ ( k ) ] ,
j ˜ ( k , t ) = k [ b ˜ + ( k , t ) e ^ + ( k ) b ˜ ( k , t ) e ^ ( k ) ] .
Thus, the helical ± components of vorticity ω ˜ ± ( k , t ) = ± k u ˜ ± ( 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, though numerical simulation is done using the non-helical representation (19) and (20).

4.3. A Dynamical System

The Fourier-transformed 3D MHD equations are found by placing expansions for ω ( x , t ) and b ( x , t ) of the form (13) into (10) and (11). The result is a set of coupled, nonlinear ordinary differential equations, i.e., the dynamical system
d ω ˜ ( k , t ) d t = S ˜ ( u , ω ; k , t ) + S ˜ ( j , b ; k , t ) + 2 i ( k · Ω o ) u ˜ ( k , t ) ν k 2 ω ˜ ( k , t ) ,
d b ˜ ( k , t ) d t = S ˜ ( u , b ; k , t ) η k 2 b ˜ ( k , t ) .
The nonlinear terms denoted by S ˜ are vector convolutions:
S ˜ ( u , b ; k , t ) = i N 3 / 2 k × p + q = k u ˜ ( p , t ) × b ˜ ( q , t ) .
The double summation in (32) is over all wavevectors p and q inside the truncation volume in k-space that satisfy p + q = k . In our 3D numerical simulations, it is the Equations (30) and (31) that are integrated forward in time, with nonlinear terms such as (32) being evaluated by fast Fourier transform (FFT) methods, as needed.
In what follows, we will, again, set ν = η = 0 . The motivation for this comes from [9], where it was seen to a good approximation that quasi-stationary, forced, dissipative MHD turbulence behaves ideally at low-k where large-scale, coherent structure is manifested; the ideal case allows us to avoid any idiosyncrasies due to the particular choice of forcing procedure. When ν = η = 0 , the dynamical systems (30) and (31) conserve energy and magnetic helicity, and, if Ω o = 0 , also cross helicity; these invariants will be discussed in Section 5. However, Ω o 0 for the Earth and most other planets, so that energy and magnetic helicity are of primary interest.
If Equation (30) is linearized (with ν = 0 ), and remembering that ω ˜ ± ( k , t ) = ± k u ˜ ± ( k , t ) , we have
d u ˜ ± ( k , t ) d t = ± i Ω k u ˜ ± ( k , t ) , Ω k = 2 k k · Ω o .
These linearized modes represent inertial waves with frequencies Ω k , so that u ˜ ± ( k , t ) exp ( ± i Ω k t ) . The corresponding periods are T k = 2 π / | Ω k | π k ¯ / Ω o ; typically, these periods will be much less than total time over which the dynamical system (30) and (31) evolves and is studied.
In a Fourier representation, x-space is modeled as a periodic box, i.e., a 3-torus. However, a representation in terms of spherical Bessel functions and spherical harmonics can also be formulated, if x-space is modeled as either a sphere [31] or as a more Earth-like spherical shell [7]. This spherical shell representation provides a dynamical system similar to (30) and (31), one with the same invariants as the Fourier case and the same statistical mechanics. Thus, the Fourier model can serve as a surrogate for a spherical shell model, enabling us to study the physics of the MHD turbulence, at least qualitatively, that exists in the Earth’s outer core. At present, this is a necessity, since a transform method numerical simulation based on such a spherical shell model is not available, although a non-transform code intrinsically limited to an N 9 grid has been developed and used to study MHD flows in a sphere [31,32]. However, this grid size is too small for our purposes, so we await the development of a transform method code.

5. Global Quantities

There are various important global quantities that can be expressed as averages over either x-space or, equivalently, k-space. We define the volume average of a quantity Φ ( x , t ) multiplied by a quantity Ψ ( x , t ) in a periodic box of side length 2 π as Φ Ψ , where
Φ Ψ ( 2 π ) 3 Φ ( x , t ) Ψ ( x , t ) d 3 x = 1 N 3 k Φ ˜ * ( k , t ) Ψ ˜ ( k , t ) .
Then, the volume-averaged energy E, enstrophy Ω , 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 a , where × a = b ) are
E K = 1 2 u 2 , E M = 1 2 b 2 , Ω = 1 2 ω 2 , J = 1 2 j 2 ,
E = E K + E M , H C = 1 2 u · b , H M = 1 2 a · b , A = 1 2 a 2 .
Since all functions are periodic in x-space, we can use Equations (10) and (11), along with integration by parts to derive the following relations [33]:
d E d t = 2 ν Ω + η J ,
d H C d t = Ω o · b × u 1 2 ( ν + η ) j · ω ,
d H M d t = η j · b .
When ν = η = 0 and Ω o = 0 , the quantities E, H C , and H M are the traditional ideal invariants of MHD turbulence [20,34,35]. If Ω o 0 , then (38) indicates that H C is no longer an ideal invariant. If an external mean magnetic field B o was imposed, then H M would also no longer be an ideal invariant [33]; here, we will always have B o 0 . We note that ‘generalized helicities’ G C and G M related to H C and H M can be defined [36] and are ideal invariants even in the presence of nonzero B o and Ω o . However, we will only need to refer to H C and H M in this paper.
Although kinetic helicity H K is not an ideal MHD invariant, it is an invariant of ideal fluid turbulence [37]. Manipulating (10), we find that, for a periodic box, we have
d H K d t = ω · j × b ν i u · i ω , H K = 1 2 u · ω .
Note that H K and E are both conserved if b ( x , t ) 0 and ν = 0 , i.e., in the case of ideal hydrodynamic turbulence.
Since b = × a , then b ˜ ( k , t ) = i k × a ˜ ( k , t ) . Using (19), (23) and (25), we see that
a ˜ ( k , t ) = i k 2 k × b ˜ ( k , t ) = 1 k b ˜ + ( k , t ) e ^ + ( k ) b ˜ ( k , t ) e ^ ( k ) .
Using this, energy E, cross helicity H C and magnetic helicity H M are represented in k-space by the quadratic sums:
E = 1 2 N 3 k | u ˜ ( k , t ) | 2 + | b ˜ ( k , t ) | 2 ,
H C = 1 2 N 3 k u ˜ + ( k , t ) b ˜ + * ( k , t ) + u ˜ ( k , t ) b ˜ * ( k , t ) ,
H M = 1 2 N 3 k 1 k | b ˜ + ( k , t ) | 2 | b ˜ ( k , t ) | 2 .
In the ideal MHD case with either periodic or homogeneous b.c.s, E and H M are invariants, as is H C when Ω o = 0 , but other quantities, such as E K , E M , H K , Ω and J, as well other functions of u ˜ ( k , t ) and b ˜ ( k , t ) , will generally be time-dependent, particularly in numerical simulations during the transition from initial conditions to an equilibrium state.

6. Statistical Mechanics

Here, we briefly discuss the statistical mechanics of ideal MHD turbulence. The Equations (30) and (31) form a finite dynamical system with a phase space defined by the independent real and imaginary components of u ˜ ( k ) and b ˜ ( k ) , k ̲ k k ¯ ; we denote that these are coordinates in phase space and not dynamical variables by omitting t from the arguments. If N = 64 and the number of independent k is M = 57,656, then the phase space has 8 M = 461,248 dimensions. As pointed out by [20], when ν = η = 0 , the system has a Liouville theorem and whose phase space Γ represents a canonical ensemble where statistical properties depend on the constants of the motion. As just discussed, these constants, also known as ideal invariants, are the energy E, the magnetic helicity H M and, if Ω o = 0 , the cross helicity H C .
The statistical mechanics of ideal MHD turbulence are based on a probability density function of the form
D = Z 1 exp ( α E β H C γ H M ) ,
Here, E, H C , and H M are given by (42), (43) and (44), respectively, and Z is the partition function; in addition, β = 0 if Ω o 0 . [20] initiated this approach and explicitly considered the dynamical system to be ergodic, an assumption that was unchallenged in the early work on ideal MHD turbulence [19,38]. It was finally challenged by [22], when apparent non-ergodicity was first noticed and reported, and confirmed later [39]. As already mentioned, this non-ergodicity is actually ‘broken ergodicity’ [21]; again, a review of broken ergodicity for ideal MHD turbulence is given by [6].
In the case of ideal MHD turbulence, expectation values of the various global quantities in (35) and (36), as well as any other phase functions, can be determined with respect to the probability density function (45). 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 u ˜ ± ( k ) and b ˜ ± ( k ) (again, the t is omitted from the arguments because these are coordinates and not dynamical variables). Given a quantity Q, the expectation value Q is defined by
Q Q D d Γ .
As an example, the velocity and magnetic field coefficients are expected to have zero mean values:
u ˜ ( k ) = b ˜ ( k ) = 0 .
In the ideal case, the integral invariants E, H M , and possibly H C , should have time-independent values E , H M , and H C that are equal to their expectation values:
E = E , H M = H M , H C = H C .
In fact, we require that (48) be true and use this to determine the ‘inverse temperatures’ α , β and γ in (45). Whereas (47) is an ‘ergodic hypothesis,’ (48) 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 Section 10.

7. Eigenvariables and Entropy

Placing the k-space representation of E, H C and H M , as given in (42)–(44), into the PDF (45) gives an expression that contains modal 4×4 Hermitian covariance matrices in the argument of the exponential:
D = k D ( k ) , D ( k ) = exp y ˜ ( k ) M k y ˜ ( k ) Z ( k ) .
In k , the notation k means that only independent modes k are included, i.e., if k is included, then k is not. Here, y ˜ = y ˜ * T is the Hermitian adjoint ( T means transpose) of the column vector y ˜ , where
y ˜ ( k ) = [ u ˜ + ( k ) u ˜ ( k ) b ˜ + ( k ) b ˜ ( k ) ] T
The Hermitian (here, real, and symmetric) 4×4 covariance matrix M k is
M k = α ^ 0 β ^ / 2 0 0 α ^ 0 β ^ / 2 β ^ / 2 0 α ^ + γ ^ / k 0 0 β ^ / 2 0 α ^ γ ^ / k .
The circumflex indicates division by N 3 : α ^ = α / N 3 , β ^ = β / N 3 and γ ^ = γ / N 3 .
Although the M k in (51) can also be expressed as 8 × 8 real symmetric matrices and the y ˜ ( k ) as 8 × 1 real arrays [19], finding eigenvalues and eigenvariables is facilitated by using the 4 × 4 matrices M k and 4 × 1 complex arrays y ˜ , along with the properties of block matrices given by [40].
The real, symmetric matrices M k can be diagonalized (and more easily than the Hermitian matrices used previously [6,30,41]) to yield the modal PDFs
D ( k ) = n = 1 4 D n ( k ) , D n ( k ) = 1 Z n ( k ) exp λ ^ k ( n ) | v ˜ n ( k ) | 2 , Z n ( k ) = π λ ^ k ( n ) .
The eigenvalues λ ^ k ( n ) = λ k ( n ) / N 3 are also written with a circumflex to indicate division by N 3 , just as for α ^ , β ^ and γ ^ . Implicit in the form of D n ( k ) given above is the transformation y ˜ ( k ) = U k v ˜ ( k ) , where U k SU(4) is a unitary transformation matrix (see below). Explicitly, v ˜ ( k ) is
v ˜ ( k ) = [ v ˜ 1 ( k ) v ˜ 2 ( k ) v ˜ 3 ( k ) v ˜ 4 ( k ) ] T .
The energy expectation values for the complex eigenvariables v ˜ n ( k ) is
E n ( k ) = v ˜ n ( k ) / N 3 = 1 / λ k ( n )
This energy contains equal contributions from the real and imaginary parts of v ˜ n ( k ) . The exact forms of the λ ^ k ( n ) and v ˜ n ( k ) in terms of α ^ , β ^ and γ ^ will be presented next.

7.1. Eigenvariables

The eigenvariables v ˜ n ( k ) in (52) can be determined for ideal MHD turbulence through a modal unitary transformation [6,7,30]. In the general case (nonrotating with zero mean magnetic field), the form of the transformation U is
v ˜ 1 ( k ) = + β ¯ ζ k u ˜ ( k ) ζ k + b ˜ ( k ) ,
v ˜ 2 ( k ) = + β ¯ ζ k u ˜ + ( k ) + ζ k + b ˜ + ( k ) ,
v ˜ 3 ( k ) = + β ¯ ζ k + u ˜ ( k ) + ζ k b ˜ ( k ) ,
v ˜ 4 ( k ) = β ¯ ζ k + u ˜ + ( k ) + ζ k b ˜ + ( k ) .
Above, β ¯ = sgn β ^ with β ¯ = 1 for β = 0 ; the functions ζ k + ( β ^ , γ ^ ) and ζ k ( β ^ , γ ^ ) are, in terms of a third function η ^ k ( β ^ , γ ^ ) ,
ζ k ± = 1 2 1 ± γ ^ k η ^ k ; η ^ k = β ^ 2 + γ ^ 2 k 2 .
In terms of η ^ k , as defined above, the eigenvalues λ ^ k ( n ) ( n = 1 , 2 , 3 , 4 ) are
λ ^ k ( 1 ) = α ^ 1 2 ( η ^ k + γ ^ / k ) , λ ^ k ( 2 ) = α ^ + 1 2 ( η ^ k + γ ^ / k ) ,
λ ^ k ( 3 ) = α ^ + 1 2 ( η ^ k γ ^ / k ) , λ ^ k ( 4 ) = α ^ 1 2 ( η ^ k γ ^ / k ) .
Although it appears that the eigenvalues given above are functions of the undetermined quantities α ^ , β ^ , and γ ^ , there is only one unknown to be determined: φ E M . Theory [6,30] tells us that
α ^ = 2 ϱ φ φ ( E φ ) H C 2 , β ^ = 2 H C φ α ^ , γ ^ = 2 φ E H M α ^ .
Here, ϱ = M / N 3 0.2194 and, again, α ^ = α / N 3 , β ^ = β / N 3 and γ ^ = γ / N 3 . Basic results of ideal MHD turbulence theory [6] are that α ^ > 0 and φ E / 2 ; thus, in the expression for γ ^ , we have 2 φ E 0 . Noting that H C and H M are pseudoscalars, we see that β ^ and γ ^ are also pseudoscalars and that β ^ H C 0 and γ ^ H M 0 .

7.2. Entropy

We use (52) to find the entropy functional σ ( φ ) = ln D / N 3 :
σ ( φ ) = 4 ρ ( 1 + ln π ) 1 N 3 k ln α ^ 2 β ^ 2 / 4 2 α ^ 2 γ ^ 2 / k 2
Above, the sum over k means, again, that only independent modes k are included (if k , then not k ). The fact that there is only one unknown quantity φ in (62) means that the entropy functional of the system is a function of only one variable. As discussed by [42], finding the (single) minimum of σ ( φ ) gives us the value φ = φ o that sets the values of α ^ , β ^ and γ ^ , as well as the system entropy s = σ ( φ o ) .
Both the non-rotating ( Ω o = 0 ) and rotating cases ( Ω o 0 ), have γ ^ 0 , and the first derivative of the entropy functional (63) with respect to φ is
d σ ( φ ) d φ = σ ( φ ) = 2 N 3 F ( φ ) G + ( φ ) + G ( φ )
F ( φ ) = φ 3 H C 2 3 φ E φ 2 φ E φ H C 2 > 0
G ± ( φ ) = k k | H M | ± φ k | H M | 1 H C 2 / φ 2 ± 2 φ E
Theory tells us that the denominators for the terms in G ± ( φ ) are positive and also that φ / | H M | k ̲ . For the Fourier case we are discussing here, k ̲ = q m i n , and for the spherical shell model of the outer core developed by [7], k ̲ 1.8638 . We define a wavenumber k M = φ / | H M | , so that k M k ̲ . For purposes of discussion, let us draw from examples found in [9]: for run NM06, E a v g = 1.0183 , H M a v g = 0.6179 and φ o = 0.8181 , so that k M = 1.324 , while, for run NM06c, E a v g = 1.0462 , H M a v g = 0.8491 and φ o = 0.9476 , so that k M = 1.1160 ; either could be used as a representative value. (The difference between runs NM06 and NM06c is that the former had steady forcing and the latter had cyclic forcing).
The important point here is that dissipative, driven magnetofluids tend to have k ̲ < k M < 2 k ̲ . Defining m ̲ as the number of independent Fourier modes with smallest wavenumber k = k ̲ , the summation G ( φ ) can be broken up into the following:
G ( φ ) = m ̲ | H M | k M k ̲ k ̲ | H M | 1 H C 2 / φ 2 2 φ E + | k | k ̲ | H M | k k M k | H M | 1 H C 2 / φ 2 2 φ E
Above, the first term on the right is negative, while all the rest are positive because k > k M = φ / | H M | for k 2 k ̲ . (Even if k M > 2 k ̲ , so that there were a few more negative terms, the following development would still be valid.) In addition, all the terms in G + ( φ ) are positive. In the limit that M ,
lim M G + ( φ ) = M 1 H C 2 / φ 2 lim M G ( φ ) = M 1 H C 2 / φ 2 m ̲ φ k ̲ | H M | k ̲ | H M | 1 H C 2 / φ 2 2 φ E
Requiring σ ( φ ) = 0 is equivalent to requiring that G + ( φ ) + G ( φ ) = 0 ; from the relations given above, we see that a small number m ̲ of negative terms (the “dipole” part, corresponding to the smallest wavenumber, k = | k | = k ̲ ) must balance a very large number 2 M m ̲ of positive terms. For a cubical periodic box (or a spherically symmetric shell), m ̲ = 3 since the independent modes with k = k ̲ are k ̲ = x ^ , y ^ , z ^ .
Putting the expressions in (68) into G + ( φ ) + G ( φ ) = 0 leads to
m ̲ φ k ̲ | H M | k ̲ | H M | 1 H C 2 / φ 2 2 φ E 2 M 1 H C 2 / φ 2 .
Defining the small quantity ϵ = m ̲ / 2 M , we get the cubic equation
( 2 + ϵ ) φ 3 ( 1 + ϵ ) k ̲ | H M | + E φ 2 ϵ H C 2 φ + ( 1 + ϵ ) k ̲ | H M | H C 2 0 .
We always have φ E / 2 , but, in the non-rotating case ( Ω o = 0 ), we also have 0 H C 2 φ ( E φ ) , so that (approximately) E / 2 φ 1 2 E + k ̲ | H M | if k ̲ | H M | < E / 2 ; or k ̲ | H M | φ 1 2 E + k ̲ | H M | if k ̲ | H M | E / 2 .
However, the rotating case ( Ω o 0 ) case applies to essentially all planets (and stars) and, in this case, we have H C = 0 . Setting H C = 0 in (70) leads, to first order in ϵ = m / 2 M ,
φ o 1 2 E + k ̲ | H M | 1 4 ϵ E k ̲ | H M | .
This approximation is needed for theoretical development, but, when exactness is required, φ o is determined by numerically finding the minimum of σ ( φ ) corresponding to E and H M for a given run, as well as H C if Ω o = 0 .
From the expression (71) for φ o E M , we can also determine the expectation value of the kinetic energy, E K = E E M = E φ o , as well as of the difference E M E K = 2 E M E = 2 φ o E :
E φ o 1 2 1 + 1 2 ϵ E k ̲ | H M | .
2 φ o E k ̲ | H M | 1 2 ϵ E k ̲ | H M | .
We will now use these results, assuming H M > 0 , to show how the k = k ̲ positive magnetic helicity eigenvariable v ˜ 4 ( k ̲ ) has an energy expectation value of | v ˜ 4 ( k ̲ ) | 2 / N 3 k ̲ | H M | / m , while all of the other eigenvariables have expected energies | v ˜ n ( k ) | 2 / N 3 M 1 . This will allow us to explain the large-scale coherent magnet structures (i.e., quasi-stationary dipole fields) that spontaneously arise within a turbulent magnetofluid such as is found in the Earth’s outer core.

8. Energy Expectation Values

In a rotating frame of reference, again, H C = 0 so that β ^ = 0 , for which β ¯ 1 . Assuming H M > 0 , so that γ ^ < 1 and thus ζ k + = 0 and ζ k = 1 , (55)–(58) become:
v ˜ 1 ( k ) = u ˜ ( k ) , v ˜ 2 ( k ) = u ˜ + ( k ) , v ˜ 3 ( k ) = b ˜ ( k ) , v ˜ 4 ( k ) = b ˜ + ( k ) .
Remember that the dynamical variables u ˜ ( k , t ) and u ˜ + ( k , t ) carry negative and positive kinetic helicity, respectively, while b ˜ ( k , t ) and b ˜ + ( k , t ) carry negative and positive magnetic helicity, respectively.
In the limit that β ^ 0 , η ^ k = | γ ^ | / k and the eigenvariables are as given in (74), while the eigenvalues (60)–(61) become
λ ^ k ( 1 ) = λ ^ k ( 2 ) = α ^ , λ ^ k ( 3 ) = α ^ + | γ ^ | / k , λ ^ k ( 4 ) = α ^ | γ ^ | / k .
In the rotating case, α ^ and γ ^ are determined by putting φ = φ o from (71) into their respective expression as given in (62) with H C = 0 ; the result is
α ^ = 2 ϱ E φ o , γ ^ = 2 φ o E H M α ^ .
Using these expressions, along with (72), (73) and (75), gives us the unnormalized eigenvalues λ k ( n ) , up to leading order:
λ k ( 1 ) = λ k ( 2 ) = 4 M E k ̲ | H M | , λ k ( 3 ) = k + k ̲ k 4 M E k ̲ | H M | , k 1 ;
λ 1 ( 4 ) = m ̲ k ̲ | H M | , λ k ( 4 ) = k k ̲ k 4 M E k ̲ | H M | , k > 1 .
The eigenvariables have real (R) and imaginary (I) parts, i.e., v ˜ n ( k ) = v ˜ n R ( k ) + i v ˜ n I ( k ) , n = 1 , 2 , 3 , 4 , each have the same eigenvalue λ k ( n ) . The associated energies of the real (R) and imaginary (I) parts are
E n R , I ( k ) = | v ˜ n R , I ( k ) | 2 / N 3 = 1 2 λ k ( n ) ,
E n ( k ) = E n R ( k ) + E n I ( k ) = 1 λ k ( n ) .
As defined in (74), 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 (77) and (78) tell us that the expected energies with respect to helicity are
E K ± ( k ) = E 1 , 2 ( k ) = E k ̲ | H M | 4 M , k k ̲ ,
E M ( k ) = E 3 ( k ) = k k + k ̲ E k ̲ | H M | 4 M , k k ̲ ,
E M + ( k ) = E 4 ( k ) = k k k ̲ E k ̲ | H M | 4 M , k > k ̲ ,
E M d ( k ) = E 4 ( k ) = k ̲ | H M | m ̲ , k = k ̲ .
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, m ̲ = 3 , since all lowest-wavenumber modes are expected to have the same energy. However, for the non-rotating 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.
We have already seen how well the prediction (84) works by considering Figure 1 and Figure 2, as well as Table 1 and Table 2, as discussed in Section 3. To see the efficacy of (81), consider Figure 3, while, for (82) and (83), consider Figure 4 and Figure 5, respectively. Please note that Figure 3c represents the general behavior originally expected of all coefficients [19,20], i.e., that of a zero-mean random variable. Figure 1 and Figure 2 demonstrate that this expected ergodicity is clearly broken, while Figure 3, Figure 4 and Figure 5 show that expected behavior may be slow in occurring.

9. The Internal Dipole Magnetic Field

9.1. Periodic Box

Let us define the magnetic field associated with the n = 4 and k = k ̲ Fourier modes as b d ( x ) (d for ‘dipole’) and the remainder as h ( x ) , where
b d ( x ) = N 3 / 2 k ̲ b ˜ + ( k ̲ ) exp ( i k ̲ · x ) , h ( x ) = b ( x ) b d ( x ) .
Note that j d ( x ) = × b d ( x ) = b d ( x ) .
The expected energy contained in b d ( x ) is E M d and in the remainder h ( x ) is E M h = E M E M d , while (72) gives us E K , to leading order:
E M d = 1 N 3 k ̲ | b ˜ + ( k ̲ ) | 2 = k ̲ | H M | ,
E M h = E K = 1 2 E k ̲ | H M | .
(Again, the sum over k ̲ means that only independent modes k ̲ are included). These may be considered ‘thermodynamic relations’ for ideal, rotating MHD turbulence. Note that (86) above, in conjunction with (44) implies that, again to leading order,
| H M | = 1 N 3 k ̲ | b ˜ + ( k ̲ ) | 2 k ̲ .
Thus, essentially all of the contribution to | H M | comes from the b ˜ + ( k ̲ ) . In the periodic box, which is a 3-torus, the magnetic field is completely internal, since there is no external volume, topologically speaking.

9.2. Spherical Shell

In the spherical shell case, however, there is an external volume and the spherical analogue of (88) allows us to relate the external dipole field to the internal poloidal field and therefore finally to know something about the internal, toroidal part of the Earth’s magnetic field. In the spherical shell model of [7], which uses homogeneous b.c.s, the statistical mechanics of the ideal magneto-fluid is essentially identical to that of the Fourier periodic box case. The only differences are superficial, in that indexing by the components of k is replaced by indexing with degree l and order m of the spherical harmonic Y l m ( θ , φ ) and the number n of the zero k l n of the radial functions g ^ l ( k l n r ) at r = 1 ; in the non-dimensional model, r = 1 corresponds to the inner core radius R I C = 1220 km. (In dealing with an ellipsoidal shell, the order m of Y l m ( θ , φ ) is also needed as k l n k l m n [9]. Here, we focus only on a spherically symmetric shell.)
In a rotating spherical shell, MHD turbulence has energy E (where E = E K + E M ) and magnetic helicity H M as ideal invariants, just as in the rotating periodic box. Thus, E = E and H M = H M , where E and H M are constants. In [7], the expressions for E M and H M are given as
E M = 1 2 l , m , n | b l m n | 2 + k l n 2 | a l m n | 2 ,
H M = l , m , n a l m n * b l m n .
Above, the coefficients a l m n are associated with the poloidal part of the internal magnetic field and the coefficients b l m n are associated with the toroidal part of the internal magnetic field.
Traditional wisdom states that we cannot know the b l m n . However, if we apply our periodic box results to the spherical shell, we see this belief may need revision. First, we recall that the smallest spherical shell dimensionless wavenumber for an Earth-like model is k ̲ = k 11 = 1.8638 (the next smallest wavenumber is k 21 = 2.1497 ). The spherical shell version of (86) is
E M d = k 11 | H M | .
Following (88), we expect that essentially all of the magnetic helicity will come from the dipole part of E M and truncate the summations (89) and (90) to the l = 1 and n = 1 terms. Setting | H M | = + H M , for H M > 0 , we have
m = 1 1 b 1 m 1 k 11 a 1 m 1 2 = 0 .
Thus, we know the dipole part of the toroidal magnetic field in terms of the poloidal field: b 1 m 1 = k 11 a 1 m 1 . This leads us to make the estimates discussed earlier in Section 3.2 and which may be reviewed there.

10. Broken Ergodicity and Broken Symmetry

Here, we discuss the differences between the k = k ̲ eigenvariables v ˜ 4 ( k ̲ ) and all the other eigenvariables v ˜ n ( k ) in regard to their dynamical behavior. It will be assumed that q x = q y = q z = 1 as differential oblateness is obviated by rotation and most planets and stars, do, in fact, rotate.

10.1. Broken Ergodicity

The expectation values (81)–(84) yield rms values | v ˜ n ( k ) | r m s | v ˜ n ( k ) | 2 1 2 , so that
( a ) | v ˜ 4 ( k ̲ ) | r m s N 3 / 2 = k ̲ | H M | 3 1 / 2 , n = 4 , k = 1 ; ( b ) | v ˜ n ( k ) | r m s N 3 / 2 E k ̲ | H M | 1 / 2 2 M 1 / 2 , all others .
Again, from (74), v ˜ 1 , 2 ( k ) = u ˜ , + ( k ) and v ˜ 3 , 4 ( k ) = b ˜ , + ( k ) . In (b), the expected magnitude is just the standard deviation because the associated mean of v ˜ n ( k ) may be taken as zero. In (a), however, the expected magnitude may represent the magnitude of the mean of b ˜ + ( k ̲ ) , rather than its standard deviation, for the following reasons.
Consider the modal dynamic Equations (30) with ν = 0 and (31) with η = 0 . Furthermore, in (30), because j d ( x ) × b d ( x ) = 0 , we have S ˜ ( j d , b d ; k ) = 0 and can then use (85) to write
S ˜ ( j , b ; k ) = S ˜ ( j d , b h ; k ) + S ˜ ( j h , b d ; k ) + S ˜ ( j h , b h ; k )
Using (93a,b) as size estimates, we see that the rms values of the first two terms on the right are ∼ M 1 / 2 larger than the third term S ˜ ( j h , b h ; k ) , which is the same size as S ˜ ( u , ω ; k ) in (30).
In (31), S ˜ ( u , b ; k ) can be written as
S ˜ ( u , b ; k ) = S ˜ ( u , b d ; k ) + S ˜ ( u , b h ; k )
Again using (93a,b), we also see that the rms value of the first term on the right is ∼ M 1 / 2 larger than the second term.
Using these estimates, estimate the rms sizes of the right sides of (30) and (31) as
| S ˜ ( j , b ; k ) | r m s N 3 / 2 | S ˜ ( u , b ; k ) | r m s N 3 / 2 1 M 1 / 2 .
From these and (93), we get
d ln | v ˜ 4 ( k ̲ , t ) | r m s d t 1 M 1 / 2 , k = 1 , n = 4 ,
d ln | v ˜ n ( k , t ) | r m s d t 1 , all others .
What this implies dynamically is that, in equilibrium, the ‘dipole’ eigenvariables v ˜ 4 ( k ̲ , t ) = b ˜ + ( k ̲ , t ) , which may be large, have, on average, fluctuations in magnitude comparable in size to the other v ˜ n ( k , t ) , which are all very small. In particular, the fluctuations of these other v ˜ n ( k , t ) are of the same size as their rms magnitudes and so they behave like zero-mean random variables, as expected. However, the rms values of one or more of the b ˜ + ( k ̲ , t ) are so large compared to their fluctuations that they exhibit nonergodic behavior, i.e., they have relatively large mean values over very long times, i.e., the exhibit ‘broken ergodicity.’ This phenomenon will be made clearer in the next subsection, where we also discuss ‘broken symmetry.’

10.2. Broken Symmetry

We have seen that, dynamically, the magnitudes of a ‘dipole’ eigenvariable v ˜ 4 ( k ̲ , t ) = b ˜ + ( k ̲ , t ) , once large enough, can become effectively constant over a long time because its fluctuation in magnitude are very small. Often, one of the b ˜ + ( k ̲ , t ) , for k ̲ = x ̲ , y ̲ or z ̲ does not become as large as predicted by (83). These predictions are just average values over the ensemble and to see what is really going on we must consider the sum of the expectation values (86). As in Section 10.1, we will use q x = q y = q z = 0 so that k ̲ = 1 and k ̲ = x ^ , y ^ or z ^ . What (83) implies is that the six components of the three b ˜ + ( k ̲ , t ) define a six component vector in a 6D real space or a three component vector in a complex 3D space; for compactness, we define a vector v ˜ d and dot product | v ˜ d | 2 = v ˜ d v ˜ d in a 3D complex space:
v ˜ d = 1 N 3 / 2 b ˜ + ( x ^ ) b ˜ + ( y ^ ) b ˜ + ( z ^ ) , | v ˜ d | 2 = | H M | .
The endpoint of v ˜ d is a quasi-stationary point on the surface of the hypersphere of radius | H M | in a 6D subspace of the 8 M -D phase space Γ . The coordinates of the 6D space can undergo an arbitrary orthogonal rotation (or a unitary transformation in the complex 3D space) that puts the point of the vector on any desired point on the surface. Thus, although (93a) 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 v ˜ d pointing in any direction its 6D space that we choose, at least for the non-rotating case.
The fact that v ˜ d will exhibit broken ergodicity (i.e., it becomes quasi-stationary while undergoing only small fluctuations. However, initial conditions will also impose a particular direction on it, leading to a broken symmetry amongst the b ˜ + ( k ̲ , t ) with regard to the expectation values (83). The phenomenon of broken symmetry has been noted before [23,43], as we have also noted many times the appearance of broken ergodicity [6,22,23,39]. Here, we see how these aspects of MHD turbulence are connected.
In equilibrium, the magnitudes of the b ˜ + ( k ̲ , t ) , for k ̲ = x ̲ , y ̲ , and z ̲ are often not that different in the non-rotating ideal case but may vary appreciably, as Figure 1a shows. However, when rotation is imposed in the ideal case, the eigenfunction b ˜ + ( k ̲ , t ) with k ̲ parallel to Ω o has essentially all the dipole energy (regardless of the values of the q j ), as Figure 2b–d demonstrate. In the real case of forced, dissipative MHD turbulence, this phenomenon is also usually observed numerically, although some forcing mechanism may disrupt this process, as seen in Figure 7 of Reference [9].
However, alignment with the rotation axis seems fairly ubiquitous in numerical simulations, as well as in planets and stars. We will now discuss possible reasons for this ‘dipole alignment’.

11. Dipole Alignment

One of the mysteries of magnetism is just how the dipole moment of the Earth and other planets align themselves with their rotation axis. Ref. [44] proposed that magnetic moment and angular momentum were directly proportional, though this conjecture has not held up to scrutiny; here we have shown that the magnetic moment is essentially proportional to | H M | . [9] suggested that a difference in the ellipticities of the Earth’s inner and outer cores could lead to alignment if the smallest wavenumber was along the rotation axis.

11.1. Differential Oblateness Effects

To test the hypothesis that differential oblateness was important with a Fourier method code, wavevector components were adjusted, as detailed in Section 4.2, to k j q j k j and x-space position components to x j x j / q j , j = x , y , z , with q x q y q z = 1 . This adjustment insured that k · x was invariant, so that the Fourier transforms (13) and (14), as well as the numerical transform method, were unaffected. When the grids were stretched with Ω o = 0 , the coefficients with the smallest wavenumber noticeably increased in magnitude compared to their unstretched values, as Figure 1 shows. However, the effect of differential oblateness was completely overshadowed when rotation Ω o 0 was introduced, as is shown in Figure 2. The oblateness of the inner and outer cores is two orders of magnitude smaller than that modeled in our runs, so, for purposes of understanding the cause of dipole alignment, we need only use q x = q y = q z = 1 in the Fourier case, and keep the spherical shell model of [7] spherical.

11.2. Dynamical Effects

We know that the vector v ˜ d , whose three complex components are the b ˜ + ( k ^ ) , where k ^ = x ^ , y ^ , z ^ , and is defined in (99), has almost constant magnitude q m i n | H M | . In the Fourier model, q m i n = 1 and in the spherical shell model, q m i n = 1.8638 . Here, we work with the Fourier model. If we look at (31) and retain only those terms that contain one of the b ˜ + ( ± k ^ ) , we get
d b ˜ + ( x ^ ) d t = i N 3 / 2 x ^ × k ^ = ± y ^ ± z ^ u ˜ ( x ^ k ^ ) × b ˜ + ( k ^ ) ,
d b ˜ + ( y ^ ) d t = i N 3 / 2 y ^ × k ^ = ± z ^ ± x ^ u ˜ ( y ^ k ^ ) × b ˜ + ( k ^ ) ,
d b ˜ + ( z ^ ) d t = i N 3 / 2 z ^ × k ^ = ± x ^ ± y ^ u ˜ ( z ^ k ^ ) × b ˜ + ( k ^ ) .
Now, (21) and (33) tell us that
u ˜ ( k , t ) = u ¯ + ( k , t ) exp ( i Ω k t ) e ^ + ( k ) + u ¯ ( k , t ) exp ( i Ω k t ) e ^ ( k ) .
Therefore, for Ω o = Ω o z ^ , on the right sides of the equations given above, all of the terms containing u ˜ ( z ^ ± x ^ ) or u ˜ ( z ^ ± y ^ ) or their complex conjugates have oscillations with angular frequencies ± Ω k = ± 2 Ω o . If we watch the dynamical system over a time T > > 2 π / Ω o , then we may expect these terms to average out to zero. In that case, (100)–(102) become
d b ˜ + ( x ^ ) d t i N 3 / 2 x ^ × k ^ = ± y ^ u ˜ ( x ^ k ^ ) × b ˜ + ( k ^ ) ,
d b ˜ + ( y ^ ) d t i N 3 / 2 y ^ × k ^ = ± x ^ u ˜ ( y ^ k ^ ) × b ˜ + ( k ^ ) ,
d b ˜ + ( z ^ ) d t 0 .
The expression (106) is true once b ˜ + ( z ^ ) grows to be very large, while (104) and (105) tell us that b ˜ + ( x ^ ) and b ˜ + ( y ^ ) are still subject to fluctuations (again, we assume Ω o = Ω o z ^ ). Presumably, b ˜ + ( x ^ ) and b ˜ + ( y ^ ) ‘regress to the mean,’ which is zero because they are subject to relatively large fluctuations, while b ˜ + ( z ^ ) has relatively insignificant fluctuations to drive it from dominance. Further analysis is beyond the scope of this paper and may require a deeper mathematical analysis along the lines of [45].

12. Conclusions

Here, we have presented theoretical and computational results in magnetohydrodynamic turbulence that we feel are essential to understanding the geodynamo. Again, these results are based on a mathematical model that focuses on MHD turbulence, but ignores compressibility and thermal effects, as well as imposing model-dependent boundary conditions. Our conclusions are: (1) In the Earth’s outer core, the energy of the magnetic dipole is equal to the magnetic helicity multiplied by the dipole wavenumber. (2) The connection between magnetic helicity and dipole field in the Earth’s liquid core gives us the toroidal part of the internal dipole field. (3) Knowing this allows the estimate the mean internal dipole field strength is 3 mT. (4) Differential oblateness can cause dipole alignment, but rotation appears to be far more important. (5) Dipole alignment appears related to hydrodynamic inertial waves, but exactly how is still an open question.

Funding

This research received no external funding.

Data Availability Statement

Data is contained within the article and references cited.

Acknowledgments

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Larmor, J. How Could a Rotating Body Such as the Sun Become a Magnet? Report British of the Association for the Advancenment of Science; John Murray: London, UK, 1919; pp. 159–160. [Google Scholar]
  2. Glatzmaier, G.A.; Roberts, P.H. A three-dimensional self-consistent computer simulation of a geomagnetic field reversal. Nature 1995, 377, 203–209. [Google Scholar] [CrossRef]
  3. Glatzmaier, G.A.; Roberts, P.H. A three-dimensional convective dynamo solution with rotating and finitely conducting inner core and mantle. Phys. Earth Planet. Int. 1995, 91, 63–75. [Google Scholar] [CrossRef]
  4. Kuang, W.; Bloxham, J. An Earth-like numerical dynamo model. Nature 1997, 389, 371–374. [Google Scholar] [CrossRef]
  5. Davidson, P.A. Turbulence in Rotating and Electrically Conducting Fluids; Cambridge University Press: Cambridge, UK, 2013; p. 532. [Google Scholar]
  6. Shebalin, J.V. Broken ergodicity in magnetohydrodynamic turbulence. Geophys. Astrophys. Fluid Dyn. 2013, 107, 411–466. [Google Scholar] [CrossRef]
  7. Shebalin, J.V. Broken ergodicity, magnetic helicity, and the MHD dynamo. Geophys. Astrophys. Fluid Dyn. 2013, 107, 353–375. [Google Scholar] [CrossRef]
  8. Shebalin, J.V. Dynamo action in dissipative, forced, rotating MHD turbulence. Phys. Plasmas 2016, 23, 062318. [Google Scholar] [CrossRef]
  9. Shebalin, J.V. Magnetohydrodynamic turbulence and the geodynamo. Phys. Earth Planet. Inter. 2018, 285, 59–75. [Google Scholar] [CrossRef] [Green Version]
  10. Hughes, D.W. Mean field electrodynamics: Triumphs and tribulations. J. Plasma Phys. 2018, 84, 735840407. [Google Scholar] [CrossRef] [Green Version]
  11. Shebalin, J.V. Magnetic Helicity and the Solar Dynamo. Entropy 2019, 21, 811. [Google Scholar] [CrossRef] [Green Version]
  12. Teissier, J.-M.; Müller, W.-C. In verse transfer of magnetic helicity in supersonic magnetohydrodynamic turbulence. J. Phys. Conf. Ser. 2020, 1623, 012011. [Google Scholar] [CrossRef]
  13. Orszag, S.A.; Patterson, G.S. Numerical simulation of three-dimensional homogeneous isotropic turbulence. Phys. Rev. Lett. 1972, 28, 76–79. [Google Scholar] [CrossRef]
  14. Gazdag, J. Time-differencing schemes and transform methods. J. Comp. Phys. 1976, 20, 196–207. [Google Scholar] [CrossRef]
  15. Chandrasekhar, S.; Kendall, P.C. On Force-Free Magnetic Fields. Astrophys. J. 1957, 12, 457–460. [Google Scholar] [CrossRef]
  16. Morrish, J.S.; Moon, W.M. Velocity field and thickness of a laminar boundary layer at the core-mantle boundary. J. Geophys. Res. 1991, 96, 14569–14575. [Google Scholar] [CrossRef]
  17. Thébault, E.; Finlay, C.C.; Beggan, C.D.; Alken, P.; Aubert, J.; Barrois, O.; Bertrand, F.; Bondar, T.; Boness, A.; Brocco, L.; et al. International Geomagnetic Reference Field: The 12th generation. Earth Planets Space 2015, 67, 79. [Google Scholar] [CrossRef]
  18. Buffett, B.A. Tidal dissipation and the strength of the Earth’s internal magnetic field. Nature 2000, 468, 952–955. [Google Scholar] [CrossRef]
  19. Frisch, U.; Pouquet, A.; Leorat, J.; Mazure, A. Possibility of an inverse cascade of magnetic helicity in magnetohydrodynamic turbulence. J. Fluid Mech. 1975, 68, 769–778. [Google Scholar] [CrossRef] [Green Version]
  20. Lee, T.D. On some statistical properties of Hydrodynamical and magneto-hydrodynamical fields. Q. Appl. Math. 1952, 10, 69–74. [Google Scholar] [CrossRef] [Green Version]
  21. Palmer, R.G. Broken ergodicity. Adv. Phys. 1982, 31, 669–735. [Google Scholar] [CrossRef]
  22. Shebalin, J.V. Anisotropy in MHD Turbulence Due to a Mean Magnetic Field. Ph.D. Thesis, College of William and Mary, Williamsburg, VA, USA, 1982. [Google Scholar]
  23. Shebalin, J.V. Broken symmetry in ideal magnetohydrodynamic turbulence. Phys. Plasmas 1994, 1, 541–547. [Google Scholar] [CrossRef] [Green Version]
  24. Bardsley, O.P.; Davidson, P.A. Inertial-Alfven waves as columnar helices in planetary cores. J. Fluid Mech. 2016, 805, R2. [Google Scholar] [CrossRef] [Green Version]
  25. Ranjan, A.; Davidson, P.A. Columnar heat transport via advection induced by inertial waves. Int. J. Heat Fluid Flow 2021, 87, 108703. [Google Scholar] [CrossRef]
  26. Soward, A.M. The Earth’s Dynamo. Geophys. Astrophys. Fluid Dyn. 1991, 62, 191–209. [Google Scholar] [CrossRef]
  27. Biskamp, D. Magnetohydrodynamic Turbulence; Cambridge University Press: Cambridge, UK, 2003; Chapter 2. [Google Scholar]
  28. Patterson, G.S.; Orszag, S.A. Spectral calculation of isotropic turbulence: Efficient removal of aliasing interaction. Phys. Fluids 1971, 14, 2538–2541. [Google Scholar] [CrossRef]
  29. Favier, B.F.N.; Godeferd, F.S.; Cambon, C. On the effect of rotation on magnetohydrodynamic turbulence at high magnetic Reynolds number. Geophys. Astrophys. Fluid Dyn. 2012, 106, 89–111. [Google Scholar] [CrossRef] [Green Version]
  30. Shebalin, J.V. Plasma relaxation and the turbulent dynamo. Phys. Plasmas 2009, 16, 072301. [Google Scholar] [CrossRef]
  31. Mininni, P.D.; Montgomery, D. Magnetohydrodynamic activity inside a sphere. Phys. Fluids 2006, 18, 116602. [Google Scholar] [CrossRef] [Green Version]
  32. Mininni, P.D.; Montgomery, D.; Turner, L. Magnetohydrodynamic activity inside a sphere. New J. Phys. 2007, 9, 303–330. [Google Scholar] [CrossRef]
  33. Shebalin, J.V. Ideal homogeneous magnetohydrodynamic turbulence in the presence of rotation and a mean magnetic field. J. Plasma Phys. 2006, 72, 507–524. [Google Scholar] [CrossRef]
  34. Elsässer, W.M. Hydromagnetic dynamo theory. Rev. Mod. Phys. 1956, 28, 135–163. [Google Scholar] [CrossRef]
  35. Woltjer, L. A theorem on force-free magnetic fields. Proc. Natl. Acad. Sci. USA 1958, 44, 489–491. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Shebalin, J.V. Global invariants in ideal magnetohydrodynamic turbulence. Phys. Plasmas 2013, 20, 102305. [Google Scholar] [CrossRef]
  37. Kraichnan, R.H. Helical turbulence and absolute equilibrium. J. Fluid Mech. 1973, 59, 745–752. [Google Scholar] [CrossRef]
  38. Fyfe, D.; Montgomery, D. High beta turbulence in two-dimensional magneto-hydrodynamics. J. Plasma Phys. 1976, 16, 181–191. [Google Scholar] [CrossRef] [Green Version]
  39. Shebalin, J.V. Broken ergodicity and coherent structure in homogeneous turbulence. Physica D 1989, 37, 173–191. [Google Scholar] [CrossRef]
  40. Silvester, J.R. Determinants of block matrices. Math. Gaz. 2000, 84, 460–467. [Google Scholar] [CrossRef] [Green Version]
  41. Shebalin, J.V. The homogeneous turbulent dynamo. Phys. Plasmas 2008, 15, 022305. [Google Scholar] [CrossRef]
  42. Khinchin, A.I. Mathematical Foundations of Statistical Mechanics; Dover: New York, NY, USA, 1949; pp. 137–145. [Google Scholar]
  43. Shebalin, J.V. Broken symmetries and magnetic dynamos. Phys. Plasmas 2007, 14, 102301. [Google Scholar] [CrossRef] [Green Version]
  44. Blackett, P.M.S., XI. The magnetic field of massive rotating bodies, London, Edinburgh, and Dublin. Phil. Mag. J. Sci. 1949, 40, 125–150. [Google Scholar] [CrossRef]
  45. Cheng, B.; Mahalov, A. General results on zonation in rotating systems with a β-effect and the electromagnetic force. In Zonal Jets: Phenomenology, Genesis, and Physics; Galperin, B., Read, P.L., Eds.; Cambridge University Press: Cambridge, UK, 2019; pp. 209–219. [Google Scholar]
Figure 1. Wavenumber k 1 , positive helicity magnetic variables, normalized: b ^ + ( k ) = b ˜ + ( k ) / N 3 | H M | , for runs (a) C0, (b) S10, (c) S20, (d) S30 of Table 1. These trajectories all began close to the origin. The behavior of the b ˜ + ( k ) , k = k z z ^ in these graphs illustrates the influence of initial conditions and the presence of broken symmetry, e.g., in (a), where all the q j = 1 , the expectation that all of the b ˜ + ( k ) have the same magnitudes, which is represented by the dashed circular arc, is not met.
Figure 1. Wavenumber k 1 , positive helicity magnetic variables, normalized: b ^ + ( k ) = b ˜ + ( k ) / N 3 | H M | , for runs (a) C0, (b) S10, (c) S20, (d) S30 of Table 1. These trajectories all began close to the origin. The behavior of the b ˜ + ( k ) , k = k z z ^ in these graphs illustrates the influence of initial conditions and the presence of broken symmetry, e.g., in (a), where all the q j = 1 , the expectation that all of the b ˜ + ( k ) have the same magnitudes, which is represented by the dashed circular arc, is not met.
Fluids 06 00099 g001
Figure 2. Wavenumber k 1 , positive helicity magnetic variables, normalized: b ^ + ( k ) = b ˜ + ( k ) / N 3 | H M | , for runs (a) S10, (b) S1X, (c) S1Y, (d) S1Z of Table 1. Comparing these figures with those in Figure 1, it is clear that rotation is more important than oblation (at least for the large value Ω o = 10 ).
Figure 2. Wavenumber k 1 , positive helicity magnetic variables, normalized: b ^ + ( k ) = b ˜ + ( k ) / N 3 | H M | , for runs (a) S10, (b) S1X, (c) S1Y, (d) S1Z of Table 1. Comparing these figures with those in Figure 1, it is clear that rotation is more important than oblation (at least for the large value Ω o = 10 ).
Fluids 06 00099 g002
Figure 3. Wavenumber k 2 2 , positive helicity kinetic variables, normalized: u ^ + ( k ) = exp ( i Ω k t ) u ˜ + ( k ) / | u ˜ + ( k ) | 2 , for run S3Z of Table 1. The factor exp ( i Ω k t ) removes the cyclic behavior of inertial waves; see Equation (33). The unit circles indicate | u ^ + ( k ) | 2 = 1 , i.e., the expected variance of unity. Initial values ( t = 0 ) of these trajectories are marked with a •. These graphs indicate that the non-dipole variables behave as expected.
Figure 3. Wavenumber k 2 2 , positive helicity kinetic variables, normalized: u ^ + ( k ) = exp ( i Ω k t ) u ˜ + ( k ) / | u ˜ + ( k ) | 2 , for run S3Z of Table 1. The factor exp ( i Ω k t ) removes the cyclic behavior of inertial waves; see Equation (33). The unit circles indicate | u ^ + ( k ) | 2 = 1 , i.e., the expected variance of unity. Initial values ( t = 0 ) of these trajectories are marked with a •. These graphs indicate that the non-dipole variables behave as expected.
Fluids 06 00099 g003
Figure 4. Run S10: wavenumber k 2 2 , negative helicity magnetic variables, normalized: b ^ ( k ) = b ˜ ( k ) / | b ˜ ( k ) | 2 . Initial values ( t = 0 ) of these trajectories are marked with a • and the darker part of the trajectory shows the last half of the run. The unit circles indicate | b ^ ( k ) | 2 = 1 , i.e., the expected variance of unity. These graphs indicate that the non-dipole variables behave as expected. In fact, all coefficients were expected to be ergodic, well-behaved zero-mean variables; however, consider Figure 1 and Figure 2, which show that ergodicity is strongly broken for some of the lowest wavenumber coefficients.
Figure 4. Run S10: wavenumber k 2 2 , negative helicity magnetic variables, normalized: b ^ ( k ) = b ˜ ( k ) / | b ˜ ( k ) | 2 . Initial values ( t = 0 ) of these trajectories are marked with a • and the darker part of the trajectory shows the last half of the run. The unit circles indicate | b ^ ( k ) | 2 = 1 , i.e., the expected variance of unity. These graphs indicate that the non-dipole variables behave as expected. In fact, all coefficients were expected to be ergodic, well-behaved zero-mean variables; however, consider Figure 1 and Figure 2, which show that ergodicity is strongly broken for some of the lowest wavenumber coefficients.
Fluids 06 00099 g004
Figure 5. Run S1X: wavenumber k 2 2 , negative helicity magnetic variables, normalized: b ^ ( k ) = b ˜ ( k ) / | b ˜ ( k ) | 2 , for run S1X of Table 1. Initial ( t = 0 ) points are marked with a • and the darker part of the trajectory shows the last half of the run. The unit circles indicate | b ^ ( k ) | 2 = 1 , i.e., the expected variance of unity. These graphs show that the non-dipole variables are tending to expected behavior, though rotation slows the process down, compared to Figure 4.
Figure 5. Run S1X: wavenumber k 2 2 , negative helicity magnetic variables, normalized: b ^ ( k ) = b ˜ ( k ) / | b ˜ ( k ) | 2 , for run S1X of Table 1. Initial ( t = 0 ) points are marked with a • and the darker part of the trajectory shows the last half of the run. The unit circles indicate | b ^ ( k ) | 2 = 1 , i.e., the expected variance of unity. These graphs show that the non-dipole variables are tending to expected behavior, though rotation slows the process down, compared to Figure 4.
Fluids 06 00099 g005
Table 1. Parameters for new ideal MHD turbulence runs, which all had the same constant energy E = 1 , with less than 0.1% change during a run. The magnetic helicity H M and H C are given below and also changed less than 0.1%, except in the runs with Ω o 0 , where H C fell quickly to ∼0 and fluctuated only slightly thereabout. The runs with cubical grids ( q x = q y = q z = 1 ) begin with the letter C and those with stretched grids with the letter S. Simulation times were t = 0 to t e n d , with  Δ t = 0.001 . The angular rotation vector of each run is Ω o , while the stretching factors are q x , q y and q z for k-space (and the inverse of these for x-space). R ̲ is the average of the ratios R = E M d / k ̲ | H M | , which are themselves averages over the last 20% of each run.
Table 1. Parameters for new ideal MHD turbulence runs, which all had the same constant energy E = 1 , with less than 0.1% change during a run. The magnetic helicity H M and H C are given below and also changed less than 0.1%, except in the runs with Ω o 0 , where H C fell quickly to ∼0 and fluctuated only slightly thereabout. The runs with cubical grids ( q x = q y = q z = 1 ) begin with the letter C and those with stretched grids with the letter S. Simulation times were t = 0 to t e n d , with  Δ t = 0.001 . The angular rotation vector of each run is Ω o , while the stretching factors are q x , q y and q z for k-space (and the inverse of these for x-space). R ̲ is the average of the ratios R = E M d / k ̲ | H M | , which are themselves averages over the last 20% of each run.
Run:C0CXCYCZS10S20S30S1XS1YS1Z
t e n d :500400400400500500500400400400
q x :11110.951 1 0.95 0.950.950.95
q y :1111 1 0.95 0.951 1 0.95 1 0.95 1 0.95
q z :11111 1 0.95 0.95111
Ω o :0 10 x ^ 10 y ^ 10 z ^ 000 10 x ^ 10 y ^ 10 z ^
H M 0.12700.13070.13070.12690.12710.12660.12710.12710.12710.1271
H C 0.05900000.05890.05880.0592000
R ¯ :0.9890.9550.9670.9751.0011.0001.0460.9821.0721.017
Table 2. Forced, dissipative MHD turbulence runs of [9]. Time step was Δ t = 0.001 and simulation times were t = 0 to 2500 for the first seven, and t = 2500 to 3500 for the last three. The angular rotation vector of each run is Ω o , while the stretching factors were all q x = q y = q z = 1 . R ̲ is the average of R = E M d / k ̲ | H M | , which are themselves averages over the last 20% of each run.
Table 2. Forced, dissipative MHD turbulence runs of [9]. Time step was Δ t = 0.001 and simulation times were t = 0 to 2500 for the first seven, and t = 2500 to 3500 for the last three. The angular rotation vector of each run is Ω o , while the stretching factors were all q x = q y = q z = 1 . R ̲ is the average of R = E M d / k ̲ | H M | , which are themselves averages over the last 20% of each run.
Run:NM01NM02NM03NM04NM05NM06NM07NM02cNM06cNM08c
E a v g :1.01511.06441.04861.07821.04041.01830.99821.04271.04620.9992
H M a v g :0.62510.69470.80490.71190.65130.61790.71180.84610.84910.7079
H C a v g :0.01630.0011−0.0068−0.0020−0.0000−0.0055−0.0001−0.02020.00590.0133
Ω o :0000 10 z ^ 10 z ^ 10 z ^ 0 10 z ^ 10 z ^
R ¯ :1.0200.9970.9770.9500.9181.0590.9881.0001.0151.038
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Shebalin, J.V. Magnetic Helicity and the Geodynamo. Fluids 2021, 6, 99. https://doi.org/10.3390/fluids6030099

AMA Style

Shebalin JV. Magnetic Helicity and the Geodynamo. Fluids. 2021; 6(3):99. https://doi.org/10.3390/fluids6030099

Chicago/Turabian Style

Shebalin, John V. 2021. "Magnetic Helicity and the Geodynamo" Fluids 6, no. 3: 99. https://doi.org/10.3390/fluids6030099

APA Style

Shebalin, J. V. (2021). Magnetic Helicity and the Geodynamo. Fluids, 6(3), 99. https://doi.org/10.3390/fluids6030099

Article Metrics

Back to TopTop