Using symbolic calculations to calculate the eigenmodes of the free damping of a geomagnetic ﬁeld

. The method for calculating the eigenmodes of free damped oscillations of the geomagnetic ﬁeld in the Earth’s core using symbolic computations is described.


Introduction
When simulating of geodynamo by spectral methods, the eigenmodes of free damping oscillations of the geomagnetic field in the core of the Earth are often useful. These modes form a complete orthogonal system in the core and can be used as a basis for magnetic field. However, the calculation of these modes themselves for large values of spherical indices is a rather complicated computational problem. This paper proposes a method for calculating modes using systems of symbolic computations (computer algebra systems -CAS) that can be easily automated.

Modes of free decay
The spectral problem of free damped oscillations of the geomagnetic field has the form [1]: where η -is eigenvalue. It is closed by the conditions of boundedness of the field in the center of the core B(0) < ∞, and the continuous transition of a potential external field for outer boundary of kernel r = 1.
It is well-known that solutions of problem (1) finite for r = 0 consist of toroidal and poloidal eigenmodes [1,2]: where Y m n (θ, ϕ) -are spherical harmonics, and X kn (r) = a kn j n λ kn r , In the Eqs. (4) j n (·) -are spherical Bessel functions first kinds, a kn and b kn -are normalization factors, λ kn and µ kn -are eigenvalues of toroidal and poloidal modes, respectively. Recall that the spherical Bessel functions are defined for each integer n by the recurrence relations [3]: We will assume that the rms norm of Y m n (θ, ϕ) on a sphere is equal to one, i.e.
and eigenmodes have a unit rms norm in the volume of the kernel: Further calculation of the eigenmodes is reduced to finding the normalizing coefficients and eigenvalues for the functions (4).

Function X kn (r) parameter calculation
For toroidal eigenmodes, the conditions of continuous transition for r = 1 lead to the condition X kn (1) = 0. We obtain the equation for eigenvalues λ: For every n, this equation has a countable set of solutions λ kn . The solutions is easy to localize, but their calculation using standard double-precision floating-point numbers results in computational instability at n 30. Most likely this is due to deep recursion in the (5) and the accumulation of computational errors. One possible waythe use of arbitrary-precision arithmetic. Such capabilities are provided by CAS.
After calculating the eigenvalue λ k n, you can determine the normalization factor a kn . Substitute expression (3) in to normalization condition (7) for the toroidal field and integrate over the sphere. Then get the normalization condition for the function X kn (r): n(n + 1) 1 0 a kn j n λ kn r 2 r 2 dr = 1, where we can determine a kn .
However, as can be seen from Eqs. (4,5), this is an improper integral of the 2nd kind with a singularity at r = 0. Given a very complex integrand, the computation of such an integral by numerical methods is almost impossible. The calculations lead to the occurrence of an overflow.
At the same time, the integral can be calculated analytically, but this is a very cumbersome procedure, the execution of which manually will inevitably lead to errors. You can avoid these errors and complex calculations if you use CAS. The system must be entrusted with calculating the integral with undefined parameters. Then the numerical values of the parameters are substituted into the algebraic expression formed by CAS.
It can be seen that the number of zeros of the function X kn (r) on the interval (0; 1) coincides with the value of the index k. This indirectly confirms the correctness of our calculations.

Function Z kn (r) parameters calculation
We work in a similar way with poloidal modes, although here the calculations are more complicated.
The conditions of continuous transition for r = 1 lead to the condition dZ(r) dr Then, using the well-known relation for spherical Bessel functions [3]: we obtain the equation for eigenvalues µ: This equation has a countable set of solutions µ kn for every n.
As in the case of toroidal modes, for large n, calculations of arbitrary-precision arithmetic are needed. After calculating the eigenvalue µ k n, you can determine the normalization factor b kn . Substitute expression (3) in to normalization condition (7) for the poloidal field and integrate over the sphere. Then get the normalization condition for the function Z kn (r): It is clearly seen that this integral is even more complicated than (9). Calculate it using CAS.
According to this scheme, we calculated the values of the parameters for n = 1, 2, . . . , 50 and k = 0, 1, . . . , 49. Fig. 2 shows the graphs of functions X kn (r) for n = 10 and k = 0, . . . , 5, and shows the continuously differentiable transition into radial functions of a potential external field. This transition of radial functions provides continuous transition of magnetic modes. It can be seen that the number of zeros of the function Z kn (r) on the interval (0; 1) coincides with the value of the index k. This indirectly confirms the correctness of our calculations.
We note that Eqs. (3) and (11) imply the coincidence of the eigenvalues of the modes T k,n and P k,n+1 . These eigenvalues were calculated by us independently of each other separately for each group of modes. The results gave the required match. This confirms the correctness of the calculations.

Conclusion
A method has been developed for calculating the parameters (eigenvalues and normalization coefficients) of eigenmodes for the free decay of the geomagnetic field for large values of spherical and radial indices. Calculations are carried out using analytical integration and arbitrary-precision calculations in a computer algebra system. This technique is implemented as programs in the CAS Maple system, which allows you to automate calculations.