# Chaotic dynamics of magnetic fields generated by thermomagnetic instability in a nonuniformly rotating electrically conductive fluid

M. Kopp, A. Tur, V. Yanovsky

► To cite this version:

M. Kopp, A. Tur, V. Yanovsky. Chaotic dynamics of magnetic fields generated by thermomagnetic instability in a nonuniformly rotating electrically conductive fluid. Journal of Physical Studies, 2023, 27 (2), pp.2403. 10.30970/jps.27.2403. hal-04894016

HAL Id: hal-04894016

https://hal.science/hal-04894016v1

Submitted on 20 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
CHAOTIC DYNAMICS OF MAGNETIC FIELDS GENERATED BY
THERMOMAGNETIC INSTABILITY IN A NONUNIFORMLY ROTATING
ELECTRICALLY CONDUCTIVE FLUID

M. I. Kopp<fnref n="1" /> ID, A. V. Tur<fnref n="3" /> ID, V. V. Yanovsky<fnref n="1" />,<fnref n="2" /> ID

(Received 10 Undecimber 2023; in final form 18 May 2023; accepted 26 May 2023; published online 15 June 2023)

The chaotic behavior of thermal convection in a nonuniformly rotating electrically conductive fluid under the action of a constant vertical magnetic field B₀ is studied. In the presence of vertical temperature gradients $\nabla T_0$ and the thermo-electromotive force coefficient $\nabla\alpha$, thermomagnetic (TM) instability arises, leading to the generation of magnetic fields. The magnetic field B₁ excited by the effect of the Biermann battery is directed perpendicular to the plane of the vectors $\nabla T_0$, $\nabla\alpha$, and the gradient of temperature disturbances $\nabla T_1$. This magnetic field changes the heat transfer regime, and due to the effects of convective heat transfer and Righi-Leduc, a positive feedback is established, which leads to an increase in magnetic field disturbances. Using the truncated Galerkin method, a nonlinear dynamic system of equations is obtained, which describes the processes of generation and regeneration of the magnetic field. Numerical analysis of these equations showed the existence of a regular, quasi-periodic, and chaotic behavior of magnetic field disturbances, accompanied by its inversions. Applying the method of perturbation theory to the nonlinear dynamic system of equations, we obtained the Ginzburg-Landau equation for the weakly nonlinear stage of nonuniformly rotating magnetic convection, taking into account TM effects. The solution of this equation showed that the stationary level of the generated magnetic fields increases with allowance for the influence of the TM instability.

**Key words:** generation of magnetic fields, thermomagnetic instability, chaotic behavior, Ginzburg-Landau equation.

DOI: https://doi.org/10.30970/jps.27.2403

I. INTRODUCTION

Recently, issues related to the physical nature of reversals and variations of the geomagnetic field have been actively studied. Reversals or polarity reversals of the Earth's magnetic field are confirmed by paleomagnetic and archeomagnetic data [1-3]. Braginsky [4] noted that the average geomagnetic field fluctuates over a period of the order of $10^3$ years. Higher frequencies in the spectra of the geomagnetic field have periods of the order of $10^2$ and shorter. Braginsky associated the main reason for the manifestation of a discrete spectrum of magnetic field variations with the excitation of MAC waves due to the action of magnetic, Archimedean, and Coriolis forces.

At present, the problem of the random nature of inversions is covered from two points of view. The first is the presence of an internal mechanism of chaotic inversion that does not require the involvement of external factors. One of such models is the convective dynamo theory, in which the magnetic field is generated due to the convective currents of the conducting medium [5]. Numerical calculations carried out in [6] perfectly reflected the dipole morphology of the Earth's magnetic field and its chaotic inversions. The electromechanical model of terrestrial magnetism proposed by Rikitake [7] can be attributed to the same category. The study of Rikitake's dynamic system of equations was also used to explain the chaotic inversion of the geomagnetic field [8-11]. In the papers [12, 13], a modified system of Rikitake

equations was studied, taking friction into account, and the possibility of magnetic field reversals was shown. At the same time, [12] states that after chaotic behavior, the system goes into a stable mode. According to the authors of [12], such a regime can describe superchrons in the geomagnetic field reversal.

The second point of view is related to the introduction of a random factor, for example, by adding random noise to the $\alpha$-effect. As shown by numerical calculations (see [14] and references therein), such a model is also capable of reproducing a chaotic reversal of the polarity of the geomagnetic field. In our opinion, the introduction of external random factors for modeling the chaotic reversal of the magnetic field gives the description a phenomenological character. Therefore, in our work, we will adhere to the first concept. The object of our research is a nonlinear system of differential equations that describes non-uniformly rotating magnetic convection in an electrically conductive liquid, taking into account thermomagnetic phenomena.

Interest in thermomagnetic phenomena in electrically conductive media arose in connection with the problem of the origin of “seed” (or initial) magnetic fields [15]. In this case, it is not always right to assume that only weak seed fields arise, which are necessary to turn on the dynamo. However, as noted in [16], magnetic fields of the Earth and planets can be created by thermoelectric currents that flow in a highly conductive region inside the planet. In this case, by analogy with the Biermann-Schlüter battery effect [17], due to the

<figure><img src="image_2.png" /></figure>

<footnote n="1">¹Institute for Single Crystals, NAS Ukraine, 60, ave. Nauky, Kharkov, UA-61072, Ukraine,</footnote>

<footnote n="2">²V. N. Karazin Kharkiv National University, 4, Svobody Sq., Kharkov, UA-61022, Ukraine</footnote>

<footnote n="3">³Université de Toulouse [UPS], CNRS, Institut de Recherche en Astrophysique et Planétologie,
9, ave. du Colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France</footnote>
non-parallelism of the temperature gradients $\nabla T$ and the thermo-electromotive force coefficient $\nabla \alpha$, a magnetic field $\partial \mathbf{B}_1/\partial t \approx [\nabla T \times \nabla \alpha]$. On the other hand, as shown in [18], and under the condition $[\nabla T_0 \times \nabla \alpha] = 0$, it is also possible to generate a magnetic field due to the development of thermomagnetic (TM) instability. Initially, TM instability was proposed to explain the spontaneous generation of strong magnetic fields ($\sim 10^6$G) in an inhomogeneous laser plasma [19–21] in negligibly short times $\sim 10^{-9}$sec. A necessary condition for the development of TM instability is the inhomogeneity and nonisothermality of the plasma. The physical mechanism of this instability is as follows. Temperature perturbations $T_1$, acting in a direction different from the initial plasma inhomogeneity, lead to the excitation of a magnetic field due to the “battery” effect $\partial \mathbf{B}_1/\partial t \approx [\nabla T_1 \times \nabla n_0]$. The magnetic field, in turn, affects the electronic thermal conductivity and changes the heat transfer mode. The emerging heat flux supplies energy to the region with elevated temperatures, thus contributing to the growth of initial perturbations.

<figure><img src="image_3.png" /><figcaption>Fig. 1. a) Scheme of shear flow in rotating flows; b) Cartesian approximation of the problem for nonuniformly rotating magnetoconvection. The nonuniform rotation in a local Cartesian coordinate system consists of a rotation with a constant angular velocity $\Omega_0$ and a shear velocity $\mathbf{U}_0 = -q\Omega_0 X e_y$ ($q = -d \ln \Omega / d \ln R$)</figcaption></figure>

Astrophysical applications of the TM instability are covered in detail in the works of Dolginov [22], Urpin et al. [23], Urpin [24, 25], where an explanation is given for the occurrence of strong magnetic fields in the nuclei of white dwarfs, binary systems, and neutron stars. Thermomagnetic instabilities in a magnetized rotating plasma were studied by Montani et al. [26], where it was shown that thermomagnetic and magnetorotational phenomena in accretion discs contribute to the transfer of angular momentum.

The purpose of this work is to study the chaotic behavior of magnetic fields generated by TM instability in an inhomogeneously rotating layer of an electrically conductive fluid (plasma) in the presence of an external constant magnetic field $\mathbf{B}_0 \parallel OZ$, collinear gradients of temperature $\nabla T_0 = -\mathbf{e}_z(dT_0/dz)$, and a thermoelectromotive force coefficient $\nabla \alpha = \mathbf{e}_z(d\alpha/dz)$. The chaotic behavior of convection has been intensively studied in rotating fluid layers [27–28], and also in conducting media rotating with a magnetic field [29]. However, these papers did not consider the dynamics of the magnetic field itself, which corresponds to the non-inductive approximation. Such problems are more important for technological applications such as crystal growth, solidification chemical processes, centrifugal casting of metals, etc. than for astrophysical and geophysical problems.

## II. BASIC EQUATIONS AND STATEMENT OF THE PROBLEM

Let us consider the physical mechanism of the generation of magnetic fields by TM-instability involving the effects associated with the inhomogeneity of the specific thermoelectric power and “magnetization” of the heat flux. The geometry of the problem is shown in Fig. 1. A higher temperature, $T_d$, is maintained on the lower plane of the layer than on the upper plane, $T_u: T_d > T_u$ — heating from below. The thermo-electromotive force coefficient on the lower (“hot”) $\alpha_d$ plane is less than on the upper (“cold”) $\alpha_u: \alpha_d < \alpha_u$. Such a situation is quite possible if we take into account the temperature dependence of the thermo-electromotive force coefficient $\alpha \sim \psi/T_0$ ($\psi$ – chemical potential) [30]. The spatially inhomogeneous distribution of $T_0(z)$ and $\alpha(z)$ inside the layer can be represented as a linear dependence on $z$:

$$T_0(z) = T_d - \frac{\Delta T}{h} \cdot z, \quad \Delta T = T_d - T_u, \\ \alpha(z) = \alpha_d + \frac{\Delta \alpha}{h} \cdot z, \quad \Delta \alpha = \alpha_u - \alpha_d.$$

Let the temperature difference at the layer boundaries in the gravity field $g$ lead to a violation of mechanical equilibrium in the system, in which convective instability develops and convective cells are formed. Temperature perturbations acting in the transverse direction on the initial inhomogeneity lead to the appearance of an eddy thermal current due to the difference in the values of the thermo-electromotive force coefficient $\Delta\alpha(z)$ at the layer boundaries. This current creates perturbations in the Y-direction of magnetic field $\mathbf{B}_1$, which will affect the heat transfer mode. The excited magnetic field $\mathbf{B}_1$ creates heat fluxes directed perpendicular to the field itself and the temperature gradient. Thus, positive feedback is established: the newly emerging heat fluxes create a vortex of thermoelectric power, which enhances the perturbations of the magnetic field $\mathbf{B}_1$. The direction of the external magnetic field $\mathbf{B}_0$ coincides with the axis of rotation of the fluid $\Omega \parallel OZ$ (see Fig. 1).
In our earlier research [18], the evolution equations for perturbations in the local Cartesian coordinate system

were discovered. These equations are expressed in dimensionless variables and read as follows:

$$ 
\begin{aligned}
& \left(\frac{\partial}{\partial t} - \nabla^2\right) \nabla^2 \psi + \sqrt{\text{Ta}} \frac{\partial v}{\partial z} - \text{Pr} \text{Pm}^{-1} \text{Q} \frac{\partial}{\partial z} \nabla^2 \phi - \text{Ra} \frac{\partial \theta}{\partial x} = \text{Pr} \text{Pm}^{-1} \text{Q} \cdot J(\phi, \nabla^2 \phi) - \text{Pr}^{-1} \cdot J(\psi, \nabla^2 \psi), \\
& \left(\frac{\partial}{\partial t} - \nabla^2\right) v - \sqrt{\text{Ta}}(1 + \text{Ro}) \frac{\partial \psi}{\partial z} - \text{Pr} \text{Pm}^{-1} \text{Q} \frac{\partial \tilde{v}}{\partial z} = \text{Pr} \text{Pm}^{-1} \text{Q} \cdot J(\phi, \tilde{v}) - \text{Pr}^{-1} \cdot J(\psi, v), \\
& \left(\frac{\partial}{\partial t} - \text{Pm}^{-1} \nabla^2\right) \phi - \text{Pr}^{-1} \frac{\partial \psi}{\partial z} = -\text{Pr}^{-1} J(\psi, \phi), \\
& \left(\frac{\partial}{\partial t} - \text{Pm}^{-1} \nabla^2\right) \tilde{v} - \text{Pr}^{-1} \frac{\partial v}{\partial z} + \text{Ro} \sqrt{\text{Ta}} \frac{\partial \phi}{\partial z} - \text{Pm}^{-1} R_{\alpha} \frac{\partial \theta}{\partial x} = \text{Pr}^{-1} (J(\phi, v) - J(\psi, \tilde{v})), \\
& \left(\text{Pr} \frac{\partial}{\partial t} - \nabla^2\right) \theta - \frac{\partial \psi}{\partial x} + q_{\alpha} \frac{\partial \tilde{v}}{\partial x} = -J(\psi, \theta) - q_{\alpha}^{(1)} \theta \cdot \frac{\partial \tilde{v}}{\partial x} + q_{\alpha}^{(2)} J(\theta, \tilde{v}),
\end{aligned}
\quad (1) $$

where

$$ J(a, b) = \frac{\partial a}{\partial X} \frac{\partial b}{\partial Z} - \frac{\partial a}{\partial Z} \frac{\partial b}{\partial X} $$

are the Jacobian or Poisson brackets $J(a, b) \equiv \{a, b\}$, $\psi$ is the hydrodynamic stream function, $\phi$ is the magnetic stream function, $\theta$ is a temperature perturbation, $v$ and $\tilde{v}$ are velocity and magnetic field perturbations for the Y-component. $\text{Pr} = \nu/\chi$ is the Prandtl number, $\text{Pm} = \nu/\eta$ is the Prandtl magnetic number, $\text{Ta} = \frac{4\Omega_0^2 h^4}{\sqrt{\rho_0 \mu \nu \eta}}$ is the Taylor number, $\text{Ha} = \frac{B_0 h}{\sqrt{\rho_0 \mu \nu \eta}}$ is the Hartman number, $\text{Ra} = \frac{g \beta_T (\Delta T) h^3}{\nu \chi}$ is the Rayleigh number, $R_\alpha = \frac{\Delta \alpha \Delta T}{\eta B_0}$ is the thermoelectromotive force number. The dimensionless parameters

$$ 
\begin{aligned}
q_{\alpha} &= \frac{\alpha T_0 B_0}{\rho_0 c_p \mu \chi (\Delta T)} \left[ \left( \frac{\mu \mathcal{L}}{\alpha} + 1 \right) \frac{\Delta T}{T_0} + \frac{\Delta \alpha}{\alpha} \right], \\
q_{\alpha}^{(1)} &= \frac{\Delta \alpha B_0}{\rho_0 c_p \mu \chi},
\end{aligned}
$$

$$ q_{\alpha}^{(2)} = \frac{\alpha B_0}{\rho_0 c_p \mu \chi} \left( \frac{\mu \mathcal{L}}{\alpha} + 1 \right) $$

are related to the influence of the thermopower effect and the Righi-Leduc effect on the heat transfer process. The system of equations (1) is supplemented with the following boundary conditions:

$$ 
\begin{aligned}
\psi\bigg|_{z=0,h} &= \nabla^2 \psi\bigg|_{z=0,h} = \frac{dv}{dz}\bigg|_{z=0,h} = \tilde{v}\bigg|_{z=0,h} = 0, \\
\frac{d\phi}{dz}\bigg|_{z=0,h} &= \theta\bigg|_{z=0,h} = 0.
\end{aligned}
\quad (2) $$

Without taking TM effects into account, the system of equations (1) was used to study weakly nonlinear and chaotic convection regimes in an inhomogeneously rotating plasma in an axial (vertical) magnetic field [31–33].

### III. A NONLINEAR SYSTEM OF DIFFERENTIAL EQUATIONS FOR A CONVECTIVE TM DYNAMO

One of the most effective methods for studying hydrodynamic instability is the Galerkin method (see, for example, [34]). The application of this method to the equations of hydrodynamics makes it possible to obtain a wide variety of different nonlinear dynamical systems of varying degrees of complexity [35]. Using the truncated expansion of Galerkin, Lorenz [36] obtained a dynamical system of equations with a small number of internal parameters for modeling free convection in the atmosphere (these are known as the Lorenz equations). Following [36], we apply the truncated expansion of Galerkin to the study of the weakly nonlinear stage of the development of convective instability described by the equations (1). Then we can represent all perturbations in the equations (1) in the form of a Fourier expansion of the following form:

$$ \psi(x, z, t) = A_1(t) \sin(kx) \sin(\pi z), $$

$$ v(x, z, t) = V_1(t) \sin(kx) \cos(\pi z), $$

$$ \phi(x, z, t) = B_1(t) \sin(kx) \cos(\pi z), \quad (3) $$

$$ \tilde{v}(x, z, t) = W_1(t) \sin(kx) \sin(\pi z), $$

$$ \theta(x, y, t) = C_1(t) \cos(kx) \sin(\pi z) + C_2(t) \sin(2\pi z), $$

where $k$ is a dimensionless wave number; $A_1, V_1, B_1, W_1, C_1$, and $C_2$ are the disturbance amplitudes. Substituting the expansions (3) into the equations (1), after orthogonalization, we obtain the evolution equations for
the disturbance amplitudes:

$$ \frac{\partial A_1}{\partial \tilde{t}} = -A_1 - \frac{\pi \sqrt{\text{Ta}}}{a^4} \cdot V_1 - \frac{\pi \text{QPr}}{a^2 \text{Pm}} \cdot B_1 + \frac{k \text{Ra}}{a^4} \cdot C_1, $$

$$ \frac{\partial V_1}{\partial \tilde{t}} = -V_1 + \frac{\pi \sqrt{\text{Ta}}}{a^2} (1 + \text{Ro}) \cdot A_1 + \frac{\pi \text{QPr}}{a^2 \text{Pm}} \cdot W_1, $$

$$ \text{Pm} \frac{\partial B_1}{\partial \tilde{t}} = -B_1 + \frac{\pi \text{Pm}}{a^2 \text{Pr}} \cdot A_1, \quad (4) $$

$$ \text{Pm} \frac{\partial W_1}{\partial \tilde{t}} = -W_1 - \frac{\pi \text{Pm}}{a^2 \text{Pr}} \cdot V_1 + \frac{\pi \text{PmRo}\sqrt{\text{Ta}}}{a^2} \cdot B_1 - \frac{k R_\alpha}{a^2} \cdot C_1, $$

$$ \text{Pr} \frac{\partial C_1}{\partial \tilde{t}} = -C_1 + \frac{k}{a^2} \cdot A_1 + \frac{\pi k}{a^2} \cdot A_1 C_2 - \frac{k q_\alpha}{a^2} \cdot W_1 + \frac{\pi k q_\alpha^{(2)}}{a^2} \cdot W_1 C_2, $$

$$ \text{Pr} \frac{\partial C_2}{\partial \tilde{t}} = -\frac{4\pi^2}{a^2} \cdot C_2 - \frac{\pi k}{2a^2} \cdot A_1 C_1 - \frac{\pi k q_\alpha^{(2)}}{2a^2} \cdot W_1 C_1, $$

where $a = \sqrt{k^2 + \pi^2}$ is the total wavenumber and $\tilde{t} = a^2 t$ is the reduced time. The system of ordinary differential equations (4) is a low-order Fourier expansion model, but it can reproduce convective processes qualitatively in a nonlinear system of equations (1). For convenience, we introduce the following notation:

$$ R = \frac{k^2 R \alpha}{a^6}, \quad T = \frac{\pi^2 \sqrt{\text{Ta}}}{a^6}, \quad H = \frac{\pi^2 \text{QPr}}{a^4 \text{Pm}}, \quad G = \frac{k^2 R_\alpha}{\pi^2 \text{Pm}}, $$

$$ q_1 = \frac{\pi^2 q_\alpha}{a^4}, \quad q_2 = \frac{\pi^2 q_\alpha^{(2)}}{a^4}, \quad b = \frac{4\pi^2}{a^2}, \quad (5) $$

and rescale the amplitudes $A_1, V_1, B_1, W_1, C_1, C_2$ in the form

$$ X(\tilde{t}) = \frac{k\pi}{a^2\sqrt{2}}A_1(\tilde{t}), \quad V(\tilde{t}) = \frac{kV_1(\tilde{t})}{\sqrt{2}}, $$

$$ U(\tilde{t}) = \frac{k B_1(\tilde{t})}{\sqrt{2}}, \quad W(\tilde{t}) = \frac{a^2 k}{\pi \sqrt{2}} W_1(\tilde{t}), $$

$$ Y(\tilde{t}) = \frac{\pi C_1(\tilde{t})}{\sqrt{2}}, \quad Z(\tilde{t}) = -\pi C_2(\tilde{t}). $$

Then the equations (4) take the form of a nonlinear dynamic system of equations:

$$ \left\{ \begin{aligned} \dot{X} &= -X + RY - TV - HU \\ \dot{V} &= -V + HW + \sqrt{\text{Ta}}(1 + \text{Ro})X \\ \dot{U} &= -\text{Pm}^{-1}U + \text{Pr}^{-1}X \\ \dot{W} &= -\text{Pm}^{-1}W - \text{Pr}^{-1}V + \text{Ro}\sqrt{\text{Ta}}U - \text{GY} \\ \dot{Y} &= \text{Pr}^{-1}(-Y + X - XZ - q_1W - q_2WZ) \\ \dot{Z} &= \text{Pr}^{-1}(-\gamma Z + XY + q_2WY), \end{aligned} \right. \quad (6) $$

where the dot (·) on the symbol indicates the differentiation with respect to time $\tilde{t}$. In the limiting case, when there are no TM effects ($q_1 = q_2 = G = 0$), the equations (6) go over to the Lorenz equations for the six-dimensional (6D) phase space, which were studied numerically in [31]. The system of nonlinear equations obtained by us (6) depends on a large number of dimensionless parameters (11 parameters). For this reason, it must have a huge variety of behavior modes, and all possible transitions to chaos can be realized in it, depending on the range of change of various dimensionless parameters. In contrast to [31], the last three equations in the system (6) describe the generation and regeneration of a W-component (or toroidal component) magnetic field. A similar process of self-excitation of the magnetic field is observed in Rikitake's well-known electromechanical model [7], which is used to explain the chaotic inversion of the Earth's magnetic field [8]. Therefore, the system of equations (6) can be attributed to the dynamic model of a nonlinear convective TM dynamo.

Note that, in contrast to the Lorenz equations and the classical Rikitake dynamo system, we have obtained nonlinear dynamic equations of a 6D dimension (6) with four quadratic nonlinearities.

## IV. STABILITY ANALYSIS

Qualitative and numerical analyses of the dynamic system of equations (6) make it possible to determine the type of stationary points and the conditions for the emergence of a chaotic regime.

### A. Homogeneity and symmetry

The trivial solution (6), corresponding to the absence of convection, leads to the appearance of a special fixed point:

$$ O(X, V, U, W, Y, Z) = O(0, 0, 0, 0, 0, 0), $$

which does not depend on the values of the parameters R, T, Ta, H, Pm, Pr, b. Given a coordinate transformation T as follows

$$ \mathbf{T}(X, V, U, W, Y, Z) \rightarrow (-X, -V, -U, -W, -Y, Z). $$

It is clear that each trajectory is symmetrical with respect to the Z-axis. That means system (6) is invariant for the given transformation T.

### B. Dissipativity

The divergence of system (6) can be calculated as

$$ \begin{aligned} \text{div} \mathbf{\Phi} &= \frac{\partial \dot{X}}{\partial X} + \frac{\partial \dot{V}}{\partial V} + \frac{\partial \dot{U}}{\partial U} + \frac{\partial \dot{W}}{\partial W} + \frac{\partial \dot{Y}}{\partial Y} + \frac{\partial \dot{Z}}{\partial Z} \\ &= -2(1 + \text{Pm}^{-1}) - \text{Pr}^{-1}(1 + \gamma) < 0. \end{aligned} $$
As a result, it follows that the system (6) is dissipative, since the divergence of the vector field is negative. Due to dissipation, the phase volume shrinks:

$$ \Phi(t) = \Phi(0) \exp \left[ \left( -2(1 + \mathrm{Pm}^{-1}) - \mathrm{Pr}^{-1}(1 + \gamma) \right) t \right]. $$

Consequently, in the phase space of dissipative systems, attracting sets can arise — attractors.

### C. Equilibrium points

Equating the left sides of equations (6) to zero, we get three equilibrium states:

$$ O_1(X_1, V_1, U_1, W_1, Y_1, Z_1), \quad O_2(X_2, V_2, U_2, W_2, Y_2, Z_2), \quad O_3(X_3, V_3, U_3, W_3, Y_3, Z_3). \quad (7) $$

After computation, we obtain the following equilibrium points:

$$ (X_1, V_1, U_1, W_1, Y_1, Z_1) = (0, 0, 0, 0, 0, 0), $$

$$ (X_2, X_3) = \pm \frac{R+r_1}{r\tilde{r}} \cdot \sqrt{\gamma \left( \frac{R+r_1}{r} \cdot (1+q_1r_2) + q_1r_3 - 1 \right)}, $$

$$ (V_2, V_3) = \pm \left[ \frac{\sqrt{\mathrm{Ta}}(1+\mathrm{Ro}+\mathrm{HRoPm}^2\mathrm{Pr}^{-1})(R+r_1)-r\mathrm{HGPm}}{r\tilde{r}(1+\mathrm{HPmPr}^{-1})} \right] \sqrt{\gamma \left( \frac{R+r_1}{r} \cdot (1+q_1r_2)+q_1r_3-1 \right)}, $$

$$ (W_2, W_3) = \pm \left[ \frac{\sqrt{\mathrm{TaPmPr}^{-1}(\mathrm{RoPm} - \mathrm{Ro} - 1)(R+r_1) - r\mathrm{GPm}}}{r\tilde{r}(1+\mathrm{HPmPr}^{-1})} \right] \sqrt{\gamma \left( \frac{R+r_1}{r} \cdot (1+q_1r_2)+q_1r_3-1 \right)}, $$

$$ (U_2, U_3) = \pm \frac{\mathrm{Pm}(R+r_1)}{r\tilde{r}\mathrm{Pr}} \cdot \sqrt{\gamma \left( \frac{R+r_1}{r} \cdot (1+q_1r_2)+q_1r_3-1 \right)}, $$

$$ (Y_2, Y_3) = \pm \frac{1}{\tilde{r}} \sqrt{\gamma \left( \frac{R+r_1}{r} \cdot (1+q_1r_2)+q_1r_3-1 \right)}, $$

$$ (Z_2, Z_3) = \frac{1}{\tilde{r}} \left( \frac{R+r_1}{r} \cdot (1+q_1r_2)+q_1r_3-1 \right), $$

where

$$ r = 1 + \frac{\mathrm{Pm}}{\mathrm{Pr}}\mathrm{H} + T\sqrt{\mathrm{Ta}} \cdot \frac{1 + \mathrm{Ro} \left( 1 + \frac{\mathrm{Pm}^2}{\mathrm{Pr}}\mathrm{H} \right)}{1 + \frac{\mathrm{Pm}}{\mathrm{Pr}}\mathrm{H}}, $$

$$ r_1 = \frac{\mathrm{THGPm}}{1 + \mathrm{HPmPr}^{-1}}, $$

$$ r_2 = \frac{\sqrt{\mathrm{TaPmPr}^{-1}(1 + \mathrm{Ro}(1 - \mathrm{Pm}))}}{1 + \mathrm{HPmPr}^{-1}}, $$

$$ r_3 = \frac{\mathrm{GPm}}{1 + \mathrm{HPmPr}^{-1}}, \quad \tilde{r} = \frac{R+r_1}{r} \cdot (1-q_2r_2) - q_2r_3. $$

If the Rayleigh parameter

$$ R = \frac{r(1 - q_1 r_3)}{1 + q_1 r_2} - r_1 = R_{cr}, $$

there is only one fixed point $O_1$ in a system. For Rayleigh parameters $R > R_{cr}$ that are two more fixed points $O_2$ and $O_3$, which are symmetric to each other. For Rayleigh parameters $R < R_{cr}$ the coordinates $O_2$ and $O_3$ are imaginary.

### D. Stability of the equilibrium points

We proceed to the study of stability, found the equilibrium states (7). To do this, we linearize the system of equations (6) in a small area of fixed points. Representing all variables in the form
$$ (X, V, U, W, Y, Z)^{\mathrm{T}} = (X_0, V_0, U_0, W_0, Y_0, Z_0)^{\mathrm{T}} + (X', V', U', W', Y', Z')^{\mathrm{T}} \cdot e^{\lambda t} $$

we obtain a linear algebraic system:

$$
\left\{
\begin{array}{l}
(\lambda + 1)X' = RY' - TV' - HU' \\
(\lambda + 1)V' = HW' + \sqrt{Ta}(1 + Ro)X' \\
(\lambda + Pm^{-1})U' = \mathrm{Pr}^{-1}X' \\
(\lambda + Pm^{-1})W' = -\mathrm{Pr}^{-1}V' + Ro\sqrt{Ta}U' - GY' \\
(\lambda + \mathrm{Pr}^{-1})Y' = \mathrm{Pr}^{-1}(X' - X_0Z' - X'Z_0 - q_1W' - q_2W_0Z' - q_2Z_0W') \\
(\lambda + \gamma\mathrm{Pr}^{-1})Z' = \mathrm{Pr}^{-1}(X_0Y' + X'Y_0 + q_2W_0Y' + q_2Y_0W')
\end{array}
\right.
\qquad (8)
$$

<table><caption>Table 1. Eigenvalues $\lambda_{1,2,3,4,5,6}$ for the fixed point $O_1$, calculated for different values of the parameter R.</caption><thead><tr><th>R</th><th>λ<sub>1</sub></th><th>λ<sub>2</sub></th><th>λ<sub>3</sub></th><th>λ<sub>4</sub></th><th>λ<sub>5</sub></th><th>λ<sub>6</sub></th></tr></thead><tbody><tr><td>1</td><td>-1.013776<br>+i1.213922</td><td>-0.029541</td><td>-0.266666</td><td>-0.487147</td><td>-1.555757</td><td>-1.013776<br>-i1.213922</td></tr><tr><td>8</td><td>0.334085</td><td>-1.098758<br>+i1.094230</td><td>-0.266666</td><td>-0.601432</td><td>-1.635135</td><td>-1.098758<br>-i1.094230</td></tr><tr><td>34.004</td><td>1.205252</td><td>-1.204008<br>+i0.765665</td><td>-0.266666</td><td>-0.744572</td><td>-2.152663</td><td>-1.204008<br>-i0.765665</td></tr><tr><td>34.005</td><td>1.205279</td><td>-1.204007<br>+i0.765656</td><td>-0.266666</td><td>-0.744575</td><td>-2.152688</td><td>-1.204007<br>-i0.765656</td></tr><tr><td>53.93</td><td>1.683372</td><td>-1.172021<br>+i0.648537</td><td>-0.266666</td><td>-0.794802</td><td>-2.644527</td><td>-1.172021<br>-i0.648537</td></tr><tr><td>58</td><td>1.768589</td><td>-1.165177<br>+i0.634200</td><td>-0.266666</td><td>-0.802489</td><td>-2.735746</td><td>-1.165177<br>-i0.634200</td></tr><tr><td>290</td><td>4.775845</td><td>-1.046055<br>+i0.486160</td><td>-0.266666</td><td>-0.938911</td><td>-5.844824</td><td>-1.046055<br>-i0.486160</td></tr></tbody></table>

Here the index 0 denotes fixed points, the prime ' denotes small perturbations, and $\lambda$ is the increase increment.

The modes of stability or instability are determined by the signs of the real parts of the eigenvalues of the determinant of the system of equations (8). The eigenvalues of the system of equations (8) $\lambda$ are found from the solution of the characteristic equation obtained by equating the determinant to zero, i.e.

$$ M_0 \cdot (M_0 \cdot M_2 + R\mathrm{Pr}^{-1}(\lambda + Pm^{-1}) \cdot M_3) + T\sqrt{Ta} \cdot M_1 \cdot M_2 \\ - G\mathrm{Pr}^{-1} \cdot M_4 - R\mathrm{Pr}^{-2} \cdot M_5 + (\lambda + Pm^{-1})\mathrm{GHT}\mathrm{Pr}^{-1} \cdot M_3 = 0, \qquad (9) $$

where the following notation is entered

$$ M_0 = (\lambda + 1)(\lambda + Pm^{-1}) + H\mathrm{Pr}^{-1}, \quad M_1 = (1 + Ro)(\lambda + Pm^{-1})^2 + H\mathrm{Ro}\mathrm{Pr}^{-1}, $$

$$ M_2 = (\lambda + \mathrm{Pr}^{-1})(\lambda + \gamma\mathrm{Pr}^{-1}) + (X_0 + q_2W_0)^2\mathrm{Pr}^{-2}, $$

$$ M_3 = (X_0 + q_2W_0)Y_0\mathrm{Pr}^{-1} - (\lambda + \gamma\mathrm{Pr}^{-1})(1 - Z_0), $$

$$ M_4 = [(\lambda + 1)^2(\lambda + Pm^{-1}) + T\sqrt{Ta}(1 + Ro)(\lambda + Pm^{-1}) + H\mathrm{Pr}^{-1}(\lambda + 1)] \cdot \tilde{M}, $$

$$ M_5 = [\sqrt{Ta}(1 + Ro)(\lambda + Pm^{-1}) - (\lambda + 1)Ro] \cdot \tilde{M}, $$

$$ \tilde{M} = \mathrm{Pr}^{-1}(X_0 + q_2W_0)q_2Y_0 + (q_1 + q_2Z_0)(\lambda + \gamma\mathrm{Pr}^{-1}). $$
Here $W_0 = (W_1, W_2, W_3)$, $X_0 = (X_1, X_2, X_3)$, $Y_0 = (Y_1, Y_2, Y_3)$, $Z_0 = (Z_1, Z_2, Z_3)$ are the coordinates of fixed points. If we substitute the values of the three equilibrium states (7) into equation (9), we obtain the characteristic equations for the eigenvalues $\lambda$ in each of these states. Moreover, for the points $O_2$ and $O_3$ the characteristic equations coincide due to the symmetry.

The characteristic equation (9) is reduced to the algebraic equation of the sixth degree:

$$P(\lambda) \equiv a_0\lambda^6 + a_1\lambda^5 + a_2\lambda^4 + a_3\lambda^3 + a_4\lambda^2 + a_5\lambda + a_6 = 0,$$

where $a_0 = 1 > 0$. We do not give the explicit form of the real coefficients $a_1, a_2, a_3, a_4, a_5, a_6$ due to a very cumbersome form. There exist methods that enable the assessment of the system stability without explicitly solving its characteristic equation. One of such methods is the Routh-Hurwitz or Liénard-Schipar [37] criterion, which provides the necessary and sufficient conditions for the stability of the system. In the last criterion, the number of determinant inequalities is approximately half as large as in the Routh-Hurwitz conditions, so it is reasonable to use it. The criterion makes it possible to determine whether all the roots of the characteristic equation have negative real parts, indicating the stability of the system. According to the Liénard-Schipar criterion, for the polynomial $P(\lambda)$ to have all roots with negative real parts, it is necessary and sufficient that the following conditions are satisfied:

<table><caption>Table 2. Eigenvalues $\lambda_{1,2,3,4,5,6}$ for the fixed point $O_{2,3}$, calculated for different values of the parameter R.</caption><thead><tr><th>R</th><th>λ<sub>1</sub></th><th>λ<sub>2</sub></th><th>λ<sub>3</sub></th><th>λ<sub>4</sub></th><th>λ<sub>5</sub></th><th>λ<sub>6</sub></th></tr></thead><tbody><tr><td>1</td><td>0.048075</td><td>-1.020790<br>+i1.207059</td><td>-0.326499</td><td>-0.487219</td><td>-1.559443</td><td>-1.020790<br>-i1.207059</td></tr><tr><td>8</td><td>0.545952</td><td>-1.138779<br>+i1.047646</td><td>-0.344958</td><td>-0.614024</td><td>-1.676077</td><td>-1.138779<br>-i1.047646</td></tr><tr><td>34.004</td><td>1.644357</td><td>-1.195875<br>+i0.673213</td><td>-0.349129</td><td>-0.767817</td><td>-2.502325</td><td>-1.195875<br>-i0.673213</td></tr><tr><td>34.005</td><td>1.644390</td><td>-1.195873<br>+i0.673206</td><td>-0.349129</td><td>-0.767820</td><td>-2.502360</td><td>-1.195873<br>-i0.673206</td></tr><tr><td>53.93</td><td>2.226098</td><td>-1.149037<br>+i0.585517</td><td>-0.349653</td><td>-0.818773</td><td>-3.126262</td><td>-1.149037<br>-i0.585517</td></tr><tr><td>58</td><td>2.331271</td><td>-1.141455<br>+i0.575978</td><td>-0.349600</td><td>-0.827007</td><td>-3.238417</td><td>-1.141455<br>-i0.575978</td></tr><tr><td>290</td><td>5.966096</td><td>-1.036574<br>+i0.473629</td><td>-0.350281</td><td>-0.951500</td><td>-6.957833</td><td>-1.036574<br>-i0.473629</td></tr></tbody></table>

To perform a numerical analysis of the equation (9), we choose the values of the parameter $b$ and the Prandtl number $Pr$ as is customary in the Lorenz equations [38]: $b = 8/3$, $Pr = 10$. The rotation parameter is assumed to be equal to $T = 0.1$ as in [28], where the emergence of a chaotic regime in a uniformly rotating fluid layer was studied. The parameter value $T = 0.1$ gives us the Taylor number $Ta = 1080$. We consider the magnetic Prandtl number to be equal to $Pm = 1$. We set the parameter for the external axial magnetic field equal to $H = 2$, which corresponds to the Chandrasekhar number

1. The coefficients of the polynomial must all be real numbers.

2. All coefficients of the polynomial $P(\lambda)$ are positive: $a_n > 0, n = 1, 2, \dots, 6$. The Hurwitz determinant is constructed from the coefficients of the polynomial $P(\lambda)$ [37]:

$$\Delta_n = \begin{vmatrix} a_1 & a_3 & a_5 & \cdots & 0 \\ a_0 & a_2 & a_4 & \cdots & 0 \\ 0 & a_1 & a_3 & \cdots & 0 \\ 0 & a_0 & a_2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \cdots \cdots \cdots \cdots \cdots & a_n \end{vmatrix}$$

3. All the principal minors of the Hurwitz determinant must be positive: $\Delta_{n-1} > 0, \Delta_{n-3} > 0, \dots$

If all these conditions are met, then the polynomial has all its roots with negative real parts, indicating the stability of the system. Using the Liénard-Schipar algorithm, we get the necessary and sufficient conditions for stability:

$$a_n > 0, n = 1, 2, \dots, 6, \Delta_3 > 0, \Delta_5 > 0.$$

Obviously, when the Liénard-Schipar criterion is satisfied, the fixed points are stable, and their equilibrium positions are classified as stable nodes.

$Q \approx 4.45$. This order of magnitude of the magnetic field was used in [29] to study the effect of the magnetic field on the chaotic behavior of convection. The generation parameters $G, q_1, q_2$ taking into account $b = 8/3$, and $Pm = 1$ take the following form:

$$G = 0.5R_{\alpha}, \quad q_{1} \approx 0.045q_{\alpha}, \quad q_{2} \approx 0.045q_{\alpha}^{(2)} \quad (10)$$

Considering the case of a low-temperature electrically conductive medium, the Righi-Leduc effect in the expressions for $q_\alpha$ and $q_\alpha^{(2)}$ can be neglected $\mu L \ll \alpha$.
As a result, $q_\alpha$ and $q_\alpha^{(2)}$ are respectively equal:

$$q_{\alpha} = q_{\alpha}^{(2)} \left( 1 + \frac{\Delta\alpha/\alpha}{\Delta T/T_0} \right), \quad q_{\alpha}^{(2)} \approx \frac{\alpha B_0}{\rho_0 c_p \mu \chi} = \frac{\alpha B_0}{\mu \kappa}. \quad (11)$$

For numerical estimates of the difference between the thermo-electromotive force coefficient $\Delta\alpha \approx 10^{-6}V/K$ ($\alpha \approx 10^{-4}V/K$ and the temperature $\Delta T \approx 2000$ K ($T_0 \approx 4000$)K we can assume $q_\alpha \approx q_\alpha^{(2)}$. The definition of the G parameter (see (III)) shows that G depends on $\Delta\alpha$ and $\Delta T$. Therefore, fixing the difference in the thermo-electromotive force coefficient $\Delta\alpha$, we will analyze the dispersion equation (9) and the system of equations (6) for different heating conditions, i.e. for different values of the Rayleigh parameter R. For a numerical estimate of $q_\alpha^{(2)}$, we use the values of physical quantities for molten iron [39], i.e. we take the thermal conductivity coefficient $\kappa = 39 \text{ W/m}\cdot\text{K}$ ($\chi = \kappa/(\rho_0 c_p)$), the coefficient permeability $\mu = 4\pi \cdot 10^{-7} \text{ V}\cdot\text{s/A}\cdot\text{m}$, thermo-electromotive force coefficient $\alpha \approx 2.2 \cdot 10^{-5} \text{ V/K}$ and the magnitude of

the external magnetic field $B_0 = 10^{-1} \text{T}$. As a result of substituting these quantities into the formula (11), we find $q_\alpha^{(2)} \approx 0.044$ or $q_1 \approx q_2 \approx 0.002$. The generation parameter is estimated to be $G \approx 0.01 \text{R}$ if

$$\frac{27\pi^4 \Delta\alpha}{8 \alpha} \frac{c_p \chi}{g\beta_T \eta^2 h^2} \approx 10^{-1}.$$

As a result, selecting the parameter values $Pm = 1$, $Pr = 10$, $H = 2$, $T = 0.1$, and $Ta = 1080$, $b = 8/3$, $q_1 = q_2 = 0.002$, $G \approx 0.01 \text{R}$, we calculate the eigenvalues $\lambda_i$ for the fixed point $O_1$ in the case of the Keplerian rotation profile $Ro = -3/4$ depending on the changes in the Rayleigh parameter R. These results are shown in Tab. 1, in the case of the Keplerian rotation profile $Ro = -3/4$ in Tab. 2. This shows that for negative $Re\lambda < 0$ the trajectories enter the point $O_1$, i.e. correspond to stable eigen directions, and for positive $Re\lambda > 0$ the trajectories leave the point $O_1$, and hence correspond to unstable eigen directions. The steady state of convection ($\lambda = 0$) corresponds to the critical value of the parameter $R_{2cr}$, which turns out to be equal to the first critical value:

$$R_{1cr} = \frac{1 + \text{HPmPr}^{-1} - q_1 \text{GPa}}{1 + \text{HPmPr}^{-1} + q_1 \text{PmPr}^{-1} \sqrt{Ta}(1 + \text{Ro}(1 - \text{Pm}))} \cdot r - r_1. \quad (12)$$

We note that the Rayleigh number depends on the temperature difference at the boundaries of the electrically conductive fluid layer. By changing the heating conditions at the boundaries of the fluid layer (the Rayleigh number), it is possible to study various modes of the generation of magnetic fields. Next, numerical solutions of the equations (6) are presented: phase portraits and time diagrams for the Y-component of the generated magnetic field for various values of R given in Table. 1–2.

# V. DISCUSSION OF NUMERICAL RESULTS

## A. Analysis of phase trajectories of regular and chaotic behavior

Let us write system (6) in a form more convenient for modeling:

$$\left\{
\begin{aligned}
\dot{x}_1 &= -x_1 + Rx_2 - 2x_4 - 0.1x_5 \\
\dot{x}_2 &= \frac{1}{10}(-x_2 + x_1 - x_1x_3 - 0.002x_6 - 0.002x_6x_3) \\
\dot{x}_3 &= \frac{1}{10}(-\frac{8}{3}x_3 + x_1x_2 + 0.002x_6x_2) \\
\dot{x}_4 &= -x_4 + \frac{1}{10}x_1 \\
\dot{x}_5 &= -x_5 + 8.21x_1 + 2x_6 \\
\dot{x}_6 &= -x_6 - 0.01Rx_2 - 24.65x_4 - \frac{1}{10}x_5,
\end{aligned}
\right.
\quad (13)$$

where new notation for variables was introduced: $x_1 = X$, $x_2 = Y$, $x_3 = Z$, $x_4 = U$, $x_5 = V$, $x_6 = W$. We used the

values of fixed parameters

$$H = 2, \quad Ta = 1080, \quad T = 0.1, \quad Pm = 1, \quad Pr = 10,$$

$$\gamma = 8/3, \quad Ro = -3/4, \quad q_1 = q_2 = 0.002, \quad G = 0.01 \text{R}$$

in the equations (6). In this section, we analyze numerical studies of the nonlinear system of equations (13) depending on the variations of the Rayleigh parameter R in the Mathematica computer environment.

Figure 2,a shows the case for $R = 1 < R_{1cr} \approx 1.463$. We see that the initial perturbation of the magnetic field decays (Fig. 2, b). As it is shown in Fig. 2,a–2,b the initial system state reaches the origin (the point $O_1$), which is the local and global attractor simultaneously. For number $R_{1cr} = 1.463 < R = 8$, there is loss of stability, and we observe the appearance of spiral trajectories (Fig. 2,c–2,d) in the phase plane $x_1x_6$, which will wind up around fixed point $O_2$ while the parameter R is increasing. It is noticeable already for the value $R = 34.004$ (see Fig. 2,e). Thus, the perturbed magnetic field performs oscillations with damped amplitude (see Fig. 2,f). In this case, eigenvalues $\lambda_i$ are complex-valued with a negative real part, and we classify the stationary point as a stable focus. The coordinates of the fixed point $O_2$ in the $x_1x_6$ plane can be calculated analytically using (7). For the Rayleigh parameter R = 34.004, we get the focus coordinates $O_2(x_1, x_6) = (7.78, -21.39)$, which coincide with the focus coordinates found numerically (see Fig. 2,e).
<figure><img src="image_4.png" /><figcaption>Fig. 2. a), c), e), g), i), k) projections of phase trajectories in the $x_1x_6$ plane upon the variation of parameter R; b), d), f), h), j), l) time dependences of magnetic component amplitude variations $U(t)$</figcaption></figure>

An insignificant increase in the Rayleigh parameter from $R = 34.004$ to $R = 34.005$ (Fig. 2,g) leads to a change in the sign (or direction) of the oscillating perturbed magnetic field, which also fades out (see Fig. 2h). Here the phase trajectories wind in a spiral

around the fixed point $O_3$, which we also classify as a stable focus. The focus coordinates for the Rayleigh parameter $R = 34.005$ obtained analytically $O_3(x_1, x_6) = (-7.78, 21.39)$ coincide with the focus coordinates found numerically (see Fig. 2,g).
<figure><img src="image_5.png" /><figcaption>Fig. 3. Projections of the phase trajectory in the plane $x_1x_6$ and the time diagram of variations in the amplitude of the magnetic component $x_6(t)$ at R = 290</figcaption></figure>

We can see the pre-chaotic state for the parameter R = 53.98 in figures 2,i and 2,j. Here the solution seems chaotic, but then it stabilizes to oscillatory and finally tends to a stationary state. This behavior of the system is called metastable chaos. As can be seen from Fig. 2j, chaotic behavior lasts approximately up to t = 660, and then a transition to damped oscillations occurs. Figures 2,k and 2,l show the case of irregular oscillations with aperiodic changes in the amplitude and direction (inversion) of the generated magnetic field at R = 58. A further increase in parameter R facilitates the evolution of chaotic behavior since the positive characteristic number $\lambda_1$ grows (see Tab. 1 and Tab. 2). As can be seen in the phase portrait and time diagram in Fig. 3 with R = 290, the amplitude of the chaotic magnetic field will also increase. The trajectories in Fig. 3 demonstrate highly developed chaos compared to the chaotic behavior at R = 58. Strange attractors for R = 58 and R = 290 belong to the class of self-excited attractors since the equilibrium points $O_1$ and $Q_{2,3}$ are unstable.

## B. Bifurcation diagrams, Lyapunov exponents and Kaplan-Yorke dimension

Bifurcation diagrams for the $x_1, x_2, x_3, x_4, x_5, x_6$ components of the system of equations (13) are shown in Fig. 4, from which one can see the appearance of a periodic, quasi-periodic, and chaotic regime depending on the values of the Rayleigh number R. One of the important criteria characterizing the chaotic behavior of a nonlinear dynamical system is the spectrum of Lyapunov exponents. With the help of the Lyapunov exponents, the rate of convergence or divergence of trajectories in the phase space is determined. The presence of at least one positive value in the spectrum of Lyapunov exponents indicates the presence of chaotic oscillations in the system. The number of Lyapunov exponents corresponds to the dimension of the phase

space of the nonlinear dynamical system. For our system (13), the number of such indicators is six. We employ the Benetin algorithm for calculating Lyapunov exponents [40, 41]. Following [42], we calculated the maximum Lyapunov exponent for the system of equations (13) at R = 290 and then used the Gram-Schmidt orthogonalization to determine all Lyapunov exponents more precisely:

$$L_1 = 0.359648, \quad L_2 = 0.11627, \quad L_3 = -0.678572,$$

$$L_4 = -0.835604, \quad L_5 = -1.45011, \quad L_6 = -1.87832. \qquad (14)$$

If the spectrum of Lyapunov exponents has two positive Lyapunov exponents, then the system (13) exhibits hyperchaotic behavior. The maximum Lyapunov exponent of the new hyperchaotic dynamo system (13) corresponds to the value $L_{1max} = 0.359648$. Since $L_1 + L_2 + L_3 + L_4 + L_5 + L_6 = -4.366688 < 0$, the new hyperchaotic dynamo system (13) is dissipative. In addition, the Kaplan-York dimension of the new hyperchaotic dynamo system (13) is calculated using the formula

$$D_{KY} = M + \frac{1}{|L_{M+1}|} \sum_{i=1}^{M} L_i, \qquad (15)$$

where M is the largest integer, for which

$$\sum_{i=1}^{M} L_i > 0, \quad \sum_{i=1}^{M+1} L_i < 0.$$

The Lyapunov dimension of this hyperchaotic attractor using (15) is $D_{KY} \approx 2.7013$. Thus, we conclude that the Lyapunov dimension of the new hyperchaotic system (13) is fractional.
<figure><img src="image_6.png" /><figcaption>Fig. 4. Bifurcation diagrams for the $x_1, x_2, x_3, x_4, x_5, x_6$ components of system (13) depending on the parameter R</figcaption></figure>

### C. Autocorrelation functions

An effective characteristic of a strongly chaotic behavior is the autocorrelation function, which is defined as the average over a certain time interval $T$ (for $T \to \infty$) of

In Fig. 5, plots of autocorrelation functions (diagonal elements) for components $x_i, i = 1, 2, 3, 4, 5, 6$ are shown, where we see that autocorrelation functions tend exponentially to zero. This is one of the criteria for dynamic chaos. On a logarithmic scale, the exponential decay area for the autocorrelation functions $K_{11}, K_{22}, K_{33}, K_{44}, K_{55}, K_{66}$, is approximated by a straight line (see Fig. 5,b,d,f,h,j,l).

## VI. EQUATION OF FINITE AMPLITUDE FOR STATIONARY CONVECTION

We obtain the nonlinear Ginzburg-Landau equation (see Appendix) for the amplitude $A_1$ by applying the

the products $x_i$ taken at times $t$ and $t + \tau$:

$$K_{ij}(\tau) = \lim_{T \to \infty} \frac{1}{T} \int_{0}^{\infty} x_i(t)x_j(t+\tau)d\tau, \quad i = 1 \div 6, \quad j = 1 \div 6.$$

asymptotic method of perturbation theory to the system of equations (6) with respect to a small parameter $\epsilon^2 = (Ra - Ra_c)/Ra_c \ll 1$:

$$A_1 \frac{\partial A_1}{\partial \tau} - A_2 A_1 + A_3 A_1^3 = 0. \quad (16)$$

This equation describes a weakly nonlinear stage of stationary magnetic convection in a nonuniformly rotating electrically conductive fluid, taking into account TM effects. In the limiting case, when there are no TM effects ($q_\alpha = 0, R_\alpha = 0$), the equation (16) corresponds to the known result [32]. The analytical solution of the equation (16) with a known initial condition $A_0 = A_1(0)$ has
the form:

$$A_1(\tau) = \frac{A_0}{\sqrt{\frac{A_3}{A_2} A_0^2 + \left(1 - A_0^2 \frac{A_3}{A_2}\right) \exp\left(-\frac{2\tau A_2}{A_1}\right)}} \quad (17)$$

The amplitude $A_1$ can be used to calculate the value of the generated Y-component (or toroidal) of the magnetic field $\tilde{v} \cong A_1 \cdot \Pi_\alpha$. For numerical calculations $\tilde{v}$ we use the following values of magnetic convection: $b = 8/3$, $Pm = 1$, $Pr = 10$, $Q = 4.45$, $T = 0.1$, $Ta = 1080$, $R_\alpha = 0.02R$, $q_\alpha = 0.044$. The initial amplitude is assumed to be $A_0 = 0.7$, and $Ra_2 \approx Ra_c$, which corresponds to the supercriticality parameter $\epsilon$ smallness. Plots of the magnetic field disturbance $\tilde{v}$ versus time $\tau$ in the absence

(dashed line) and presence (solid line) of TM effects are shown in Fig. 6. Here, we can see that the growth of magnetic field disturbances $\tilde{v}$ continues until a certain time $\tau_0$, after which the magnetic field's stationary level $\tilde{v}_{max}$ is established. The amplitude of the stationary magnetic field with TM effects is greater in this case than without TM effects. Figure 6a-6b shows that with an increase in the Rayleigh parameter $R = 260 \rightarrow R = 320$ for the Keplerian rotation profile $Ro = -3/4$ the level of the stationary magnetic field increases. We observe a similar situation for the Rayleigh rotation profile $Ro = -1$ (see Fig. 6,c,d).

Thus, with an increase in temperature at the lower boundary of the layer, an increase in the stationary level of the generated magnetic field occurs at a fixed difference in the thermo-electromotive force coefficient.

<figure><img src="image_7.png" /><figcaption>Fig. 5. a), c), e), g), i), k) are plots of autocorrelation functions $K_{11}$, $K_{22}$, $K_{33}$, $K_{44}$, $K_{55}$, $K_{66}$ from time $\tau$ for $R = 290$; b), d), f), h), j), l) are a linear dependence of autocorrelation functions $K_{11}$, $K_{22}$, $K_{33}$, $K_{44}$, $K_{55}$, $K_{66}$ on logarithmic scales of the time interval $\tau$ for a strongly chaotic motion</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 6. Dependence of the amplitude $\tilde{v}(\tau)$ of the generated magnetic field on the time $\tau$. With TM effects taken into account (solid line), the steady-state amplitude $\tilde{v}(\tau)$ is larger than in the case without TM effects (dashed line). Plots a)-b) correspond to the Rayleigh parameters $R = 260$ and $R = 320$ at $Ro = -3/4$, and plots c)-d) correspond to Rayleigh parameters $R = 260$ and $R = 320$ for $Ro = -1$</figcaption></figure>

## VII. CONCLUSIONS

In this research, we have derived a novel 6D hyperchaotic system of nonlinear dynamic equations. These equations characterize the generation of a magnetic field through thermomagnetic effects in a nonuniformly rotating layer of an electrically conductive fluid in the presence of a constant vertical magnetic field $B_0$. The derivation of these equations follows an approach similar to the Lorenz equations [36], employing a Fourier series expansion of minimum order to represent the physical fields. However, unlike the Lorenz equations and the classical Rikitake dynamo system [7], our equations involve four quadratic nonlinear terms in a 6D dimensional space. Furthermore, the application of the electromechanical model to problems related to terrestrial magnetism has raised concerns among experts in dynamo theory, as Rikitake's model is more similar to an electric machine than a tool for describing the MHD dynamo process [43]. The results obtained are in good agreement with our previous works [31]-[33] on the chaotic dynamics of nonuniformly rotating magnetic convection, in which thermomagnetic instability was not considered. In this work, for the first time, a theoretical model is constructed to describe the chaotic behavior of the generated magnetic field based on the combined action of the effects of convection, nonuniform rotation, and inhomogeneity of the chemical composition of the electrically conductive medium. The main results of this work are as follows:

* A dynamic analysis of the stability of a new 6D hyperchaotic system is carried out depending on changes in the Rayleigh parameter $R$. For this system, phase portraits are constructed, which show the regimes of regular and chaotic behavior

of magnetic field disturbances for different values of the parameter $R$.

* It was demonstrated that chaotic variations might emerge in the magnitude and direction (inversion) of magnetic field disturbances. Therefore, the resulting new 6D hyperchaotic system can be used as an alternative to the Rikitake equations in describing the inversions of the generated magnetic fields.

* An increase in the temperature difference (parameter $R$) at the boundaries of a layer of a nonuniformly rotating medium leads to an increase in the amplitude of stationary perturbations of the magnetic field.

* The stationary level of magnetic field disturbances increases when TM phenomena are taken into account for both Rayleigh ($Ro = -1$) and Keplerian ($Ro = -3/4$) rotation profiles.

Thus, various modes of regular and chaotic behavior of excited magnetic fields were discovered under certain ratios of dimensionless parameters $b, R, Ro, Ta, T, H, Pm, Pr, G, q_1, q_2$. Given the large number of parameters, it can be assumed that most of the currently known scenarios of the transition to chaos are realized in this system of equations.

The results obtained in this work can be used not only for the problem of terrestrial magnetism but also for other astrophysical objects.

## ACKNOWLEDGMENTS

We thank the anonymous reviewer for their interest in the article and valuable comments.
APPENDIX

Derivation of the Ginzburg–Landau equation from the dynamic system of equations (6)

The Ginzburg–Landau equation is obtained from a nonlinear dynamic system of equations (6) in this section. All perturbed quantities in the equations (6) can be represented as a series expansion in the small supercritical parameter $\epsilon$:

$$ 
\begin{gathered}
\mathbf{X}(\tilde{t}) = \epsilon \mathbf{X}_1 + \epsilon^2 \mathbf{X}_2 + \epsilon^3 \mathbf{X}_3 + \dots, \quad \mathbf{X}(\tilde{t}) = [X, V, U, W, Y, Z]^T \\
R = R_0 + \epsilon^2 R_2 + \dots \tag{18}
\end{gathered}
$$

The amplitudes of the perturbed quantities depend only on the slow time $\tilde{\tau} = \epsilon^2 \tilde{t}$. For the first order in $\epsilon$, after substituting the expansion (VII) into (6), we obtain a linear system of equations:

$$ \mathcal{L} \mathbf{X}_1 = 0, \quad \mathbf{X}_1 = [X_1, V_1, U_1, W_1, Y_1, Z_1]^T, \tag{19} $$

where matrix $\mathcal{L}$ has the form

$$ \mathcal{L} = \begin{bmatrix} -1 & -T & -H & 0 & R_0 & 0 \\ \sqrt{\mathrm{Ta}}(1+\mathrm{Ro}) & -1 & 0 & H & 0 & 0 \\ \mathrm{Pr}^{-1} & 0 & -\mathrm{Pm}^{-1} & 0 & 0 & 0 \\ 0 & -\mathrm{Pr}^{-1} & \mathrm{Ro}\sqrt{\mathrm{Ta}} & -\mathrm{Pm}^{-1} & -\mathrm{G} & 0 \\ 1 & 0 & 0 & -q_1 & -1 & 0 \\ 0 & 0 & 0 & 0 & 0 & -b \end{bmatrix}. $$

The solutions of the equations (19) with respect to the variable $X_1$ are respectively equal to

$$ V_1 = \Xi \cdot X_1, \quad \Xi = \frac{\sqrt{\mathrm{Ta}}(1+\mathrm{Ro})(1-q_1\mathrm{GPm}) + \mathrm{H}\mathrm{Ro}\sqrt{\mathrm{Ta}}\frac{\mathrm{Pm}^2}{\mathrm{Pr}} - \mathrm{GHPm}}{1 + \mathrm{HPm}\mathrm{Pr}^{-1} - q_1\mathrm{GPm}}, $$

$$ U_1 = \frac{\mathrm{Pm}}{\mathrm{Pr}} \cdot X_1, \quad W_1 = -\Pi \cdot X_1, \quad \Pi = \frac{\mathrm{Pm}\mathrm{Pr}^{-1}\sqrt{\mathrm{Ta}}(1+\mathrm{Ro}(1-\mathrm{Pm})) + \mathrm{GPm}}{1 + \mathrm{HPm}\mathrm{Pr}^{-1} - q_1\mathrm{GPm}}, $$

$$ Y_1 = \Lambda \cdot X_1, \quad \Lambda = \frac{1 + \mathrm{HPm}\mathrm{Pr}^{-1} + q_1 \mathrm{Pm}\mathrm{Pr}^{-1} \sqrt{\mathrm{Ta}}(1 + \mathrm{Ro}(1 - \mathrm{Pm}))}{1 + \mathrm{HPm}\mathrm{Pr}^{-1} - q_1 \mathrm{GPm}}, \quad Z_1 = 0. $$

For the second order in $\epsilon$, we have the following equation:

$$ \mathcal{L} \mathbf{X}_2 = [\mathcal{R}_{21}, \mathcal{R}_{22}, \mathcal{R}_{23}, \mathcal{R}_{24}, \mathcal{R}_{25}, \mathcal{R}_{26}]^T, \quad \mathbf{X}_2 = [X_2, V_2, U_2, W_2, Y_2, Z_2]^T, \tag{20} $$

where there are non-linear terms

$$ \mathcal{R}_{21} = 0, \quad \mathcal{R}_{22} = 0, \quad \mathcal{R}_{23} = 0, \quad \mathcal{R}_{24} = 0, \quad \mathcal{R}_{25} = X_1 Z_1 + q_2 W_1 Z_1, \quad \mathcal{R}_{26} = -X_1 Y_1 - q_2 W_1 Y_1. $$

Solutions of equations (20) look like:

$$ V_2 = \Xi \cdot X_2, \quad U_2 = \frac{\mathrm{Pm}}{\mathrm{Pr}} \cdot X_2, \quad W_2 = -\Pi \cdot X_2, $$

$$ Y_2 = X_2 - q_1 W_2, \quad Z_2 = \Lambda(1 - q_2 \cdot \Pi) \frac{X_1^2}{b}. $$

Next, we proceed to equations of the third order in $\epsilon$:

$$ \mathcal{L} \mathbf{X}_3 = [\mathcal{R}_{31}, \mathcal{R}_{32}, \mathcal{R}_{33}, \mathcal{R}_{34}, \mathcal{R}_{35}, \mathcal{R}_{36}]^T, \quad \mathbf{X}_3 = [X_3, V_3, U_3, W_3, Y_3, Z_3]^T, \tag{21} $$
where

$$ \mathcal{R}_{31} = \frac{\partial X_1}{\partial \tilde{\tau}} - R_2 Y_1, \quad \mathcal{R}_{32} = \frac{\partial V_1}{\partial \tilde{\tau}}, \quad \mathcal{R}_{33} = \frac{\partial U_1}{\partial \tilde{\tau}}, \quad \mathcal{R}_{34} = \frac{\partial W_1}{\partial \tilde{\tau}}, $$

$$ \mathcal{R}_{35} = \mathrm{Pr} \frac{\partial Y_1}{\partial \tilde{\tau}} + X_1 Z_2 + X_2 Z_1 + q_2 W_1 Z_2, $$

$$ \mathcal{R}_{36} = \frac{\partial Z_1}{\partial \tilde{\tau}} - \mathrm{Pr}^{-1}(X_1 Y_2 + X_2 Y_1 + q_2 W_1 Y_2 + q_2 W_2 Y_1). $$

The solvability condition (Fredholm's alternative) for nonlinear equations in the third order in $\epsilon$ is as follows:

$$ \sum_{j=1}^{5} \mathcal{R}_{3j} \mathbf{X}_{1}^{\dagger} = 0, \qquad (22) $$

where

$$ \mathcal{R}_{31} = \mathcal{M} \mathcal{K} \mathcal{R}_{31}, \quad \mathcal{R}_{32} = \mathcal{K} R_0 ((\mathcal{R}_{35} - q_1 \mathrm{Pm} \mathcal{R}_{34} - q_1 \mathrm{Pm}^2 \mathrm{Ro} \sqrt{\mathrm{Ta}} \mathcal{R}_{33})(1 + \mathrm{HPmPr}^{-1}) \\ + q_1 \mathrm{PmPr}^{-1} (\mathcal{R}_{32} + \mathrm{HPm} \mathcal{R}_{34} + \mathrm{HPm}^2 \mathrm{Ro} \sqrt{\mathrm{Ta}} \mathcal{R}_{33})), $$

$$ \mathcal{R}_{33} = \mathcal{R}_{36}, \quad \mathcal{R}_{34} = -\mathcal{M} \mathcal{K} \mathcal{H} \mathcal{R}_{36}, $$

$$ \mathcal{R}_{35} = \mathcal{M} \mathrm{T} ((\mathcal{R}_{32} + \mathrm{HPm}\mathcal{R}_{34} + \mathrm{HPm}^2\mathrm{Ro}\sqrt{\mathrm{Ta}}\mathcal{R}_{33})(1 - q_1 \mathrm{GPm}) \\ - \mathrm{GHPm}(\mathcal{R}_{35} - q_1 \mathrm{Pm}\mathcal{R}_{34} - q_1 \mathrm{Pm}^2\mathrm{Ro}\sqrt{\mathrm{Ta}}\mathcal{R}_{33})), $$

$$ \mathcal{M} = (1 - q_1 \mathrm{Pm}^2 \mathrm{Pr}^{-1} \mathrm{Ro} \sqrt{\mathrm{Ta}})(1 + \mathrm{HPmPr}^{-1}) + q_1 \mathrm{PmPr}^{-1} \cdot \mathcal{P}, $$

$$ \mathcal{P} = \sqrt{\mathrm{Ta}} \left( 1 + \mathrm{Ro} + \mathrm{HRo} \frac{\mathrm{Pm}^2}{\mathrm{Pr}} \right), $$

$$ \mathcal{K} = (1 - q_1 \mathrm{Pm}^2 \mathrm{Pr}^{-1} \mathrm{Ro} \sqrt{\mathrm{Ta}}) \mathrm{GHPm} - (1 - q_1 \mathrm{GPm}) \cdot \mathcal{P}. $$

Elements of the matrix $\mathbf{X}_1^\dagger = [X_1^\dagger, Y_1^\dagger, Z_1^\dagger, U_1^\dagger, V_1^\dagger]^T$ are solutions of the linear self-adjoint problem $\mathcal{L}^\dagger \mathbf{X}_1^\dagger = 0$, where the self-adjoint matrix $\mathcal{L}^\dagger$ is defined as

$$ \mathcal{L}^{\dagger} = \begin{bmatrix} -\mathcal{M}\mathcal{K} & R_0\mathcal{M}\mathcal{K} & 0 & -\mathcal{H}\mathcal{M}\mathcal{K} & -\mathcal{T}\mathcal{M}\mathcal{K} \\ R_0\mathcal{M}\mathcal{K} & R_0\mathcal{K}\mathcal{N} & 0 & 0 & 0 \\ 0 & 0 & -\gamma & 0 & 0 \\ -\mathcal{H}\mathcal{M}\mathcal{K} & 0 & 0 & -\mathcal{HPm}^{-1}\mathrm{Pr}\mathcal{K}\mathcal{M} & 0 \\ -\mathcal{T}\mathcal{M}\mathcal{K} & 0 & 0 & 0 & -\mathcal{T}\mathcal{L} \end{bmatrix} $$

where

$$ \mathcal{N} = -(1 + \mathrm{HPmPr}^{-1} - q_1 \mathrm{GPm}), \quad \mathcal{L} = -\mathcal{N}. $$

Finally, we derive the nonlinear Ginzburg-Landau equation for the amplitude $A_1(\tau)$ from the equation (22):

$$ A_1 \frac{\partial A_1}{\partial \tau} - A_2 A_1 + A_3 A_1^3 = 0 \qquad (23) $$
Using the definition for R, T, H, q₁, q₂, G, the coefficients A₁,₂,₃ take the following form:

$$A_1 = \frac{a^2}{\mathrm{Pr}} + \frac{k_c^2 \mathrm{Ra}_c}{a^4} \cdot \frac{\left(1 + q_\alpha \Pi_\alpha \left(1 + \frac{\mathrm{Pm}}{\mathrm{Pr}}\right) - q_\alpha \frac{\pi^2 \mathrm{Pm}^3}{a^4 \mathrm{Pr}^2} \mathrm{Ro} \sqrt{\mathrm{Ta}}\right) (a^4 + \pi^2 \mathrm{Q}) + q_\alpha \frac{\pi \mathrm{Pm}}{a^2 \mathrm{Pr}^2} \cdot \delta}{a^4 + \pi^2 \mathrm{Q} - q_\alpha R_\alpha k_c^2} \\ - \frac{\pi^2 \mathrm{QPm}}{a^2 \mathrm{Pr}} - \frac{\pi \sqrt{\mathrm{Ta}}}{a^6 (a^4 + \pi^2 \mathrm{Q} - k_c^2 q_\alpha R_\alpha) \mathrm{Pr}} \\ \times \left[ \delta \cdot (a^4 - k_c^2 q_\alpha R_\alpha) - \pi k_c^2 a^2 \tilde{\mathrm{Q}} \mathrm{Pr} R_\alpha \left(1 + q_\alpha \Pi_\alpha \left(1 + \frac{\mathrm{Pm}}{\mathrm{Pr}}\right) - q_\alpha \frac{\pi^2 \mathrm{Pm}^3}{a^4 \mathrm{Pr}^2} \mathrm{Ro} \sqrt{\mathrm{Ta}}\right) \right],$$

$$\delta = \frac{\pi a^2 \sqrt{\mathrm{Ta}}}{a^4 + \pi^2 \mathrm{Q} - k_c^2 q_\alpha R_\alpha} \cdot [(1 + \mathrm{Ro})(a^4 - k_c^2 q_\alpha R_\alpha) + \pi^2 \mathrm{QPm}(\mathrm{PmRo} - 1)] \\ - \frac{\pi a^2 k_c^2 R_\alpha \tilde{\mathrm{Q}}(1 + \mathrm{Pm})}{a^4 + \pi^2 \mathrm{Q} - q_\alpha R_\alpha k_c^2} + \frac{\pi^3 \mathrm{Pm}^2}{a^2} \mathrm{QRo}\sqrt{\mathrm{Ta}}, \quad \tilde{\mathrm{Q}} = \mathrm{QPm}^{-1} \mathrm{Pr},$$

$$A_2 = \frac{k_c^2 \mathrm{Ra}_2}{a^2 \mathrm{Pr}} (1 + q_\alpha \Pi_\alpha),$$

$$A_3 = \frac{k_c^4 \mathrm{Ra}_c}{8a^4 \mathrm{Pr}} \cdot \frac{(1 + q_\alpha \Pi_\alpha)(1 - q_\alpha^{(2)} \Pi_\alpha)^2 (a^4 + \pi^2 \mathrm{Q})}{a^4 + \pi^2 \mathrm{Q} - k_c^2 q_\alpha R_\alpha} + R_\alpha \cdot \frac{\pi^2 k_c^4 \mathrm{QPm}^{-1} \sqrt{\mathrm{Ta}} (1 + q_\alpha \Pi_\alpha)(1 - q_\alpha^{(2)} \Pi_\alpha)^2}{8a^4 (a^4 + \pi^2 \mathrm{Q} - k_c^2 q_\alpha R_\alpha)},$$

where

$$\Pi_{\alpha} = \frac{\pi^2 \mathrm{Pm} \mathrm{Pr}^{-1} \sqrt{\mathrm{Ta}} (1 + \mathrm{Ro}(1 - \mathrm{Pm})) + k_c^2 R_{\alpha}}{a^4 + \pi^2 \mathrm{Q} - k_c^2 q_{\alpha} R_{\alpha}}.$$

When deriving the equation (23), we used the rescaled derivative with respect to slow time $\tilde{\tau}$: $\frac{\partial}{\partial \tilde{\tau}} = \frac{1}{a^2} \frac{\partial}{\partial \tau}$.

[1] T. H. Jordan, Proc. Natl. Acad. Sci. USA 76, 4192 (1979); https://doi.org/10.1073/pnas.76.9.4192.

[2] T. Rikitake, Electromagnetism and the Earth's Interior (Elsevier, Amsterdam-London-New York, 1966).

[3] B. M. Yanovsky, Earth Magnetism (Leningrad University Publishing, Leningrad, 1978).

[4] S. I. Braginsky, Geomagn. Aeron. 4, 698 (1964).

[5] R. Hollerbach, Phys. Earth Planet. Inter. 98, 163 (1996); https://doi.org/10.1016/S0031-9201(96)03185-8.

[6] P. H. Roberts, G. A. Glatzmaier, Geophys. Astrophys. Fluid Dynam. 94, 47 (2001); https://doi.org/10.1080/03091920108204131.

[7] T. Rikitake, Proc. Cambridge Philos. Soc. 54, 89 (1958); https://doi.org/10.1017/S0305004100033223.

[8] A. E. Cook, P. H. Roberts, Math. Proc. Camb. Philos. Soc. 68, 547 (1970); https://doi.org/10.1017/S0305004100046338.

[9] Z. Wez, B. Zhu, J. Yang, M. Perc, M. Slavinec, Appl. Math. Comput. 347, 265 (2019); https://doi.org/10.1016/j.amc.2018.10.090.

[10] L. Cristian, T. Binzar, Int. J. Bifurc. Chaos 22, 1250274 (2012); https://doi.org/10.1142/S0218127412502744.

[11] T. Binzar, L. Cristian, Discret. Contin. Dyn. Syst. Ser. B 18, 1755 (2013); https://doi.org/10.3934/dcdsb.2013.18.1755.

[12] I. A. Ilyin, D. S. Noshchenko, A. S. Perezhogin, Bull. KRASEC. Phys. Math. Sci. 2, 43 (2013).

[13] V. I. Potapov, Rus. J. Nonlin. Dyn. 6, 255 (2010); https://doi.org/10.20537/nd1002002.

[14] M. Yu. Reshetnyak, Russ. J. Earth Sci. 15, ES4001 (2015); https://doi.org/10.2205/2015ES000558.

[15] F. Krause, K.H. Rädler, Mean-field Magnetohydrodynamics and Dynamo Theory (Pergamon Press, Oxford, 1980).

[16] A. I. Laptukhov, Geomagn. Aeron. 20, 530 (1980).

[17] A. Schlüter and L. Biermann, Z. Naturforsch. A 5, 65 (1950); https://doi.org/10.1515/zna-1950-0201.

[18] M. I. Kopp, K. N. Kulik, A. V. Tur, V. V. Yanovsky, J. Phys. Stud. 25, 2401 (2021); https://doi.org/10.3097/jps.25.2401.

[19] L. A. Bol'shov, Yu. A. Dreizin, A. M. Dykhne, JETP Lett. 19, 168 (1974).

[20] B. A. Al'terkop, E. V. Mishin, A. A. Rukhadze, JETP Lett. 19, 170 (1974).
[21] D. A. Tidman, R. A. Shanny, Phys. Fluids 17, 1207 (1974).

[22] A. Z. Dolginov, Phys.-Uspekhi 30, 475 (1987); https://doi.org/10.1070/PU1987v030n06ABEH002851.

[23] V. A. Urpin, S. A. Levshakov, D. G. Yakovlev, Mon. Not. Roy. Astron. Soc. 219, 703 (1986); https://doi.org/10.1093/mnras/219.3.703.

[24] V. Urpin, Mon. Not. Roy. Astron. Soc. 472, L5 (2017); https://doi.org/10.1093/mnrasl/slx127.

[25] V. Urpin, Plasma Phys. Rep. 45, 366 (2019); https://doi.org/10.1134/S1063780X19030103.

[26] G. Montani, R. Benini, N. Carlevaro, A. Franco, Mon. Not. Roy. Astron. Soc. 436, 327 (2013); https://doi.org/10.1093/mnras/stt1568.

[27] P. Vadasz, S. Olek, Int. J. Heat Mass Transf. 41, 1417 (1999); 10.1016/S0017-9310(97)00265-2

[28] Vinod K. Gupta, B. S. Bhadauria, I. Hasim, J. Jawdat, A. K. Singh, Alex. Eng. J. 54, 981 (2015); https://doi.org/10.1016/j.aej.2015.09.002.

[29] R. Prasad, A. K. Singh, J. Appl. Fluid Mech. 9, 2887 (2016); https://doi.org/10.29252/jafm.09.06.24811.

[30] L. D. Landau, L. P. Pitaevskii, E. M. Lifshitz, Electrodynamics of Continuous Media (Butterworth-Heinemann, 1984).

[31] M. I. Kopp, A. V. Tour, V. V. Yanovsky, JETP 127, 1173 (2018); https://doi.org/10.1134/S106377611812018x.

[32] M. Kopp, A. Tur, V. Yanovsky, East Eur. J. Phys. 1, 4 (2019); https://doi.org/10.26565/2312-4334-2019-1-01.

[33] M. I. Kopp, A. V. Tur, V. V. Yanovsky, Probl. At. Sci. Technol. 4, 116 (2018); https://arxiv.org/pdf/1805.11894.pdf.

[34] C. A. J. Fletcher, Computational Galerkin Methods (Springer, Berlin-Heidelberg, 1984).

[35] F. V. Dolzhansky, Phys.-Uspekhi 48, 1205 (2005); https://doi.org/10.1070/PU2005v048n12ABEH002375.

[36] E. N. Lorenz, J. Atmos. Sci. 20, 130 (1963); https://doi.org/10.1175/1520-0469(1963)020<0130:DNF>2.0.C0;2.

[37] F. Gantmacher, Lectures in Analytical Mechanics (Mir Publishers, Moscow, 1975).

[38] C. Sparrow, The Lorenz Equations: Bifurcations, Chaos and Strange Attractors (Springer-Verlag, New York, 1982).

[39] V. N. Zharkov, Internal Structure of the Earth and Planets (Gordon and Breach, New York, 1983).

[40] G. Benettin, L. Galgani, A. Giorgilli, J. M. Strelcyn, Meccanica 15, 21 (1980); https://doi.org/10.1007/BF02128237.

[41] A. Wolf, J. B. Swift, H. L. Swinney, J. A. Vastano, Physica D 16, 285 (1985); https://doi.org/10.1016/0167-2789(85)90011-9.

[42] M. Sandri, Mathematica J. 6, 78 (1996); http://www.mathematica-journal.com/issue/v6i3/article/sandri/contents/63sandri.pdf.

[43] H. K. Moffatt, Magnetic Field Generation in Electrically Conducting Fluids (Cambridge University Press, New York, 1978).

ХАОТИЧНА ДИНАМІКА МАГНІТНИХ ПОЛІВ, ІЯКІ ГЕНЕРУЮТЬСЯ ТЕРМОМАГНІТНОЮ НЕСТІЙКІСТЮ В ЕЛЕКТРОПРОВІДНІЙ РІДИНІ, ЩО НЕОДНОРІДНО ОБЕРТАЄТЬСЯ

М. Й. Копш<fnref n="1" />, А. В. Тур<fnref n="3" />, В. В. Яновський<fnref n="1" />,<fnref n="2" />

Досліджено хаотичну поведінку теплової конвекції в електропровідній рідині, що неоднорідно обертається, під дією постійного вертикального магнітного поля В₀. За наявності вертикальних градієнтів температури ∇T₀ і коефіцієнта термо-ЕРС ∇α виникає термомагнітна (ТМ) нестійкість, що приводить до генерації магнітних полів. Магнітне поле В₁, що збуджується дією батареї Бірма-на, спрямоване перпендикулярно до площини векторів ∇T₀, ∇α і градієнта температурних збурень ∇T₁. Це магнітне поле змінює режим теплообміну, а за рахунок ефектів конвективного тепло-обміну та Ригі-Ледюка встановлюється позитивний зворотний зв'язок, що спричиняє збільшення збурень магнітного поля. За допомогою усіченого методу Гальоркіна отримано нелінійну динамі-чну систему рівнянь, що описує генерацію та регенерацію магнітного поля. Чисельний аналіз цих рівнянь показав наявність рег'улярної, квазіперіодичної та хаотичної поведінки збурень магнітно-го поля, що супроводжуються його інверсіями. Застосовуючи метод теорії збурень до нелінійної динамічної системи рівнянь, ми отримали рівняння Гінзбурга-Ландау для слабонелінійної стадії магнітної конвекції, що нерівномірно обертається з урахуванням ефектів ТМ. Розв'язок цього рів-няння показав: стаціонарний рівень магнітних полів, що генеруються, збільшується з урахуванням впливу ТМ-нестійкості.

**Ключові слова:** генерація магнітних полів, термомагнітна нестійкість, хаотична поведінка, рівняння Гінзбурга-Ландау.

<footnote n="1">¹Інститут монокристалів, Національна академія наук України, просп. Науки, 60, Харків, 61072, Україна</footnote>

<footnote n="2">²Харківський національний університет імені В. Н. Каразіна, майдан Свободи, 4, Харків, 61022, Україна</footnote>

<footnote n="3">³Університет Тулузи, CNRS, Інститут досліджень астрофізики та планології</footnote>