Hereditary Oscillator Associated with the Model of a Large-Scale αω -Dynamo

: Cosmic magnetic fields possess complex time dynamics. They are characterized by abrupt polarity changes (reversals), fluctuations of fixed polarity, bursts and attenuations. These dynamic conditions can replace each other, including both regular and chaotic components. Memory in dynamo systems manifests itself in a feedback mechanism when a strong magnetic field begins to change the properties of turbulent flows. A hereditary oscillator can be the simplest model of such complex oscillatory systems with memory. The article suggests the construction of such oscillator by means of two-mode approximation of magnetic field components in the αω -dynamo model. The hereditary member describes the suppression of a field turbulent generator by magnetic helicity and determines the shape of oscillator potential. The article describes the implicit difference scheme for numerical research of oscillator. It also describes the results of numerical simulation for two cases—instantaneous feedback and delay in feedback. The results of simulation are interpreted in terms of oscillator theory. It is shown that the observed dynamic regimes in the model go well with the change of potential shape.


Introduction
The existence of large-scale magnetic fields of planets, stars and galaxies is usually associated with the action of hydromagnetic dynamo [1]. The dynamo theory proves the process of magnetic field generation by the movements of conductive liquids and gases. The source of such movements in celestial bodies is convection in the liquid kernels of planets and convective zones of stars. The existing models of such convection reproduce hydrodynamic flows that can generate magnetic fields close to real [2,3].
The productivity of modern computing systems does not allow for the conduction of the direct numerical simulation of the three-dimensional tasks of planetary and stellar dynamo on the time scales of celestial bodies existence. In this respect, the numerical models either reproduce convective movements and magnetic fields on fairly shallow spatial grid on small time scales (∼10 5 years) or allow for the calculation of a long evolution (∼10 9 years) of only large-scale spatial structures.
The models of the second type include low-mode dynamic systems based usually on the αω-dynamo mechanism. This mechanism was suggested by Parker [4]. This kind of dynamo is typical for astrophysical objects (planets, stars, galactic disks) and assumes differential rotation of the object and the turbulent nature of conductive medium movement in this object.
The essence of such dynamo is as follows. The initial moment implies the existence of a dipole type poloidal field. With differential rotation, the magnetic field lines of well-conducting medium swirl around the axis of rotation, which leads to the appearance of toroidal field in the convective zone. To close the cycle, it is necessary to get a new poloidal field from this toroidal field. It is assumed that this is related to the mirror symmetry violation of currents in the convective zone. The turbulent mirror-asymmetric flow generates an effective electromotive force in the direction of toroidal field (α-effect), which leads to the excitation of new poloidal field. The theory of α-effect was developed in the works of [5,6].
The equations of magnetic hydrodynamics are symmetrical in relation to the change of magnetic field sign. This makes the existence of reversals possible-rapid shifts in magnetic field polarity without convection structure restructuring. These reversals are observed in real cosmic dynamo systems [7,8]. The reversals series can contain both regular and chaotic components.
For instance, the reversals of the Sun's magnetic field occur approximately every 11 years, although they contain the chaotic component [8]. The information about reversals of the geomagnetic field is obtained from paleomagnetic records, which lay the basis for construction of geomagnetic polarity time scale. The sequence of geomagnetic field reversal moments is a nonperiodic random sequence [7]. Different scales of geomagnetic polarity form self-similar fractal structures [9][10][11]. The intervals between reversals (polarity intervals) differ by several orders of magnitude; there are long intervals without reversals (superchrones) [7,12].
In real dynamo systems, oscillations are also observed around certain direction of the field (vacillations), bursts, excursions (short-term changes in polarity after which it is restored). Therefore, it is reasonable to suggest that space dynamo systems are complex oscillatory systems.
In self-consistent nonlinear models αω-dynamo, a large-scale magnetic field impacts upon turbulent generator, providing generation for finite field. In the simplest cases, this feedback is considered to be instantaneous in time and local in space. However, the correct description of turbulent transfer may include both spatial nonlocality [13] and memory effect [14][15][16]. The consideration of memory in low-mode dynamo models can affect essentially the dynamic regimes. For instance, the introduction of memory into Rickitake two-disc dynamo leads to the emergence of hyperchaotic oscillations [17]; non-Markovian fluctuations of α-effect intensity lead to power distributions of polarity intervals and realistic fractal dimension of polarity scales [18].
Thus, it can be said that real dynamo systems are complex oscillatory systems with memory effect. Therefore, it is possible to make the attempt to describe their basic properties with the help of hereditary oscillators. This work is devoted to the development and research of one oscillator of the kind.
The popular method of introduction of heredity into mathematical models is the use of fractional derivatives [19,20]. In the equations of classical models, the whole order derivatives are often replaced by fractional order derivatives [21,22]. However, the formal replacement of conventional differentiation operators in models by fractional operators is physically difficult to interpret. In this work we will not formally introduce fractional derivatives but derive the oscillator equation from the dynamo equations.
The hereditary member appears in the feedback mechanism in the form of an integral operator with the kernel of sufficiently arbitrary shape. This hereditary member qualitatively determines the shape of oscillator potential.
The work also describes the implicit difference scheme for numerical study of the integrodifferential equation of oscillator and presents some results of numerical simulations and their interpretation in terms of oscillator theory. It is shown that the dynamic modes observed in numerical simulation are well coordinated with the change of potential shape.
The work bears purely theoretical character and does not involve the comparison of model dynamic regimes with the real regimes of stars or planets dynamo.

Hereditary Two-Modes αω-Dynamo Model
We assume that some area Ω (planet liquid core or stellar convective shell) is filled with an electrically conductive medium. The generation of axis-symmetric magnetic field B(r, t) within the class of αω-dynamo models can be described by the Equations [1]: where B T (r, t) and B P (r, t) are the toroidal and poloidal components of the field B respectively, v T (r, t) is a predetermined axis-symmetric toroidal velocity field (differential rotation), η denotes the turbulent magnetic diffusivity, and the functionα determines the helicity of small-scale turbulence. It is assumed that the surrounding of the Ω is electrically nonconducting. Therefore, vacuum boundary conditions are specified for the magnetic field: We consider that the spatial structure of the field is simple and can be described by one-poloidal and one-toroidal modes. Thus, we envisage only the largest-scale structure of the field. This approximation has the following form: where the dimensionless modes b T (r) and b P (r) satisfy the boundary conditions (2). Then we substitute approximation (3) into Equation (1) and use the Galerkin method. As a result, we obtain the following equations: The Galerkin coefficients: ω-is the measure of generation toroidal mode by large-scale differential rotation, α-is the measured generation poloidal mode by helicity of small-scale turbulence, 1/η T and 1/η P -are the character times of modes dissipation. These coefficients are determined by the formulas: Equations (4) are the simplest model of a two-mode dynamo. This model is linear, so the generation of a finite magnetic field is impossible. It is also clear that generation will occur in case D = αω η T η P > 1. Exactly in this case the zero point will be unstable and the small initial values of the field will increase. Therefore, the dimensionless parameter D is a dynamo-number.
In order to ensure the stable generation of a finite field in the model it is necessary to introduce the feedback into the system (4). In the real physical dynamo system this feedback is the change of the turbulent flow characteristics by the magnetic field in the result of the Lorentz force, in particular, the change of small-scale turbulence helicity. The Lorentz force is quadratic in the magnetic field, so the feedback can be introduced into the model in the form of where α 0 -is the helicity in the absence of a strong magnetic field, and Q(·, ·)-is some quadratic form. Among the various quadratic forms the following two are of physical interest. The first is related to the energy of the field: The second is related to the magnetic helicity: In this paper we shall observe the second case. It is related to the fact that this form has an alternating sign and a wide variety of dynamic regimes can be expected. If we use the sign-positive form of energy, the suppression of generation could be too strong.
In the simplest case, the functional dependencies are introduced, where Q(·, ·)-is the function of two scalar variables. Such type models are known as algebraic quenching models and the helicity value depends on the field current value, i.e., a response to the changes in the field of turbulence is instant. The different variants of algebraic quenching were studied, for example, in [1,23,24].
In the introduction we have already mentioned that the inclusion of nonlocality and memory in the dynamo model can have a significant effect. We also mention the results of [25]. This paper studies the model multiscale dynamo. It integrates the equations of mean field dynamo and the equations of shell model of magnetohydrodynamic turbulence. The authors of this work found that the simultaneous values of field and helicity are uncorrelated. Moreover, if the response of field on the helicity is fast, the inverse response occurs with a noticeable delay, and correlation decay is slow. The authors [25] came to the conclusion that the response of helicity to changes of the field 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 our model too. Let the quadratic form Q be not a function but the t-parametric functional of B T B P : where K(s) ≥ 0, s ∈ R-is some dimensionless kernel, with the property K(+∞) = 0, t F -is the time scale of kernel, B 0 -some characteristic value of field. This predetermined expression specifies the model of hereditary quenching. Obviously, if we change the kernel K by a constant factor c > 0, this is equivalent to changing B 0 by a factor √ c. Therefore, in what follows, we always consider the kernel K in a normalized suitable way, for example sup K = 1.
From Equations (4)-(6) we obtain the model of a two-mode αω-dynamo with memory: The model is closed by the initial conditions B T (0) = B T 0 , B P (0) = B P 0 . For planetary and stellar dynamo systems, it is reasonable to assume that B T 0 = 0. It is related to the fact that a small external field, which is poloidal, is required to start the dynamo system at the initial moment [1]. Now let us make the system (7) dimensionless. We will use the diffusion time of the poloidal field 1/η P as the time scale, and B 0 √ α 0 t F as the magnetic field scale. Let us now turn to new dimensionless variables (keep the notation for time) and to the new dimensionless parameters: These parameters: D-is the dynamo-number; σ −1 -is the dimensionless time of toroidal field decay; p-is the dimensionless time scale of the kernel. Then the dimensionless model of a two-modes dynamo takes the form: It is known that if we take one toroidal and one poloidal mode of free decay of a magnetic field with the same spatial scales, the poloidal mode has an eigenvalue smaller [26]. Therefore, in the future, we consider that σ > 1. For the classic Parker's dynamo [4], for example, the work [26] results give the values σ ≈ 3.37. We shall also study the case of a «working dynamo», i.e., D > 1.

Hereditary Oscillator
Now, we transform the Equations (8) to the form of an oscillator with memory. Let us differentiate with respect to time t the first equation of (8) and use the second equation. As a result, we obtain We substitute from the first equation of (8) into the (9). Then Further, from the first and third equations of (8) we obtain We apply integration by parts in the first integral in (11): We substitute this expression in (11) and obtain Now we introduce the new (σ, p)-parametric kernel J σp (s) using equality: Using (12) and (13) we can write Equation (10) in the following form: So, the two-mode dynamo model (8) with the initial conditions x(0) = x 0 , y(0) = y 0 are equivalent to the integrodifferential Equation (14) with the initial conditions x(0) = x 0 , x (0) = σ(y 0 − x 0 ). We have already mentioned above that from the point of view of dynamo theory, it is reasonable to assume that the initial conditions for a toroidal field are zero, i.e., x 0 = 0.
The condition x 0 = 0 does not impose strong restrictions. Since K(+∞) = 0, the term with x 0 in (14) affects only the transient process. In the steady state, we can assume that x 0 = 0. Therefore, further we always assume that x 0 = 0.
The Equation (14) can be considered (formally) as a dissipative oscillator with hereditary potential: This point of view on the two-mode dynamo model allows for the determination of the possible dynamic regimes without a numerical solution of the integrodifferential equation. However, for the detailed study of dynamics, numerical modeling is necessary. For this, it is necessary to compose difference schemes and develop software.
It should be noted that for some classes of kernels K(s) and J σp (s), integrodifferential systems (8) and (14) can be reduced to "purely" differential systems with additional initial conditions. The Appendix proves the following. If the kernel K(s) or J σp (s) is a solution to the homogeneous linear differential equation with constant coefficients of the n-th order, then the integrodifferential system can be reduced to the differential system of order n + 2.
We especially note the special case of the kernel K(s) = e −s , i.e., J σp (s) = e −s σ − 1 2p . Using the scheme described in the Appendix A, we obtain that the (8) is equivalent to This is the classical Lorenz system [27], the dynamics of which at σ = 10 and p = 3/8 is well studied [28,29]. Note that the Lorenz system as the simplest model of the Solar αω-dynamo was proposed in the work [30] and as the model of disk dynamo was used in the work [31].
This makes it possible to verify the numerical schemes and programs for numerical simulation of the dynamic of the oscillator (14).

Difference Scheme and Software for Numeric Simulation
To construct a difference scheme, we write the oscillator in the form of the following system: In numerical simulation we shall use nonlocal initial conditions for x(t) in the form x(≤ 0) = 0. We introduce the time step h and the corresponding values of the phase variables x n , v n and w n .
We also denote J σp nh p through J n .
For the integral part of the system (16) (third equation), we will use the Simpson formula: where for even n: and for odd n: For the differential part of the system, we use the implicit second-order Runge-Kutta method [32]: It is convenient to write (20) in the following form, taking into account that X(v) = v: Now we shall combine the schemes for integral and differential parts. We substitute in the second Equation (21) the expression for w n+1 from (17). Then we substitute the expression for x n+1 from the first Equation (21) into the second equation. As a result, we obtain the scheme: Then arises the following algorithm for the numerical solution of the oscillator equation on uniform grid t n , n = 0, 1, . . . , N: 1. We set the initial conditions x 0 = 0, v 0 = x (0), w 0 = 0. 2. Let at n-th time step, the values x n , v n , w n , J n are known.
4. We calculate L n by Formulas (18) or (19). 5. We find v n+1 from the first Equation (22). This is a cubic equation with respect to v n+1 . We solve it by Newton's method, using v n as the starting value for v n+1 . 6. We calculate x n+1 (second Equation (22)). 7. We calculate w n+1 (third Equation (22)). 8. Increase the time index n by 1 and go to step 2.
This algorithm was programmed in Java code. To verify the code, the numerical simulation of the oscillator was carried out in case of the kernel K(s) = e −s , i.e., J σp (s) = e −s σ − 1 2p , and parameter values σ = 10, p = 3/8. Earlier we noted that such choice of the kernel and parameters corresponds to the classical case of the Lorenz system. The simulation results are shown in Figure 1. Three characteristic regimes for the Lorenz system are shown: asymptotically stationary regime (Figure 1a), chaotic regime (Figure 1c,d) and periodic regime (Figure 1b). It can be concluded that the difference scheme and the program code are working correctly.

Simulation Results
Now let us look at some results of numerical simulation of the oscillator (14). They will be interpreted in terms of the general theory of oscillators. In numerical simulation we will use the value σ = 3.37, which is typical for the Parker's dynamo. The case of the power kernel K(s) with the asymptotic 1/s α , 0 < α < 1 will be considered. The kernels of this kind are of a special interest in terms of hereditary systems. In case the kernel is integrated on the interval [0; +∞), an effective memory time scale can be entered. Then the memory in the system is finite and the system will not be "truly" hereditary.
The shape of the potential graph (15) is determined by the signs of the coefficients A and B. If the coefficient A is fixed, then the coefficient B is determined by the averaged current and previous values of the phase variable x(t). Therefore, it could be expected that the phase variable will respond to changes in the coefficient B with delay or omission of rapid changes.
We note that for the energy E of the oscillator Consequently, in case of sharp decrease in the value of the coefficient B, i.e., dB dt 0, the energy of the oscillator may increase. An important thing is that it is not related somehow to the sign of the coefficient B.
The value K(0) determines the degree of potential U. In case of K(0) = 0, we can speak about instantaneous feedback. The potential in this case is of the 4th degree. In case of K(0) = 0, we can speak about feedback delay and quadratic potential. Further, let us look at the examples for these two cases.

Instantaneous Feedback: K(0) = 0
In this case A = K(0)/8 > 0. Figure 2 schematically presents the potential graph. For small values of t the coefficient equals to B ≈ σ(D − 1) > 0 and the movement occurs in the potential with two wells (Figure 2a).
In this case for s > max {0, 1 − α/ (2pσ)} the kernel J σp (s) > 0, and for the coefficient B both positive and negative values are possible. It means that either potential forms shown in Figure 2 are possible. If in case of B > 0 the phase variable x(t) is at the bottom of the potential wells I and II (Figure 2a), it corresponds to the asymptotically stable generation of a fixed-sign field. However, in case of B < 0 the variable x(t) cannot be at the bottom of potential well III (Figure 2b). This would correspond to a zero solution but for D > 1 zero point is unstable. Two cases are possible here. In the first case B remains negative but decreases rapidly. Then the energy of the oscillator increases and the variable x(t) rises from the bottom of well III. In the second case, B changes its sign to positive and x(t) leaves the unstable zero point. Both case combinations are certainly possible.
Numerical simulations show that for p = 0.1, the solutions are asymptotic stationary even for D = 1000. We can assume that for a kernel small-time scale, the solutions will be only asymptotically stationary. It is interesting to note in relation to this fact that even in the Lorenz system the nonzero rest points are stable for any values of D if p < 1/(σ − 1). It is probable that for the small kernel scale p and instantaneous feedback only asymptotically stationary solutions exist.
Then let us look at the numerical simulation results shown in Figure 3a. It demonstrates damped oscillations in potential well II in case of the same potential shape with two potential wells. Figure 3b shows the chaotic regime accompanied by damped oscillations in the potential well. This is the regime of metastable chaos. It can be seen that the change of the sign x(t) goes with the sharp transition of B into negative band. In this case the potential has one minimum and then returns to the shape with two potential wells. It proves the change of x(t) sign. However, sometimes the negative excursus B happens to be too short and x(t) does not have time to change the sign. This is probably how the memory effect in the oscillator manifests itself.  Figure 3c shows the chaotic regime with the change of x(t) sign. Here B keeps positive value almost all the time, i.e., the potential has two wells and one maximum. It is seen that sharp decrease of B is seen in case the oscillator receives energy keeps x(t) in the extreme position. The energy reserve is sufficient to move from the maximum level to the value with a different sign. Consequently, the fluctuations with a change of sign occur above the maximum potential level. Figure 3d shows the burst regime with the transition to stationary regime. Here it could also be seen that the change of the sign x(t) is associated with rather long stay of B below zero level. Figure 3e shows the burst regime. The mechanism of x(t) sign change is similar to the previous one. Figure 3f shows the chaotic regime similar to the transition regime in Figure 3b.
In general, we can say that the dynamics of the phase variable x(t) are well interpreted in terms of the theory of oscillators.

Feedback Delay: K(0) = 0
In this case A = K(0)/8 = 0. Figure 4 shows schematically the potential shape. We take the kernel K(s) as K(s) = s/(1 + s) 1+α . For numerical simulation we used the value As in the previous case (25) shows that for 2pσ < α kernel J σp (s) < 0 for small s > 0 and positive for large s. If 2pσ > α, then J σp (s) > 0 for all s ≥ 0. In any case for s > max {0, 1 − α/ (2pσ)} kernel J σp (s) > 0 and for the coefficient B both positive and negative values are possible.
First, consider the simulation results shown in Figure 5a,b. Figure 5a shows the damped oscillations near stationary nonzero value. Figure 5b shows the chaotic regime with variable x(t) sign changes and vacillations around two symmetrical levels. We consider the coordination of these regimes with the shape of potential. For different signs of the coefficient B, the branches of the parabola are directed up or down. Figuratively, we can say that the potential flaps its wings. When the wing is raised x(t) tends to zero. When the wing descends x(t), on the contrary, goes away from zero. Therefore, if the direction of parabola branches changes fast enough, then x(t) fluctuates around some average nonzero value. This is how vacillations occur. It can be seen that the sign of x(t) changes when moving from B < 0 to B > 0, i.e., when the potential well turns into its peak. Then the variable x(t), sliding to zero, may have time to slide into the band of another sign. This is how reversals occur. The asymptotically stationary regime (Figure 5a) corresponds to the stabilization of the coefficient B at zero level, i.e., the horizontal position of the branches.
The scheme for changing the sign of x(t) described above also describes the mode shown in Figure 5c. Here large values of the dynamo-number D and the sharp decrease of B provide high energy. This leads to the fact that x(t) always moves to the area of another sign and no vacillations occur.
Figures 5d-f show the burst regimes. When B falls sharply below zero level, i.e., the branches of the potential go up sharply, x(t) quickly turns out to be near zero. Then the branches of the potential slowly descend. When they go down hard enough x(t) will quickly move away from zero in a random direction. The following sharp rise of branches starts the process to be repeated. This is how these observed regimes can be explained.
In general, we can say that the dynamics of the phase variable x(t) are well interpreted in terms of the theory of oscillators again.

Discussion
In this article: 1. A two-mode dynamo model with hereditary quenching of the alpha-effect by magnetic helicity has been obtained. A rigorous mathematical derivation of the model is carried out. The general idea of building such a model was previously described in the work of the author [33]. 2. The model is converte to the form of an oscillator with a potential hereditary term. The phase variable of the oscillator is the amplitude of the toroidal component of the magnetic field. The potential of the oscillator is determined, among other things, by the current and previous values of the square of the phase variable.
3. It is proved that if the kernel of the hereditary term is a solution of a linear homogeneous differential equation with constant coefficients, then the integrodifferential equation of the oscillator can be written as a differential system. 4. A difference scheme for the numerical simulation of the oscillator is obtained. The correctness of the circuit and verification of the program code were carried out for the special case of the oscillator, when it is equivalent to the Lorenz system. 5. Numerical simulation of dynamic regimes for two power nonintegrable kernels is carried out.
It is shown that the calculation results are well interpreted from the point of view of the theory of oscillators.
In the author's opinion, the description of the dynamo system as an oscillator will make it possible to better comprehend the processes in this system from a general physical point of view. In addition, it provides a way to identify possible dynamic modes without performing numerical simulations. It is possible that the class of oscillators described in the work can be applied to other physical systems.
As already mentioned in the introduction, this work bears purely theoretical character and does not involve the comparison of model dynamic regimes with the real regimes of stars or planets dynamo. For such complex systems it is necessary to compare not the solutions and real regimes themselves but their statistical characteristics-distributions of polarity intervals, fractal dimensions of polarity scales, Hurst exponents, etc. This implies the pursuance of a separate study whose content differs essentially from the described study. The conduct of such work in the future could evidently be advisable. Then, we multiply the k-th Equation (A2) by a k p k and sum over all k = 0 . . . n. We obtain the following differential equation a n p n d n w dt n + a n−1 p n−1 d n−1 dt n−1 + · · · + a 1 p dw dt The initial conditions for the w(t) and its derivatives are determined from Equations (A2) by the substitution t = 0.
So, integral relation is formally equivalent to differential Equation (A3) with suitable initial conditions. Therefore, the integrodifferential Equation (14) can be specified using the (n + 2)-th order differential system. Obviously, similar transformations can be performed for the last equation of a two-mode dynamo system (8) if the kernel K(s) satisfies the Equation (A1).