Global fluctuations in magnetohydrodynamic dynamos

The spectrum of temporal fluctuations of total magnetic energy for several dynamo models is different from white noise at frequencies smaller than the inverse of the turnover time of the underlying turbulent velocity field. Examples for this phenomenon are known from previous work and we add in this paper simulations of the G.O. Roberts dynamo and of convectively driven dynamos in rotating spherical shells. The appearance of colored noise in the magnetic energy is explained by simple phenomenological models. The Kolmogorov theory of turbulence is used to predict the spectrum of kinetic and magnetic energy fluctuations in the inertial range.


Introduction
Several experiments have been carried out in the last decade in liquid sodium at high magnetic Reynolds numbers and in highly turbulent flows. Measurements of magnetic field fluctuations, either due to an externally imposed magnetic field or due to magnetic field generated through a dynamo effect by the sodium flow itself, generally reveal spectra with a 1/f -noise at low frequency [1,2,3]. Numerical simulations of magnetohydrodynamic (MHD) flows showed the same phenomenon [4,5] although sometimes the precise form of the low frequency noise varies [6]. The precise value of the exponent in the algebraic decay in a spectrum is not so important, but a behavior other than white noise, i.e. colored noise, always invites closer inspection [7]. For example, a 1/f -noise, frequently called flicker noise, cannot extend down to zero frequency if the power integral is to stay finite. The low frequency cutoff to the 1/f -spectrum should carry some informations about the dynamics of an MHD dynamo.
Measurements of velocity in a fixed point of a turbulent flow frequently find white noise at low frequency (see for example the data compilation in fig. 6.14 of [8]), but there are also examples of colored noise [9]. The simultaneous presence of colored noise in both magnetic and velocity spectra looks like a plausible combination: If flow velocities increase in a dynamo, the magnetic field is amplified more rapidly and the amplitude of the magnetic field increases. If the velocity has no white noise, so should the magnetic field have a nontrivial spectrum. In the Karlsruhe experiment on the other hand [1] the 1/f -fluctuations in magnetic signals had no correspondence in pump pressures or volumetric flow rates. This prompts us to investigate in more detail the mechanisms leading to colored noise in magnetic spectra. This paper will deal with the fluctuations of total magnetic energy, which is a global measure of the field amplitude. A theory for the fluctuations of a global quantity is of interest even though laboratory experiments usually report local measurements. Numerical simulations always compute total magnetic and kinetic energies. Their statistics are less noisy than those of field amplitudes at a given point because of the spatial averaging inherent to the computation of total energies. Spectra of the geomagnetic field are frequently plotted as variations of another global quantity, the dipole moment [10], because the geomagnetic field is dominated by the dipole component and low frequency variations reflect variations of the dipole moment. Secular variations of the Earth's magnetic field again exhibit colored noise [11,12]. In astrophysics, we do have local measurements of solar magnetic fields, but we only know a global amplitude for most stellar magnetic fields.
The goal of this paper is to explore conditions leading to colored noise in the total magnetic energy and to show that this type of noise can appear even if the kinetic energy has white noise. Three simplified models of MHD flows are presented in section 2. The first two model the slow fluctuations per se, the last one by contrast computes the spectrum of inertial range fluctuations in homogeneous and isotropic turbulence. Section 3 presents numerical simulations of the G.O. Roberts dynamo [13] and compares the results with section 2, whereas section 4 does the same for simulations of convection driven dynamos in rotating spherical shells.

Phenomenological models
The basic problem in this paper is to find the temporal spectrum of a quantity derived from a magnetic field B(r, t), such as the magnetic energy. The induction equation governs the evolution of B in a liquid conductor with magnetic diffusivity λ whose movement is given by the velocity field v(r, t): Exact solutions of this equation are difficult to obtain, so that we resort to phenomenological models.

A single magnetic mode
Mean field magnetohydrodynamics have proven a most fruitful simplification of the induction equation [14]. In this approach, the effect of small scale fluctuations on the large scales are not computed exactly but are modeled, in the simplest case as an α−effect. The number of magnetic degrees of freedom which need to be retained is thus reduced and in an extreme simplification, only one mode remains. If we call B the amplitude of that mode,α(t) and β the coefficients describing the α−effect and its quenching, respectively, and µ a coefficient related to the magnetic dissipation, the simplest model reproducing the main features of the induction equation is: α is allowed to be time dependent in order to reflect a time dependent velocity field. This time dependence will be essentially random for a turbulent velocity field. The reduction of the α−effect by the term βB 3 models the retroaction of the magnetic field on the velocity field via the Lorentz force (which is quadratic in the magnetic field) in the Navier-Stokes equation. We now consider α(t) =α(t) − µ to be a random process with mean square α 2 and remove the dimensions from eq. (2) by expressing time in multiples of α 2 −1 and the magnetic field amplitude in multiples of ( α 2 /β) 1/2 . The adimensional quantities t ′ , α ′ and B ′ are given by t ′ = t α 2 , α ′ = α/ α 2 and B ′ = B β/ α 2 . In the remainder of this section, all quantities are understood to be nondimensional and the primes are omitted for convenience. The nondimensional variables then obey the equation: in which α(t) is a random variable with α 2 = 1. As long as B is small, the solution to this equation is For times small enough so that the exponent can be considered small, one has Taking the Fourier transform of this equation, it follows that the spectrum of B is, apart from frequency independent prefactors, the same as the spectrum of α divided by the square of the angular frequency, ω 2 . For example, if the spectrum of α is a white noise, the spectrum of B behaves as ω −2 . For large times t, eqs. (4) and (3) become a poor approximation, which means that the ω −2 will not be observable below some cutoff-frequency. If the mean of α, α , is different from zero, B will be large enough for the nonlinear term in eq. (3) to become dominant after a time on the order of α −1 .
In that regime, and concentrating on slow fluctuations, eq. (3) reduces to B 2 = α.
Considering again the example of α(t) with a white noise, one finds a spectrum of B which is a white noise, too. The transition in the spectrum of B from ω 0 to ω −2 occurs at a frequency which increases with increasing α , because eq. (5) fails at earlier times t.
In order to test these ideas, we solved eq. (3) numerically. The random α(t) was generated by sending the output of a gaussian deviate random number generator through a Butterworth filter [15,16]. The filter was adjusted such that its output had a spectrum as a function of frequency f = ω/(2π) in 1/(1 + (f /f 1 ) 4 ) with f 1 = 50. The spectrum of B is shown in fig. 1 for different α . As expected, this spectrum decays as ω −6 for ω > ω 1 = 2πf 1 , and as ω −2 for ω c < ω < ω 1 . Below the cut-off ω c , the spectrum is independent of ω, and ω c ∝ α .
Even though the spectra in fig. 1 differ only in the value of ω c , the qualitative appearence of the underlying time series is quite variable: The time series consists of intermittent bursts for small α and of random fluctuations around a well defined mean for large α . A study of the probability distribution function of the solutions of eq. (2) together with a few graphs of representative time series is presented in [17].

Several magnetic modes
We will now investigate under which conditions the single mode model is applicable to more general systems and extend the discussion of the previous section to include several modes. It will be shown that the predictions of the single mode model are recovered in the limit of small fluctuations.
The precise form of the dynamical system used as a model does not matter much for the following analysis, but a specific system has to be chosen for the numerical examples. In order to stay as close as possible to the previous section, let us assume ∇ · v = 0 and rewrite the left hand side of eq. (1) as We then proceed through the same steps as before and replace the combination of velocity and derivation by a random variable in which we absorb the dissipative term and the right hand side, remove dimensions, and model saturation through a cubic term. This leads to the following system: in which the α i (t) are random variables. This system bears only a metaphorical relation with the original induction equation and will be used to exemplify two limiting cases: If the fluctuations of the α i (t) are small compared with the mean of the α i (t), the solution of (6) will be close to the solution of the time independent system in which each α i (t) in (6) is replaced by its mean α i . Let us assume α 1 = α 2 = α 3 < 0. An eigenvalue analysis of the left hand side of (6) then reveals one neutral mode and two modes with equal and positive growth rate. In the presence of small fluctuations, the neutral mode will not contribute significantly to the dynamics. If the initial conditions and the nonlinear term select an arbitrary direction in the space spanned by the two degenerate growing modes, we expect (6) to behave the same as the single mode model. This is a fortiori true for a dynamical system with a single non-degenerate growing mode.
If on the other hand the fluctuations of the α i (t) are large compared with their means, the dynamics is not dominated by a single mode anymore and the analysis of the previous section breaks down. The spectrum of the fluctuations of B 2 i can now be different and must be found from numerical computation. Fig. 2 shows some examples of solutions of eq. (6) in which the spectrum of the fluctuations of the α i (t) is in 1/f [18] and α 2 i = 1. Let us first consider the case in which the fluctuations of the α i are small compared with their mean. The spectrum of B 2 i shown in fig. 2 should then behave as predicted by the single mode model: At the smallest frequencies, the spectrum of B 2 i must decay the same as the spectrum of the α i , i.e. as 1/f in the present example. Above a frequency on the order of α i , there must be a factor f 2 between the power laws followed by the spectra of B 2 i and the α i , which implies a spectrum in 1/f 3 for B 2 i in the example considered here. All these predictions fit well the spectrum shown in fig. 2 for α 2 i = 1 and α i = −5. In order to further support the applicability of the single mode model, we also computed the angle between the instantaneous vector B(t) = (B 1 (t), B 2 (t), B 3 (t)) and its mean B . The cosine of that angle, B · B(t)/ | B | 2 |B(t)| 2 in the statistically stationary state stays larger than 0.99 for α i = −5 but is scattered over a large interval for α i = −0.01. In the latter case the fluctuations of the α i are large compared with their mean and different exponents unrelated to the single mode model become possible. A decay in 1/f 2 appears in fig. 2 which will become important again in connection with the convectively driven dynamos in spherical shells discussed below.

Homogeneous and isotropic turbulence
The main focus of this paper is on slow fluctuations, but it is also worthwhile to have a look at the fast fluctuations in order to see where the differences are. We also would like to check whether it was reasonable in section 2.1 to assume a spectrum for α which is flat at low frequencies and steeper at large frequencies. We will specifically look at fluctuations of total kinetic and magnetic energy. The spectra of these fluctuations must be different from flow to flow, but we can expect a unique behavior in flows to which the Kolmogorov phenomenology (K41 in short after [19]) of homogeneous and isotropic turbulence applies. We will employ this phenomenology in a cartesian domain with periodic boundaries, so that it is useful to expand all fields in Fourier series. Much effort has been spent in the past to theoretically deduce the wavenumber dependence of the Fourier coefficients, the best known theory being of course K41. But the Fourier coefficients also depend on time, or after an additional transformation, on frequency. Theories such as K41 are concerned with the temporal mean of those coefficients. The spectrum of temporal fluctuations of those Fourier coefficients has received much less attention.
This section will present an extension of the Kolmogorov phenomenology which predicts the spectrum of fluctuations of the total kinetic energy of a flow. There is no magnetic field involved, but the result is directly applicable to magnetic field spectra under certain assumptions.
The mean kinetic energy E kin of a turbulent flow is decomposed into contributions made by wavevectors of modulus k in the form E kin = ∞ 0 E k dk. According to K41, E k ∝ k −5/3 . The number of modes in a thin shell of radius k in spectral space is proportional to k 2 , so that the contribution of a single mode with wavevector k to the total mean kinetic energy is proportional to k −11/3 . We now consider the root mean square of the fluctuations of the total kinetic energy: σ tot = (E 2 kin − E kin 2 ) 1/2 . The angular brackets denote average over time. The total root mean square is again decomposed into contributions of different wavevector shells: σ tot = ∞ 0 σ k dk. In order to compute σ k , we assume that the fluctuations of the k 2 modes in the wavevector shell of radius k are all uncorrelated so that the mean squares of the fluctuations of each single mode simply add to give σ 2 k . We next invoke the concept of self-similarity, which is a central tenant of the K41 theory: All structures (or eddies or modes) in the inertial range are statistically indistinguishable from each other after a rescaling of length and time. This implies that the histogram of the fluctuations in every mode has the same shape if the fluctuations are given as multiples of the mean amplitude. We then have to conclude that the root mean square of the fluctuations scales the same as the mean amplitude, i.e. is proportional to k −11/3 . It follows that σ k ∝ √ k 2 k −11/3 = k −8/3 . Note that it is not possible to deduce the scaling of σ k from the usual dimensional arguments of K41. The ratio σ k / E k tends to zero if the number of modes in a wavevector shell tends to infinity because the fluctuations of the modes in that shell average out to zero. The rms of the total fluctuation depends therefore on the number of degrees of freedom, which in turn depends on the ratio of integral to dissipative length scale. These two length scales are assumed to be irrelevant in the K41 theory, however.
The energy of the modes in the wavevector shell of radius k is now written in the form in which h k (t) has to obey h k = 0 and h 2 k = 1. The Fourier transform of the fluctuations of the total energy is We now use again the hypothesis of self-similarity and note that the typical time scale for a mode with wavenumber k is its turn-over time which is proportional to k −2/3 [20,8]. Restricting ourselves to frequencies ω and wavenumbers k in the inertial range, we expect with an unknown but universal function g. The form of the argument reflects that allĥ k (ω) should have the same form once ω is expressed in multiples of the turn-over frequency, and the amplitude factor k −1/3 follows from the requirement that h 2 k = 1, which implies that |ĥ k (ω)| 2 dω is a constant independent of k. Substituting z = ωk −2/3 into |ĥ k (ω)| 2 dω = k −2/3 |g(ωk −2/3 | 2 )dω shows that this is indeed fulfilled.
We can now proceed to finally compute the spectrum of the energy fluctuations, using the same substitution z = ωk −2/3 , to find: The spectrum of the energy fluctuations, i.e. the square of the Fourier spectrum above, must thus decay as ω −6 for ω in the inertial range. This prediction will be verified below.
Spectra in MHD turbulence can have a variety of shapes, depending on the presence of Alfvén waves, an applied magnetic field, the magnetic Prandtl number etc. [21,22] The above calculation can be directly reproduced for the magnetic energy as long as the magnetic energy density follows a Kolmogorov spectrum in k −5/3 . This happens in homogenous and isotropic flows at high magnetic Reynolds numbers without applied external field [23]. In that case, magnetic field fluctuations should decay as ω −6 , too.

The G.O. Roberts dynamo
Dynamos based on periodic flows first investigated by G.O. Roberts [13] have been useful for a number of basic studies of the dynamo effect [24,25,26,27,28] and have inspired the Karlsruhe dynamo experiment [29]. In order to verify the ideas developed in the previous section, we wish to numerically simulate the Navier-Stokes and induction equations in the following non-dimensional form: Re and Pm are two control parameters standing for Reynolds and magnetic Prandtl numbers, respectively. The forcing F will take the form One solution to the above equations is then B = 0, u = u 0 with This solution is unstable for large enough Re, either because the flow is hydrodynamically unstable or because the dynamo effect leads to B = 0. The velocity field (16) is a cartesian arrangement of helical eddies with their axes along z. Periodic boundary conditions will be applied in all cartesian directions x, y, z with periodicity lengths l x = 1, l y = 1, and three different choices of l z , namely l z = 1, 1.5 and 2.
If we call in dimensional units L x the periodicity length in x−direction, the laminar velocity v 0 = V 0 u 0 , and the kinematic viscosity and magnetic diffusivity ν and λ, respectively, the nondimensional parameters are given by Pm = ν/λ and Re = V 0 L x /ν.
The above problem is conveniently discretized with a spectral method. However, we need very long time series if we are interested in small frequencies and the availability of computer hardware for long runs becomes an important consideration in the choice of the numerical method. The most readily available platform happened to be a set of machines equipped with CUDA (compute unified device architecture) capable graphic processing units. These architectures are not well adapted to spectral methods, nor to Poisson solvers. That is why a fully explicit finite difference scheme was implemented to solve the following set of equation: These equations describe a compressible fluid with the equation of state p = c 2 ρ and a linearized continuity equation. Solutions of the incompressible Navier-Stokes equation will be recovered in the limit of the sound speed c tending to infinity. The sound speed was chosen large enough so that the peak Mach number never exceeded 0.12 in the computations presented below. The typical Mach number was more like 0.05. The finite difference code was validated against a spectral code [30] by comparing the results of short runs, but all the heavy computation was done with the finite difference code. We will first look at non-magnetic flows. Figure 3 shows spectra of kinetic energy density E kin for l z = 2, 1.5 and 1. E kin is defined as where V is the computational volume. In all cases, the spectra of E kin follow an ω −6 over at least a decade in ω. This power law identifies the inertial range of a Kolmogorov cascade according to section 2.3 (It was not checked directly whether the spatial spectrum of the velocity fields follows the K41 law because the simulations were not done with a spectral code and the evaluation of the spatial spectra would have been cumbersome). An even steeper decay follows at higher frequencies before the spectrum dives under the noise level introduced by roundoff error. On the low frequency side, the inertial range ends in a conspicuous hump. The frequency of the local maximum in the spectrum, f max , scales with what one may identify   Table 1. l z , Re, Re 2 and Re 2 f max for the simulations in fig. 3.
as the energy injection scale: Re is a control parameter in eq. (14), but the Reynolds number determined a posteriori is Re 2 = √ E kin . The values of Re 2 given in table 1 vary from 340 to 2100 indicating that all flows are turbulent. The product Re 2 f max varies from 27 to 41. Despite this variation, we identify f max with the injection scale because the definition of Re 2 does not take into account that the flow is anisotropic so that we do not expect Re 2 f max to be strictly constant.
The spectra at frequencies smaller than the injection frequency are different for l z = 1 and l z = 2 or 1.5: For l z = 1, the spectrum is purely white noise, whereas a section of 1/f -noise is visible for l z = 2 and 1.5. The section of 1/f -noise shrinks with increasing Re.
We next turn to dynamos. For Pm slightly above the critical value given in table 1, the spectra of the magnetic energy density E B , defined as E B = 1 2V B 2 dV , show a segment of spectrum in ω −3 for l z = 2 and 1.5 and in ω −2 for l z = 1. A white spectrum is always found at the smallest frequencies ( figure 4). This matches the prediction of section 2.1 in as far as the velocity spectrum without magnetic field contains an interval with a spectrum close to ω −1 for l z = 2 and 1.5 and nothing but white noise below the injection frequency for l z = 1. According to section 2.1, there should be a factor ω 2 between magnetic and velocity spectra at low frequencies, which is compatible with the results in figures 3 and 4.
The transition from white noise to either ω −2 or ω −3 occurs at a higher frequency for larger Pm in figure 4. This is again in agreement with section 2.1 because a larger Pm at constant Re corresponds to a flow farther above onset, which corresponds to a larger α in the model of section 2.1.
Experiments use local probes to characterize the magnetic field. In the present simulations, the magnitude |B| of the field at position x = y = 0, z = l z /2 was recorded as a function of time. This location corresponds to a position at which the magnetic field has been measured in the Karlsruhe experiment [1]. The spectra of the magnetic field amplitude in one point show a power law at low frequencies over a larger interval than the spectra of E B , but with an exponent which bears no obvious relation with the exponent of the energy spectra ( figure 5). There does not seem to be any theoretical tool for predicting local fluctuations, so that we simply note the spectra in figure 5 as an empirical fact. However, they also approach a white noise spectrum as Pm increases, just as the magnetic energy spectra do in fig. 4.

Dynamos in spherical shells
Simulations of convectively driven dynamos in spherical shells reported in [6] yielded spectra with power laws which are not compatible with the single mode model. According to the phenomenology developed in section 2.2, this simply means that the fluctuations in these dynamos are large. It is expected that at sufficiently low Rayleigh numbers and in sufficiently quiet convection, one finds again spectra in agreement with the single mode model. In order to test this idea, we also performed simulations of dynamos in spherical shells. Because these computations are more expensive, the spectra are more noisy and not as extended in frequency as for the G.O. Roberts dynamos, but they demonstrate the main effect.
We consider the same physical model as in [6] of a rotating spherical shell with its gap filled with fluid driven by convection. The variable quantifying buoyancy will be called temperature here, but the employed boundary conditions better fit compositional convection: Fixed temperature at the inner boundary and zero heat flux through the outer boundary. Both boundaries are assumed no slip and electrically insulating. The following dimensionless equations ∇ · u = 0 , ∇ · B = 0 (25) are solved in a spherical shell for the velocity, magnetic and temperature fields u, B and T . The control parameters are the Ekman number E, the Rayleigh number Ra, the Prandtl number Pr, and the magnetic Prandtl number Pm. The ratio of inner and outer radii is fixed at 0.35 with the dimensionless outer radius being r 0 = 1/0.65, and ǫ, the variable modeling the buoyancy source [6] is fixed at ǫ = 1. The equations have been solved using the spectral method described in [31] with a resolution of 33 Chebychev polynomials in radius and spherical harmonics of degree up to 64. Fig. 6 shows two calculations in which all parameters are held constant except for the Rayleigh number. In going from the low to the high Rayleigh number, the spectrum of kinetic energy changes from white noise to approximately f −1 . The spectrum of magnetic energy shows a decay in f −2 in both cases. This means that at low Ra, the exponents in the power laws for the kinetic and magnetic energies differ by 2, just as they should according to the single mode model. At higher Ra, i.e. in the flow with larger fluctuations, the exponents are different as expected from section 2.2 and differ by 1. Incidentally, the combination of f −1 and f −2 for kinetic and magnetic energies is identical with a combination of exponents found in the low dimensional dynamic system of section 2.2.

Conclusion
Colored noise has been noted in fluctuations of magnetic field strength in both experiments and simulations. According to common experience with turbulent velocity fields, the fluctuations of mechanical quantities have on the contrary a white noise at low frequencies, a low frequency being a frequency smaller than the inverse of the integral time scale. It has been shown in the present paper that even the simplest phenomenology based on a single mode model predicts a spectrum in 1/f 2 for the total magnetic energy if the velocity field is characterized by white noise. The single mode model is justified as long as the fluctuations in the velocity field are small enough in amplitude. These simple models are corroborated by numerical solutions of the dynamo equations for 2D periodic dynamos and for convectively driven dynamos in spherical shells.
Colored noise must have a low frequency cut-off for the power integral to converge. This low frequency is independent of diffusive time scales according to the models studied here. Instead, it depends on how strongly a dynamo is driven. The more supercritical a dynamo is, the larger is the cut-off frequency. It is also possible to deduce the spectrum of turbulent magnetic fluctuations from the Kolmogorov phenomenology. Within the inertial range, fluctuations of the total kinetic and magnetic energy should decay as f −6 , which is verified by numerical simulation.
It is left to future work to deal with fluctuations of the magnetic field at a given point in space. If the magnetic field is dominated by a large scale component, local fluctuations are of course similar to global fluctuations. At large magnetic Reynolds numbers, small scale fluctuations exceed in magnitude the fluctuations of the large scale field and the spectrum of the fluctuations of local magnetic fields can be different from the spectrum of total magnetic energy. The spectra computed for the 2D periodic dynamos at high magnetic Reynolds numbers reveal power laws at low frequencies in the range f −1.3 to f −0.5 .