Reversals in the large-scale αΩ-dynamo with memory

Inversion of the magnetic field in a large-scale model of αΩ-dynamo with α-effect with stochastic memory is under the investigation. The model allows us to reproduce the main features of the geomagnetic field reversals. It was established that the polarity intervals in the model are distributed according to the power law. Model magnetic polarity time scale is fractal. Its dimension is consistent with the dimension of the real geomagnetic polarity time scale.


Introduction
The existence of large-scale magnetic fields of planets, stars and galaxies is usually attributed to the action of the dynamo mechanism (Zeldovich et al., 1983).Magnetohydrodynamic equations are symmetric with respect to the change of sign of the magnetic field, which leads to a potential reversal in the dynamo system.These reversals are observed in cosmic dynamo systems.For example, the reversal of the magnetic field of the Sun occurs approximately every 11 years (Stix (1989).We get the information on the geomagnetic field reversals from palaeomagnetic records, on the basis of which the geomagnetic polarity timescale is constructed.The sequence of moments of geomagnetic field reversals is a nonperiodic random sequence (Merril et al., 1996).Thus, the statistical reversals of magnetic fields of the Sun and Earth are very different.However, concerning geomagnetic reversals, we mean the transition between stable states of a geomagnetic dipole, averaged over a few thousands of years (Merril et al., 1996).Therefore, the difference in the reversals of the magnetic fields of the Sun and the Earth is the difference in the processes at absolutely different timescales.
It is known that different scales of geomagnetic polarity form self-similar fractal structures (Ermushev et al., 1992;Ivanov, 1993;Pechersky et al., 1997).Intervals between the reversals (polarity intervals) differ by several orders of magnitude; there are long intervals without reversals, superchrons (Gaffin, 1989;Merril et al., 1996).A large scatter of the interval lengths does not allow us to use such characteristics as mean or variance correctly.It is known that the random variables with the properties such as self-similarity of the set of realizations, the range, and the infinity of points may be well described by power law distributions (Sornette, 2006).
Of course, one can not rely on the construction of a geodynamo model, which would fully reproduce the real palaeomagnetic scale.It is only possible to get similar statistical characteristics.Different models of geodynamo allow us to obtain random sequence reversals, the properties of which are very different.In some models the solutions are periodic or quasiperiodic (Hejda and Anufriev, 2003;Rikitake, 1965); in others they exhibit fractal properties (Anufriev and Sokoloff, 1994;Hollerbach et al., 1992).
In mean-field theory the α-effect and turbulent diffusion are usually assumed to be instantaneous in time and local in space.However, a proper description of turbulent transport involves a convolution of an integral kernels with the mean field (Moffat, 1978;Krause and Rädler, 1980).Hori and Yoshida (1992) showed that the memory effect for magnetic Reynolds number R m 4 strongly affects the dynamo action.Brandenburg (2009) used the formalism of response functions and showed that the effect of the integral kernels can be significant for anisotropic flows.
In this paper, we consider the simple model of a large-scale α -dynamo in which there is the perturbation of α-effect of a non-Markovian random pulse process.Physically, this process may be interpreted as the effect of rejected modes of mean field.According to the authors, non-Markov character of the process is very important, because it can describe the L. K. Feschenko and G. M. Vodinchar: Reversals in the large-scale α -dynamo with memory "memory" (hereditarity) in the model in the simplest form (without convolutions).
The mechanism of α -dynamo was proposed by Parker (1955).This kind of dynamo is typical for astrophysical objects (planets, stars, galactic disks) and suggests differential rotation of the object and turbulence in the character of motion of a conducting medium in this object.
The essence of such a dynamo is as follows.During the initial moment, the existence of a poloidal field of dipole type is supposed.During the differential rotation, the magnetic field lines of a highly conducting medium curl around the axis of rotation; this leads to the appearance of the toroidal field in the convective zone of a star or a planet liquid core.
To close the cycle, it is necessary to get a new poloidal field from this toroidal one.It is assumed that this is due to the breaking of mirror symmetry flows in the convection zone.Turbulent mirror-asymmetric flow generates effective the electromotive force in the direction of the toroidal field (α-effect), which leads to the excitation of a new poloidal field.The theory of α-effect was developed in Steenbek et al. (1966) and Steenbek and Krause (1969).A detailed description of the mean-field theory is given in the works of Krause and Rädler (1980), Moffat (1978), and Zeldovich et al. (1983).
Induction equation for the magnetic field in a conducting medium is the following: where v is the velocity field of the medium, and ν m is the magnetic viscosity, which is assumed constant.
If the velocity field is defined, then Eq. ( 1) is linear and defines the kinematic dynamo problem.However, the magnetic field affects the flow of a medium by the Lorentz force.The effect of this force in the equations of motion of the medium is quadratic in the magnetic field, so in the case of small magnetic fields, we can be restricted to kinematic approximation.The formal criterion of the non-applicability of the kinematic approximation is the satisfaction of the ratio E K E B , where E K and E B are the kinetic energy of the moving medium and the energy of the magnetic field, respectively.
In this case, it is necessary either to solve Eq. ( 1) together with the equations of motion or to enter a modelling approach, where v is the given functional of B. In any case the solved equations become nonlinear.
In the mean-field theory the expansion of fields v and B in the mean U and B and fluctuations u and b are introduced.We do not assume the smallness of fluctuations.Then from Eq. ( 1) we obtain the equation for the mean-field generation (Zeldovich et al., 1983): Here α and β are, in general, second-rank tensors, depending on the velocity and magnetic field.The determination of the form of these curves is the main task of mean-field theory.Convolution of αB determines the turbulent EMF (α-effect), and β B gives the diffusion of the magnetic field, which consists of molecular and turbulent diffusions.
We will further consider the isotropic case of scalar α and β; β is assumed to be the constant.
Let the B 0 and the u 0 be the characteristic value of the field and the velocity, respectively.α 0 is the maximum value of the α-effect, and L is the linear dimension.We choose for length scale, timescale, velocity scale, and field scale the value L, L 2 /β, u 0 , and B 0 respectively.Then Eq. ( 2) takes the dimensionless form where R m = Lu 0 /β is the magnetic Reynolds number, and R α = Lα 0 /β is the measure of the α-effect.

The simplest equations of the large-scale α -dynamo
In this paper we consider the simplest model of α -dynamo.The generation of the magnetic field of the planet can be described using the α 2 -dynamo when the toroidal and poloidal field are generated with α-effect.This dynamo may work in a situation where the differential rotation is weak or non-existent.It is known that α 2 -dynamo can lead to generation of the axisymmetric field with a dominant dipole structure.However, the possibility of oscillating solutions associated with inversions is problematic (Krause and Rädler, 1980).
In this dynamo the toroidal field is generated from poloidal by joined α-effect and -effect.If the intensity of the αeffect has a much smaller intensity of -effect, the role of αeffect in the generation of the toroidal field can be neglected.
For the core of the Earth, there is no single opinion about the ratio of α-effect and -effect (Glatzmaier and Roberts, 1995;Anufriev et al., 1997;Schrinner et al., 2007).However, following Ruzmaikin and Starchenko (1991), we believe that a strong dominance in the axisymmetric dipole in the magnetic field is the result of the predominant role of differential rotation.In this regard, we completely have neglected the αeffect in the generation of a toroidal field.
We believe that the spatial structure of the field is simple and can be described by one-poloidal and one-toroidal modes.Thus we consider only the largest-scale structure of the mean field.This approximation has the form and does not necessarily imply an axial symmetry.
We also assume that the mean flow U is the toroidal shear flow.In the simplest case, U is the differential rotation.Then velocity scale u 0 = LG, where G > 0 is the measure of the differential rotation.If r is the distance to the rotation axis and the is the angular velocity, then G ∼ |r∂ r |.Then R m is a measure of -effect; therefore, it is denoted R .
Then, the evolution of field is determined by the behaviour of scalar amplitudes B T (t) and B P (t) for which the following equations are valid: As mentioned earlier, we have neglected in Eq. ( 5) the αeffect.These equations can be obtained from Eqs. ( 3) and ( 4) using the Galerkin method and rescaling t, R α , and R .Here for t with a rescaling factor ∼ 1, therefore, timescale remains the diffusion.
The first part of Eq. ( 5) describes the -stage, and the second is for the α-stage cycle.
In the assumption of the constancy of R and R α , field generation (i.e. the growth of small fluctuations of B = |B T | 2 + |B P | 2 ) occurs at R α > R −1 .The field increases indefinitely at an exponential rate.If R α < R −1 , the field is attenuated.Limited-largest non-vanishing solution can occur only if R α R = 1.Thus, D = R α R is the dynamo-number.When D = 1, except for the zero steady-state solution, a lot of stationary regimes of the form B T = R B P appear in Eq. ( 5), forming a straight line in an asymptotically stable phase plane.
Limited nonvanishing solutions of Eq. ( 2) are obtained by taking into account the feedback that is the change of the turbulent flow characteristics by the magnetic field in the result of the Lorentz force.In the models of Eq. ( 5) type, this mechanism is implemented in the form of the prescribed dependence R α on B. In the simplest case, functional dependencies of the form R α = f (B(t)) are introduced.Such type models are known as algebraic quenching models and the αeffect value depends on the field current value, i.e. a response to the changes in the field of turbulence instant.The simplest version of this dependence is given in Zeldovich et al. (1983).More complex variants, based on the representation of α as the differences between the kinetic helicity and current helicity, were studied, for example, in Field and Blackman (2002) and Brandenburg and Sandin (2004).
In the Introduction we have already mentioned that the inclusion of non-locality and memory in the dynamo model can have a significant effect.We mention also the results of Frick et al. (2006).In this paper the model multi-scale dynamo was researched.The equations of mean field dynamo and the equations of shell model of MHD turbulence were integrated.In the large-scale part of the model, the authors used the α 2 -dynamo.The α-effect values were calculated by the variables of shell model.As we have already mentioned, the behaviour of large-scale field was described of two scalar time-dependent functions.Frick et al. (2006) found that simultaneous values of B and α are uncorrelated.Moreover, if the response of B on the α is fast, the inverse response occurs with a noticeable delay, and correlation decay is slow.The authors came to the conclusion that the response of R α to B has dynamic character and may not be described in terms of algebraic quenching.In our opinion, the slow decay of the correlation is an indication of memory in this model.
We introduce memory in the Eq. ( 5) model too.This can be done in two ways.In the first case, R α is not a function, but a functional of B (convolution of integral kernel and B).In the second case, R α is a function of B and a non-Markovian randomly process ξ(t).Physically, this process may be interpreted as the effect of rejected modes of fields U and B. The dependence of R α on previous values B will be implemented through the stochastic memory of the process ξ(t).These two variants of hereditarity will be further called the dynamic and random hereditarity, respectively.Of course, combination of these two types of memory is also possible.
Further, the simplest variant of the algebraic quenching will be used as the original form of the feedback where ε > 0 is the model parameter, which determines the efficiency of the feedback.A similar form of the dependence was considered in Zeldovich et al. (1983).
For the model in Eqs. ( 5) and ( 7), there are three stationary points.First and foremost is the zero point, which is unstable, which provides generation of the field.In addition, there are rest points of the form B T = ±R 1 + R 2 −1/2 , and . It is easy to show that these points are asymptotically stable.Thus, the model in Eqs. ( 5) and ( 7) gives field generation with the output to the characteristic value of B = 1.In this case, R determines the ratio of characteristic values of the toroidal and poloidal components.Therefore, during model calculations we will always assume that R = 1.
Introduction of hereditarity of dynamic type requires the following modification of Eq. ( 7): where h(t) ≥ 0 defines the "memory" of the system.The model in Eqs. ( 5) and ( 8) has the same equilibrium points as Eqs.( 5) and ( 7), and the computational experiments have shown that the nature of their stability remains the same.
The asymptotic stability of the points for this hereditarity model gives no possibility to reversals.
We return to Eq. ( 5) with the constants R α and R and consider their solutions more carefully.When R α R > 1, the solution increases indefinitely without oscillations, and the characteristic time of the increase is Supposing now that R α is a variable, we can say that negative peaks of this value are required for the occurrence of reversals, since during negative R α in the linear case oscillations appear.They must be strong enough so that the oscillation period is less than the characteristic decay time, and rare enough so that a feedback mechanism recovers the field, decreasing during the reversal.

Model with random hereditarity
Peaks in the value of R α , necessary for the formation of reversals, may be obtained by the introduction of a random hereditarity into the model in the following form: where ξ(t) is some non-Markov random pulse process with zero mean value.We define the random process ξ(t) by the following formula: Here θ k is the increasing sequence of random instants of exponential pulses, η k is random pulse amplitude, and the constant λ −1 > 0 determines the pulse width.
We assume that the time intervals between pulses τ k = θ k − θ k−1 are independent and identically distributed with the probability density function (pdf) p τ (t).Amplitudes η k are assumed to be Gaussian random variables that are independent between each other and with the times of the pulses having zero mean and variance σ 2 .
The important element of this model is the law of distribution of p τ (t).If it is exponential, the pulse sequence forms a Poisson processes of events, and the process ξ(t) turns out to be a Markov one.Any other kind of law p τ (t) leads to the fact that the waiting time of the next pulse will depend on the time from the previous pulse.Thus, ξ(t) will turn out to be a non-Markov process.
We assume that the law p τ (t) has a power asymptotic dependence ∼ 1/t γ , γ > 1.We will give a number of arguments in favour of this assumption.
Random intervals τ k may be considered as the result of the joint effect of a large number of independent factors.If we suppose the additive character of the joint effect, then, according to the generalized central limit theorem, p τ (t) should refer to the class of stable laws (Samorodnitsky and Taqqu, 1994).All such laws, except for the Gaussian one, have the power asymptotic dependence.Note that for stable power laws 1 < γ < 3.
In addition, just the power distributions have the property of self-similarity, manifesting themselves in the reversals of geomagnetic field.Finally, the power of statistics is generally characteristic of turbulent phenomena, which include αeffect.
The explicit form of the pdf for the unilateral power stable laws is unknown, except for the distribution of Levi-Smirnov (γ = 3/2).It causes difficulties in obtaining their computer implementations.
Therefore, in the calculations we used the following expression for the pdf: This form of the distribution law allows us to obtain the random variables τ k easily.
The accepted distribution coincides with the stable one only asymptotically; therefore for the distribution of polarity intervals, we will further be interested only in its asymptotics.

Simulation results
In order to analyse the effect of pulses in ξ(t) on the field, we first made some calculations in the model of Eqs. ( 5), ( 9) and (10) for a case of non-random and regular process ξ(t), which is a sequence of pulses with alternating signs of the type ±Ae −t , A ≥ 0. The interval between the pulses was 50 time units; the value ε = 0.5.The initial conditions were given as B T (0) = 0 and B P (0) = 10 −2 .The results of these calculations are shown in Fig. 1.
It is seen that positive pulses cause a sharp rise in the field, but they are not accompanied by reversals.Field response on the negative peaks depends on the magnitude of these peaks.We see that for the small pulses the poloidal component does not change the sign (A = 2); B P (t) changes the sign for a short time and returns to its original value (A = 4.2).Such behaviour of the field is well known in the palaeomagnetic data and is called excursion (Merril et al., 1996).Then there is the reversal (A = 10).During the subsequent growth, the reversal is replaced by field excursion (A = 15) again.Then excursion combination appears with the subsequent reversal (A = 30), followed by two consecutive excursion (A = 50).The trend shown in Fig. 1 continues further.For example, when A = 100 the combination of two excursions and a reversal appears.However, such sharp peaks in R α are difficult to admit in a real system.It is also clear that there are critical values of the amplitude A, separating the different types of field reversal.In particular, the critical value of A, separating the cases A = 4.2 and A = 10 from Fig. 1, is 4.455 ± 0.005.We also see that for the chosen value of ε = 0.5 the time of t B field transfer to a steady state is about 30 units.In general, as the numerical experiments showed, this dependence has a power law of t B ∼ ε −0.9 .Now consider the results of simulation in the model of Eqs. ( 5), ( 9) and ( 11) with random process ξ(t).
In the calculations, the values of the parameters R = 1 and λ = 1 were applied.Standard deviation of random amplitudes η k of pulses in Eq. ( 9) is σ = 6.6.For this value of σ the reversal in the result of negative pulse in ξ(t) occurs with the probability of 0.5.The initial values for the field components were chosen as B T (0) = 0 and B P (0) = 10 −2 .
We suppose the characteristic size of the Earth is L = 3.48 × 10 6 m (the radius of the liquid core) and the turbulent magnetic diffusion is β = 10 m 2 /s.Then our dimensionless time 5 × 10 4 corresponds to the length of the longest scale of geomagnetic polarity (Pechersky, 1997) in 1700 Myr.Therefore, calculations in the model were carried out up to t = 5 × 10 4 .Figure 2 shows an example of a segment of one of the magnetic field realizations and the toroidal and poloidal components for ε = 5.
We calculated the different values of the parameters ε and γ .For ε parameter, the values 0.1, 0.5, 1.0, and 5.0 were used; for γ , the values from 1.1 to 2.9 with 0.2 step were used.For each parameter combination, a histogram of the interval lengths of polarity and the fractal dimension of polarity model scale were estimated.First, the obtained distribution of interval lengths of polarity ζ is considered.They are illustrated in Fig. 3.
It is clear that we may speak about the power asymptotic dependence of distribution of these intervals ∼ 1/ζ δ .More specifically, the power type for ε = 0.1 and ε = 0.5 appears from ζ = 10, and for ε = 1 and ε = 5, from ζ = 30.Also there is a deviation from the power law at large ζ (more then 10 3 ).They correspond to a single event.Therefore, these deviations may be explained by insufficient data.
We have calculated the value of the index δ on the straight section of the chart shown in Fig. 3 The obtained values and the correlation coefficients corresponding to the straight sections are shown in Table 1.
According to these values, it is easy to show that, for different ε, the correlation coefficient between γ and δ is more than 0.92.It means that γ and δ are linear, the coefficients of which depend on the parameter ε.
Integrability conditions of pdf for ζ imply that δ > 1.These values are obtained from the above-mentioned linear relations for γ > 2.1.
Also note that Fig. 3 does not show the distributions for γ < 1.This is due to the fact that for such values of γ the rectilinear sections, corresponding to the power laws, do not occur on the graphs.
It may be concluded that the degree distribution of polarity interval occurs in the model at γ > 2.1.Now consider the fractal dimension of the derived polarity scales.In the calculation, we followed the procedure proposed in Pechersky et al. (1997) for real geomagnetic polarity timescale.
The technique is as follows.On the scale of T length, some interval of length is distinguished.N ( ) is the number of intervals of length on this scale, on which at least one reversal occurs.If T and the reversal are distributed uniformly, then N ( ) ∼ −1 .If T and reversal are distributed unevenly, then we may expect the dependence of the form N ( ) ∼ −d .In this case, d is the Hausdorff dimension of the scale reversals and for 0 < d < 1 the reversal series is fractal.
We made calculations in the model for the above mentioned values of γ and ε.The value of decreased in the geometric progression from 5000 ( 5 × 10 4 ) to ∼ 10.Graphs of the obtained dependencies are illustrated in  dorff dimension d.It is clear that, in all the cases, 0 < d < 1, and the reversal series is fractal, although there is a tendency to achieve the boundary of the fractal region when γ increases.

Conclusions
The geomagnetic polarity scale is characterized by the absence of standard reversal wait time and by the self-similarity on various timescales.The values of interval cover several magnitude orders; the polarity scale is chaotic (Gaffin, 1989;Ermushev et al., 1992;Ivanov, 1993;Merril et al., 1996;Pechersky et al., 1997).One of the problems of the geodynamo theory is the explanation of this phenomenon and re-L.K. Feschenko and G. M. Vodinchar: Reversals in the large-scale α -dynamo with memory producing a series with similar statistical properties in numerical simulation.
Direct simulation has shown that the structure of the magnetic field and its reversals may vary depending on the values of the relevant dimensionless numbers (Kutzner and Christensen, 2002;Busse and Simitev, 2006).It can be concluded that each inversion of the simulation has its own unique character, which can differ greatly in various aspects to the others (Coe et al., 2000).
The above-mentioned models greatly differ in total simulation time and reproduce reversals with exponentially distributed of polarity intervals.In the new preprint Wicht and Meduri (2015) described results of an extremely long simulation in different models to study the statistical behaviour of dipole variations and reversals.The distribution of polarity intervals is also obtained by the exponential.They also mentioned slower decreasing tails, but their statistical significance is very low.
Our main task was to simulate power law distribution of the polarity intervals and self-similarity polarity scale.
We use a simple model of a dynamo, where the behaviour of the toroidal and poloidal modes are described by two scalar amplitudes.This allowed for multiple simulation of the geomagnetic polarity scale for different values of the parameters.The total simulation time was always 5×10 4 times of turbulent diffusion, which corresponds ∼ 1.5 × 10 9 years.
The model uses an algebraic quenching and random perturbations in α-effect.The pulse process with a random intensity of the pulses and a random waiting time realize these perturbations.This perturbation is supposed to be interpreted as the influence of rejected modes of mean fields.
The power law was applied as the distribution law of the pulse waiting time.The reason for this was the power-law character of stable distributions, limiting distributions in the scheme of summation of independent random variables with slowly decaying pdf.
The power law distribution of waiting time of the pulse leads to the fact that perturbing processes is non-Markovian.Therefore, the model has simplest form of memory, through the memory of non-Markov perturbations.
It was discovered that the power law of polarity interval distribution is asymptotically realized in this model, if the exponent γ in the distribution of the pulse waiting time is not less than 2.1.The exponent δ in the distribution of polarity interval is linearly related to the γ .The coefficients of this linear relation depend on the effectiveness of the feedback in α-effect.
It is shown that the model scale of geomagnetic polarity is a fractal set with Hausdorff dimension of 0.7.It is consistent with the actual Hausdorff dimension of geomagnetic scale according to the paper by Pechersky et al. (1997).
Thus, it was established that the proposed model largescale dynamo allows us to reproduce the main features of the process of geomagnetic field reversals.
the solution decays without oscillations for the characteristic time of ∼ 1.If R α R < 0, the solution oscillates with the frequency √ |R α R | and decays for the characteristic time of ∼ 1.

Figure 1 .
Figure 1.The response of the poloidal component B P (t) on the regular sequence of alternating pulses with different amplitudes A.

Figure 2 .
Figure 2. The segment of the magnetic field realization for ε = 5, γ = 2.5: the toroidal component of B T , the poloidal component of B P , field value of B = |B T | 2 + |B P | 2 (log scale).

Figure 4 .
Figure 4. Number of N( ) intervals of length, which contain at least one inversion.