Possible Role of Non-Stationarity of Magnetohydrodynamic Turbulence in Understanding of Geomagnetic Excursions

: The natural simplifying assumptions often put forward in the theoretical investigations of the magnetohydrodynamic turbulence are that the turbulent ﬂow is statistically isotropic, homogeneous and stationary. Of course, the natural turbulence in the planetary interiors, such as the liquid core of the Earth is neither, which has important consequences for the dynamics of the planetary magnetic ﬁelds generated via the hydromagnetic dynamo mechanism operating in the interiors of the planets. Here we concentrate on the relaxation of the assumption of statistical stationarity of the turbulent ﬂow and study the effect of turbulent wave ﬁelds in the Earth’s core, which induces non-stationarity, on the turbulent resistivity in the non-reﬂectionally symmetric ﬂow and the geodynamo effect. It is shown that the electromotive force, including the so-called α -effect and the turbulent magnetic diffusivity ¯ η , induced by non-stationary turbulence, evolves slowly in time. However, the turbulent ¯ α coefﬁcient, responsible for the dynamo action and ¯ η evolve differently in time, thus creating periods of enhanced and suppressed turbulent diffusion and dynamo action somewhat independently. In particular, periods of enhanced ¯ α may coincide with periods of suppressed diffusion, leading to a stable and strong ﬁeld period. On the other hand, it is shown that when enhanced diffusion occurs simultaneously with suppression of the α -effect, this leads to a sharp drop in the intensity of the large-scale ﬁeld, corresponding to a geomagnetic excursion.


Introduction
The terrestrial magnetic field is generated by the hydromagnetic dynamo action in the Earth's liquid core driven thermally and compositionally [1][2][3]. The turbulent mechanism of dynamo action, which is typically invoked as responsible for the generation of large-scale fields, is the so-called α-effect, based on nonlinear interactions of the small-scale fluctuating components of the non-reflectionally symmetric turbulent state, which generate the largescale electromotive force (EMF). For dynamo action to occur, the wave field must exhibit chirality, i.e., lack of reflectional symmetry [4][5][6]. The turbulent flow in the core of the Earth contains a rich wave field composed of nonlinearly interacting magnetohydrodynamic waves such as, e.g., the inertial waves [7], the so-called MAC waves [8,9] or magnetic Rossby waves [10]; for a more complete review see [11] on the dynamics of the core. These interactions could have a profound effect on the dynamics of the core, in particular the dynamo process as they practically rule out statistical stationarity of the core turbulence-a feature often invoked in theoretical investigations of magnetohydrodynamic turbulence, as it greatly simplifies the mathematical approach.
In fact, non-stationarity has been recently shown by [12][13][14][15], as an effective mechanism of generation of the large-scale electromotive force through interactions of waves with distinct but close frequencies of oscillations (the beating effect). Such an EMF slowly evolves on timescales comparable or larger than the typical time scales of variation of the large-scale field. It is shown here that statistically non-stationary, helical turbulence induces a slowly varying in time α-effect and turbulent magnetic diffusion, and their variations are out of phase.
Some well-known features of the evolution of the geomagnetic field are the so-called geomagnetic excursions and polarity reversals, the former being the periods of a significant drop in the large-scale field's intensity. In fact, from the point of view of the possible impact that such a phenomenon could have on the high-tech human civilization, it is the decrease in the intensity of the geomagnetic dipole, which is crucial, since it is the dipolar field that greatly prevents the solar wind from entering the atmosphere. Whether or not the actual polarity reversal takes place is of much less importance and hence in here, for short, under the term 'excursion', we will contain both the actual excursions and reversals. The aim of this paper is to comment on the phenomenon of geomagnetic excursions in light of the non-stationarity of the core turbulence induced by wave-interactions and the slow time evolution of the α-effect and turbulent diffusion. This evolution is shown to create periods of enhanced diffusion, which may be correlated with suppression of the α-effect, thus creating favourable conditions for the excursions.
Numerical simulations do not reveal any significant large-scale alterations in the flow of the conducting liquid during a reversal [16][17][18], although it must be said that the currently available computer power does not allow to reach the strongly constrained Earth-like parameter regime. The first stage of a polarity reversal in simulations by [18] is associated with local intensification of turbulence, where the vigorous flow twists and bends the field lines to locally reverse the magnetic field direction. The enhanced turbulence intensifies the reversed field, which spreads into surrounding regions until the polarity in the core becomes mixed and the dipole moment is weakened. This suggests that the magnetic excursions and polarity reversals are manifestations of a chaotic turbulent behaviour of the liquid core, which is in qualitative agreement with the temporal variations of turbulent transport coefficients conjectured here.

Mathematical Formulation of the Mean Field Dynamo Problem
Evolution of the large-scale magnetic fields induced by the complex flow of an incompressible conducting fluid is governed by the following dynamical equations where the velocity field of the fluid flow is denoted by U(t, x), the magnetic field by B(t, x) and the total pressure Without loss of generality, we may assume solenoidal forcing and for the purpose of simplicity, we rescale the magnetic field in the following way so that the prefactor 1/µ 0 ρ in the Lorentz-force term in the Navier-Stokes (1a) equation is lost. Next, denoting by angular brackets the ensemble mean, · − ensemble mean let us introduce the standard Reynolds decomposition and write down separately the equations for the mean fields U and B and the turbulent fluctuations u and b; this yields is the large-scale electromotive force (EMF) and Furthermore, we adopt the "first-order smoothing approximation" [6] in which squares and products of fluctuating quantities are ignored Introducing taking the Fourier transforms of (8a)-(8c) and eliminating the pressure from the Fourier transformed velocity equation with the use of the projection operator one obtains where we have already used (11c) in projecting the velocity equation on the plane perpendicular to the wave vector k. Next we rearrange the terms to put all the terms involving gradients of means on the right hand sidê where and assume scale separation between the means and the fluctuations, so that the gradients of means can be assumed small and treated in a perturbational manner. The large scale EMF, by the use of (6) and (12b), can be expressed in the following way where we have used a short notationû j =û j (ω , k ). On substituting once again forb p from (12b) we obtain at the leading order where higher-order terms in the gradients G and G have been neglected.

Simple Kinematic Theory for Homogeneous Turbulence
It is a common practice to assume that the background small scale turbulence has statistical properties that are given beforehand, and thus the situation becomes kinematic. For example, we may assume that the turbulence is homogeneous and isotropic and thus the turbulence correlation tensor takes the general form where the H(ω, ω , k) is responsible for the lack of reflectional symmetry required for the large-scale dynamo process [6]. Since the mean velocity U only creates a shift of the frequency ω → ω − k · U , we simply absorb it into the frequency. The electromotive force takes the form where we have used (15), (16), (A8a) and (A8b) in Appendix B and we have also introduced a short notation q = (ω, k), Next we consider the case of non-stationary turbulence and assume the following where the 'non-stationarity' function with Γ = const, is chosen in such a way so that the correlations in the real space have a simple sinusoidal time dependence Let us set for simplicity U = 0, though it is straightforward to include it, and concentrate on the α-effect only. Consequently, for the EMF in the real space, we obtain where we have made use of the fact that for any even function of ω, say f e (ω), the integral ∞ −∞ ω f e (ω)dω = 0 vanishes and the integrals I j , for j = 1, 2, 3 are given in the Appendix A. In other words, the full EMF, including the molecular diffusion, takes the following form Clearly, there is a significant phase shift between theᾱ coefficient and the turbulent diffusionη.
The excursions and in particular polarity reversals have been observed in numerical simulations, although none of the numerical models were able to reach the very demanding, extreme parameter regime of the core. It was reported, however, that increasing the vigour of convection tends to increase the frequency of polarity reversals [16,17]. In simulations of [18], the first stage of a polarity reversal is characterised by intensification of turbulence in some region of the core where the vigorous flow twists and bends the field lines to locally reverse the magnetic field direction. Since the magnitude of the turbulent flow in that region is particularly vigorous, so is the local magnetic Reynolds number, thus favouring conditions for amplification of the reversed field, which spreads into surrounding regions until the polarity in the core becomes mixed. This implies a sharp decay of the dipole moment. Consequently, we may expect that the magnetic excursions and polarity reversals are manifestations of a chaotic turbulent behaviour of the liquid core and no significant large-scale phenomena seem to take place during excursions.
The real core turbulence is obviously bound to be non-stationary, including propagation of well-known waves such as the inertial waves [7], MAC-waves [8,9] or the magnetic Rossby waves [10], likely forming wave packets. These wave fields interact nonlinearly to form the turbulent electromotive force (EMF), and interactions of waves with distinct but close frequencies were shown to be effective in creation of the EMF by [12][13][14]. For simplicity, we will neglect higher-order interactions and include two-wave interactions only, so that in the simplest case when only two waves with distinct frequencies but the same wave vector are present (homogeneity but non-stationarity), the turbulent point correlations at time t can be written as c.c. means the complex conjugate, and we have assumed that the means of quantities oscillating with frequencies ω 1 , ω 2 and ω 1 + ω 2 , which correspond to mean field variations on the timescales from days to decades can be assumed to vanish (without loss of generality, we may assume ω 1 > 0, ω 2 > 0). However, since the core perturbations often possess close frequencies of oscillation, the frequency differenceω can be small and even smaller than the typical frequencies of oscillations of the large-scale field; in such a case, the first two terms in the second line of (25) cannot be neglected and lead to long time variations of the turbulence correlation tensor. This is in keeping with the chosen form of the 'nonstationarity' function (20), which leads to correlations in the form (21), corresponding to the form based on two-wave interactions.
In reality even such two-wave interactions lead to more than one (say N)ω-mode, because the turbulent wave field consists of a number of distinct waves so that more generally where the phase shift ψ n accounts for both the sine-type and cosine-type time dependence of the correlation tensor.

Dynamic Theory
Instead of assuming some given form of the turbulence correlation tensor, we can use Equation (12a) to express it in the following form where higher-order terms in the gradients G and G have been neglected, and we have used the short notation q = (ω, k). The EMF can then be calculated with the use of (15) and (28) and on assuming the force correlations in a statistically homogeneous form we have derived the formula for the mean EMF in Appendix B. This involves the mean velocity gradients, which within the considered model correspond to effects such as the cross-helicity dynamo or the shear-current effect [19][20][21]. Although it is straightforward to include them in the calculations, for clarity, we drop all the mean velocity terms (we set U = 0) and concentrate only on the α-effect and the effects of turbulent diffusion. For simplicity, we will also consider the case when ∇ B 2 × B is negligible, hence the mean EMF is reduced to (cf. Appendix B) where the integrals I 1 , I 2 and I 3 are provided in (A10a)-(A10c) in Appendix B. Note that lack of reflectional symmetry is introduced here by the term proportional to F 1 (ω, ω , k) in the force corrections (29), hence it is crucial for calculation of the α-effect, which is not generated by reflectionally symmetric flows [6]. In natural systems, the reflectional symmetry is broken by the presence of background rotation, which introduces the Coriolis force into the dynamics. However, rapidly rotating flows are naturally anisotropic; thus for simplicity, we have chosen an isotropic model.

Stationary Turbulence
First let us explicitly demonstrate the simplest result obtained for stationary turbulence, when i.e., the force correlations are given by In such a case, the integrals in (30) can be calculated to yield where X = cos θ and θ is the polar angle in the spherical coordinates (k, θ, φ) in the wave-vector space, hence the induction equation takes the following form where we have made use of the fact that for any even function of ω, say f e (ω) the integral ∞ −∞ ω f e (ω)dω = 0 vanishes.

Non-Stationary Turbulence
The dynamics of the Earth's core is strongly influenced by the wave field, composed of the fundamental modes such as the inertial waves, the MAC waves and the magnetic Rossby waves. These waves are described by different dispersion relations, which imply the existence of a large number of waves with distinct frequencies of oscillation. As argued at the end of Section 3, their interactions lead to a time-dependent turbulent correlation tensor. We consider the simplest case of only two-wave interactions and introduce a simple model by postulating that the turbulence is forced by the following isotropic, homogeneous but non-stationary forcing (cf. (20)) i.e., In the above, the function ∆(ω, ω ;ω, Γ) is the 'non-stationarity' function given in (20), so that indeed the two-wave interaction form is preserved. Of course, equally well, one could choose this to be a cosine-type time dependence. The parameterω has the physical interpretation of frequency shift between distinct modes, i.e., MAC or magnetic Rossby waves. The calculation of the mean EMF in the non-stationary turbulence is postponed until Appendix C. Although the limit is not geophysically relevant because of very small viscosity of the core, to make analytical progress, we have consideredω νk 2 , ηk 2 , where E is the Ekman number based on the length scale corresponding to most energetic turbulent eddies in the core. The induction equation then reads where, again,η ns is the full effective magnetic diffusivity, including molecular effects and where the coefficientsᾱ st andη st are given in (35a) and (35b), respectively, whereas formulae forᾱ ∆ω andη ∆ω can be found in (A24), (A26), (A27a) and (A27b) in Appendix C. The latter expressions can be rearranged to yield where

The Mean EMF in Non-Stationary, Low-Pm Turbulence
Let us introduce the Hartmann number based on characteristic fluctuational length scale M and a new variable , We can evaluate the integrals expressing the turbulent diffusivity and the α-effect in non-stationary turbulence (cf. Appendix C) in the geophysically relevant limit of a strong magnetic field however, on top of this assumption we will also need in order to make analytical progress. This means we assume large Hartmann numbers (strong field), but low magnetic Prandtl numbers Pm = ν/η 1. Since in the Earth's core the Hartmann number based on the core depth L is of the order 10 7 , and Pm ≈ 5 × 10 −7 , the latter assumption PmM 2 (k) 1 for the Hartmann number M (k) based on fluctuational wavelengths is not necessarily satisfied. Still, it is reasonable at least in a large part of the short-wavelength spectrum, and thus, it is put forward in order to simplify the calculations. Note that since we consider the strong field limit, it is not applicable to a linear dynamo regime, i.e., quenching of the electromotive force is significant. Let us take the following simple statistical model (cf. e.g., [22,23], etc.) with D 0 = const > 0, D 1 = const ≤ k min D 0 [6], and introduce explicitly the scale separation, i.e., we introduce the scale of the largest/most energetic turbulent eddies = 2π/k min , where k min is the smallest fluctuational wave number. This allows to calculate the relations between coefficientsᾱ st ,η st andᾱ ∆ω ,η ∆ω (see Appendix C for general formulae), which yieldsᾱ and we recall here the assumptions that led to the final form of the transport coefficients Note, however, that these assumptions were only necessary to clearly demonstrate how the turbulent diffusivity and the α coefficient can be calculated in non-stationary turbulence, and in particular, we have shown that the significant phase shift in the slow time dependence between the turbulent diffusivity and the α coefficient is a rather general feature of non-stationary turbulence. Under the current assumptions (51) the phase shifts (44) can be calculated to yield where is the Hartmann number based on the maximal fluctuational lengthscale, which in the Earth's core can be expected to be at the order of 10 4 or even 10 5 . It follows that φ α ≈ 0 is negligibly small and the phase shift between the α-effect and the turbulent magnetic diffusivity is entirely determined by φ η , which is at least a few orders of magnitude larger We emphasise that φ η is proportional to the frequency shiftω between the slowly evolving MAC/magnetic Rossby waves, which is not uniquely determined, and in fact, more realistically, the non-stationary turbulence should be modelled with a superposition of at least a few (say N) distinct valuesω n , n = 1, . . . , N, each generating a different phase shift φ η n between the α-effect and the turbulent diffusivity.
Toy Model-Energy of a Force-Free Mode at U = 0 The aim of the simple calculation shown in this section is to demonstrate how a phase shift between the field-amplifying α-effect and the resistive decay can lead to magnetic field excursions, through an explicit numerical solution obtained for a force-free mode defined by where κ(x) is a scalar function of position; this ensures B × (∇ × B ) = 0. Such forcefree states are known to exist and have been intensively investigated, e.g., in seminal works of [24,25] (cf. also more recent works of [26,27]). On defining the magnetic energy the induction equation in the absence of large-scale flows, U = 0, yields for the energy Based on the previous analysis we may propose the followinḡ where instead of the magnetic field dependence obtained in the strong field limit with Pm 1 considered in previous sections, we added a standard quenching factor, which models action of the Lorentz force in the simplest way (cf. [28,29]), withĀ = const. Taking gives a well-suited case study, since it implies that the maximal enhancement of the field takes place when the effect of the resistive decay is the weakest and vice versa-the strongest resistive decay of the mean field is associated with the weakest amplification by the α-effect; the latter situation leads to possible excursions (field decay), whereas the former one to 'stable' periods with a strong magnetic field. In such a case, we obtain where we have inserted Γ = 1. A selected solution of the latter equation is shown on Figure 1, where indeed regular short-lived excursions are manifested in the evolution of the magnetic energy, i.e., short periods with a suppressed magnetic field. More generally, accounting for the possibility of a more complex time dependence of the force correlations with more than one slowω-mode, the magnetic energy evolution equation can be rewritten in the form (61) Selected solutions are shown in Figure 2 for the cases of two and threeω-modes. Especially in the latter case, it is evident that longer periods of relative stability of the magnetic energy are separated by much shorter excursions.

Discussion of Relevance of the Results to the Problem of Geomagnetic Excursions
The presented toy model is obviously a great simplification. Firstly, it considers forcefree modes, but it also relies on derivations obtained under restrictive assumptions of statistic homogeneity and isotropy. The toy model, however, is only used to visualise the studied effect of non-stationarity of the turbulent coefficientsᾱ andη, which is generic and independent of the simplifying assumptions.
It needs to be stressed that in the presented analysis, the large-scale EMF is calculated with the use of fluctuational equations, but then we concentrate on the evolution of the large-scale field, leaving the dynamics of the small-scale component of the magnetic field unexplored. It follows that temporal enhancement of turbulent magnetic diffusivity does not need to suppress the small-scale dynamo, which in reality is actually rather expected to be vividly operating on much shorter timescales, largely independently of the large-scale process. Furthermore, as stated above, the problem is greatly simplified from the start by the assumption of isotropy and neglection of the strongly anisotropic effect of the background rotation and density stratification. At first, a rough estimate of the anisotropy leads to a separation between horizontal and vertical turbulentᾱ andη coefficients. It is also known that in the weak seed field limit (linear stage of growth of the mean magnetic field) the mean induction equation with anisotropicᾱ ij =ᾱ h δ ij + (ᾱ v −ᾱ h )δ i3 δ j3 andη ij = η h δ ij + (η v −η h )δ i3 δ j3 (say rotation is along the z-axis) separates the evolution of even and odd modes, thus dipolar parity evolves independently with the quadrupolar one. Moreover, one should bear in mind that the reality is even more complex-not only non-stationary but also inhomogeneous, leading to spatially dependent turbulent coefficients, which further differentiates the evolution of different spatial modes. Summarising, the small-scale dynamo is by no means excluded during the periods of enhanced turbulent diffusion, and the simple isotropic, homogeneous model obviously does not grasp the complexity of the strongly anisotropic Earth's core turbulence and the mean-field evolution; and it is the anisotropy and inhomogeneity (rotation, stratification, buoyancy, inhomogeneous wave interactions, etc.) that lead to significant differences in the evolution of the mean dipolar and quadrupolar fields. Although theoretical models, including non-stationarity, inhomogeneity and anisotropy, may be too cumbersome to investigate, it could be of interest to utilise DNS to study the long time evolution of turbulent coefficients in spherical, rapidly rotating shells and its dynamical connection to magnetic excursions.
Furthermore, there is a number of other important models, which identify and study other effects as underlying mechanisms of the geomagnetic reversals/excursions and present a more general approach to the problem. The dynamics of reversals has been modelled with sets of nonlinear ODE's involving supercritical bifurcations and describing the evolution of field modes with different parities, e.g., by [30] or [31]. Later, ref. [32] attempted to grasp the main features of evolution of the planetary field by ODE models with a stochastic noise. A saddle-node-bifurcation model capturing the dynamical features of the geomagnetic field evolution and those of the field obtained in the VKS experiment (cf. [33]) was developed by [34]. A group of nonlinear, one-dimensional evolutional models, which also reported Earth-like features were represented by [35][36][37]. In short, the described models point to nonlinearities causing chaotic behaviour of the system as the mechanism of generation of reversals/excursions in the system evolution. Of course this does not exclude non-stationarity of the mean-field coefficients as a trigger of the excursion events.

Conclusions
Although obviously no decisive conclusion about predictions of geomagnetic excursions can be made from the presented analysis, the paper points to the non-stationarity of the turbulent transport coefficientsᾱ andη, which are likely to slowly evolve in time due to the non-stationarity of the entire background turbulence in the core. Their evolution was shown to be out of phase, implying the existence of periods of enhanced and suppressed turbulent dynamo process, which may correlate or not with periods of enhanced/suppressed diffusion. In particular, enhancement of the α-effect may coincide with suppression of diffusion (stable field) and vice versa, enhancement of diffusion may coincide with suppression of the dynamo effect (field decay). It would be instructive to study such temporal characteristics of the turbulent diffusivity and the α-effect in numerical simulations, in particular the correlations between the value of the ratioη/ᾱ and the occurrences of excursions.
The magnetic diffusion in the Earth's core is thought to be η = 2 m 2 /s (cf. e.g., [38]), which implies the magnetic diffusion time of about 10 4 years. The non-stationary turbulence can, however, enhance the effective magnetic diffusion of the large-scale dipole to make it a few times larger, which decreases the magnetic diffusion time to a value of a few thousand years, which is in line with the realistic time scales of the geomagnetic excursions/reversals. As shown, in a non-stationary turbulence, such an enhancement of diffusion can be correlated with simultaneous suppression of the α-effect, in which case it is likely to result in a sharp drop of the geomagnetic field intensity, i.e., the excursion. With such a picture, the geomagnetic excursions are manifestations of the chaotic core turbulence and their occurrences are also bound to be chaotic, and no characteristic, significant alterations in the core large-scale flow are expected during excursions, in accordance with results of numerical simulations [16][17][18]. Perhaps with the development of the available computational abilities and power, in the future, it may become possible to numerically study the Earth-like parameter regime and construct reliable fits to the geomagnetic data in order to provide estimates of the turbulent magnetic diffusivity and the α-effect, which could then be monitored since a drop in the value of the ratioᾱ/η for the outer core could indicate a higher possibility for a geomagnetic excursion.
Funding: The support of the National Science Centre of Poland (grant no. 2017/26/E/ST3/00554) is gratefully acknowledged. This work was partially financed as a part of the statutory activity from the subvention of the Ministry of Science and Higher Education in Poland.

Conflicts of Interest:
The author declares no conflict of interest.
Next, using the following formulae whereΩ denotes the solid angle and the spherical coordinates (k, θ, φ) have been used (with a substitution X = cos θ) one obtains with where g 6 (X) = X 2 1 − X 2 , (A11c) The mean velocity gradients within the considered model correspond to effects such as the cross-helicity dynamo or the shear-current effect [19][20][21], but for clarity, we will exclude them here. Hence, we now set U = 0 and concentrate only on the α-effect and the effects of turbulent diffusion. For simplicity, we will also consider the case when ∇ B 2 × B is negligible, hence the mean EMF is reduced to ljk d 4 qe i(k·x−ωt) d 4 q e i(k ·x−ω t) û j (ω, k)b k ω , k =I 1 ν, η, B 2 B l − I 2 ν, η, B 2 + B 2 I 3 ν, η, B 2 [∇ × B ] l . (A12)