# Numerical technique¶

**MagIC** is a pseudo-spectral MHD code. This numerical technique was
originally developed by P. Gilman and G. Glatzmaier for the spherical geometry.
In this approach the unknowns are expanded into complete sets of functions in
radial and angular directions: **Chebyshev polynomials or Finite differences in
the radial direction** and **spherical harmonic functions in the azimuthal and
latitudinal directions**. This allows to express all partial derivatives
analytically. Employing orthogonality relations of spherical harmonic
functions and using collocation in radius then lead to algebraic equations that
are integrated in time with a **mixed implicit/explicit time stepping scheme**.
The nonlinear terms and the Coriolis force are evaluated in the physical (or
grid) space rather than in spectral space. Although this approach requires
costly numerical transformations between the two representations (from spatial
to spectral using Legendre and Fourier transforms), the resulting decoupling of
all spherical harmonic modes leads to a net gain in computational speed.
Before explaining these methods in more detail, we introduce the
poloidal/toroidal decomposition.

## Poloidal/toroidal decomposition¶

Any vector \(\vec{v}\) that fulfills \(\vec{\nabla}\cdot\vec{v}=0\), i.e.
a so-called *solenoidal field*,
can be decomposed in a poloidal and a toroidal part \(W\) and \(Z\),
respectively

Three unknown vector components are thus replaced by two scalar fields, the poloidal potential \(W\) and the toroidal potential \(Z\). This decomposition is unique, aside from an arbitrary radial function \(f(r)\) that can be added to \(W\) or \(Z\) without affecting \(\vec{v}\).

In the anelastic approximation, such a decomposition can be used for the mass flux \(\tilde{\rho}\vec{u}\) and the magnetic field \(\vec{B}\). This yields

The two scalar potentials of a divergence free vector field can be extracted from its radial component and the radial component of its curl using the fact that the toroidal field has not radial component:

where the operator \(\Delta_H\) denotes the horizontal part of the Laplacian:

The equation (1) can be expanded in spherical coordinates. The three components of \(\tilde{\rho}\vec{u}\) are given by

while the curl of \(\tilde{\rho}\vec{u}\) is expressed by

Using the horizontal part of the divergence operator

above expressions can be simplified to

and

Below we will use the fact that the horizontal components of the poloidal field depend on the radial derivative of the poloidal potential.

## Spherical harmonic representation¶

Spherical harmonic functions \(Y_\ell^m\) are a natural choice for the horizontal expansion in colatitude \(\theta\) and longitude \(\phi\):

where \(\ell\) and \(m\) denote spherical harmonic degree and order, respectively, \(P_\ell^m\) is an associated Legendre function. Different normalization are in use. Here we adopt a complete normalization so that the orthogonality relation reads

This means that

As an example, the spherical harmonic representation of the magnetic poloidal potential \(g(r,\theta,\phi)\), truncated at degree and order \(\ell_{max}\), then reads

with

The potential \(g(r,\theta,\phi)\) is a real function so that \(g_{\ell m}^\star(r)=g_{\ell,-m}(r)\), where the asterisk denotes the complex conjugate. Thus, only coefficients with \(m \ge 0\) have to be considered. The same kind of expansion is made for the toroidal magnetic potential, the mass flux potentials, pressure, entropy (or temperature) and chemical composition.

The equations (8) and (9) define a two-step transform from the longitude/latitude representation to the spherical harmonic representation \((r,\theta,\phi)\longrightarrow(r,\ell,m)\). The equation (7) formulates the inverse procedure \((r,\ell,m)\longrightarrow(r,\theta,\phi)\). Fast-Fourier transforms are employed in the longitudinal direction, requiring (at least) \(N_\phi = 2 \ell_{max}+1\) evenly spaced grid points \(\phi_i\). MagIC relies on the Gauss-Legendre quadrature for evaluating the integral (8)

where \(\theta_j\) are the \(N_{\theta}\) Gaussian quadrature points defining the latitudinal grid, and \(w_j\) are the respective weights. Pre-stored values of the associated Legendre functions at grid points \(\theta_j\) in combination with a FFT in \(\phi\) provide the inverse transform (7). Generally, \(N_\phi= 2 N_\theta\) is chosen, which provides isotropic resolution in the equatorial region. Choosing \(\ell_{max}= [ \min(2 N_\theta,N_\phi)-1]/3\) prevents aliasing errors.

See also

In MagIC, the Legendre functions are defined in the subroutine
`plm_theta`

. The Legendre transforms
from spectral to grid space are computed in the module
`legendre_spec_to_grid`

, while the backward transform (from grid
space to spectral space) are computed in the module
`legendre_grid_to_spec`

. The fast Fourier transforms are computed
in the module `fft`

.

### Special recurrence relations¶

The action of a horizontal Laplacian (3) on spherical harmonics can be analytically expressed by

They are several useful recurrence relations for the Legendre polynomials that will
be further employed to compute Coriolis forces and the \(\theta\) and \(\phi\)
derivatives of advection and Lorentz forces.
Four of these operators are used in **MagIC**. The first one is defined by

The action of this operator on a Legendre polynomials is given by

where \(c_\ell^m\) is defined by

*How is that implemented in the code?* Let’s assume we want the spherical harmonic contribution
of degree \(\ell\) and order m for the expression

In order to employ the operator \(\vartheta_1\) for the derivative, we thus define a new function

so that

Expanding \(F(\theta)\) in Legendre polynomials and using the respective orthogonality relation we can then map out the required contribution in the following way:

Here, we have assumed that the Legendre functions are completely normalized such that

See also

This operator is defined in the module `horizontal_data`

by the variables
`dTheta1S`

for the first part of the right-hand
side of (12) and `dTheta1A`

for the
second part.

The second operator used to formulate colatitude derivatives is

The action of this operator on the Legendre polynomials reads

so that

See also

This operator is defined in the module `horizontal_data`

by the variables
`dTheta2S`

for the first part of the right-hand
side of (13) and `dTheta2A`

for the
second part.

The third combined operator is defined by:

where \(-L_H/r^2=\Delta_H\).

Acting with \(\vartheta_3\) on a Legendre function gives:

which results into:

See also

This operator is defined in the module `horizontal_data`

by the variables
`dTheta3S`

for the first part of the right-hand
side of (14) and `dTheta3A`

for the
second part.

The fourth (and last) combined operator is defined by:

Acting with \(\vartheta_3\) on a Legendre function gives:

which results into:

See also

This operator is defined in the module `horizontal_data`

by the variables
`dTheta4S`

for the first part of the right-hand
side of (15) and `dTheta4A`

for the
second part.

## Radial representation¶

In MagIC, the radial dependencies are either expanded into complete sets of functions, the Chebyshev polynomials \({\cal C}(x)\); or discretized using finite differences. For the former approach, the Chebyshev polynomial of degree \(n\) is defined by

When truncating at degree \(N\), the radial expansion of the poloidal magnetic potential reads

with

The Chebyshev definition space \((-1\leq x\leq 1)\) is then linearly mapped onto a radius range \((r_i\leq r \leq r_o)\) by

In addition, nonlinear mapping can be defined to modify the radial dependence of the grid-point density.

When choosing the \(N_r\) extrema of \({\cal C}_{N_r-1}\) as radial grid points,

the values of the Chebyshev polynomials at these points are simply given by the cosine functions:

This particular choice has two advantages. For one, the grid points become denser toward the inner and outer radius and better resolve potential thermal and viscous boundary layers. In addition, type I Discrete Cosine Transforms (DCTs) can be employed to switch between grid representation (16) and Chebyshev representations (17), rendering this procedure a fast-Chebyshev transform.

See also

The Chebyshev (Gauss-Lobatto) grid is defined in the module
`chebyshev_polynoms_mod`

. The cosine transforms are computed in the
modules `cosine_transform_even`

and `fft_fac_mod`

.

## Spectral equations¶

We have now introduced the necessary tools for deriving the
spectral equations.
Taking the **radial components** of the Navier-Stokes equation
and the induction equation provides the equations
for the poloidal potentials \(W(r,\theta,\phi)\) and \(g(r,\theta,\phi)\).
The **radial component of the curl** of these equations provides
the equations for the toroidal counterparts
\(Z(r,\theta,\phi)\) and \(h(r,\theta,\phi)\).
The pressure remains an additional unknown. Hence one more equation
involving \(W_{\ell mn}\) and \(p_{\ell mn}\)
is required. It is obtained by taking the
**horizontal divergence** of the Navier-Stokes equation.

Expanding all potentials in spherical harmonics and Chebyshev polynomials, multiplying with \({Y_{\ell}^{m}}^\star\), and integrating over spherical surfaces (while making use of the orthogonality relation (6) results in equations for the coefficients \(W_{\ell mn}\), \(Z_{\ell mn}\), \(g_{\ell mn}\), \(h_{\ell mn}\), \(P_{\ell mn}\) and \(s_{\ell mn}\), respectively.

### Equation for the poloidal potential \(W\)¶

The temporal evolution of \(W\) is obtained by taking \(\vec{e_r}\cdot\) of each term entering the Navier-Stokes equation. For the time-derivative, one gets using (2):

For the viscosity term, one gets

Note

In case of a constant kinematic viscosity, the \(d\ln\nu/dr\) terms vanish. If in addition,the background density is constant, the \(d\ln\tilde{\rho}/dr\) terms also vanish. In that Boussinesq limit, this viscosity term would then be simplified as

Using Eq. (10) then allows to finally write the time-evolution equation for the poloidal potential \(W_{\ell m n}\):

Here, \(d\Omega\) is the spherical surface element. We use the summation convention for the Chebyshev index \(n\). The radial derivatives of Chebyshev polynomials are denoted by primes.

### Equation for the toroidal potential \(Z\)¶

The temporal evolution of \(Z\) is obtained by taking the radial component of the curl of the Navier-Stokes equation (i.e. \(\vec{e_r}\cdot\vec{\nabla}\times\)). For the time derivative, one gets using (2):

The pressure gradient, one has

This expression has no component along \(\vec{e_r}\), as a consequence, there is no pressure gradient contribution here. The gravity term also vanishes as \(\vec{\nabla}\times(\tilde{\rho}g(r)\vec{e_r})\) has no radial component.

Note

Once again, this viscous term can be greatly simplified in the Boussinesq limit:

Using Eq. (10) then allows to finally write the time-evolution equation for the poloidal potential \(Z_{\ell m n}\):

### Equation for pressure \(P\)¶

The evolution of equation for pressure is obtained by taking the horizontal divergence (i.e. \(\vec{\nabla}_H\cdot\)) of the Navier-Stokes equation. This operator is defined such that

This relates to the total divergence via:

The time-derivative term is thus expressed by

We note that the gravity term vanishes since \(\vec{\nabla}_H\cdot(\tilde{\rho} g(r)\vec{e_r}) = 0\). Concerning the pressure gradient, one has

The viscosity term then reads

Note

Once again, this viscous term can be greatly simplified in the Boussinesq limit:

Using Eq. (10) then allows to finally write the equation for the pressure \(P_{\ell m n}\):

Note

We note that the terms on the left hand side of (20), (21) and (22) resulting from the viscous term, the pressure gradient, the buoyancy term, and the explicit time derivative completely decouple in spherical harmonic degree and order.

The terms that do not decouple, namely Coriolis force, Lorentz force and advection of momentum, are collected on the right-hand side of (20), (21) and (22) into the forcing term \(\vec{F}\):

Resolving \(\vec{F}\) into potential functions is not required. Its numerical evaluation is discussed below.

### Equation for entropy \(s\)¶

The equation for the entropy (or temperature in the Boussinesq limit) is given by

In this expression, \(j=\vec{\nabla}\times\vec{B}\) is the current. Once again, the numerical evaluation of the right-hand-side (i.e. the non-linear terms) is discussed below.

### Equation for chemical composition \(\xi\)¶

The equation for the chemical composition is given by

Once again, the numerical evaluation of the right-hand-side (i.e. the non-linear term) is discussed below.

### Equation for the poloidal magnetic potential \(g\)¶

The equation for the poloidal magnetic potential is the radial component of the dynamo equation since

The spectral form then reads

### Equation for the toroidal magnetic potential \(h\)¶

The equation for the toroidal magnetic field coefficient reads

Note

We note that the terms on the left hand side of (26) and (27) resulting from the magnetic diffusion term and the explicit time derivative completely decouple in spherical harmonic degree and order.

The dynamo term does not decouple:

We have now derived a full set of equations (20), (21), (22), (24), (26) and (27), each describing the evolution of a single spherical harmonic mode of the six unknown fields (assuming that the terms on the right hand side are given). Each equation couples \(N+1\) Chebyshev coefficients for a given spherical harmonic mode \((\ell,m)\). Typically, a collocation method is employed to solve for the Chebyshev coefficients. This means that the equations are required to be exactly satisfied at \(N-1\) grid points defined by the equations (18) and (19). Excluded are the points \(r=r_i\) and \(r=r_o\), where the boundary conditions provide additional constraints on the set of Chebyshev coefficients.

## Time-stepping schemes¶

Implicit time stepping schemes theoretically offer increased stability and allow for larger time steps. However, fully implicit approaches have the disadvantage that the nonlinear-terms couple all spherical harmonic modes. The potential gain in computational speed is therefore lost at higher resolution, where one very large matrix has to be dealt with rather than a set of much smaller ones. Similar considerations hold for the Coriolis force, one of the dominating forces in the system and therefore a prime candidate for implicit treatment. However, the Coriolis term couples modes \((\ell,m,n)\) with \((\ell+1,m,n)\) and \((\ell-1,m,n)\) and also couples poloidal and toroidal flow potentials. An implicit treatment of the Coriolis term therefore also results in a much larger (albeit sparse) inversion matrix.

We consequently adopt in **MagIC** a mixed implicit/explicit algorithm.
The general differential equation in time can be written in the form

where \(\mathcal{I}\) denotes the terms treated in an implicit time step and \(\mathcal{E}\) the terms treated explicitly, i.e. the nonlinear and Coriolis contributions. In MagIC, two families of time-stepping schemes are supported: IMEX multistep and IMEX multistage methods.

First of all, the IMEX multistep methods correpond to time schemes where the solution results from the combination of several previous steps (such as for instance the Crank-Nicolson/Adams-Bashforth scheme). In this case, a general \(k\)-step IMEX multistep scheme reads

where \(I\) is the identity matrix. The vectors \(\vec{a}\), \(\vec{b}^\mathcal{E}\) and \(\vec{b}^\mathcal{I}\) correspond to the weighting factors of an IMEX multistep scheme. For instance, the commonly-used second-order scheme assembled from the combination of a Crank-Nicolson for the implicit terms and a second-order Adams-Bashforth for the explicit terms (hereafter CNAB2) corresponds to the following vectors: \(\vec{a}=(1,0)\), \(\vec{b}^{\mathcal{I}}=(1/2,1/2)\) and \(\vec{b}^{\mathcal{E}}=(3/2,-1/2)\) for a constant timestep size \(\delta t\).

In addition to CNAB2, MagIC supports several semi-implicit backward differentiation schemes of second, third and fourth order that are known to have good stability properties (heareafter SBDF2, SBDF3 and SBDF4), a modified CNAB2 from Ascher et al. (1995) (termed MODCNAB) and the CNLF scheme (combination of Crank-Nicolson and Leap-Frog for the explicit terms).

MagIC also supports several IMEX Runge-Kutta multistage methods, frequently
called DIRK, an acronym that stands for *Diagonally Implicit Runge Kutta*. For
such schemes, the equation (29) is time-advanced from \(t_n\) to \(t_{n+1}\) by solving \(\nu\) sub-stages

where \(y_i\) is the intermediate solution at the stage \(i\). The matrices \(a_{i,j}^\mathcal{E}\) and \(a_{i,j}^\mathcal{I}\) constitute the so-called Butcher tables that correspond to a DIRK scheme. MagIC supports several second and third order schemes: ARS222 and ARS443 from Ascher et al. (1997), LZ232 from Liu and Zou (2006), PC2 from Jameson et al. (1981) and BPR353 from Boscarino et al. (2013).

In the code the equation (29) is formulated for each unknown spectral coefficient (expect pressure) of spherical harmonic degree \(\ell\) and order \(m\) and for each radial grid point \(r_k\). Because non-linear terms and the Coriolis force are treated explicitly, the equations decouple for all spherical harmonic modes. The different radial grid points, however, couple via the Chebychev modes and form a linear algebraic system of equations that can be solved with standard methods for the different spectral contributions.

For example the respective system of equations for the poloidal magnetic potential \(g\) time advanced with a CNAB2 reads

with

and \({\cal C}_{nk}={\cal C}_n(r_k)\). \(\mathcal{A}_{kn}\) is a matrix that converts the poloidal field modes \(g_{\ell mn}\) to the radial magnetic field \(B_r(r_k,\ell,m)\) for a given spherical harmonic contribution.

Here \(k\) and \(n\) number the radial grid points and the Chebychev coefficients, respectively. Note that the Einstein sum convention is used for Chebychev modes \(n\).

\(\mathcal{I}_{kn}\) is the matrix describing the implicit contribution which is purely diffusive here. Neither \(\mathcal{A}_{kn}\) nor \(\mathcal{I}_{kn}\) depend on time but the former needs to be updated when the time step \(\delta t\) is changed. The only explicit contribution is the nonlinear dynamo term

\(\mathcal{E}_{k\ell m}\) is a one dimensional vector for all spherical harmonic combinations \(\ell m\).

**Courant’s condition** offers a guideline
concerning the value of \(\delta t\), demanding that \(\delta t\) should be smaller
than the advection time between two grid points. Strong Lorentz forces require
an additional stability criterion that is obtained by replacing the flow speed
by Alfvén’s velocity in a modified Courant criterion.
The explicit treatment of the Coriolis force requires that the time step is
limited to a fraction of the rotation period, which may be the relevant
criterion at low Ekman number when flow and magnetic field remain weak.
Non-homogeneous grids and other numerical effects generally require an
additional safety factor in the choice of \(\delta t\).

## Coriolis force and nonlinear terms¶

### Nonlinear terms entering the equation for \(W\)¶

The nonlinear term \({\cal N}^W\) that enters the equation for the poloidal potential (20) contains the radial component of the advection, the Lorentz force and Coriolis force. In spherical coordinate, the two first contributions read:

The Coriolis force can be expressed as a function of the potentials \(W\) and \(Z\) using (4)

The nonlinear terms that enter the equation for the poloidal potential (20) thus reads:

The \(\theta\)-derivative entering the radial component of the Coriolis force is thus the operator \(\vartheta_2\) defined in (12). Using the recurrence relation, one thus finally gets in spherical harmonic space:

To get this expression, we need to first compute \({\cal A}_r\) in the physical space. This
term is computed in the subroutine `get_nl`

in
the module `grid_space_arrays_mod`

. \({\cal A}_r\) is then transformed to the
spectral space by using a Legendre and a Fourier transform to produce \({{\cal A}_r}_\ell^m\).

### Nonlinear terms entering the equation for \(Z\)¶

The nonlinear term \({\cal N}^Z\) that enters the equation for the toroidal potential (21) contains the radial component of the curl of the advection and Coriolis force. The Coriolis force can be rewritten as a function of \(W\) and \(Z\):

Using the \(\vartheta\) operators defined in (12)-(15) then allows to rewrite the Coriolis force in the following way:

The contributions of nonlinear advection and Lorentz forces that enter the equation for the toroidal potential are written this way:

To make use of the recurrence relations (12)-(15), the actual strategy is to follow the following steps:

Compute the quantities \({\cal A}_\phi/r\sin\theta\) and \({\cal A}_\theta/r\sin\theta\) in the physical space. In the code, this step is computed in the subroutine

`get_nl`

in the module`grid_space_arrays_mod`

.Transform \({\cal A}_\phi/r\sin\theta\) and \({\cal A}_\theta/r\sin\theta\) to the spectral space (thanks to a Legendre and a Fourier transform). In MagIC, this step is computed in the modules

`legendre_grid_to_spec`

and`fft`

. After this step \({{\cal A}t}_{\ell}^m\) and \({{\cal A}p}_{\ell}^m\) are defined.Calculate the colatitude and theta derivatives using the recurrence relations:

(34)¶\[\vartheta_1\,{{\cal A}p}_{\ell}^m-\dfrac{\partial {{\cal A}t}_{\ell}^m}{\partial \phi}\,.\]

Using (33) and (34), one thus finally gets

### Nonlinear terms entering the equation for \(P\)¶

The nonlinear term \({\cal N}^P\) that enters the equation for the pressure (22) contains the horizontal divergence of the advection and Coriolis force. The Coriolis force can be rewritten as a function of \(W\) and \(Z\):

Using the \(\vartheta\) operators defined in (14)-(15) then allows to rewrite the Coriolis force in the following way:

The contributions of nonlinear advection and Lorentz forces that enter the equation for pressure are written this way:

To make use of the recurrence relations (12)-(15), we then follow the same three steps as for the advection term entering the equation for \(Z\).

Using (36) and (37), one thus finally gets

### Nonlinear terms entering the equation for \(s\)¶

The nonlinear terms that enter the equation for entropy/temperature (24) are twofold: (i) the advection term, (ii) the viscous and Ohmic heating terms (that vanish in the Boussinesq limit of the Navier Stokes equations).

Viscous and Ohmic heating are directly calculated in the physical space by the
subroutine `get_nl`

in
the module `grid_space_arrays_mod`

. Let’s introduce \({\cal H}\), the sum
of the viscous and Ohmic heating terms.

Expanding this term leads to:

This term is then transformed to the spectral space with a Legendre and a Fourier transform to produce \({\cal H}_\ell^m\).

The treatment of the advection term \(-\vec{u}\cdot\vec{\nabla}s\) is a bit different. It is in a first step rearranged as follows

The quantities that are calculated in the physical space are thus simply the product of
entropy/temperature \(s\) by the velocity components. This defines three variables
defined in the grid space that are computed in the subroutine `get_nl`

:

To get the actual advection term, one must then apply the divergence operator to get:

To make use of the recurrence relations (12)-(15), the actual strategy is then to follow the following steps:

Compute the quantities \(r^2\,\mathcal{US}_r\), \(\mathcal{US}_\phi/r\sin\theta\) and \(\mathcal{US}_\theta/r\sin\theta\) in the physical space. In the code, this step is computed in the subroutine

`get_nl`

in the module`grid_space_arrays_mod`

.Transform \(r^2\,\mathcal{US}_r\), \(\mathcal{US}_\phi/r\sin\theta\) and \(\mathcal{US}_\theta/r\sin\theta\) to the spectral space (thanks to a Legendre and a Fourier transform). In MagIC, this step is computed in the modules

`legendre_grid_to_spec`

and`fft`

. After this step \({\mathcal{US}r}_{\ell}^m\), \({\mathcal{US}t}_{\ell}^m\) and \({\mathcal{US}p}_{\ell}^m\) are defined.Calculate the colatitude and theta derivatives using the recurrence relations:

(40)¶\[-\dfrac{1}{\tilde{\rho}}\left[ \dfrac{1}{r^2}\dfrac{\partial\, {\mathcal{US}r}_\ell^m}{\partial r}+ \vartheta_1\,{\mathcal{US}t}_\ell^m+ \dfrac{\partial\,{\mathcal{US}p}_\ell^m}{\partial \phi}\right]\,.\]

Using (39) and (40), one thus finally gets

### Nonlinear terms entering the equation for \(\xi\)¶

The nonlinear term that enters the equation for chemical composition (25) is the advection term. This term is treated the same way as the advection term that enters the entropy equation. It is in a first step rearranged as follows

The quantities that are calculated in the physical space are thus simply the
product of composition \(\xi\) by the velocity components. This
defines three variables defined in the grid space that are computed in the
subroutine `get_nl`

:

To get the actual advection term, one must then apply the divergence operator to get:

To make use of the recurrence relations (12)-(15), the actual strategy is then to follow the following steps:

Compute the quantities \(r^2\,\mathcal{UX}_r\), \(\mathcal{UX}_\phi/r\sin\theta\) and \(\mathcal{UX}_\theta/r\sin\theta\) in the physical space. In the code, this step is computed in the subroutine

`get_nl`

in the module`grid_space_arrays_mod`

.Transform \(r^2\,\mathcal{UX}_r\), \(\mathcal{UX}_\phi/r\sin\theta\) and \(\mathcal{UX}_\theta/r\sin\theta\) to the spectral space (thanks to a Legendre and a Fourier transform). In MagIC, this step is computed in the modules

`legendre_grid_to_spec`

and`fft`

. After this step \({\mathcal{UX}r}_{\ell}^m\), \({\mathcal{UX}t}_{\ell}^m\) and \({\mathcal{UX}p}_{\ell}^m\) are defined.Calculate the colatitude and theta derivatives using the recurrence relations:

\[-\dfrac{1}{\tilde{\rho}}\left[ \dfrac{1}{r^2}\dfrac{\partial\, {\mathcal{UX}r}_\ell^m}{\partial r}+ \vartheta_1\,{\mathcal{UX}t}_\ell^m+ \dfrac{\partial\,{\mathcal{UX}p}_\ell^m}{\partial \phi}\right]\,.\]

One thus finally gets

### Nonlinear terms entering the equation for \(g\)¶

The nonlinear term that enters the equation for the poloidal potential of the magnetic field (26) is the radial component of the induction term (28). In the following we introduce the electromotive force \({\cal F} = \vec{u}\times\vec{B}\) with its three components

The radial component of the induction term then reads:

To make use of the recurrence relations (12)-(15), we then follow the usual following steps:

Compute the quantities \(r^2\,\mathcal{F}_r\), \(\mathcal{F}_\phi/r\sin\theta\) and \(\mathcal{F}_\theta/r\sin\theta\) in the physical space. In the code, this step is computed in the subroutine

`get_nl`

in the module`grid_space_arrays_mod`

.Transform \(r^2\,\mathcal{F}_r\), \(\mathcal{F}_\phi/r\sin\theta\) and \(\mathcal{F}_\theta/r\sin\theta\) to the spectral space (thanks to a Legendre and a Fourier transform). In MagIC, this step is computed in the modules

`legendre_grid_to_spec`

and`fft`

. After this step \({\mathcal{F}_r}_{\ell}^m\), \({\mathcal{F}_\theta}_{\ell}^m\) and \({\mathcal{F}_\phi}_{\ell}^m\) are defined.Calculate the colatitude and theta derivatives using the recurrence relations:

\[\vartheta_1\,{\mathcal{F}_\phi}_\ell^m- \dfrac{\partial\,{\mathcal{F}_\theta}_\ell^m}{\partial \phi}\,.\]

We thus finally get

### Nonlinear terms entering the equation for \(h\)¶

The nonlinear term that enters the equation for the toroidal potential of the magnetic field (27) is the radial component of the curl of the induction term (28):

To make use of the recurrence relations (12)-(15), we then follow the same steps than for the nonlinear terms that enter the equation for poloidal potential of the magnetic field \(g\):

We thus finally get

## Boundary conditions and inner core¶

### Mechanical boundary conditions¶

Since the system of equations is formulated on a radial grid, boundary conditions can simply be satisfied by replacing the collocation equation at grid points \(r_i\) and \(r_o\) with appropriate expressions. The condition of zero radial flow on the boundaries implies that the poloidal potential has to vanish, i.e. \(W(r_o)=0\) and \(W(r_i)=0\). In Chebychev representation this implies

Note that the summation convention with respect to
radial modes \(n\) is used again.
**The no-slip** condition further requires that the
horizontal flow components also have to vanish, provided
the two boundaries are at rest. This condition is fulfilled for

at the respective boundary. In spectral space these conditions read

and

for all spherical harmonic modes \((\ell,m)\). The conditions (45)-(47) replace the poloidal flow potential equations (20) and the pressure equation (22), respectively, at the collocation points \(r_i\) and \(r_o\).

If the inner-core and/or the mantle are allowed to react to torques, a condition based on the conservation of angular momentum replaces condition (47) for the mode \((\ell =1,m=0)\):

The tensor \(\mathsf{I}\) denotes the moment of inertia of inner core or mantle, respectively, \(\omega\) is the mantle or inner-core rotation rate relative to that of the reference frame, and \(\Gamma_{L,\nu}\) are the respective torques associated with Lorentz or viscous forces. The torques are expressed by

and

where \(\mathrm{d}S = r^2\sin\theta \mathrm{d}\theta\mathrm{d}\phi\) and \(r\in[r_i,r_o]\) in the above expressions. Using the following equality

the viscous torques can be expressed by

where the sign in front depends whether \(r=r_o\) or \(r=r_i\).

**Free-slip boundary conditions** require that the viscous stress vanishes, which
in turn implies that the non-diagonal components \(\mathsf{Sr}_{r\phi}\) and
\(\mathsf{S}_{r\theta}\) of the rate-of-strain tensor vanish.
Translated to the spectral representation this requires

We show the derivation for the somewhat simpler Boussinesq approximation which yields the condition

where the index H denotes the horizonal flow components. In terms of poloidal and toroidal components this implies

and

which can be fulfilled with

and

In spectral representation this then reads

### Thermal boundary conditions¶

For Entropy or temperature in the Boussinesq approximation either fixed value of fixed flux conditions are used. The former implies

at \(r_i\) and/or \(r_o\), while the latter means

In spectral representation for example the respective entropy condition read

Appropriate constant values need to be chosen and are instrumental in driving the dynamo when flux conditions are imposed.

### Boundary conditions for chemical composition¶

For the chemical composition, either the value or the flux is imposed at the boundaries. The former implies:

at \(r_i\) and/or \(r_o\), while the latter means

In spectral representation, this then reads

### Magnetic boundary conditions and inner core¶

Three different magnetic boundary conditions are implemented in MagIC. The most simple one is the conceptual condition at the boundary to an infinite conductor. Surface current in this conductor will prevent the internally produced magnetic field from penetrating so that the field has to vanish at the boundary. The condition are thus the same as for a rigid flow (with boundaries at rest). We only provide the spectral representation here:

Note that the summation convention with respect to
radial modes \(n\) is used again.
**The no-slip** condition further requires that the
horizontal flow components also have to vanish, provided
the two boundaries are at rest. This condition is fulfilled for

More complex are the conditions to an electrical insulator. Here we actually use matching condition to a potential field condition that are formulated like boundary conditions. Since the electrical currents have to vanish in the insulator we have \(\nabla\times\vec{B}\), which means that the magnetic field is a potential field \(\vec{B}^I=-\vec{\nabla} V\) with \(\Delta V=0\). This Laplace equation implies a coupling between radial and horizontal derivatives which is best solved in spectral space. Two potential contributions have to be considered depending whether the field is produced above the interface radius \(r_{BC}\) or below. We distinguish these contributions with upper indices I for internal or below and E for external or above. The total potential then has the form:

with the two spectral potential representations \(V_{\ell m}^I\) and \(V_{\ell m}^E\). This provides well defined radial derivative for both field contributions. For boundary \(r_o\) we have to use the first contribution and match the respective field as well as its radial derivative to the dynamo solution. The toroidal field cannot penetrate the insulator and thus simply vanishes which yields \(h=0\) or

in spectral space. The poloidal field then has to match the potential field which implies

for the horizontal components and

for the radial. In spectral space these condition can be reduce to

Combining both allows to eliminate the potential and finally leads to the spectral condition used in MagIC:

Analogous consideration lead to the respective condition at the interface to an insulating inner core:

If the inner core is modelled as an electrical conductor, a simplified dynamo equation has to be solved in which the fluid flow is replaced by the solid-body rotation of the inner core. The latter is described by a single toroidal flow mode \((\ell=1,m=0)\). The resulting nonlinear terms can be expressed by a simple spherical harmonic expansion, where the superscript \(I\) denotes values in the inner core and \(\omega_I\) its differential rotation rate:

The expensive back and forth transformations between spherical harmonic and grid representations are therefore not required for advancing the inner-core magnetic field in time.

In the inner core the magnetic potentials are again conveniently expanded into Chebyshev polynomials. The Chebyshev variable \(x\) spans the whole diameter of the inner core, so that grid points are dense near the inner-core boundary but sparse in the center. The mapping is given by:

Each point in the inner core is thus represented twice, by grid points \(({r,\theta,\phi})\) and \(({-r,\pi-\theta,\phi+\pi})\). Since both representations must be identical, this imposes a symmetry constraint that can be fulfilled when the radial expansion comprises only polynomials of even order:

An equivalent expression holds for the toroidal potential in the inner core. FFTs can again by employed efficiently for the radial transformation, using the \(M\) extrema of \({\cal C}_{2 M-1}(r)\) with \(x>0\) as grid points.

The sets of spectral magnetic field equations for the inner and the outer core
are coupled via continuity equations for the magnetic field and the
horizontal electric field.
Continuity of the magnetic field is assured by (**i**) continuity of the
toroidal potential, (**ii**) continuity of the poloidal potential, and
(**iii**) continuity of the radial derivative of
the latter. Continuity of the horizontal electric field demands (**iv**) that
the radial derivative of the toroidal potential is continuous, provided
that the horizontal flow and the electrical conductivity are continuous
at the interface.
These four conditions replace the spectral equations
(26), (27) on the outer-core side and
equations (50), (51) on the inner-core side.
Employing free-slip conditions or allowing for electrical conductivity differences
between inner and outer core leads to more complicated and even non-linear matching
conditions.