Computational Technology for the Basis and Coefficients of Geodynamo Spectral Models in the Maple System

Abstract

Spectral models are often used in the study of geodynamo problems. Physical fields in these models are presented as stationary basic modes combinations with time-dependent amplitudes. To construct a model it is necessary to calculate the modes parameters, and to calculate the model coefficients (the Galerkin coefficients). These coefficients are integrals of complex multiplicative combinations of modes and differential operators. The paper proposes computing technology for the calculation of parameters, the derivation of integrands and the calculation of the integrals themselves. The technology is based on computer algebra methods. The main elements for implementation of technology in the Maple system are described. The proposed computational technology makes it possible to quickly and accurately construct fairly wide classes of new geodynamo spectral models.

1. Introduction

Tackling mathematical problems for geophysical magnetohydrodynamics (MHD) consists in solving MHD equations in a rotating coordinate system in an area that is usually a ball or a spherical shell [1,2,3,4]. Only the equation for magnetic field generation is solved for the kinematic dynamo problem [5]. This problem is linear. If feedback is taken into account in the form of the influence of the Lorentz force on the motions of a fluid, then a nonlinear system arises that describes a self-consistent process of interaction between hydrodynamic and magnetic fields. Accounting for the thermal source of the motion of the medium leads to the equations of magneto-convection [6,7,8]. The exact solutions are unknown for such nonlinear systems, and one has to use simplifications and apply numerical methods [9,10,11,12,13,14,15,16,17]. The spectral methods and models based on them are very popular [2,6,11,18,19,20,21,22,23,24,25,26,27,28]. Solutions to non-stationary problems in these methods are sought in the form of linear combinations of a finite number of stationary basis fields (modes). The coefficients (amplitudes) in these combinations depend on time. The modes themselves are usually orthogonal in the area of solution calculation and, preferably, are selected from the complete system. For such a system, one chooses the eigenmodes of a suitable spectral problem. The dynamic system is compiled for non-stationary amplitudes using the Galerkin method. This system, together with the chosen finite set of basic modes, forms the spectral model of the MHD system.

If the basis mode systems are complete in the respective function spaces, this guarantees convergence to the exact solution. The rate of this convergence needs to be studied separately. However, completeness is not required if only large-scale spatial structures are modelled. Then, the model may have a relatively small (10–20) number of modes. In this case it is sufficient that the system of basic modes has good approximating properties. In general, from a formal point of view, any orthogonal system can be suitable for building models. However, it is desirable that the modes have an obvious physical meaning. Then, the modes should form spatial structures that are most natural for the physical system under study. According to the authors, the most appropriate from this point of view are the eigenmodes of the free oscillations or the free dissipation of the fields of linearized problems. The corresponding spectral problems are usually Hermitian. This ensures the orthogonality of such modes. For the Navier–Stokes equation, the application of the Galerkin method with similar modes was studied in detail in [29].

To apply spectral models, two problems must be solved. Firstly, it is necessary to compose a dynamical system for the amplitudes. Secondly, this system must be solved, often with many variants of the initial conditions. The dynamical system itself is quadratically nonlinear with constant coefficients, which are determined by the control parameters of the model and the set of selected modes. Its numerical solution does not cause any particular difficulties. It is much more difficult to compose the dynamic system itself for mode amplitudes, i.e., to calculate all the Galerkin coefficients. Calculation of the coefficients consists in the calculation of volume integrals over a spherical area from various multiplicative combinations of basic modes and their spatial derivatives. These combinations form very complex algebraic expressions. Of course, various methods of numerical integration can be used to calculate the coefficients, but even the error-free writing of integrands is problematic. Most likely, this represents the main difficulty in using spectral methods. The absolute value of each coefficient is a measure of the modes’ interaction in the physical process, which is described by the corresponding term of the original equation. For this reason, the exact equality of some coefficients to zero can carry significant physical information, making it possible to reveal chains of interacting modes. Therefore, it is important to be able to tell whether the coefficient approximates to zero or is exactly zero. Therefore, for the construction of low-mode models, a computational technology is very useful, with the help of which one can quickly calculate the Galerkin coefficients for various sets of basis modes and determine the zero coefficients. This was the main motivation for this research.

According to the authors, the basis of such technology can be computer algebra methods (symbolic computations) [30,31]. Using these methods, one can obtain error-free expressions for integrands and calculate the integrals themselves, partly analytically and partly numerically. Obviously, for a spherical region, all calculations should be carried out with spherical coordinates and the modes should be expressed using spherical functions. Then, when calculating the Galerkin coefficients, integration over the sphere can be carried out analytically, while numerical integration will occur in the radial direction. Computer algebra systems (CAS) software provides a floating-point computation environment. This allows flexibly control of the accuracy of computations and identification of zero coefficients even with numerical integration. CAS are intensively used in research in various fields of mathematics [32,33,34,35,36,37], physics [31,38,39,40,41,42], engineering [43], education [44], and other areas.

The authors have previously used some of the algorithms described in this article to calculate dynamo spectral models as reported in [23,25,45,46,47]. However, the details for calculating the model coefficients were not discussed in these works. The axis-symmetric case for a kinematic dynamo is described in [48]. A complete description of the developed algorithms for the three-dimensional self-consistent problem of magneto-convection is presented in this paper for the first time.

So, this paper describes the new computational technology for compiling spectral geodynamo models. The technology is based on computer algebra methods. It is clear that the last two items can be performed by standard means using any CAS, since any such system allows performance of both analytical and numerical integration. To implement the first two items, it is necessary to calculate vector analysis operations in an analytical form in spherical coordinates. Therefore, it is convenient to use CAS containing packages of vector analysis operators. The Maple system includes the VectorCalculus package, which implements all standard differential and integral operations of vector analysis in various orthogonal coordinate systems. Therefore, the Maple system was used in the work.

The process consists of the following main stages:

• derivation of explicit analytical expressions for each basic mode;
• derivation of explicit analytical expressions for integrands that are multiplicative combinations of modes and their spatial derivatives;
• analytical integration over the sphere (this integration gives a zero value or some expression depending on the radial variable);
• if the integral over a sphere is non-zero, then analytical or numerical integration are executed over a radial variable.

2. Magnetoconvection Equations

Let us consider a spherical shell of an incompressible electrically conducting fluid (liquid core), which rotates uniformly about the z-axis. In a spherical coordinate system $(r,\theta ,\phi )$, the inner core boundary (ICB) is $r=r_i$, and the core-mantle boundary (CMB) is $r=r_o$. The temperature at the ICB and CMB is constant and equal to $T_i$ and $T_o$, respectively.

The fluid velocity $\mathbf{v}$, the magnetic field $\mathbf{B}$, and the temperature T, are governed by the dimensionless equations of magnetoconvection:

$$\begin{array}{c}{\displaystyle \frac{\partial \mathbf{v}}{\partial t}+\left(\mathbf{v}\nabla \right)\mathbf{v}=-\nabla p+\Delta \mathbf{v}-2{\mathsf{E}}^{-1}{\mathbf{e}}_z\times \mathbf{v}+\mathsf{RaP}{\mathsf{r}}^{-1}T\mathbf{r}+\mathrm{rot}\mathbf{B}\times \mathbf{B},}\hfill \\ {\displaystyle \frac{\partial T}{\partial t}+\left(\mathbf{v}·\nabla \right)\left(T+T_s\right)={\mathsf{Pr}}^{-1}\Delta T,}\hfill \\ {\displaystyle \frac{\partial \mathbf{B}}{\partial t}=\mathrm{rot}\left(\mathbf{v}\times \mathbf{B}\right)+{\mathsf{R}}_{\alpha }\mathrm{rot}\left(\alpha \left(\mathbf{r}\right)\mathbf{B}\right)+{\mathsf{Pm}}^{-1}\Delta \mathbf{B},}\hfill \\ {\displaystyle \nabla \mathbf{v}=0,\nabla \mathbf{B}=0.}\hfill \end{array}$$

(1)

The fluid dynamic equations have to be understood as Boussinesq approximations. The governing dimensionless parameters of the model are the Ekman number $\mathsf{E}$, the Rayleigh number $\mathsf{Ra}$, the Prandtl number $\mathsf{Pr}$, the magnitude of the $\alpha$-effect ${\mathsf{R}}_{\alpha }$, and the magnetic Prandtl number $\mathsf{Pm}$. The liquid core outer radius $r_o$ and $\delta T=T_i-T_o>0$ serve as the length scale and the temperature scale, respectively. To maintain the correct ratio at this length scale, the radius of the inner core $r_i=0.35$.

On the right side of the Navier–Stokes equation (the first Equation (1)), there are three force sources: the Coriolis force $2{\mathsf{E}}^{-1}{\mathbf{e}}_z\times \mathbf{v}$, the buoyancy force $\mathsf{RaP}{\mathsf{r}}^{-1}T\mathbf{r}$, and the Lorentz force $\mathrm{rot}\mathbf{B}\times \mathbf{B}$. On the right side of the induction equation (third Equation (1)), the term $\mathrm{rot}\left(\mathbf{v}\times \mathbf{B}\right)$ describes the field generation by large-scale motions of the fluid, and the term ${\mathsf{R}}_{\alpha }\mathrm{rot}\left(\alpha \left(\mathbf{r}\right)\mathbf{B}\right)$ describes a small-scale turbulent generator ($\alpha$-effect).

In Equation (1), the variable T denotes the temperature deviation from the equilibrium profile

$$\begin{array}{c}{\displaystyle T_s=\frac{r_i}{1-r_i}\left(\frac{1}{r}-1\right)+T_i-1.}\hfill \end{array}$$

(2)

The temperature deviation T is zero at the ICB and CMB. For the velocity, no-slip boundary conditions at the ICB and CMB are applied.

Moreover, it is assumed that the magnetic permeabilities of the outer and inner cores are equal, and that the surrounding of the core is electrically non-conducting. Therefore, the magnetic field must be the potential for $r>1$. Then, $\mathbf{B}\left(r>1\right)={\mathbf{B}}^{\mathrm{out}}\left(\mathbf{r},t\right)$, where ${\mathbf{B}}^{\mathrm{out}}=-\nabla u\left(\mathbf{r},t\right)$ and $\Delta u=0$. The magnetic field must disappear at infinity, so

$$\begin{array}{c}{\displaystyle u=-\sum _{n=1}^{\infty }{r}^{-n-1}\sum _{m=-n}^n{g}_{nm}\left(t\right)Y_n^m(\theta ,\phi ),}\hfill \end{array}$$

(3)

where $Y_n^m(\theta ,\phi )$—spherical harmonics, and $g_{nm}\left(t\right)$—multipole amplitudes.

Turbulence is considered to be isotropic. Therefore, the $\alpha$-effect is used for scalar parametrization in the form $\alpha \left(\mathbf{r}\right)=A\left(r\right)cos\theta$, where $max\left|A\right(r\left)\right|=1$ and $A\left(r\right)=0$ for $r<r_i$ and $r>1$.

It is important to note that convection and magnetic field generation occur only in the liquid core, but that the magnetic field dissipates throughout the entire volume of the core. Therefore, the magnetic field modes must be determined for $0\le r\le 1$, and the velocity and temperature modes only for $r_i\le r\le 1$. Thus, Equation (1) is closed by the boundary conditions:

$$\begin{array}{c}{\displaystyle \mathbf{v}\left(r=r_i\right)=\mathbf{v}\left(r=1\right)=0,}\hfill \\ {\displaystyle T\left(r=r_i\right)=T\left(r=1\right)=0,}\hfill \\ {\displaystyle \mathbf{B}\left(r=0\right)<\infty ,\mathbf{B}(r=1)={\mathbf{B}}^{out}(r=1).}\hfill \end{array}$$

(4)

3. General Form of Spectral Model Equations

To obtain the general form of the spectral model equations, assume that the field components are approximately represented by linear decompositions of a finite number of modes:

$$\begin{array}{c}{\displaystyle \mathbf{v}\left(\mathbf{r},t\right)=\sum _{l=1}^{L_{max}}{\beta }_l\left(t\right){\mathbf{v}}_l\left(\mathbf{r}\right),T\left(\mathbf{r},t\right)=\sum _{s=1}^{S_{max}}{\alpha }_s\left(t\right)T_s\left(\mathbf{r}\right),\mathbf{B}\left(\mathbf{r},t\right)=\sum _{p=1}^{P_{max}}{\gamma }_p\left(t\right){\mathbf{B}}_p\left(\mathbf{r}\right).}\hfill \end{array}$$

(5)

All these modes are free decay eigenmodes of the corresponding fields, i.e., solutions of spectral problems

$$\begin{array}{c}{\displaystyle \mu \mathbf{v}-\nabla p+\Delta \mathbf{v}=0,\nabla \mathbf{v}=0,}\hfill \\ {\displaystyle \lambda T+\Delta T=0,}\hfill \\ {\displaystyle \eta \mathbf{B}+\Delta \mathbf{B}=0,\nabla \mathbf{B}=0}\hfill \end{array}$$

(6)

with boundary conditions (4).

Each of these spectral problems describes the free decay of the corresponding field when there are no sources. This decay for the velocity field is related to the viscosity of the medium, for the magnetic field, it is related to the finite conductivity, and, for the temperature, it is related to the transfer of heat due to heat conduction. The eigenvalues determine the characteristic decay times of the modes—the larger the eigenvalue, the faster the mode decays.

Explicit expressions for these modes and how their parameters and eigenvalues can be determined using CAS will be considered below.

For now, only note that all these problems are Hermitian with respect to inner products:

$$\begin{array}{c}{\displaystyle {〈{\mathbf{v}}_1,{\mathbf{v}}_2〉}_v=\underset{r_i\le r\le 1}{\int \int \int }{\mathbf{v}}_1{\mathbf{v}}_{2}dV={\int }_{r_i}^1{r}^{2}dr{\int }_0^{\pi }sin\theta d\theta {\int }_{-\pi }^{\pi }{\mathbf{v}}_1{\mathbf{v}}_{2}d\phi ,}\hfill \\ {\displaystyle {〈T_1,T_2〉}_T=\underset{r_i\le r\le 1}{\int \int \int }{T}_1{T}_{2}dV={\int }_{r_i}^1{r}^{2}dr{\int }_0^{\pi }sin\theta d\theta {\int }_{-\pi }^{\pi }{T}_1{T}_{2}d\phi ,}\hfill \\ {\displaystyle {〈{\mathbf{B}}_1,{\mathbf{B}}_2〉}_B=\underset{0\le r\le 1}{\int \int \int }{\mathbf{B}}_1{\mathbf{B}}_{2}dV={\int }_0^1{r}^{2}dr{\int }_0^{\pi }sin\theta d\theta {\int }_{-\pi }^{\pi }{\mathbf{B}}_1{\mathbf{B}}_{2}d\phi ,}\hfill \end{array}$$

(7)

and their eigenmodes systems are full [29]. Therefore, the modes are orthogonal with respect to these inner products. Moreover, they are assumed to be normalized.

The decompositions (5) are substituted into the Equation (1) and the Galerkin procedure is applied. The result is a dynamic system for amplitudes:

$$\begin{array}{c}{\displaystyle \frac{d{\beta }_l}{dt}=\sum _{i,j=1}^{L_{max}}{B}_{lij}{\beta }_i{\beta }_j-{\mu }_l{\beta }_l+{\mathsf{E}}^{-1}\sum _{i=1}^{L_{max}}{E}_{li}{\beta }_i+\mathsf{RaP}{\mathsf{r}}^{-1}\sum _{i=1}^{S_{max}}{C}_{li}{\alpha }_i+\sum _{i,j=1}^{P_{max}}{Q}_{lij}{\gamma }_i{\gamma }_j,}\hfill \\ {\displaystyle l=1,\dots ,L_{max},}\hfill \\ {\displaystyle \frac{d{\alpha }_s}{dt}=\sum _{i,j=1}^{L_{max},S_{max}}{F}_{sij}{\beta }_i{\alpha }_j+\sum _{i=1}^{L_{max}}{H}_{si}{\beta }_i-{\mathsf{Pr}}^{-1}{\lambda }_s{\alpha }_s,}\hfill \\ {\displaystyle s=1,\dots ,S_{max},}\hfill \\ {\displaystyle \frac{d{\gamma }_p}{dt}=\sum _{i,j=1}^{L_{max},P_{max}}{W}_{pij}{\beta }_i{\gamma }_j+{\mathsf{R}}_{\alpha }\sum _{i=1}^{P_{max}}{W}_{pij}^{\alpha }{\gamma }_i-{\mathsf{Pm}}^{-1}{\eta }_p{\gamma }_p,}\hfill \\ {\displaystyle p=1,\dots ,P_{max}.}\hfill \end{array}$$

(8)

In this system, ${\mu }_l>0$, ${\lambda }_s>0$ and ${\eta }_p>0$ are the eigenvalues of the corresponding modes, and the remaining coefficients (Galerkin coefficients) are determined by the inner products, which are calculated by the Formula (7):

$$\begin{array}{c}{\displaystyle B_{lij}=-{\langle \left({\mathbf{v}}_i\nabla \right){\mathbf{v}}_j,{\mathbf{v}}_l\rangle }_v,E_{li}=-2{\langle {\mathbf{e}}_z\times {\mathbf{v}}_i,{\mathbf{v}}_l\rangle }_v,C_{li}={\langle T_i\mathbf{r},{\mathbf{v}}_l\rangle }_v,}\hfill \\ {\displaystyle Q_{lij}={\langle \mathrm{rot}{\mathbf{B}}_i\times {\mathbf{B}}_j,{\mathbf{v}}_l\rangle }_v,F_{sij}=-{\langle {\mathbf{v}}_i\left(\nabla T_j\right),T_s\rangle }_T,H_{si}=\frac{r_i}{1-r_i}{\langle r^{-2}{\mathbf{v}}_i{\mathbf{e}}_r,T_s\rangle }_T,}\hfill \\ {\displaystyle W_{pij}={\langle \mathrm{rot}\left({\mathbf{v}}_i\times {\mathbf{B}}_j\right),{\mathbf{B}}_s\rangle }_v,W_{pi}^{\alpha }={\langle \mathrm{rot}\left(A\left(r\right)cos\theta {\mathbf{B}}_i\right),{\mathbf{B}}_s\rangle }_v.}\hfill \end{array}$$

(9)

It is necessary to provide an explanation for the coefficients $W_{pij}$ and $W_{pi}^{\alpha }$. All Galerkin coefficients for the third equation of the system (1) are calculated using the inner product ${\langle ·,·\rangle }_B$ where integration over the radial variable goes over the segment $[0;1]$. But the first factor in $W_{pij}$ contains the ${\mathbf{v}}_i$ velocity modes, which are zero for $r<r_i$. Therefore, in fact, integration will only be for $r_i\le r\le 1$, i.e., as in the inner product ${\langle ·,·\rangle }_v$. The situation is similar with $W_{pi}^{\alpha }$ because $A\left(r<r_i\right)=0$.

The system (8) does not have a term corresponding to the pressure field. The reason is that the Galerkin coefficients for terms with pressure are zero. This is a well-known result [29], but the derivation is described in Appendix A.1.

From the well-known triple scalar product property $\left({\mathbf{e}}_z\times {\mathbf{v}}_i\right){\mathbf{v}}_l=\left({\mathbf{e}}_z\times {\mathbf{v}}_l\right){\mathbf{v}}_i$, it follows that $E_{li}=-E_{il}$. In particular, $E_{ll}=0$.

It is necessary to note the following well-known equality [29]

$$\begin{array}{c}{\displaystyle B_{lij}+B_{lji}+B_{ilj}+B_{ijl}+B_{jil}+B_{jli}=0.}\hfill \end{array}$$

(10)

It is usually obtained from physical considerations as a consequence of the law of conservation of kinetic energy for an ideal fluid. The mathematical derivation of these relations, using only the boundary conditions for the velocity components and their divergence-freeness, is given in Appendix A.2. Note that, from (10), it follows that $B_{lll}=0$.

The coefficients $B_{lij}$ and $B_{lji}$ enter the equations of the spectral model (8) symmetrically, although they are defined as Galerkin coefficients by expressions asymmetric in i and j. However, based on the formula (A8) from Appendix A.2, they can be redefined in the following symmetrical way:

$$\begin{array}{c}{\displaystyle B_{lij}=B_{lji}=\frac{{\langle {\mathbf{v}}_i\times \mathrm{rot}{\mathbf{v}}_j+{\mathbf{v}}_j\times \mathrm{rot}{\mathbf{v}}_i,{\mathbf{v}}_l\rangle }_v}{2}.}\hfill \end{array}$$

(11)

The equality (10) will be preserved after such a redefinition.

The coefficients $Q_{lij}$ and $Q_{lji}$ also enter the equations of the spectral model (8) in a symmetrical way, although they are defined by expressions that are asymmetric in i and j. For their symmetrical redefinition, the relationship (A15) between the coefficients can be used from Appendix A.3:

$$\begin{array}{c}{\displaystyle Q_{lij}+Q_{lji}+W_{ilj}+W_{jli}=0.}\hfill \end{array}$$

(12)

Then, the Lorenz terms coefficients can be redefined in a way that is symmetric in i and j:

$$\begin{array}{c}{\displaystyle Q_{lij}=Q_{lji}=-\frac{W_{ilj}+W_{jli}}{2}.}\hfill \end{array}$$

(13)

Another advantage of this definition is that there is no need to separately calculate the integrals for the coefficients of the Lorentz terms. Note that, for spectral models of the type (8), the equality (12) expresses the total energy conservation law in ideal magnetohydrodynamics.

So, in what follows, it is assumed that the coefficients $B_{lij}$ and $Q_{lij}$ are defined by the formulas (11) and (13), respectively.

There is another important relationship between the coefficients (Appendix A.4)

$$\begin{array}{c}{\displaystyle F_{sij}=-F_{jis}\mathrm{and}{F}_{sis}=0.}\hfill \end{array}$$

(14)

4. Calculation of Basis Modes

Let us consider the calculation of the basis modes, i.e., the solutions of spectral problems (6) with boundary conditions (4). The general form of these solutions is well known, but it is necessary to obtain the numerical values of the parameters, normalization coefficients, and eigenvalues. Following the general idea of the work, these problems are solved using the CAS.

The dependence of the modes on $\theta$ and $\phi$ is expressed in terms of the spherical harmonics $Y_n^m(\theta ,\phi )$. The harmonics are assumed to be normalized on the sphere surface:

$$\begin{array}{c}{\displaystyle {\int }_0^{\pi }sin\theta d\theta {\int }_{-\pi }^{\pi }{\left(Y_n^m(\theta ,\phi )\right)}^{2}d\phi =1.}\hfill \end{array}$$

(15)

Then, it follows from the properties of the spherical functions that

$$\begin{array}{c}{\displaystyle {\int }_0^{\pi }sin\theta d\theta {\int }_{-\pi }^{\pi }\left[{\left(\frac{\partial Y_n^m}{\partial \theta }\right)}^2+\frac{1}{{sin}^2\theta }{\left(\frac{\partial Y_n^m}{\partial \phi }\right)}^2\right]d\phi =n\left(n+1\right).}\hfill \end{array}$$

(16)

These relations lead to normalization conditions for radial functions, which are presented below.

4.1. Velocity Modes

The problem (4,6) for velocity is solved separately in the class of toroidal and poloidal fields.

Toroidal eigenmodes have the form:

$$\begin{array}{c}{\displaystyle {\mathbf{v}}_{knm}^T=\mathrm{rot}\left(R_{kn}^T\left(r\right)Y_n^m(\theta ,\phi )\mathbf{r}\right)=R_{kn}^T\left(r\right)\left(\frac{1}{sin\theta }{\mathbf{e}}_{\theta }\frac{\partial }{\partial \phi }-{\mathbf{e}}_{\phi }\frac{\partial }{\partial \theta }\right)Y_n^m(\theta ,\phi ),}\hfill \\ {\displaystyle k=0,1,2,\dots ,n=1,2,3,\dots ,m=-n,\dots ,n.}\hfill \end{array}$$

(17)

The radial functions

$$R_{kn}^T\left(r\right)=A_{kn}^T{j}_n\left(\sqrt{{\mu }_{kn}^T}r\right)+B_{kn}^T{y}_n\left(\sqrt{{\mu }_{kn}^T}r\right),$$

(18)

where $j_n(·)$ and $y_n(·)$ are the spherical Bessel functions of the 1st and 2nd kind, and ${\mu }_{kn}^T$ are the eigenvalues. The coefficients $A_{kn}^T$ and $B_{kn}^T$ must provide the normalization condition

$$\begin{array}{c}{\displaystyle n(n+1){\int }_{r_i}^1{r}^2{\left[R_{kn}^T\left(r\right)\right]}^{2}dr=1,}\hfill \end{array}$$

(19)

and the boundary conditions

$$\begin{array}{c}{\displaystyle R_{kn}^T\left(r_i\right)=R_{kn}^T\left(1\right)=0.}\hfill \end{array}$$

(20)

The equalities (19) and (20) ensure that the ${\mathbf{v}}_{knm}^T$ modes are normalized with respect to the inner product ${\langle ·,·\rangle }_v$ and they satisfy the boundary conditions (4).

It can be seen from the expressions (18) that the boundary conditions (20) are a system of two homogeneous linear equations with respect to the coefficients $A_{kn}^T$ and $B_{kn}^T$. For the existence of a non-zero solution of this system, it is necessary that the determinant of its matrix be zero. The need for a zero determinant provides equations for the eigenvalues:

$$\begin{array}{c}{\displaystyle j_n\left(\sqrt{\mu }{r}_i\right)y_n\left(\sqrt{\mu }\right)-j_n\left(\sqrt{\mu }\right)y_n\left(\sqrt{\mu }{r}_i\right)=0.}\hfill \end{array}$$

(21)

For each $n=1,2,\dots$, this equation has a countable set of solutions ${\mu }_{kn}^T$, numbered by index $k=0,1,2,\dots$.

Each such eigenvalue makes the equations of the system (20) linearly dependent. Therefore, the coefficients $A_{kn}^T$ and $B_{kn}^T$ can only be determined from the second equation. It is clear from (18) that one of the solutions to this equation is $A_{kn}^T=y_n\left(\sqrt{{\mu }_{kn}^T}\right)$ and $B_{kn}^T=-j_n\left(\sqrt{{\mu }_{kn}^T}\right)$. Such values of the coefficients will ensure the fulfillment of the boundary condition. However, they will still need to be renormalized to ensure that the conditions (19) are met.

So, to calculate the modes, it is necessary to find the eigenvalues ${\mu }_{kn}^T$ and the coefficients $A_{kn}^T$ and $B_{kn}^T$. Let us consider how this problem is solved by means of Maple. To do this, the maximum values N and K for the indices n and k are chosen. Then the following calculations are organized in Maple:

• Arrays of explicit expressions for $j_n\left(z\right)$ and $y_n\left(z\right)$ are formed, based on the general recurrence relations for spherical Bessel functions $w_n\left(z\right)$ and explicit expressions for functions of the 1st and 2nd kind [49]:

  $$\begin{array}{c}{\displaystyle w_{n+1}\left(z\right)+w_{n-1}\left(z\right)=\frac{2n+1}{z}{w}_n\left(z\right),}\hfill \\ {\displaystyle j_0\left(z\right)=\frac{sinz}{z},j_1\left(z\right)=\frac{sinz}{z^2}-\frac{cosz}{z},}\hfill \\ {\displaystyle y_0\left(z\right)=-\frac{cosz}{z},y_1\left(z\right)=-\frac{sinz}{z}-\frac{cosz}{z^2}.}\hfill \end{array}$$

  (22)
  Here is the corresponding Maple code:
  j:=array(0..N,[]): y:=array(0..N,[]):
  j[0]:=sin(z)/z: j[1]:=sin(z)/z^2-cos(z)/z:
  y[0]:=-cos(z)/z: y[1]:=-sin(z)/z-cos(z)/z^2:
  for n from 1 to N-1 do
    j[n+1]:=simplify((2∗n+1)/z∗j[n]-j[n-1]):
    y[n+1]:=simplify((2∗n+1)/z∗y[n]-y[n-1]):
  end do:
• The normalization expression (19) is calculated and written to the variable $U_{\mathrm{symbol}}$ in symbolic form with undefined values of the coefficients, eigenvalues, and the integration limit $r_i$:
  assume(r_i>0):
  RT_:=AT_∗eval(j[n],z=sqrt(mu_)∗r)+BT_∗eval(y[n],z=sqrt(mu_)∗r):
  Usymbol:=n∗(n+1)∗int(r^2∗RT ^2,r=r_i..1):
• The left side of the Equation (21) is formed and its first $K+1$ roots ${\mu }_{kn}^T$ are numerically determined. The roots are first separated with a step of 1 and then refined with the required degree of accuracy. Note that, Maple allows you to control the precision of floating point calculations. The precision of calculations is determined by the Digits system variable. The Maple code:
  LPEq:=eval(j[n],z=sqrt(mu_)∗.35)∗eval(y[n],z=sqrt(mu_))-
              eval(j[n],z=sqrt(mu_))∗eval(y[n],z=sqrt(mu_)∗.35):
  mu_left:=1.e-03:
  mu_right:=mu_left+1.:
  LPEq_left:=eval(LPEq,mu_=mu_left):
  k:=0:
  while (k<=K) do
    LPEq_right:=eval(LPEq,mu_=mu_right):
    if ((LPEq_left∗LPEq_right)<=0) then
      muT[k,n]:=fsolve(LPEq=0,mu_=mu_left..mu_right):
      k:=k+1:
    end if:
    mu_left=mu_right:
    mu_right:=mu_left+1.:
    LPEq_left:=LPEq_right:
  end do:
• Assignments are made:

  $$\begin{array}{c}{\displaystyle A_{kn}^T:=y_n\left(\sqrt{{\mu }_{kn}^T}\right),B_{kn}^T:=-j_n\left(\sqrt{{\mu }_{kn}^T}\right).}\hfill \end{array}$$

  (23)
  Such values of the coefficients will ensure the fulfillment of the boundary conditions.
  for k from 0 to K do
    AT[k,n]:=eval(y[n],z=sqrt(muT[k,n])):
    BT[k,n]:=-eval(j[n],z=sqrt(muT[k,n])):
  end do:
• The numeric value $U_{\mathrm{number}}$ is calculated by substituting the numeric values $A_{kn}^T$, $B_{kn}^T$, ${\mu }_{kn}^T$ and $r_i$ into $U_{\mathrm{symbol}}$. Then renormalization is performed using assignments:

  $$\begin{array}{c}{\displaystyle A_{kn}^T:=A_{kn}^T/\sqrt{U_{\mathrm{number}}},B_{kn}^T:=B_{kn}^T/\sqrt{U_{\mathrm{number}}}.}\hfill \end{array}$$

  (24)
  for k from 0 to K do
    Unumber:=eval(Usymbol,
                  [AT_=AT[k,n],BT_=BT[k,n],mu_=muT[k,n],r_i=.35]):
    AT[k,n]:=AT[k,n]/sqrt(Unumber):
    BT[k,n]:=BT[k,n]/sqrt(Unumber):
  end do:

The calculations from points (2)–(5) must be performed in a loop for $n=1,\dots ,N$.

For readers unfamiliar with the Maple language, let us explain that the procedure simplify(expr) performs a simplification of the expression expr; the procedure eval(expr,[x1=expr1,x2=expr2,...]) performs a substitution into the expression expr instead of the variables x1, x2,... of the expressions expr1, expr2,...; the procedure fsolve(F(x)=0, x=a..b) numerically finds the root of the equation $F\left(x\right)=0$ localized on the interval $(a;b)$; and the procedure int(expr,x=a..b) performs analytical integration of the expression expr over the variable x in the range from a to b.

The found values $A_{kn}^T$, $B_{kn}^T$, ${\mu }_{kn}^T$, together with the expressions (17,18), determine the basic toroidal velocity modes.

In Figure 1, graphs of $R_{kn}^T$ functions for $n=2,5$ and $k=0,\dots ,4$ are given. It can be seen that the functions are zero at $r_i=0.35$ and $r_o=1$ in accordance with the conditions (20). Also note that the functions $R_{kn}^T$ have k zero values on the interval $(r_i;r_o)$. This indirectly confirms the correctness of the calculation of ${\mu }_{kn}^T$.

Now, let us find the poloidal velocity eigenmodes. They have a general form

$$\begin{array}{c}{\displaystyle {\mathbf{v}}_{knm}^P=\mathrm{rot}\mathrm{rot}\left(R_{kn}^P\left(r\right)Y_n^m(\theta ,\phi )\mathbf{r}\right)=n(n+1)\frac{R_{kn}^P\left(r\right)}{r}{Y}_n^m(\theta ,\phi ){\mathbf{e}}_r+}\hfill \\ {\displaystyle +\left(\frac{d{R}_{kn}^P\left(r\right)}{dr}+\right)\frac{R_{kn}^P\left(r\right)}{r}\left({\mathbf{e}}_{\theta }\frac{\partial }{\partial \theta }+\frac{1}{sin\theta }{\mathbf{e}}_{\phi }\frac{\partial }{\partial \phi }\right)Y_n^m(\theta ,\phi ),}\hfill \\ k=0,1,2,\dots , n=1,2,3,\dots , m=-n,\dots ,n,\hfill \end{array}$$

(25)

where radial functions

$${\displaystyle R_{kn}^P\left(r\right)=C_{kn}^1{j}_n\left(\sqrt{{\mu }_{kn}^P}r\right)+C_{kn}^2{y}_n\left(\sqrt{{\mu }_{kn}^P}r\right)+C_{kn}^3{r}^n+C_{kn}^4{r}^{-n-1}.}$$

(26)

Comparing the expression (25) with the expressions (15) and (16) shows that the coefficients $C_{kn}^i$ must provide the normalization condition

$$\begin{array}{c}{\displaystyle n^2{(n+1)}^2{\int }_{r_i}^1{\left[R_{kn}^P\left(r\right)\right]}^{2}dr+n(n+1){\int }_{r_i}^1{r}^2{\left[\frac{d{R}_{kn}^P\left(r\right)}{dr}+\frac{R_{kn}^P\left(r\right)}{r}\right]}^{2}dr=1.}\hfill \end{array}$$

(27)

This ensures that the ${\mathbf{v}}_{knm}^P$ eigenmodes are normalized with respect to the inner product ${\langle ·,·\rangle }_v$.

It is also necessary that such conditions be satisfied for the radial functions $R_{kn}^P\left(r\right)$:

$$\begin{array}{c}{\displaystyle R_{kn}^P\left(r_i\right)=R_{kn}^P\left(1\right)=\frac{d{R}_{kn}^P}{dr}{|}_{r=r_i}=\frac{d{R}_{kn}^P}{dr}{|}_{r=1}=0.}\hfill \end{array}$$

(28)

This will enforce the boundary conditions (4).

The conditions (28) lead to a homogeneous linear system of equations with respect to the coefficients $C_{kn}^i$. The condition for the existence of non-zero solutions of this system gives equations for the eigenvalues in the form of the determinant being equal to zero:

$${\displaystyle \left|\begin{array}{cccc}{j}_n\left(\sqrt{\mu }{r}_i\right)& y_n\left(\sqrt{\mu }{r}_i\right)& r_i^n& r_i^{-n-1}\\ j_n\left(\sqrt{\mu }\right)& y_n\left(\sqrt{\mu }\right)& 1& 1\\ {\displaystyle \sqrt{\mu }{j}_n^{\prime }\left(\sqrt{\mu }{r}_i\right)}& {\displaystyle \sqrt{\mu }{y}_n^{\prime }\left(\sqrt{\mu }{r}_i\right)}& n{r}_i^{n-1}& (-n-1)r_i^{-n-2}\\ {\displaystyle \sqrt{\mu }{j}_n^{\prime }\left(\sqrt{\mu }\right)}& {\displaystyle \sqrt{\mu }{y}_n^{\prime }\left(\sqrt{\mu }\right)}& n& -n-1\end{array}\right|=0.}$$

(29)

For each $n=1,2,\dots$, this equation has a countable set of solutions ${\mu }_{kn}^P$, numbered by index $k=0,1,2,\dots$.

Using Maple, let us find the eigenvalues ${\mu }_{kn}^P$ and the coefficients $C_{kn}^i$.

• The package of linear algebra operations is connected.
  with(LinearAlgebra):
• The normalization expression (27) is calculated and written to the variable $U_{\mathrm{symbol}}$ in symbolic form with undefined values of the coefficients, eigenvalues, and the integration limit $r_i$.
  RP_:=C1_∗eval(j[n],z=sqrt(mu_)∗r)+C2_∗eval(y[n],z=sqrt(mu_)∗r)+
      C3_∗r^n+C4_∗r^(-n-1):
  Usymbol:=n^2∗(n+1)^2∗int(RP_^2,r=r_i..1)+
           n∗(n+1)∗int((diff(RP_,r)+RP_/r)^2∗r^2,r=r_i..1):
• The matrix G is formed for the determinant from the Equation (29).
  G:=Matrix(4,4):
  G[1,1]:=eval(j[n],z=sqrt(mu_)∗.35): G[1,2]:=eval(y[n],z=sqrt(mu_)∗.35):
  G[1,3]:=.35^n: G[1,4]:=.35^(-n-1):
  G[2,1]:=eval(j[n],z=sqrt(mu_)): G[2,2]:=eval(y[n],z=sqrt(mu_)):
  G[2,3]:=1: G[2,4]:=1:
  G[3,1]:=sqrt(mu_)∗eval(diff(j[n],z),z=sqrt(mu_)∗.35):
  G[3,2]:=sqrt(mu_)∗eval(diff(y[n],z),z=sqrt(mu_)∗.35):
  G[3,3]:=n∗.35^(n-1): G[3,4]:=(-n-1)∗.35^(-n-2):
  G[4,1]:=sqrt(mu_)∗eval(diff(j[n],z),z=sqrt(mu_)):
  G[4,2]:=sqrt(mu_)∗eval(diff(y[n],z),z=sqrt(mu_)):
  G[4,3]:=n: G[4,4]:=-n-1:
• The left part of the Equation (29) is formed and its first $K+1$ roots ${\mu }_{kn}^P$ are numerically determined. The roots are first separated with a step of 1 and then refined with the required degree of accuracy.
  LPEq:=Determinant(G):
  k:=0:
  mu_left:=1e-06:
  mu_right:=mu_left+.01:
  while (k<=K) do
    if (eval(LPEq,mu_=left_mu)∗eval(LPEq,mu_=right_mu)<=0) then
      muP[k,n]:=fsolve(LPEq=0,mu_=left_mu..right_mu):
      k:=k+1:
    end if:
    left_mu=right_mu:
    right_mu:=left_mu+.01:
  end do:
• $C_{kn}^4:=1$ is assigned, and the remaining coefficients are determined as solutions of the linear system:

  $${\displaystyle \left[\begin{array}{ccc}{j}_n\left(\sqrt{{\mu }_{kn}^P}{r}_i\right)& y_n\left(\sqrt{{\mu }_{kn}^P}{r}_i\right)& r_i^n\\ j_n\left(\sqrt{{\mu }_{kn}^P}\right)& y_n\left(\sqrt{{\mu }_{kn}^P}\right)& 1\\ {\displaystyle \sqrt{\mu }{j}_n^{\prime }\left(\sqrt{\mu }\right)}& {\displaystyle \sqrt{\mu }{y}_n^{\prime }\left(\sqrt{\mu }\right)}& n\end{array}\right]·\left[\begin{array}{c}{C}_{kn}^1\\ C_{kn}^2\\ C_{kn}^3\end{array}\right]=\left[\begin{array}{c}-r_i^{-n-1}\\ -1\\ n+1\end{array}\right]}$$

  (30)
  Such values of the coefficients will ensure the fulfillment of the boundary conditions (28).
  for k from 0 to K do
    C4_[k,n]:=1.:
    L:=SubMatrix(eval(G,[mu_=muP[k,n]]),[1..2,4],[1..3]):
    Q:=SubMatrix(eval(G,[mu_=muP[k,n]]),[1..2,4],[4]):
    S:=LinearSolve(L,-Q):
    C1[k,n]:=S[1,1]: C2[k,n]:=S[2,1]:
    C3[k,n]:=S[3,1]: C4[k,n]:=S[4,1]:
  end do:
• The numeric value $U_{\mathrm{number}}$ is calculated by substituting into $U_{\mathrm{symbol}}$ the numeric values $C_{kn}^i$ and ${\mu }_{kn}^P$ and $r_i$. Then renormalization is performed using assignments:

  $$\begin{array}{c}{\displaystyle C_{kn}^i:=C_{kn}^i/\sqrt{U_{\mathrm{number}}}.}\hfill \end{array}$$

  (31)
  for k from 0 to K do
    Unumber:=eval(Usymbol,
             [C1_=C1[k,n],C2_=C2[k,n],C3_=C3[k,n],C4_=C4[k,n],
              mu_=muP[k,n],r_i=.35]):
    C1[k,n]:=C1[k,n]/sqrt(Unumber): C2[k,n]:=C2[k,n]/sqrt(Unumber):
    C3[k,n]:=C3[k,n]/sqrt(Unumber): C4[k,n]:=C4[k,n]/sqrt(Unumber):
  end do:

The calculations from points (2)–(6) must be performed in a loop for $n=1,\dots ,N$.

The found values of $C_{kn}^i$ and ${\mu }_{kn}^P$, together with the expressions (25,26), determine the basis poloidal velocity modes.

In Figure 2, graphs of $R_{kn}^P$ functions for $n=1,4$ and $k=0,\dots ,4$ are given. It can be seen that both the functions themselves and their derivatives are zero at $r_i=0.35$ and $r_o=1$ in accordance with the conditions (28). Also note that, as in the toroidal case, these functions are zero on the interval $(r_i;r_o)$ at k distinct points. This indirectly confirms the correctness of the calculation of ${\mu }_{kn}^P$.

4.2. Magnetic Modes

The problems (4) and (6) for the magnetic field are also solved separately in the class of toroidal and poloidal fields.

Toroidal solutions of the spectral problem (6) for a magnetic field have the form

$$\begin{array}{c}{\displaystyle {\mathbf{B}}_{knm}^T=\mathrm{rot}\left(X_{kn}^T\left(r\right)Y_n^m(\theta ,\phi )\mathbf{r}\right)=X_{kn}^T\left(\frac{1}{sin\theta }{\mathbf{e}}_{\theta }\frac{\partial }{\partial \phi }-{\mathbf{e}}_{\phi }\frac{\partial }{\partial \theta }\right)Y_n^m(\theta ,\phi ),}\hfill \\ {\displaystyle X_{kn}^T\left(r\right)=a_{kn}^T{j}_n\left(\sqrt{{\eta }_{kn}^T}r\right)+b_{kn}^T{y}_n\left(\sqrt{{\eta }_{kn}^T}r\right),}\hfill \\ k=0,1,2,\dots ,n=1,2,3,\dots ,m=-n,\dots ,n,\hfill \end{array}$$

(32)

where eigenvalues ${\eta }_{kn}^T>0$.

The boundary condition $\mathbf{B}(r=0)<\infty$ from (4) leads to $b_{kn}^T=0$ because $y_n\left(0\right)=\infty$. Consider now the boundary condition for $r=1$. First of all, it must be noted that the field ${\mathbf{B}}^{\mathrm{out}}=-\nabla u$ can be represented as

$$\begin{array}{c}{\displaystyle {\mathbf{B}}^{\mathrm{out}}=\sum _{n=1}^{\infty }\sum _{m=-n}^n{g}_{nm}\left(t\right){\mathbf{B}}_{nm}^{\mathrm{out}}\left(\mathbf{r}\right),\mathrm{where}{\mathbf{B}}_{nm}^{\mathrm{out}}\left(\mathbf{r}\right)=\mathrm{rot}\mathrm{rot}\left(\frac{r^{-n-1}}{n}{Y}_n^m(\theta ,\phi )\mathbf{r}\right),}\hfill \end{array}$$

(33)

so it is a poloidal field. Therefore, the toroidal modes must be zero at $r=1$, i.e., $X_{kn}^T\left(1\right)=0$. Then the eigenvalue equation

$$\begin{array}{c}{\displaystyle j_n\left(\sqrt{\eta }\right)=0.}\hfill \end{array}$$

(34)

The normalization condition has the form

$$\begin{array}{c}{\displaystyle n(n+1){\int }_0^1{r}^2{\left[X_{kn}^T\left(r\right)\right]}^{2}dr=n(n+1){\int }_0^1{r}^2{\left[a_{kn}^T{j}_n\left(\sqrt{{\eta }_{kn}^{T}r}\right)\right]}^{2}dr=1.}\hfill \end{array}$$

(35)

From these conditions appear explicit expressions for the coefficients:

$$\begin{array}{c}{\displaystyle a_{kn}^T={\left[n(n+1){\int }_0^1{r}^2{j}_n^2\left(\sqrt{{\eta }_{kn}^{T}r}\right)dr\right]}^{-1/2}.}\hfill \end{array}$$

(36)

Let us now describe the calculation scheme for ${\eta }_{kn}^T$ and $a_{kn}^T$ in Maple. The following calculations are carried out in the loop for $n=1,\dots ,N$.

• The expression (36) for the coefficient $a_{kn}^T$ is calculated and written to the variable $U_{\mathrm{symbol}}$ in symbolic form with an undefined eigenvalue.
  Usymbol:=(n∗(n+1)∗int(r^2∗eval(j[n],z=sqrt(eta_)∗r)^2,r=0..1))^(-1/2):
• The left side of the Equation (34) is formed and its first $K+1$ roots ${\eta }_{kn}^T$ are determined numerically. The roots are first separated with a step of 1 and then refined with the required degree of accuracy.
  LPEq:=eval(j[n],z=sqrt(eta_)):
  eta_left:=1e-03:
  eta_right:=eta_left+1.:
  LPEq_left:=eval(LPEq,eta_=eta_left):
  k:=0:
  while (k<=K) do
    LPEq_right:=eval(LPEq,eta_=eta_right):
    if ((LPEq_left∗LPEq_right)<=0) then
      etaT[k,n]:=fsolve(LPEq=0,eta_=eta_left..eta_right):
      k:=k+1:
    end if:
    eta_left=eta_right:
    eta_right:=eta_left+1.:
    LPEq_left:=LPEq_right:
  end do:
• The numeric value $a_{kn}^T$ is calculated by substituting the numeric value ${\eta }_{kn}^T$ found above into $U_{\mathrm{symbol}}$.
  for k from 0 to K do
    aT[k,n]:=eval(Usymbol,[eta_=eta_T[k,n]]):
  end do:

The found values of $a_{kn}^T$ and ${\eta }_{kn}^T$, together with the expressions (32), determine the toroidal modes of the magnetic field.

In Figure 3, graphs of $X_{kn}^T$ functions for $n=2,5$ and $k=0,\dots ,4$ are given. Again, one can clearly see the fulfillment of the boundary conditions and the coincidence of the number of zero values of these functions on the interval $(0;1)$ with the value of the index k.

The poloidal solutions of the spectral problem (6) for a magnetic field have the form

$${\displaystyle \begin{array}{c}{\mathbf{B}}_{knm}^P=\mathrm{rot}\mathrm{rot}\left(X_{kn}^P\left(r\right)Y_n^m(\theta ,\phi )\mathbf{r}\right)=n(n+1)\frac{X_{k,n}^P\left(r\right)}{r}{Y}_n^m(\theta ,\phi ){\mathbf{e}}_r+\hfill \\ {\displaystyle +\left(\frac{d{X}_{kn}^P\left(r\right)}{dr}+\frac{X_{kn}^P\left(r\right)}{r}\right)\left({\mathbf{e}}_{\theta }\frac{\partial }{\partial \theta }+\frac{1}{sin\theta }{\mathbf{e}}_{\phi }\frac{\partial }{\partial \phi }\right)Y_n^m(\theta ,\phi ),}\hfill \\ {\displaystyle X_{kn}^P\left(r\right)=a_{kn}^P{j}_n\left(\sqrt{{\eta }_{kn}^P}r\right)+b_{kn}^P{y}_n\left(\sqrt{{\eta }_{kn}^P}r\right),}\hfill \\ k=0,1,2,\dots , n=1,2,3,\dots , m=-n,\dots ,n,\hfill \end{array}}$$

(37)

where eigenvalues ${\eta }_{kn}^P>0$.

The boundary condition $\mathbf{B}(r=0)<\infty$ from (4), as in the case of the toroidal modes, requires $b_{kn}^P=0$.

For $r=1$, a continuous transition of the ${\mathbf{B}}_{knm}^P$ modes to the external field ${\mathbf{B}}_{nm}^{out}$ modes

$${\displaystyle {\mathbf{B}}_{nm}^{out}=\mathrm{rot}\mathrm{rot}\left(\frac{r^{-n-1}}{n}{Y}_n^m(\theta ,\phi )\mathbf{r}\right)}$$

is needed. To do this, the components of the corresponding modes must be proportional. From this proportionality, it is easy to obtain the boundary condition

$$\begin{array}{c}{\displaystyle \frac{d{X}_{kn}^P}{dr}{|}_{r=1}+\left(n+1\right)X_{kn}^P\left(1\right)=0,}\hfill \end{array}$$

(38)

and then the eigenvalues equation

$$\begin{array}{c}{\displaystyle \sqrt{\eta }{j}_n^{\prime }\left(\sqrt{\eta }\right)+(n+1)j_n\left(\sqrt{\eta }\right)=0.}\hfill \end{array}$$

(39)

The normalization condition

$$\begin{array}{c}{\displaystyle n^2{(n+1)}^2{\int }_0^1{\left[X_{kn}^P\left(r\right)\right]}^{2}dr+n(n+1){\int }_0^1{r}^2{\left[\frac{d{X}_{kn}^P\left(r\right)}{dr}+\frac{X_{kn}^P\left(r\right)}{r}\right]}^{2}dr=1}\hfill \end{array}$$

(40)

gives the explicit expressions for the coefficients

$$\begin{array}{c}{\displaystyle a_{kn}^P=\left[n^2{(n+1)}^2{\int }_0^1{j}_n^2\left(\sqrt{{\eta }_{kn}^P}r\right)dr+\right.}\hfill \\ {\displaystyle +n(n+1){\int }_0^1{r}^2{\left.{\left(\frac{j_n\left(\sqrt{{\eta }_{kn}^P}r\right)}{dr}+\frac{j_n\left(\sqrt{{\eta }_{kn}^P}r\right)}{r}\right)}^{2}dr\right]}^{-1/2}.}\hfill \end{array}$$

(41)

Let us describe the calculation scheme for ${\eta }_{kn}^P$ and $a_{kn}^P$ in Maple. The following calculations are performed in the loop for $n=1,\dots ,N$:

• The expression (41) for the coefficient $a_{kn}^P$ is calculated and written to the variable $U_{\mathrm{symbol}}$ in symbolic form with undefined eigenvalues.
  Usymbol:=(n^2∗(n+1)^2∗int(eval(j[n],z=sqrt(eta_)∗r)^2,r=0..1)+
           n∗(n+1)∗int(r^2∗(diff(eval(j[n],z=sqrt(eta_)),r)+
           eval(j[n],z=sqrt(eta_))/r)^2,r=0..1))^(-1/2):
• The left part of the Equation (39) is formed and its first $K+1$ roots ${\eta }_{kn}^P$ are determined numerically. The roots are first separated with a step of 1 and then refined with the required degree of accuracy.
   LPEq:=sqrt(eta_)∗eval(diff(j[n],z),z=sqrt(eta_))+
         (n+1)∗eval(j[n],z=sqrt(eta_)):
   eta_left:=1e-03:
   eta_right:=eta_left+1.:
   LPEq_left:=eval(LPEq,eta_=eta_left):
   k:=0:
   while (k<=K) do
     LPEq_right:=eval(LPEq,eta_=eta_right):
     if ((LPEq_left∗LPEq_right)<=0) then
       etaP[k,n]:=fsolve(LPEq=0,eta_=eta_left..eta_right):
       k:=k+1:
     end if:
     eta_left=eta_right:
     eta_right:=eta_left+1.:
     LPEq_left:=LPEq_right:
   end do:
• The numeric value $a_{kn}^P$ is calculated by substituting into $U_{\mathrm{symbol}}$ the numeric value ${\eta }_{kn}^P$ found above.
  for k from 0 to K do
    aP[k,n]:=eval(Usymbol,[eta_=eta_P[k,n]]):
  end do:

The found values of $a_{kn}^P$ and ${\eta }_{kn}^P$, together with the expressions (37), define the poloidal magnetic modes.

In Figure 4, graphs of the $X_{kn}^P$ functions for $n=1,3$ and $k=0,\dots ,4$ on the segment $[0,1]$ are shown. They are extended to $r>1$ by the external field modes’ radial functions ${\mathbf{B}}_{nm}^{\mathrm{out}}$. It is clearly seen from the figure that the boundary conditions (38) ensure the differentiability of this continuation. Then the transition of the magnetic field itself will be continuous. As for the other modes, the number of zero values of these functions on the interval $(0;1)$ coincides with the value of the index k.

4.3. Temperature Modes

The solution of the spectral problem (4,6) for temperature is given by the following family of eigenmodes

$${\displaystyle \begin{array}{c}{\displaystyle T_{knm}=Z_{kn}\left(r\right)Y_n^m(\theta ,\phi ),}\hfill \\ {\displaystyle Z_{kn}\left(r\right)=a_{kn}{j}_n\left(\sqrt{{\lambda }_{kn}}r\right)+b_{kn}{y}_n\left(\sqrt{{\lambda }_{kn}}r\right),}\hfill \\ k=0,1,2,\dots , n=0,2,3,\dots , m=-n,\dots ,n,\hfill \end{array}}$$

(42)

where eigenvalues ${\lambda }_{kn}>0$. This general solution structure follows from the general solution form for the problem of a sphere free of oscillations [50].

Note that, the index n for the temperature modes can also be zero, in contrast to the velocity and magnetic field modes.

The coefficients $a_{kn}$ and $b_{kn}$ must provide the normalization condition

$$\begin{array}{c}{\displaystyle {\int }_{r_i}^1{r}^2{\left[Z_{kn}^T\left(r\right)\right]}^{2}dr=1,}\hfill \end{array}$$

(43)

and boundary conditions

$$\begin{array}{c}{\displaystyle X_{kn}\left(r_i\right)=X_{kn}\left(1\right)=0.}\hfill \end{array}$$

(44)

The Equalities (43) and (44) ensure that the $T_{knm}$ modes are normalized with respect to the scalar product ${\langle ·,·\rangle }_T$ and that the boundary conditions (4) are satisfied.

The eigenvalues of ${\lambda }_{kn}$ are solutions to the equation

$$\begin{array}{c}{\displaystyle j_n\left(\sqrt{\lambda }{r}_i\right)y_n\left(\sqrt{\lambda }\right)-j_n\left(\sqrt{\lambda }\right)y_n\left(\sqrt{\lambda }{r}_i\right)=0.}\hfill \end{array}$$

(45)

In form, this equation coincides with the Equation (21). Therefore, ${\lambda }_{kn}={\mu }_{kn}^T$ for $n\ge 1$. Comparison of (43) and (19) also shows that $a_{kn}=\sqrt{n(n+1)A_{kn}^T}$ and $b_{kn}=\sqrt{n(n+1)B_{kn}^T}$. Therefore, coefficients of the function $Z_{kn}\left(r\right)$ and eigenvalues can be determined from (43) and (45) only for $n=0$. However, in order to systematize the computational scheme, it is advisable to find them separately using the Maple tools, without reducing them to the parameters of toroidal velocity modes. Let us carry out the following calculations in Maple:

• The normalization expression (43) is calculated and written to the variable $U_{\mathrm{symbol}}$ in symbolic form with undefined values of the coefficients, eigenvalues, and the integration limit $r_i$.
  assume(r_i>0):
  Z_:=a_∗eval(j[n],z=sqrt(lambda_)∗r)+b_∗eval(y[n],z=sqrt(lambda_)∗r):
  Usymbol:=int(r^2∗Z_^2,r=r_i..1):
• The left part of the Equation (45) is formed and its first $K+1$ roots ${\lambda }_{kn}$ are determined numerically. The roots are first separated with a step of 1 and then refined with the required degree of accuracy.
  LPEq:=eval(j[n],z=sqrt(lambda_)∗.35)∗eval(y[n],z=sqrt(lambda_))-
              eval(j[n],z=sqrt(lambda_))∗eval(y[n],z=sqrt(lambda_)∗.35):
  lambda_left:=1e-03:
  lambda_right:=lambda_left+1.:
  LPEq_left:=eval(LPEq,lambda_=lambda_left):
  k:=0:
  while (k<=K) do
    LPEq_right:=eval(LPEq,lambda_=lambda_right):
    if ((LPEq_left∗LPEq_right)<=0) then
      lambda[k,n]:=fsolve(LPEq=0,lambda_=lambda_left..lambda_right):
      k:=k+1:
    end if:
    lambda_left=lambda_right:
    lambda_right:=lambda_left+1.:
    LPEq_left:=LPEq_right:
  end do:
• Assignments are made:

  $$\begin{array}{c}{\displaystyle a_{kn}:=y_n\left(\sqrt{{\lambda }_{kn}}\right),b_{kn}:=-j_n\left(\sqrt{{\lambda }_{kn}}\right).}\hfill \end{array}$$

  (46)
  Such values of the coefficients will ensure the fulfillment of the boundary conditions.
  for k from 0 to K do
    a[k,n]:=eval(y[n],z=sqrt(lambda[k,n])):
    b[k,n]:=-eval(j[n],z=sqrt(lambda[k,n])):
  end do:
• The numeric value $U_{\mathrm{number}}$ is calculated by substituting the numeric values $a_{kn}$, $b_{kn}$, ${\lambda }_{kn}$ into $U_{\mathrm{symbol}}$ and $r_i$. Then, renormalization is performed using assignments:

  $$\begin{array}{c}{\displaystyle a_{kn}:=a_{kn}/\sqrt{U_{\mathrm{number}}},b_{kn}:=b_{kn}/\sqrt{U_{\mathrm{number}}}.}\hfill \end{array}$$

  (47)
  for k from 0 to K do
    Unumber:=eval(Usymbol,
                  [a_=a[k,n],b_=b[k,n],lambda_=lambda[k,n],r_i=.35]):
    a[k,n]:=a[k,n]/sqrt(Unumber):
    b[k,n]:=b[k,n]/sqrt(Unumber):
  end do:

The calculations from points (1)–(4) must be performed in a loop for $n=0,\dots ,N$.

The found values $a_{kn}$, $b_{kn}$, ${\lambda }_{kn}$, together with the expressions (42), determine the basic toroidal velocity modes.

In Figure 5, graphs of $Z_{kn}\left(r\right)$ functions for $n=0,8$ and $k=0,\dots ,4$ are given. One can see the fulfillment of the boundary conditions and the coincidence of the number of zero values with the index k.

5. Calculation of the Spectral Model Coefficients

When compiling a spectral model (8), the main computational problem is to calculate the coefficients (9). These coefficients are volume integrals of multiplicative combinations of basic modes and differential operators of vector analysis operators in spherical coordinates. The expressions for the modes themselves are composed of spherical harmonics, spherical Bessel functions, and the curl operator (for vector modes). It is clear that very complex expressions will be obtained as a result. It is very difficult to obtain them using manual formula transformations. In addition, it will be difficult to avoid mistakes when writing integrands in the program code. There is actually only one way out—to automate the calculation of integrands. Obviously, there is simply no alternative to using CAS to compose models like (8).

Here, it is also necessary to take into account that, after selecting a finite set of modes and compiling the model, it may turn out that it does not describe all the necessary processes and interactions that exist in the dynamo system. For example, if all $E_{li}$ coefficients are equal to zero, then the model (8) does not take into account shell rotation. If all coefficients $W_{pij}$ are equal to zero, then no field generation occurs. In such cases, it will be necessary to change modes and recalculate the coefficients.

Let us now describe the technology for calculating the coefficients of the model using the VectorCalculus package of the Maple system.

First of all, it is necessary to list the set of commands and functions of the Maple system that will be required for the calculations:

• SetCoordinates(’spherical’[r,theta,phi])—setting spherical coordinates $(r,\theta ,\phi )$.
• VectorField(<expr1,expr2,expr3>)—vector field constructor with expr1, expr2, expr3 components in a local spherical basis as an object of class VectorField. These components must depend on the coordinates $(r,\theta ,\phi )$.
• DotProduct(F,P)—calculation of the dot product of two objects of class VectorField.
• CrossProduct(F,P)—calculation of the two objects of class VectorField cross product.
• Divergence(F)—calculation of the VectorField class object divergence.
• Curl(F)—calculation of the VectorField class object curl.
• Gradient(expr)—calculation of the gradient of the scalar field, which is determined by the expression expr.

The velocity modes ${\mathbf{v}}_{knm}^T$ and ${\mathbf{v}}_{knm}^P$ are determined by the triples of indices $(k,n,m)$ and the type (toroidal or poloidal), so the velocity modes’ indices l in the decompositions (5) are multi-indices: $l\sim (k_l,n_l,m_l,typ{e}_l)$, where $typ{e}_l=0$ for toroidal modes and $typ{e}_l=1$ for poloidal modes. These multi-indices can be defined in Maple using a one-dimensional array of four-element lists indV. For example, if in the decompositions (5) ${\mathbf{v}}_3={\mathbf{v}}_{2,4,7}^T$, then the Maple code will be indV[3]: =[2,4,7,0]. Then, for example, the type value of the (toroidal) mode ${\mathbf{v}}_{2,4,7}^T$ will be indV[3][4], and its first index value will be indV[3][1].

The situation is completely similar with the ${\mathbf{B}}_p$ magnetic modes in (5) expansions.They can be identified using a one-dimensional array of four-element lists indB. This array is organized in the same way as indV.

To identify the temperature modes $T_s$, a one-dimensional array of three-element lists indT is used. Here, there are only three elements in each list, since it is not necessary to specify the mode type, only the indices $(k_s,n_s,m_s)$

Let K, N and M be the maximum values of the first three elements in the multi-indices. Let us assume that the program has already calculated or read from the files the elements of the array ($k=0,\dots ,K,n=0,\dots ,N,m=0,\dots ,M$):

–

expressions j[n] and y[n] for spherical Bessel functions $j_n\left(z\right)$ and $y_n\left(z\right)$ and expressions P[n,m] for associated Legendre functions $P_n^m(cos\theta )$;

–

muT[k,n], AT[k,n], BT[k,n] eigenvalues ${\mu }_{kn}^T$ and coefficients $A_{kn}^T$ and $B_{kn}^T$ toroidal velocity modes;

–

muP[k,n], C1[k,n], … ,C4[k,n] eigenvalues ${\mu }_{kn}^P$ and coefficients $C_{kn}^1,\dots ,C_{kn}^4$ poloidal velocity modes;

–

etaT[k,n], aT[k,n] eigenvalues ${\eta }_{kn}^T$ and coefficients $a_{kn}^T$ of toroidal magnetic modes;

–

etaP[k,n], aP[k,n] eigenvalues ${\eta }_{kn}^P$ and coefficients $a_{kn}^P$ of poloidal magnetic modes;

–

lambda[k,n], a[k,n], b[k,n] eigenvalues ${\lambda }_{kn}$ and coefficients $a_{kn}$ and $b_{kn}$ temperature modes.

The associated Legendre functions $P_n^{\left|m\right|}(cos\theta )$ must be normalized so that the spherical harmonics

$$\begin{array}{c}{\displaystyle Y_n^m\left(\theta ,\phi \right)=\left\{\begin{array}{cc}{P}_n^{-m}(cos\theta )cos\left(m\phi \right),\hfill & \mathrm{if}m\le 0;\hfill \\ P_n^m(cos\theta )sin\left(m\phi \right),\hfill & \mathrm{if}m>0\hfill \end{array}\right.}\hfill \end{array}$$

(48)

are normalized on the sphere in the sense of equality (15).

The spherical harmonics are a complete orthogonal system of functions on the sphere. Therefore, they are included in the expressions for all basic modes. However, in Maple code, it is convenient to define an array of expressions for the associated Legendre functions. Then, instead of a separate array of expressions for spherical harmonics, one can use the right-hand sides of (48).

Assume that the multi-index arrays indV, indT, indB are filled. These arrays define a set of basis modes (5).

Next, arrays of expressions for modes are compiled, but a specific type of dependence on the r variable is not specified yet. To calculate the Galerkin coefficients, it is also convenient to compose an array of curls of the velocity and magnetic field modes.

The corresponding Maple code would be like this.

As a result of executing the code shown in Listing 1, expressions will be formed for the expression arrays for the basic modes, where the dependence on the angular variables will be explicitly specified. However, explicit expressions for the radial functions have not yet been given. Instead, there will be abstract R_[l](r), X_[p](r), Z_[s](r).

Listing 1. Modes expressions calculations.

#

with(VectorCalculus):

SetCoordinates(’spherical’[r,theta,phi]):

#----velocity_modes-----------------------

V:=array(1..Lmax,[]):

rotV:=array(1..Lmax,[]):

for l from 1 to Lmax do

n:=indV[l][2]: m:=indV[l][3]:

if (m<=0)then

func_phi:=cos(-m∗phi):

else

func_phi:=sin(m∗phi):

end if:

V[l]:=simplify(

Curl(R_[l](r)∗P[n,abs(m)]∗func_phi∗VectorField(<r,0,0>))):

if (indV[l][4]=1) then

V[l]:=simplify(Curl(V[l])):

end if:

rotV[l]:=simplify(Curl(V[l])):

end do:

#----magnetic_modes-----------------------

B[p]:=array(1..Pmax,[]):

rotB:=array(1..Pmax,[]):

for p from 1 to Pmax do

n:=indB[p][2]: m:=indB[p][3]:

if (m<=0)then

func_phi:=cos(-m∗phi):

else

func_phi:=sin(m∗phi):

end if:

B[p]:=simplify(

Curl(X_[p](r)∗P[n,abs(m)]∗func_phi∗VectorField(<r,0,0>))):

if (indB[p][4]=1) then

B[p]:=simplify(Curl(B[p])):

end if:

rotB[p]:=simplify(Curl(B[p])):

end do:

#----thermal_modes-----------------------

T[s]:=array(1..Smax,[]):

for s from 1 to Smax do

n:=indT[s][2]: m:=indT[s][3]:

if (m<=0)then

func_phi:=cos(-m∗phi):

else

func_phi:=sin(m∗phi):

end if:

T[s]:=simplify(Z_[s](r)∗P[n,abs(m)]∗func_phi):

end do:

Let us consider such an example. Let ${\mathbf{v}}_5={\mathbf{v}}_{1,4,-3}^P$; then the element indV[5] of the multi-indexes array will be equal to [1,4,−3,1]. After executing the Listing 1 code, the element V[5] of the velocity modes array will be equal to the expression shown in Figure 6. The structure of (26) is clearly visible in the expression in this Figure.

The explicit expressions for the radial functions in the arrays of modes are not specified for the following reason: First, analytical integration over the angular variables $\phi$ and $\theta$ is carried out when calculating the coefficients (9).

First of all, integration over $\phi$ is performed. For this, specific expressions for the radial functions are not needed. CAS only needs to know that these are some functions of the r argument. If the integration result is zero, no further integration is necessary. If it is non-zero, then analytic integration over the variable $\theta$ is performed. And, only if this integral is also different from zero, the expression of radial functions is substituted into it and numerical integration over r is performed. Expressions for the radial functions contain parameters that were previously determined numerically. In such a situation, analytic integration over r is no longer appropriate. This is the basic scheme for calculating all the coefficients. Note that, the experience of constructing models shows that the equality of the Galerkin coefficient to zero is usually revealed precisely when integrating over a sphere.

Now let us define explicit expressions for all the radial functions for use in future substitutions.

Further, it is already possible to calculate the coefficients (9) of the spectral model. It is assumed that $r_i=0.35$.

Now, consider the Maple code for calculating the $B_{lij}$ coefficients based on symmetrical definitions (11):

Here, the outer loops go through indexes i and j because this allows only once to evaluate the expression ${\mathbf{v}}_i\times \mathrm{rot}{\mathbf{v}}_j+{\mathbf{v}}_j\times \mathrm{rot}{\mathbf{v}}_i$ and store it in a temporary variable. Then, in the cycle over the index l, the integrand for $B_{lij}$ is compiled and is integrated successively over $\phi$, $\theta$, and r. In this case, the scheme of gradual substitution of dependence on r, described above, is used.

All the other coefficients, except for $Q_{lij}$, are determined using similar calculations based on (9) definitions. The corresponding codes are given in Appendix B. The coefficients $Q_{lij}$, symmetric in i and j, are easily computed from $W_{pij}$ using the equalities (13).

It was already mentioned in the introduction that, in order to distinguish chains of interacting modes, it is important to know whether one or another coefficient of the model is exactly zero or not. This is important physical information. If the equality to zero is established at the stage of analytical integration, then everything is clear. But it is also possible that the integration over r should give this zero value. In numerical integration, you will never get an exactly zero result.

Floating-point computations in CAS are implemented using the software environment. This makes them slow, but makes it easy to control computation accuracy. For example, in Maple, the value of the Digits environment variable is responsible for this, and is easily changed (the default is 10). If some coefficient is actually zero, then its value will be a small value, the order of which is related to the value of Digits. Therefore, recalculating the coefficients with different accuracy, it is easy to single out the zero ones.

6. Example of a Spectral Model Construction

In the authors’ work [25], a spectral geodynamo model based on the hypothesis of a 6-cell large-scale convection in the Earth’s core was studied. Let us show how the proposed approaches make it possible to calculate the coefficients of the model. The physical substantiation of the 6-cell convection hypothesis, the selection rule for fundamental modes, and the results of numerical simulation will not be discussed here. The purpose of our present work is to describe the computational technology for compiling models with any given set of modes. Therefore, only the calculation of the mode parameters and coefficients of the model is of interest.

First of all, let us list the modes that are used in the model [25], using the notation of our paper. The list of modes is given in

Table 1. List of basic modes.

Index i	Velocity Modes ${\mathbf{v}}_{\mathit{i}}$	Magnetic Modes ${\mathbf{B}}_{\mathit{i}}$	Temperature Modes ${\mathit{T}}_{\mathit{i}}$
1	${\mathbf{v}}_{0,5,-2}^T$	${\mathbf{B}}_{0,1,0}^P$	$T_1=T_{0,4,-2}$
2	${\mathbf{v}}_{0,5,2}^T$	${\mathbf{B}}_{0,3,-2}^P$	$T_2=T_{0,4,2}$
3	${\mathbf{v}}_{0,4,-2}^P$	${\mathbf{B}}_{0,3,2}^P$	$T_3=T_{1,0,0}$
4	${\mathbf{v}}_{0,4,2}^P$	${\mathbf{B}}_{1,1,0}^P$	–
5	${\mathbf{v}}_{1,3,-2}^T$	${\mathbf{B}}_{0,4,-2}^T$	–
6	${\mathbf{v}}_{1,3,2}^T$	${\mathbf{B}}_{0,4,2}^T$	–
7	–	${\mathbf{B}}_{0,5,-2}^P$	–
8	–	${\mathbf{B}}_{0,5,2}^P$	–

.

The poloidal velocity modes ${\mathbf{v}}_{0,4,\pm 2}^P$ and the temperature modes $T_{0,4,\pm 2}$ and $T_{1,0,0}$ describe the 6-cell convection. The toroidal velocity modes ${\mathbf{v}}_{0.4,\pm 2}^P$ and ${\mathbf{v}}_{0.4,\pm 2}^P$ simulate the Coriolis drift of the poloidal modes. The magnetic modes ${\mathbf{B}}_{0,1,0}^P$ and ${\mathbf{B}}_{1,1,0}^P$ describe the main geomagnetic dipole. The rest of the magnetic modes are needed to transfer the convection energy to the main dipole. This is the minimum mode set that generates a dipole in the 6-cell convection [25].

The eigenvalues and coefficients of the radial functions for the index values $k=0,\dots ,8$ and $n=0,\dots ,8$ were calculated using the approaches from Section 1. The calculation results are presented in the files (attached in the Supplementary Materials):

–

function $R_{kn}^T\left(r\right)$ – file “RT.txt”;

–

function $R_{kn}^P\left(r\right)$ – file “RP.txt”;

–

function $X_{kn}^T\left(r\right)$ – file “XT.txt”;

–

function $X_{kn}^P\left(r\right)$ – file “XP.txt”;

–

function $Z_{kn}\left(r\right)$ – file “Z.txt”;

Only the eigenvalues of the modes used in the model are given here. They are presented in

Table 2. Modes eigenvalues.

Index i	${\mathit{\mu }}_{\mathit{i}}$	${\mathit{\eta }}_{\mathit{i}}$	${\mathit{\lambda }}_{\mathit{i}}$
1	87.900199617096118967	9.8696044010893586188	67.86985576121043418
2	87.900199617096118967	33.21746191426836886	67.86985576121043418
3	100.15439579986503383	33.21746191426836886	93.440041667118188107
4	100.15439579986503383	39.478417604357434475	–
5	125.97590077869094824	66.954311925104805326	–
6	125.97590077869094824	66.954311925104805326	–
7	–	66.954311925104805326	–
8	–	66.954311925104805326	–

.

It can be seen that some eigenvalues coincide. The reason is that eigenmodes that differ only in the order m value of the spherical harmonics have the same eigenvalues. This explains the following equalities:

–

${\mu }_1={\mu }_2$: eigenvalues for ${\mathbf{v}}_{0,5,-2}^T$ and ${\mathbf{v}}_{0,5,2}^T$;

–

${\mu }_3={\mu }_4$: eigenvalues for ${\mathbf{v}}_{0,4,-2}^P$ and ${\mathbf{v}}_{0,4,2}^P$;

–

${\mu }_5={\mu }_6$: eigenvalues for ${\mathbf{v}}_{1,3,-2}^T$ and ${\mathbf{v}}_{1,3,2}^T$;

–

${\eta }_2={\eta }_3$: eigenvalues for ${\mathbf{B}}_{0,3,-2}^P$ and ${\mathbf{B}}_{0,3,2}^P$;

–

${\eta }_5={\eta }_6$: eigenvalues for ${\mathbf{B}}_{0,4,-2}^T$ and ${\mathbf{B}}_{0,4,2}^T$;

–

${\eta }_7={\eta }_8$: eigenvalues for ${\mathbf{B}}_{0,5,-2}^P$ and ${\mathbf{B}}_{0,5,2}^P$;

–

${\lambda }_1={\lambda }_2$: eigenvalues for $T_{0,4,-2}^P$ and $T_{0,4,2}^P$.

The coincidence of the last four eigenvalues for magnetic modes in Table 2 is explained as follows: For any integer l, the following equality is true [49]:

$${\displaystyle \frac{2l+1}{z}{j}_l\left(z\right)-j_{l-1}\left(z\right)=-z^l\frac{d}{dz}\left[z^{-l}{j}_l\left(z\right)\right].}$$

(49)

From this equality it is easy to obtain that

$${\displaystyle z{j}_l^{\prime }\left(z\right)+(l+1)j_l\left(z\right)=z{j}_{l-1}\left(z\right).}$$

(50)

Let us now compare (50) with the eigenvalue Equations (34) and (39). It can be seen that the positive eigenvalues ${\eta }_{k,n}^P$ and ${\eta }_{k,n-1}^T$ are solutions of equivalent equations. Therefore, these eigenvalues coincide. Then, it turns out that ${\eta }_5={\eta }_6={\eta }_{0,4}^T={\eta }_{0,5}^P={\eta }_7={\eta }_8$.

For modes from Table 1, it is necessary to form the corresponding arrays of multi-indexes. Here is a Maple code snippet.

Lmax:=6: Pmax:=8: Smax:=3:

indV:=array(1..Lmax,[]):

indB:=array(1..Pmax,[]):

indT:=array(1..Smax,[]):

#--------velocity modes indices-----------------

indV[1]:=[0,5,-2,0];

indV[2]:=[0,5,2,0];

indV[3]:=[0,4,-2,1];

indV[4]:=[0,4,2,1];

indV[5]:=[1,3,-2,0];

indV[6]:=[1,3,2,0];

#--------magnetic modes indices-----------------

indB[1]:=[0,1,0,1];

indB[2]:=[0,3,-2,1];

indB[3]:=[0,3,2,1];

indB[4]:=[1,1,0,1];

indB[5]:=[0,4,-2,0];

indB[6]:=[0,4,2,0];

indB[7]:=[0,5,-2,1];

indB[8]:=[0,5,2,1];

#--------temperature modes indices-----------------

indT[1]:=[0,4,-2];

indT[2]:=[0,4,2];

indT[3]:=[1,0,0];

After that, arrays of expressions are calculated for the modes V[l], B[p] and T[s] with undefined radial functions R_[l](r), X_[p](r) and T_[s](r) (Listing 1). The Maple output for the velocity modes’ expressions is shown in Figure 7.

Next, the arrays of explicit expressions of the radial functions R[l], X[p] and T[s] are calculated (Listing 2). For the first three velocity modes, the expressions for the radial functions are shown in Figure 8.

Listing 2. Radial function expressions calculation.

#

#----velocity_modes_radial_functions--------------

for l from 1 to Lmax do

k:=indV[l][1]: n:=indV[l][2]: type_:=indV[l][4]

if (type_=1) then

R[l]:=C1[k,n]∗eval(j[n],[z=sqrt(muP[k,n])∗r])+

C2[k,n]∗eval(y[n],[z=sqrt(muP[k,n])∗r])+

C3[k,n]∗r^n+C4[k,n]∗r^(-n-1):

else

R[l]:=AT[k,n]∗eval(j[n],[z=sqrt(muT[k,n])∗r])+

BT[k,n]∗eval(y[n],[z=sqrt(muT[k,n])∗r]):

end if:

end do:

#----magnetic_modes_radial_functions-------------

for p from 1 to Pmax do

k:=indB[p][1]: n:=indB[p][2]: type_:=indB[p][4]

if (type_=1) then

X[p]:=aP[k,n]∗eval(j[n],[z=sqrt(etaP[k,n])∗r]):

else

X[p]:=aT[k,n]∗eval(j[n],[z=sqrt(etaT[k,n])∗r]):

end if:

end do:

#----thermal_modes_radial_functions-------------

for s from 1 to Smax do

k:=indT[s][1]: n:=indT[s][2]:

Z[s]:=a[k,n]∗eval(j[n],[z=sqrt(lambda[k,n])∗r])+

b[k,n]∗eval(y[n],[z=sqrt(lambda[k,n])∗r]):

end do:

Next, all the coefficients of the model can be calculated. To illustrate the symbolic computations, the values of the coefficients $E_{2,4}$ and $E_{3,4}$ after integration over the variable $\phi$ are shown in Figure 9. These are only non-zero coefficients. Next comes the integration over the variable $\theta$. The results are shown in Figure 10. The final step in the code is substitution of the radial function expressions and integration over the radius. The final values of the coefficients are $E_{2,4}=-0.859376$ and $E_{2,4}=0.2$.

Note that, all model coefficients were calculated using Maple (Listings 3 and A1–A6). Calculations $W_{pi}^{\alpha }$ were carried out at $A\left(r\right)=1$ and $A\left(r\right)=r$, as in work [25]. All the $B_{lij}$ coefficients turned out to be zero. The values of all non-zero coefficients of the model are presented in the files (attached in the Supplementary Materials):

–

$E_{li}$ – file “E_coeff.txt”;

–

$C_{li}$ – file “C_coeff.txt”;

–

$Q_{lij}$ – file “Q_coeff.txt”;

–

$F_{sij}$ – file “F_coeff.txt”;

–

$H_{si}$ – file “H_coeff.txt”;

–

$W_{sij}$ – file “W_coeff.txt”;

–

$W_{sij}^{\alpha }$ for $A\left(r\right)=1$ – file “Walpha_coeff_1.txt”;

–

$W_{sij}^{\alpha }$ for $A\left(r\right)=r$ – file “Walpha_coeff_r.txt”.

Listing 3. calculation.]Coefficients B[l,i,j] calculation.

#----------------------------------------------------------------------

r_i:=.35:

B_:=array(1..Lmax,1..Lmax,1..Lmax,[]):

for i1 from 1 to Lmax do

for i2 from 1 to i1 do

tmp1:=CrossProduct(V[i],rotV[j])+CrossProduct(V[j],rotV[i]):

for l from 1 to L do

if (l=i1 and l=i2) then

B_[l,i1,i2]:=0:

else

tmp2:=DotProduct(tmp1,V[l]):

B_[l,i1,i2]:=int(tmp2,phi=-Pi..Pi):

end if:

if (B_[l,i,j]<>0) then

B_[l,i,j]:=int(sin(theta)∗B_[l,i,j],theta=0..Pi):

end if:

if (B_[l,i,j]<>0) then

B_[l,i,j]:=int(r^2∗eval(B_[l,i,j],[R_[l](r)=R[l],

R_[i](r)=R[i],R_[j](r)=R[j]]),r=r_i..1)/2:

end if:

B_[l,j,i]:=B_[l,i,j]:

end do:

end do:

end do:

7. Conclusions

A new computational technology was developed for calculating the parameters of the basic modes in spectral geodynamo models and calculating the coefficients of these models. The technology is based on a combination of symbolic computation and numerical computation in the Maple system.

The main difficulty in studying the geodynamo using spectral models is calculation of the basic modes and the Galerkin coefficients. The proposed technology makes it easy to solve these problems in automatic mode. Algorithms and programs were developed for:

—

derivation of explicit analytical expressions for the modes of free decay of the fields of speed, temperature and magnetic induction;

—

numerical calculation of eigenvalues and parameters of eigenmodes;

—

derivation of explicit analytical expressions for the integrands of the Galerkin coefficients;

—

combined analytical and numerical calculation of the Galerkin coefficients.

In the future, it will be possible to develop the results of the work for the cases of a compressible fluid, of the cylindrical geometry of the domain, and of a free boundary. This will make it possible to generalize the results of the work to some other problems of geophysical hydrodynamics.

The proposed computational technology makes it possible to quickly and accurately construct fairly wide classes of new spectral models of geodynamo. This will make it possible to advance the study of the geodynamo mechanism.