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

\[\vec{v} = \vec{\nabla}\times\left(\vec{\nabla}\times W\,\vec{e_r}\right) + \vec{\nabla}\times Z\,\vec{e_r}.\]

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

(1)\[\begin{split}\boxed{ \begin{aligned} \tilde{\rho} \vec{u} & = \vec{\nabla}\times(\vec{\nabla}\times W\,\vec{e_r}) + \vec{\nabla}\times Z\,\vec{e_r}\,, \\ \vec{B} & = \vec{\nabla}\times(\vec{\nabla}\times g\,\vec{e_r}) + \vec{\nabla}\times h\,\vec{e_r}\,. \end{aligned}}\end{split}\]

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:

(2)\[\begin{split}\begin{aligned} \vec{e_r}\cdot \tilde{\rho}\vec{u} &= - \Delta_H\,W, \\ \vec{e_r}\cdot\left(\vec{\nabla}\times\vec{u}\right) & = - \Delta_H\,Z, \end{aligned}\end{split}\]

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

(3)\[\Delta_H= \frac{1}{r^{2}\sin{\theta}} \frac{\partial}{\partial\theta}\left(\sin{\theta}\frac{\partial}{\partial\theta}\right) + \frac{1}{r^{2}\sin^2{\theta}} \frac{\partial^{2}}{\partial^{2}\phi}.\]

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

(4)\[\tilde{\rho}\vec{u} = -(\Delta_H W)\,\vec{e_r} + \left( \dfrac{1}{r} \dfrac{\partial^2 W}{\partial r \partial \theta} + \dfrac{1}{r\sin\theta}\dfrac{\partial Z}{\partial \phi}\right)\,\vec{e_\theta} +\left(\dfrac{1}{r\sin\theta}\dfrac{\partial^2 W}{\partial r\partial \phi}- \dfrac{1}{r}\dfrac{\partial Z}{\partial\theta} \right)\,\vec{e_\phi},\]

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

(5)\[\begin{split}\begin{aligned} \vec{\nabla}\times\tilde{\rho}\vec{u} = &-(\Delta_H Z)\,\vec{e_r}+ \left[-\dfrac{1}{r\sin\theta}\dfrac{\partial}{\partial\phi}\left( \dfrac{\partial^2}{\partial r^2}+\Delta_H \right) W + \dfrac{1}{r} \dfrac{\partial^2 Z}{\partial r\partial\theta}\right]\vec{e_\theta} \\ &+\left[\dfrac{1}{r}\dfrac{\partial}{\partial\theta}\left( \dfrac{\partial^2}{\partial r^2}+\Delta_H\right) W + \dfrac{1}{r \sin\theta} \dfrac{\partial^2 Z}{\partial r\partial\phi}\right]\vec{e_\phi}, \end{aligned}\end{split}\]

Using the horizontal part of the divergence operator

\[\vec{\nabla}_H = \dfrac{1}{r\sin\theta} \dfrac{\partial}{\partial \theta}\sin\theta\;\vec{e}_\theta + \dfrac{1}{r\sin\theta} \dfrac{\partial}{\partial \phi}\;\vec{e}_\phi\]

above expressions can be simplified to

\[\tilde{\rho}\vec{u} = -\Delta_H\;\vec{e_r}\; W + \vec{\nabla}_H \dfrac{\partial}{\partial r}\;W + \vec{\nabla}_H\times\vec{e}_r\;Z\]

and

\[\nabla\times\tilde{\rho}\vec{u} = -\Delta_H\;\vec{e}_r\;Z + \vec{\nabla}_H \dfrac{\partial}{\partial r}\;Z - \vec{\nabla}_H\times\Delta_H\vec{e}_r\;W\;\;.\]

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\):

\[Y_\ell^m(\theta,\phi) = P_{\ell}^m(\cos{\theta})\,e^{i m \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

(6)\[\int_{0}^{2\pi} d\,\phi \int_{0}^{\pi} \sin{\theta}\, d\theta\; Y_\ell^m(\theta,\phi)\,Y_{\ell^\prime}^{m^\prime} (\theta,\phi) \; = \; \delta_{\ell \ell^\prime}\delta^{m m^\prime}.\]

This means that

\[Y_{\ell}^{m}(\theta,\phi) = \left(\dfrac{(2\ell+1)}{4\pi}\dfrac{(\ell-|m|)!}{(\ell+|m|)!}\right)^{1/2} P_\ell^m(\cos{\theta})\,e^{i m \phi},\]

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

(7)\[g(r,\theta,\phi) = \sum_{\ell=0}^{\ell_{max}}\sum_{m=-\ell}^{\ell} g_{\ell m}(r)\,Y_{\ell}^{m}(\theta,\phi),\]

with

(8)\[g_{\ell m}(r) = \frac{1}{\pi}\,\int_{0}^{\pi} d \theta \sin{\theta}\; g_m(r,\theta)\; P_\ell^m(\cos{\theta}),\]

(9)\[g_{m}(r,\theta) = \frac{1}{2\pi}\,\int_{0}^{2\pi} d \phi\; g(r,\theta,\phi)\; e^{- i m \phi} .\]

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)

\[g_{\ell m}(r) = \frac{1}{N_{\theta}} \sum_{j=1}^{N_{\theta}}\,w_j\,g_m(r,\theta_j)\; P_\ell^m(\cos{\theta_j}),\]

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

(10)\[\boxed{ \Delta_H Y_{\ell}^{m} = -\dfrac{\ell(\ell+1)}{r^2}\,Y_{\ell}^{m}\,. }\]

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

\[\vartheta_1 = \dfrac{1}{\sin\theta}\dfrac{\partial}{\partial\theta}\sin^2\theta =\sin\theta\dfrac{\partial}{\partial\theta}+2\cos\theta\,.\]

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

\[\vartheta_1 = (\ell+2)\,c_{\ell+1}^m\,P_{\ell+1}^m(\cos\theta) -(\ell-1)\,c_\ell^m\,P_{\ell-1}^m(\cos\theta),\]

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

(11)\[c_\ell^m = \sqrt{\dfrac{(\ell+m)(\ell-m)}{(2\ell-1)(2\ell+1)}}\,.\]

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

\[\dfrac{1}{\sin\theta}\dfrac{\partial}{\partial\theta}(\sin\theta\,f(\theta))\,.\]

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

\[F(\theta)=f(\theta)/\sin\theta\,,\]

so that

\[\dfrac{1}{\sin\theta}\dfrac{\partial}{\partial\theta}[\sin\theta\,f(\theta)] =\vartheta_1 F(\theta)\,.\]

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

(12)\[\boxed{ \int_0^\pi d\theta\,\sin\theta\,P_\ell^m\vartheta_1\sum_{\ell'}F_{\ell'}^m P_{\ell'}^m =(\ell+1)\,c_{\ell}^m\,F_{\ell-1}^m-\ell\,c_{\ell+1}^m\,F_{\ell+1}^m\,.}\]

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

\[\int_0^\pi d\theta\,\sin\theta\,P_\ell^m P_{\ell'}^m = \delta_{\ell \ell'}\,.\]

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

\[\vartheta_2 = \sin\theta\dfrac{\partial}{\partial\theta}\,.\]

The action of this operator on the Legendre polynomials reads

\[\vartheta_2 P_\ell^m(\cos\theta)=\ell\,c_{\ell+1}^m\,P_{\ell+1}^m(\cos\theta) -(\ell+1)\,c_\ell^m\,P_{\ell-1}^m(\cos\theta)\,,\]

so that

(13)\[\boxed{ \int_0^\pi d\theta\,\sin\theta \,P_\ell^m\vartheta_2\sum_{\ell'}f_{\ell'}^m P_{\ell'}^m =(\ell-1)\,c_{\ell}^m\,f_{\ell-1}^m-(\ell+2)\,c_{\ell+1}^m\,f_{\ell+1}^m\,.}\]

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:

\[\vartheta_3 = \sin\theta\dfrac{\partial}{\partial\theta}+\cos\theta\,L_H\,,\]

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

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

\[\vartheta_3 P_\ell^m(\cos\theta)=\ell(\ell+1)\,c_{\ell+1}^m\,P_{\ell+1}^m(\cos\theta) +(\ell-1)(\ell+1)\,c_\ell^m\,P_{\ell-1}^m(\cos\theta)\,,\]

which results into:

(14)\[\boxed{ \int_0^\pi d\theta\,\sin\theta\,P_\ell^m\vartheta_3\sum_{\ell'}f_{\ell'}^m P_{\ell'}^m =(\ell-1)(\ell+1)\,c_{\ell}^m\,f_{\ell-1}^m+\ell(\ell+2)\,c_{\ell+1}^m\,f_{\ell+1}^m\,.}\]

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:

\[\vartheta_4 = \dfrac{1}{\sin\theta}\dfrac{\partial}{\partial\theta}\sin^2\theta\,L_H =\vartheta1\,L_H\,,\]

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

\[\vartheta_4 P_\ell^m(\cos\theta)=\ell(\ell+1)(\ell+2)\,c_{\ell+1}^m\,P_{\ell+1}^m(\cos\theta) -\ell(\ell-1)(\ell+1)\,c_\ell^m\,P_{\ell-1}^m(\cos\theta)\,,\]

which results into:

(15)\[\boxed{ \int_0^\pi d\theta\,\sin\theta\,P_\ell^m\vartheta_4\sum_{\ell'}f_{\ell'}^m P_{\ell'}^m =\ell(\ell-1)(\ell+1)\,c_{\ell}^m\,f_{\ell-1}^m-\ell(\ell+1)(\ell+2)\,c_{\ell+1}^m\,f_{\ell+1}^m\,.}\]

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

\[{\cal C}_n(x)\approx\cos\left[n\,\arccos(x)\right]\quad -1\leq x \leq 1\,.\]

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

(16)\[g_{\ell m}(r) = \sum_{n=0}^{N} g_{\ell mn}\;{\cal C}_n(r) ,\]

with

(17)\[g_{\ell mn} = \frac{2-\delta_{n0}}{\pi}\int_{-1}^{1} \frac{d x\, g_{\ell m}(r(x))\;{\cal C}_n(x)}{\sqrt{1-x^2}} .\]

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

(18)\[x(r)= 2 \frac{r-r_i}{r_o-r_i} - 1 .\]

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,

(19)\[x_k=\cos{\left[\pi \frac{(k-1)}{N_r-1}\right]}\;\;\;,\;\;\; k=1,2,\ldots,N_r ,\]

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

\[{\cal C}_{nk} = {\cal C}_n(x_k)=\cos{\left[\pi \frac{ n (k-1)}{N_r-1}\right]} .\]

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):

\[\tilde{\rho}\vec{e_r}\cdot\left(\dfrac{\partial \vec{u}}{\partial t}\right) = \dfrac{\partial}{\partial t}(\vec{e_r}\cdot\tilde{\rho}\vec{u})=-\Delta_H\dfrac{\partial W}{\partial t}.\]

For the viscosity term, one gets

\[\begin{split}\begin{aligned} \vec{e_r}\cdot\vec{\nabla}\cdot \mathsf{S} = & -\nu\,\Delta_H\left[\dfrac{\partial^2 W} {\partial r^2} +\left\lbrace 2\dfrac{d\ln\nu}{dr}-\dfrac{1}{3}\dfrac{d\ln\tilde{\rho}}{dr}\right\rbrace \dfrac{\partial W}{\partial r} \right. \\ & -\left. \left\lbrace -\Delta_H + \dfrac{4}{3}\left(\dfrac{d^2\ln\tilde{\rho}}{dr^2} +\dfrac{d\ln\nu}{dr} \dfrac{d\ln\tilde{\rho}}{dr} + \dfrac{1}{r}\left[3\dfrac{d\ln\nu}{dr}+ \dfrac{d\ln\tilde{\rho}}{dr}\right] \right) \right\rbrace W\right], \end{aligned}\end{split}\]

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

\[\vec{e_r}\cdot\Delta \vec{u} = -\Delta_H\left[\dfrac{\partial^2 W}{\partial r^2} +\Delta_H\,W\right]\,.\]

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

(20)\[\begin{split}\boxed{ \begin{aligned} E\,\dfrac{\ell(\ell+1)}{r^2}\left[\left\lbrace\dfrac{\partial}{\partial t} + \nu\,\dfrac{\ell(\ell+1)}{r^2} + \dfrac{4}{3}\,\nu\,\left(\dfrac{d^2\ln\tilde{\rho}}{dr^2} +\dfrac{d\ln\nu}{dr} \dfrac{d\ln\tilde{\rho}}{dr} + \dfrac{1}{r}\left[3\dfrac{d\ln\nu}{dr}+ \dfrac{d\ln\tilde{\rho}}{dr}\right] \right)\right\rbrace\right. & \,{\cal C}_n & \\ -\nu\,\left\lbrace 2\dfrac{d\ln\nu}{dr}-\dfrac{1}{3}\dfrac{d\ln\tilde{\rho}}{dr}\right\rbrace &\,{\cal C}'_n & \\ -\nu & \,{\cal C}''_n \left. \phantom{\dfrac{d\nu}{dr}}\right]& W_{\ell m n} \\ + \left[{\cal C}'_n -\dfrac{d\ln\tilde{\rho}}{dr}{\cal C}_n\right] & & P_{\ell m n} \\ - \left[\dfrac{Ra\,E}{Pr}\,\tilde{\rho}\,g(r)\right] & \,{\cal C}_n & s_{\ell m n} \\ - \left[\dfrac{Ra_\xi\,E}{Sc}\,\tilde{\rho}\,g(r)\right] & \,{\cal C}_n & \xi_{\ell m n} \\ = {\cal N}^W_{\ell m} = \int d\Omega\,{Y_{\ell}^{m}}^\star\,{\cal N}^W =\int d\Omega\,{Y_{\ell}^{m}}^\star\,\vec{e_r}\cdot \vec{F}\,. & & \end{aligned}}\end{split}\]

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.

See also

The exact computation of the linear terms of (20) are coded in the subroutines get_wpMat

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):

\[\vec{e_r}\cdot\vec\nabla\times\left(\dfrac{\partial\tilde{\rho}\vec{u}}{\partial t}\right)= \dfrac{\partial}{\partial t}(\vec{e_r}\cdot\vec{\nabla}\times\tilde{\rho} \vec{u})=-\dfrac{\partial}{\partial t}(\Delta_H Z) = -\Delta_H\dfrac{\partial Z}{\partial t}\,.\]

The pressure gradient, one has

\[\vec{\nabla}\times \left[\tilde{\rho}\vec{\nabla}\left(\dfrac{p'}{\tilde{\rho}}\right)\right] = \vec{\nabla} \tilde{\rho} \times \vec{\nabla}\left(\dfrac{p'}{\tilde{\rho}}\right) + \underbrace{\tilde{\rho} \vec{\nabla} \times \left[\vec{\nabla}\left( \dfrac{p'}{\tilde{\rho}} \right)\right]}_{=0}\,.\]

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.

\[\begin{split}\begin{aligned} \vec{e_r}\cdot\vec{\nabla}\times\left[\vec{\nabla}\cdot\mathsf{S}\right] = & -\nu\,\Delta_H\left[\dfrac{\partial^2 Z}{\partial r^2} +\left(\dfrac{d\ln\nu}{dr}-\dfrac{d\ln\tilde{\rho}}{dr}\right)\,\dfrac{\partial Z}{\partial r} \right.\\ & \left. - \left(\dfrac{d\ln\nu}{dr}\dfrac{d\ln\tilde{\rho}}{dr}+ \dfrac{2}{r}\dfrac{d\ln\nu}{dr}+ \dfrac{d^2\ln\tilde{\rho}}{dr^2}+\dfrac{2}{r} \dfrac{d\ln\tilde{\rho}}{dr}-\Delta_H\right) Z \right]. \end{aligned}\end{split}\]

Note

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

\[\vec{e_r}\cdot\vec{\nabla}\times\left(\Delta \vec{u}\right) = -\Delta_H\left[\dfrac{\partial^2 Z}{\partial r^2} +\Delta_H\,Z\right]\,.\]

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

(21)\[\begin{split}\boxed{ \begin{aligned} E\,\dfrac{\ell(\ell+1)}{r^2}\left[\left\lbrace\dfrac{\partial}{\partial t} + \nu\,\dfrac{\ell(\ell+1)}{r^2} + \nu\,\left(\dfrac{d\ln\nu}{dr}\dfrac{d\ln\tilde{\rho}}{dr}+ \dfrac{2}{r}\dfrac{d\ln\nu}{dr}+ \dfrac{d^2\ln\tilde{\rho}}{dr^2}+\dfrac{2}{r} \dfrac{d\ln\tilde{\rho}}{dr}\right)\right\rbrace\right. & \,{\cal C}_n & \\ -\nu\,\left(\dfrac{d\ln\nu}{dr}-\dfrac{d\ln\tilde{\rho}}{dr}\right) &\,{\cal C}'_n & \\ -\nu & \,{\cal C}''_n \left. \phantom{\dfrac{d\nu}{dr}}\right]& Z_{\ell m n} \\ = {\cal N}^Z_{\ell m} = \int d\Omega\,{Y_{\ell}^{m}}^\star\,{\cal N}^Z = \int d\Omega\,{Y_{\ell}^{m}}^\star\,\vec{e_r}\cdot \left(\vec{\nabla} \times\vec{F}\right)\,. & & \end{aligned}}\end{split}\]

See also

The exact computation of the linear terms of (21) are coded in the subroutines get_zMat

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

\[\vec{\nabla}_H\cdot\vec{a} = r\sin \dfrac{\partial (\sin\theta\,a_\theta)}{\partial \theta} +r\sin \dfrac{\partial a_\phi}{\partial \phi}.\]

This relates to the total divergence via:

\[\vec{\nabla}\cdot\vec{a}= \dfrac{1}{r^2}\dfrac{\partial(r^2 a_r)}{\partial r}+ \vec{\nabla}_H\cdot\vec{a}.\]

The time-derivative term is thus expressed by

\[\begin{split}\begin{aligned} \vec{\nabla}_H\cdot\left(\tilde{\rho}\dfrac{\partial \vec{u}}{\partial t}\right) &= \dfrac{\partial}{\partial t}\left[\vec{\nabla}_H\cdot(\tilde{\rho}\vec{u} )\right], \\ & = \dfrac{\partial}{\partial t}\left[\vec{\nabla}\cdot(\tilde{\rho}\vec{u}) -\dfrac{1}{r^2}\dfrac{\partial(r^2\tilde{\rho} u_r)}{\partial r}\right], \\ & = -\dfrac{\partial}{\partial t}\left[\dfrac{\partial (\tilde{\rho} u_r)}{\partial r} +\dfrac{2\tilde{\rho} u_r}{r}\right], \\ & = \dfrac{\partial}{\partial t}\left[\dfrac{\partial (\Delta_H W)}{\partial r} +\dfrac{2}{r}\Delta_H W\right], \\ & = \Delta_H\dfrac{\partial}{\partial t}\left(\dfrac{\partial W}{\partial r}\right). \end{aligned}\end{split}\]

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

\[-\vec{\nabla}_H\cdot\left[\tilde{\rho} \vec{\nabla}\left(\dfrac{p'}{\tilde{\rho}} \right)\right] = -\left\lbrace\vec{\nabla}\cdot\left[\tilde{\rho} \vec{\nabla} \left(\dfrac{p'}{\tilde{\rho}}\right)\right]- \dfrac{1}{r^2}\dfrac{\partial}{\partial r}\left[ r^2 \tilde{\rho} \dfrac{\partial}{\partial r}\left(\dfrac{p'}{\tilde{\rho}}\right)\right] \right\rbrace = -\Delta_H \, p'.\]

The viscosity term then reads

\[\begin{split}\begin{aligned} \vec{\nabla}_H\cdot \left( \vec{\nabla}\cdot\mathsf{S} \right) = & \nu\,\Delta_H\left[ \dfrac{\partial^3 W}{\partial r^3} + \left(\dfrac{d\ln\nu}{dr}- \dfrac{d\ln\tilde{\rho}}{dr}\right) \dfrac{\partial^2 W}{\partial r^2} \right. \\ & - \left[\dfrac{d^2\ln\tilde{\rho}}{dr^2} + \dfrac{d\ln\nu}{dr}\dfrac{d\ln\tilde{\rho}}{dr}+ \dfrac{2}{r}\left(\dfrac{d\ln\nu}{dr}+\dfrac{d\ln\tilde{\rho}}{dr}\right) -\Delta_H \right]\dfrac{\partial W}{\partial r} \\ & \left. -\left( \dfrac{2}{3}\dfrac{d\ln\tilde{\rho}}{dr}+\dfrac{2}{r}+\dfrac{d\ln\nu}{dr} \right)\Delta_H\,W \right]. \end{aligned}\end{split}\]

Note

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

\[\vec{\nabla}_H\cdot\left(\Delta \vec{u}\right) = -\Delta_H\left[\dfrac{\partial^3 W}{\partial r^3} +\Delta_H\,\dfrac{\partial W}{\partial r}-\dfrac{2}{r}\Delta_H\,W\right]\,.\]

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

(22)\[\begin{split}\boxed{ \begin{aligned} E\,\dfrac{\ell(\ell+1)}{r^2}\left[ -\nu\,\left( \dfrac{2}{3}\dfrac{d\ln\tilde{\rho}}{dr}+\dfrac{2}{r}+\dfrac{d\ln\nu}{dr} \right)\dfrac{\ell(\ell+1)}{r^2} \right. & \,{\cal C}_n & \\ \left\lbrace\dfrac{\partial}{\partial t} + \nu\,\dfrac{\ell(\ell+1)}{r^2} + \nu\,\left[\dfrac{d^2\ln\tilde{\rho}}{dr^2}+ \dfrac{d\ln\nu}{dr}\dfrac{d\ln\tilde{\rho}}{dr}+ \dfrac{2}{r}\left(\dfrac{d\ln\nu}{dr}+\dfrac{d\ln\tilde{\rho}}{dr}\right)\right]\right\rbrace & \,{\cal C}'_n & \\ -\nu\,\left( \dfrac{d\ln\nu}{dr}-\dfrac{d\ln\tilde{\rho}}{dr} \right) &\,{\cal C}''_n & \\ -\nu & \,{\cal C}'''_n \left. \phantom{\dfrac{d\nu}{dr}}\right]& W_{\ell m n} \\ + \left[\dfrac{\ell(\ell+1)}{r^2}\right] & \,{\cal C}_n & P_{\ell m n} \\ = {\cal N}^P_{\ell m} = -\int d\Omega\,{Y_{\ell}^{m}}^\star\,{\cal N}^P=-\int d\Omega\,{Y_{\ell}^{m}}^\star\,\vec{\nabla}_H\cdot\vec{F}\,. & & \end{aligned}}\end{split}\]

See also

The exact computation of the linear terms of (22) are coded in the subroutines get_wpMat

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}\):

(23)\[\vec{F}=-2\,\tilde{\rho}\,\vec{e_z}\times\vec{u} - E\,\tilde{\rho}\, \vec{u}\cdot\vec{\nabla}\,\vec{u} +\frac{1}{Pm}\left(\vec{\nabla}\times\vec{B}\right)\times\vec{B}\,.\]

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

(24)\[\begin{split}\boxed{ \begin{aligned} \dfrac{1}{Pr}\left[\left(Pr\dfrac{\partial}{\partial t} + \kappa\,\dfrac{\ell(\ell+1)}{r^2} \right)\right. & \,{\cal C}_n & \\ -\kappa\,\left(\dfrac{d\ln\kappa}{dr}+\dfrac{d\ln\tilde{\rho}}{dr}+ +\dfrac{dln\tilde{T}}{dr}+\dfrac{2}{r}\right) &\,{\cal C}'_n & \\ -\kappa & \,{\cal C}''_n \left. \phantom{\dfrac{d\nu}{dr}}\right]& s_{\ell m n} \\ = {\cal N}^S_{\ell m} = \int d\Omega\,{Y_{\ell}^{m}}^\star\,{\cal N}^S = \int d\Omega\,{Y_{\ell}^{m}}^\star\,\left[-\vec{u}\cdot\vec{\nabla}s+ \dfrac{Pr\,Di}{Ra}\dfrac{1}{\tilde{\rho}\tilde{T}}\left(\Phi_\nu+ \dfrac{\lambda}{Pm^2\,E}\,j^2\right) \right]\,. & & \end{aligned}}\end{split}\]

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.

See also

The exact computation of the linear terms of (24) are coded in the subroutines get_sMat

Equation for chemical composition \(\xi\)

The equation for the chemical composition is given by

(25)\[\begin{split}\boxed{ \begin{aligned} \dfrac{1}{Sc}\left[\left(Sc\dfrac{\partial}{\partial t} + \kappa_\xi\,\dfrac{\ell(\ell+1)}{r^2} \right)\right. & \,{\cal C}_n & \\ -\kappa_\xi\,\left(\dfrac{d\ln\kappa_\xi}{dr}+\dfrac{d\ln\tilde{\rho}}{dr}+ +\dfrac{2}{r}\right) &\,{\cal C}'_n & \\ -\kappa_\xi & \,{\cal C}''_n \left. \phantom{\dfrac{d\nu}{dr}}\right]& \xi_{\ell m n} \\ = {\cal N}^\xi_{\ell m} = \int d\Omega\,{Y_{\ell}^{m}}^\star\,{\cal N}^\xi = \int d\Omega\,{Y_{\ell}^{m}}^\star\,\left[-\vec{u}\cdot\vec{\nabla}\xi \right]\,. & & \end{aligned}}\end{split}\]

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

See also

The exact computation of the linear terms of (25) are coded in the subroutines get_xiMat

Equation for the poloidal magnetic potential \(g\)

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

\[\vec{e_r}\cdot\left(\dfrac{\partial \vec{B}}{\partial t}\right) = \dfrac{\partial}{\partial t}(\vec{e_r}\cdot\vec{B})=-\Delta_H\dfrac{\partial g}{\partial t}.\]

The spectral form then reads

(26)\[\begin{split}\boxed{ \begin{aligned} \dfrac{\ell(\ell+1)}{r^2}\left[\left(\dfrac{\partial}{\partial t} + \dfrac{1}{Pm}\lambda\,\dfrac{\ell(\ell+1)}{r^2} \right)\right. & \,{\cal C}_n & \\ -\dfrac{1}{Pm}\,\lambda & \,{\cal C}''_n \left. \phantom{\dfrac{d\nu}{dr}}\right]& g_{\ell m n} \\ = {\cal N}^g_{\ell m} = \int d\Omega\,{Y_{\ell}^{m}}^\star\,{\cal N}^g=\int d\Omega\,{Y_{\ell}^{m}}^\star\,\vec{e_r}\cdot \vec{D}\,. & & \end{aligned}}\end{split}\]

See also

The exact computation of the linear terms of (26) are coded in the subroutines get_bMat

Equation for the toroidal magnetic potential \(h\)

The equation for the toroidal magnetic field coefficient reads

(27)\[\begin{split}\boxed{ \begin{aligned} \dfrac{\ell(\ell+1)}{r^2}\left[\left(\dfrac{\partial}{\partial t} + \dfrac{1}{Pm}\lambda\,\dfrac{\ell(\ell+1)}{r^2} \right)\right. & \,{\cal C}_n & \\ -\dfrac{1}{Pm}\,\dfrac{d\lambda}{dr} &\,{\cal C}'_n & \\ -\dfrac{1}{Pm}\,\lambda & \,{\cal C}''_n \left. \phantom{\dfrac{d\nu}{dr}}\right]& h_{\ell m n} \\ = {\cal N}^h_{\ell m}= \int d\Omega\,{Y_{\ell}^{m}}^\star\,{\cal N}^h = \int d\Omega\,{Y_{\ell}^{m}}^\star\,\vec{e_r}\cdot \left(\vec{\nabla}\times \vec{D}\right)\,. & & \end{aligned}}\end{split}\]

See also

The exact computation of the linear terms of (27) are coded in the subroutines get_bMat

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:

(28)\[\vec{D}=\vec{\nabla}\times\left(\vec{u}\times\vec{B}\right)\,.\]

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

(29)\[\dfrac{\partial }{\partial t} y = \mathcal{I}(y,t) + \mathcal{E}(y,t),\quad y(t_o)=y_o\;.\]

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 correspond 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

\[\left(I-b_o^{\mathcal I} \delta t\,\mathcal{I}\right)y_{n+1}=\sum_{j=1}^k a_j y_{n+1-j}+\delta t\sum_{j=1}^k \left(b_j^\mathcal{E} \mathcal{E}_{n+1-j}+b_{j}^\mathcal{I}\mathcal{I}_{n+1-j}\right)\,,\]

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

\[(I-a_{ii}^\mathcal{I}\delta t \mathcal{I})y_{i}=y_n+\delta t \sum_{j=1}^{i-1}\left(a_{i,j}^{\mathcal{E}}\mathcal{E}_j+a_{i,j}^\mathcal{I}\mathcal{I}_j\right), \quad 1\leq i\leq \nu,\]

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 Chebyshev 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

(30)\[\left( \mathcal{A}_{kn} - \dfrac{1}{2}\,\delta t\,\mathcal{I}_{kn}\right)\;g_{\ell mn}(t+\delta t) = \left( \mathcal{A}_{kn} + \dfrac{1}{2}\,\delta t\,\mathcal{I}_{kn} \right)\;g_{\ell mn}(t) + \frac{3}{2}\,\delta t\,\mathcal{E}_{k\ell m}(t) - \frac{1}{2}\,\delta t\,\mathcal{E}_{k\ell m}(t-\delta t)\]

with

\[\mathcal{A}_{kn} = \dfrac{\ell (\ell+1)}{r_k^2}\,{\cal C}_{nk}\;,\]

\[\mathcal{I}_{kn}=\dfrac{\ell(\ell+1)}{r_k^2}\,\dfrac{1}{Pm}\left({\cal C}''_{nk} - \dfrac{\ell(\ell+1)}{r_k^2}\; {\cal C}_{nk}\right)\;,\]

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 Chebyshev coefficients, respectively. Note that the Einstein sum convention is used for Chebyshev 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}(t)= {\cal N}_{k\ell m}^g = \int d\Omega\; {Y_{\ell}^{m}}^\star\; \vec{e_r} \cdot \vec{D}(t,r_k,\theta,\phi)\;\; .\]

\(\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:

(31)\[\begin{split}\tilde{\rho}\left(\vec{u}\cdot\vec{\nabla}\vec{u}\right)= \left\lbrace \begin{aligned} {\cal A}_r \\ {\cal A}_\theta \\ {\cal A}_\phi \end{aligned} \right\rbrace = \left\lbrace \begin{aligned} -\tilde{\rho}\,E\,\left( u_r\dfrac{\partial u_r}{\partial r}+ \dfrac{u_\theta}{r}\dfrac{\partial u_r}{\partial \theta}+ \dfrac{u_\phi}{r\sin\theta}\dfrac{\partial u_r}{\partial \phi} -\dfrac{u_\theta^2+u_\phi^2}{r}\right)+ \dfrac{1}{Pm}\left(j_\theta\,B_\phi-j_\phi\,B_\theta\right)\, , \\ -\tilde{\rho}\,E\,\left( u_r\dfrac{\partial u_\theta}{\partial r}+ \dfrac{u_\theta}{r}\dfrac{\partial u_\theta}{\partial \theta} + \dfrac{u_\phi}{r\sin\theta}\dfrac{\partial u_\theta}{\partial \phi}+ \dfrac{u_r u_\theta}{r}-\dfrac{\cos\theta}{r\sin\theta}u_\phi^2\right)+ \dfrac{1}{Pm}\left(j_\phi\,B_r-j_r\,B_\phi\right)\, ,\\ -\tilde{\rho}\,E\,\left( u_r\dfrac{\partial u_\phi}{\partial r}+ \dfrac{u_\theta}{r}\dfrac{\partial u_\phi}{\partial \theta} + \dfrac{u_\phi}{r\sin\theta}\dfrac{\partial u_\phi}{\partial \phi}+ \dfrac{u_r u_\phi}{r} +\dfrac{\cos\theta}{r\sin\theta}u_\theta u_\phi\right)+ \dfrac{1}{Pm}\left(j_r\,B_\theta-j_\theta\,B_r\right)\, , \end{aligned} \right\rbrace\end{split}\]

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

\[2\tilde{\rho} \vec{e_r}\cdot(\vec{u}\times\vec{e_z})=2\sin\theta\,\tilde{\rho} u_\phi=\dfrac{2}{r}\left(\dfrac{\partial^2 W}{\partial r\partial \phi}-\sin\theta \dfrac{\partial Z}{\partial \theta}\right)\,.\]

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

\[{\cal N}^W = \dfrac{2}{r}\left(\dfrac{\partial^2 W}{\partial r\partial \phi}-\sin\theta \dfrac{\partial Z}{\partial \theta}\right)+{\cal A}_r\,.\]

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:

(32)\[\boxed{ {\cal N}^W_{\ell m} = \dfrac{2}{r}\left[i m \dfrac{\partial W_\ell^m}{\partial r}-(\ell-1)c_\ell^m Z_{\ell-1}^m+(\ell+2)c_{\ell+1}^m Z_{\ell+1}^m\right] +{{\cal A}_r}_\ell^m\, . }\]

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\).

See also

The final calculations of (32) are done in the subroutine get_td.

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\):

\[\begin{split}\begin{aligned} \vec{e_r}\cdot\vec{\nabla}\times\left[(2\tilde{\rho}\vec{u})\times \vec{e_z}\right] & =2\vec{e_r}\cdot\left[(\vec{e_z}\cdot\vec{\nabla})(\tilde{\rho} \vec{u})\right], \\ & = 2\left[\cos\theta\dfrac{\partial (\tilde{\rho} u_r)}{\partial r} -\dfrac{\sin\theta}{r}\dfrac{\partial (\tilde{\rho} u_r)}{\partial \theta}+\dfrac{\tilde{\rho} u_\theta\sin\theta}{r}\right], \\ & = 2\left[-\cos\theta\dfrac{\partial}{\partial r}(\Delta_H W)+ \dfrac{\sin\theta}{r}\dfrac{\partial}{\partial \theta}(\Delta_H W)+\dfrac{\sin\theta}{r^2}\dfrac{\partial^2 W}{\partial r\partial \theta}+ \dfrac{1}{r^2}\dfrac{\partial Z}{\partial \phi}\right]. \end{aligned}\end{split}\]

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

(33)\[\vec{e_r}\cdot\vec{\nabla}\times\left[(2\tilde{\rho}\vec{u})\times \vec{e_z}\right]=\dfrac{2}{r^2}\left(\vartheta_3\,\dfrac{\partial W}{\partial r} -\dfrac{1}{r}\,\vartheta_4\,W+ \dfrac{\partial Z}{\partial \phi} \right)\,.\]

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

\[\dfrac{1}{r\sin\theta}\left[ \dfrac{\partial (\sin\theta{\cal A}_\phi)}{\partial \theta} - \dfrac{\partial {\cal A}_\theta}{ \partial\phi}\right]\,.\]

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

1. 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.

2. 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.

3. 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

(35)\[\begin{split}\boxed{ \begin{aligned} {\cal N}^Z_{\ell m} = & \dfrac{2}{r^2}\left[(\ell-1)(\ell+1)\,c_\ell^m\, \dfrac{\partial W_{\ell-1}^m}{\partial r}+\ell(\ell+2)\,c_{\ell+1}^m\, \dfrac{\partial W_{\ell+1}^m}{\partial r} \right. \\ & \left. -\dfrac{\ell(\ell-1)(\ell+1)}{r}\,c_\ell^m\,W_{\ell-1}^m+ \dfrac{\ell(\ell+1)(\ell+2)}{r}\,c_{\ell+1}^m\,W_{\ell+1}^m+ im\,Z_\ell^m\right] \\ & + (\ell+1)\,c_\ell^m\,{{\cal A}p}_{\ell-1}^m- \ell\,c_{\ell+1}^m\,{{\cal A}p}_{\ell+1}^m -im\,{{\cal A}t}_{\ell}^m\,. \end{aligned} }\end{split}\]

See also

The final calculations of (35) are done in the subroutine get_td.

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\):

\[\begin{split}\begin{aligned} \vec{\nabla}_H\cdot\left[(2\tilde{\rho}\vec{u})\times \vec{e_z}\right] & =2\vec{e_z}\cdot\left[\vec{\nabla}\times(\tilde{\rho} \vec{u})\right] -\left(\dfrac{\partial}{\partial r}+\dfrac{2}{r}\right) \left[\vec{e_r}\cdot(2\tilde{\rho}\vec{u}\times\vec{e_z})\right],\\ & = -2\cos\theta\,\Delta_H Z-2\sin\theta\left[-\dfrac{1}{r\sin\theta} \dfrac{\partial}{\partial\phi}\left( \dfrac{\partial^2}{\partial r^2}+\Delta_H \right) W + \dfrac{1}{r}\dfrac{\partial^2 Z}{\partial r\partial\theta}\right] \\ & \phantom{=\cos\theta} -\left(\dfrac{\partial}{\partial r}+\dfrac{2}{r}\right) \left[2\sin\theta\tilde{\rho}u_\phi\right], \\ & = 2\left[\dfrac{1}{r}\left(\Delta_H+\dfrac{\partial^2}{\partial r^2}\right) \dfrac{\partial W}{\partial \phi}-\cos\theta\Delta_H Z -\dfrac{\sin\theta}{r} \dfrac{\partial^2 Z}{\partial r \partial \theta}\right] \\ & \phantom{=\cos\theta} -\left(\dfrac{\partial}{\partial r}+\dfrac{2}{r}\right) \left[\dfrac{2}{r}\left(\dfrac{\partial^2 W}{\partial r\partial\phi}-\sin\theta \dfrac{\partial Z}{\partial \theta}\right)\right], \\ & = 2\left(\dfrac{\Delta_H}{r}\dfrac{\partial W}{\partial \phi}-\dfrac{1}{r^2} \dfrac{\partial^2 W}{\partial\phi\partial r} -\cos\theta\Delta_H\,Z +\dfrac{\sin\theta}{r^2}\dfrac{\partial Z}{\partial \theta}\right). \end{aligned}\end{split}\]

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

(36)\[\vec{\nabla}_H\cdot\left[(2\tilde{\rho}\vec{u})\times \vec{e_z}\right]=\dfrac{2}{r^2}\left(-\dfrac{L_H}{r}\,\dfrac{\partial W}{\partial \phi} -\dfrac{\partial^2 W}{\partial\phi\partial r}+\vartheta_3\, Z \right)\,.\]

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

\[\dfrac{1}{r\sin\theta}\left[ \dfrac{\partial (\sin\theta{\cal A}_\theta)}{\partial \theta} + \dfrac{\partial {\cal A}_\phi}{ \partial\phi}\right]\,.\]

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\).

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

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

(38)\[\begin{split}\boxed{ \begin{aligned} {\cal N}^P_{\ell m} = & \dfrac{2}{r^2}\left[-im\,\dfrac{\ell(\ell+1)}{r}\,W_\ell^m -im\,\dfrac{\partial W_\ell^m}{\partial r}+(\ell-1)(\ell+1)\,c_\ell^m\, Z_{\ell-1}^m+\ell(\ell+2)\,c_{\ell+1}^m\, Z_{\ell+1}^m \right] \\ & + (\ell+1)\,c_\ell^m\,{{\cal A}t}_{\ell-1}^m- \ell\,c_{\ell+1}^m\,{{\cal A}t}_{\ell+1}^m +im\,{{\cal A}p}_{\ell}^m\,. \end{aligned} }\end{split}\]

See also

The final calculations of (38) are done in the subroutine get_td.

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.

\[{\cal H} = \dfrac{Pr\,Di}{Ra}\dfrac{1}{\tilde{\rho}\tilde{T}}\left(\Phi_\nu+ \dfrac{\lambda}{Pm^2\,E}\,j^2\right)\,.\]

Expanding this term leads to:

(39)\[\begin{split}\begin{aligned} {\cal H}=& \dfrac{Pr\,Di}{Ra}\dfrac{1}{\tilde{\rho}\tilde{T}}\left[ \tilde{\rho}\nu\left\lbrace 2\left(\dfrac{\partial u_r}{ \partial r}\right)^2 +2\left(\dfrac{1}{r}\dfrac{\partial u_\theta}{\partial\theta}+\dfrac{u_r}{r} \right)^2+2\left( \dfrac{1}{r\sin\theta}\dfrac{\partial u_\phi}{\partial\phi} + \dfrac{u_r}{r}+\dfrac{\cos\theta}{r\sin\theta}u_\theta \right)^2\right.\right. \\ & \phantom{\dfrac{Pr\,Di}{Ra}\dfrac{1}{\tilde{\rho}\tilde{T}}} +\left(r\dfrac{\partial}{\partial r}\left(\dfrac{u_\theta}{r} \right)+\dfrac{1}{r}\dfrac{\partial u_r}{\partial\theta}\right)^2+ \left(r\dfrac{\partial}{\partial r}\left(\dfrac{u_\phi}{r}\right)+ \dfrac{1}{r\sin\theta}\dfrac{\partial u_r}{\partial\phi} \right)^2 \\ & \phantom{\dfrac{Pr\,Di}{Ra}\dfrac{1}{\tilde{\rho}\tilde{T}}}\left. + \left(\dfrac{\sin\theta}{r}\dfrac{\partial}{\partial\theta}\left( \dfrac{u_\phi}{\sin\theta}\right)+\dfrac{1}{r\sin\theta} \dfrac{\partial u_\theta}{\partial\phi}\right)^2 -\dfrac{2}{3}\,\left(\dfrac{d\ln\tilde{\rho}}{dr}\,u_r\right)^2 \right\rbrace \\ & \phantom{\dfrac{Pr\,Di}{Ra}\dfrac{1}{\tilde{\rho}\tilde{T}}}\left. + \dfrac{\lambda}{Pm^2\,E}\,\left\lbrace j_r^2+j_\theta^2+j_\phi^2\right\rbrace\right]\,. \end{aligned}\end{split}\]

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

\[-\vec{u}\cdot\vec{\nabla}s = -\dfrac{1}{\tilde{\rho}}\left[ \vec{\nabla}\cdot\left(\tilde{\rho}s\vec{u} \right)- s\underbrace{\vec{\nabla}\cdot\left(\tilde{\rho}\vec{u} \right)}_{=0}\right]\,.\]

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:

\[\mathcal{US}_r = \tilde{\rho}s u_r,\quad \mathcal{US}_\theta = \tilde{\rho}s u_\theta, \quad \mathcal{US}_\phi = \tilde{\rho}s u_\phi,\]

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

\[-\vec{u}\cdot\vec{\nabla}s = -\dfrac{1}{\tilde{\rho}}\left[ \dfrac{1}{r^2}\dfrac{\partial}{\partial r}\left(r^2\,\mathcal{US}_r\right)+ \dfrac{1}{r\sin\theta}\dfrac{\partial}{\partial\theta}\left(\sin\theta\,\mathcal{US}_\theta \right)+\dfrac{1}{r\sin\theta}\dfrac{\partial\,\mathcal{US}_\phi}{\partial\phi}\right]\,.\]

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

1. 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.

2. 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.

3. 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

(41)\[\boxed{ {\cal N}^S_{\ell m} = -\dfrac{1}{\tilde{\rho}}\left[ \dfrac{1}{r^2}\dfrac{\partial\, {\mathcal{US}r}_\ell^m}{\partial r} + (\ell+1)\,c_\ell^m\,{\mathcal{US}t}_{\ell-1}^m- \ell\,c_{\ell+1}^m\,{\mathcal{US}t}_{\ell+1}^m +im\,{\mathcal{US}p}_\ell^m\right]+{\cal H}_\ell^m\,. }\]

See also

The \(\theta\) and \(\phi\) derivatives that enter (41) are done in the subroutine get_td. The radial derivative is computed afterwards at the very beginning of updateS.

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

\[-\vec{u}\cdot\vec{\nabla}\xi = -\dfrac{1}{\tilde{\rho}}\left[ \vec{\nabla}\cdot\left(\tilde{\rho}\xi\vec{u} \right)- \xi\underbrace{\vec{\nabla}\cdot\left(\tilde{\rho}\vec{u} \right)}_{=0}\right]\,.\]

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:

\[\mathcal{UX}_r = \tilde{\rho}\xi u_r,\quad \mathcal{US}_\theta = \tilde{\rho}\xi u_\theta, \quad \mathcal{UX}_\phi = \tilde{\rho}\xi u_\phi,\]

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

\[-\vec{u}\cdot\vec{\nabla}\xi = -\dfrac{1}{\tilde{\rho}}\left[ \dfrac{1}{r^2}\dfrac{\partial}{\partial r}\left(r^2\,\mathcal{UX}_r\right)+ \dfrac{1}{r\sin\theta}\dfrac{\partial}{\partial\theta}\left(\sin\theta\,\mathcal{UX}_\theta \right)+\dfrac{1}{r\sin\theta}\dfrac{\partial\,\mathcal{UX}_\phi}{\partial\phi}\right]\,.\]

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

1. 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.

2. 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.

3. 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

(42)\[\boxed{ {\cal N}^\xi_{\ell m} = -\dfrac{1}{\tilde{\rho}}\left[ \dfrac{1}{r^2}\dfrac{\partial\, {\mathcal{UX}r}_\ell^m}{\partial r} + (\ell+1)\,c_\ell^m\,{\mathcal{UX}t}_{\ell-1}^m- \ell\,c_{\ell+1}^m\,{\mathcal{UX}t}_{\ell+1}^m +im\,{\mathcal{UX}p}_\ell^m\right]\,. }\]

See also

The \(\theta\) and \(\phi\) derivatives that enter (42) are done in the subroutine get_td. The radial derivative is computed afterwards at the very beginning of updateXi.

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

\[{\cal F}_r=u_\theta B_\phi-u_\phi B_\theta,\quad {\cal F}_\theta=u_\phi B_r-u_r B_\phi,\quad {\cal F}_\phi=u_r B_\theta-u_\theta B_r\,.\]

The radial component of the induction term then reads:

\[{\cal N}^g = \vec{e_r}\cdot\left[\vec{\nabla}\times\left(\vec{u}\times\vec{B}\right)\right] =\dfrac{1}{r\sin\theta}\left[\dfrac{\partial\,(\sin\theta {\cal F}_\phi)}{\partial\theta} -\dfrac{\partial {\cal F}_\theta}{\partial \phi}\right]\,.\]

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

1. 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.

2. 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.

3. 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

(43)\[\boxed{ {\cal N}^g_{\ell m} = (\ell+1)\,c_\ell^m\,{\mathcal{F}_\phi}_{\ell-1}^m-\ell\,c_{\ell+1}^m\, {\mathcal{F}_\phi}_{\ell+1}^m -im\,{\mathcal{F}_\theta}_{\ell}^m\,. }\]

See also

The final calculations of (43) are done in the subroutine get_td.

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):

\[\begin{split}\begin{aligned} {\cal N}^h = \vec{e_r}\cdot\left[\vec{\nabla}\times\vec{\nabla}\times\left(\vec{u}\times\vec{B}\right) \right] & =\vec{e_r}\cdot\left[\vec{\nabla}\left(\vec{\nabla}\cdot\vec{\mathcal{F}}\right) -\Delta\vec{\mathcal{F}}\right], \\ & = \dfrac{\partial}{\partial r}\left[\dfrac{1}{r^2} \dfrac{\partial(r^2 {\mathcal{F}}_r)}{\partial r} + \dfrac{1}{r\sin\theta} \dfrac{\partial(\sin\theta\,{\mathcal{F}}_\theta)}{\partial\theta}+\dfrac{1}{r\sin\theta} \dfrac{\partial{\mathcal{F}}_\phi}{\partial\phi} \right] \\ & \phantom{=\ }- \Delta {\mathcal{F}}_r+\dfrac{2}{r^2}\left[{\mathcal{F}}_r +\dfrac{1}{\sin\theta} \dfrac{\partial(\sin\theta\,{\mathcal{F}}_\theta)}{\partial\theta}+ \dfrac{1}{\sin\theta}\dfrac{\partial {\mathcal{F}}_\phi}{\partial \phi}\right], \\ & = \dfrac{1}{r^2}\dfrac{\partial}{\partial r}\left[\dfrac{r}{\sin\theta}\left( \dfrac{\partial(\sin\theta\,{\mathcal{F}}_\theta)}{\partial\theta}+ \dfrac{\partial{\mathcal{F}}_\phi}{\partial\phi} \right)\right]-\Delta_H\,{\mathcal{F}}_r\,. \end{aligned}\end{split}\]

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\):

\[\dfrac{1}{r^2}\dfrac{\partial }{\partial r}\left[r^2\left(\vartheta_1\, {\mathcal{F}t}_\ell^m+\dfrac{\partial\,{\mathcal{F}p}_\ell^m}{\partial \phi}\right)\right] +L_H\, {\mathcal{F}r}_\ell^m\,.\]

We thus finally get

(44)\[\boxed{ {\cal N}^h_{\ell m} =\ell(\ell+1)\,{\mathcal{F}r}_{\ell}^m+ \dfrac{1}{r^2}\dfrac{\partial}{\partial r}\left[r^2\left\lbrace (\ell+1)\,c_\ell^m\,{\mathcal{F}t}_{\ell-1}^m-\ell\,c_{\ell+1}^m\, {\mathcal{F}t}_{\ell+1}^m +im\,{\mathcal{F}p}_{\ell}^m\right\rbrace \right]\,. }\]

See also

The \(\theta\) and \(\phi\) derivatives that enter (44) are computed in the subroutine get_td. The remaining radial derivative is computed afterwards at the very beginning of updateB.

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

(45)\[{\cal C}_n(r) W_{\ell mn} = 0 \;\;\mbox{at}\;\; r=r_i,r_o\;\;.\]

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

\[\dfrac{\partial W}{\partial r}=0\;\;\mbox{and}\;\; Z=0,\]

at the respective boundary. In spectral space these conditions read

(46)\[{\cal C}'_n(r) W_{\ell mn} = 0\;\;\mbox{at}\;\; r=r_i,r_o\,,\]

and

(47)\[{\cal C}_n(r) Z_{\ell mn} = 0\;\;\mbox{at}\;\; r=r_i,r_o\,,\]

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)\):

\[\mathsf{I} \dfrac{\partial\omega}{\partial t}= \Gamma_L+\Gamma_\nu+\Gamma_g\,.\]

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

\[\Gamma_L = \dfrac{1}{E\,Pm}\oint B_r B_\phi\,r\sin\theta\,\mathrm{d}S\,,\]

and

\[\Gamma_\nu = \oint \tilde{\rho} \tilde{\nu} r\dfrac{\partial}{\partial r}\left(\dfrac{u_\phi}{r}\right) r\sin\theta\,\mathrm{d}S\,,\]

and

\[\Gamma_g = -\Gamma\tau(\omega_{ic}-\omega_{ma})\,,\]

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

\[\oint \tilde{\rho} r\sin\theta u_\phi\,\mathrm{d} S=4\sqrt{\dfrac{\pi}{3}}Z_{10}r^2,\]

the viscous torques can be expressed by

\[\Gamma_\nu = \pm4\sqrt{\dfrac{\pi}{3}}\tilde{\nu}r^2\left[\dfrac{\partial Z_{10}}{\partial r}-\left(\dfrac{2}{r}+\beta\right)Z_{10}\right]\,,\]

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

\[\left[{\cal C}''_n(r) -\left(\frac{2}{r}+\dfrac{d\ln\tilde{\rho}}{dr}\right)\,{\cal C}'_n(r) \right] W_{\ell mn} = 0 \;\;\mbox{and}\;\; \left[{\cal C}'_n(r) -\left(\frac{2}{r}+\dfrac{d\ln\tilde{\rho}}{dr}\right)\,{\cal C}_n(r) \right] Z_{\ell mn} = 0\;.\]

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

\[\dfrac{\partial}{\partial r} \dfrac{\vec{u}_H}{r} = 0\]

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

\[\dfrac{\partial}{\partial r} \dfrac{1}{r} \left( \vec{\nabla}_H \dfrac{\partial W}{\partial r}\right) = \vec{\nabla}_H \dfrac{1}{r} \left( \dfrac{\partial^2}{\partial r^2} - \dfrac{2}{r} \dfrac{\partial}{\partial r} \right) W = 0\]

and

\[\dfrac{\partial}{\partial r} \dfrac{1}{r} \nabla\times \vec{e}_r Z = \nabla\times \vec{e}_r \dfrac{1}{r} \left( \dfrac{\partial}{\partial r} - \dfrac{2}{r} \right) Z = 0\]

which can be fulfilled with

\[\left( \dfrac{\partial^2}{\partial r^2} - \dfrac{2}{r} \dfrac{\partial}{\partial r} \right) W = 0\]

and

\[\left( \dfrac{\partial}{\partial r} - \dfrac{2}{r} \right) Z = 0\;.\]

In spectral representation this then reads

\[\left({\cal C}''_n - \dfrac{2}{r}{\cal C}'_n \right) W_{\ell mn} = 0 \;\;\mbox{and}\;\; \left({\cal C}'_n - \frac{2}{r}{\cal C}_n \right) Z_{\ell mn} = 0\;.\]

Thermal boundary conditions

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

\[s=\mbox{const.}\;\;\mbox{or}\;\;T=\mbox{const.}\]

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

\[\dfrac{\partial}{\partial r} s=\mbox{const.}\;\;\mbox{or}\;\;\dfrac{\partial}{\partial r} T=\mbox{const.}\]

In spectral representation for example the respective entropy condition read

\[{\cal C}_n s_{\ell mn}=\mbox{const.}\;\;\mbox{or}\;\;{\cal C'}_n s_{\ell mn}=\mbox{const.}\]

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:

\[\xi=\mbox{const.}\]

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

\[\dfrac{\partial}{\partial r} \xi=\mbox{const.}\]

In spectral representation, this then reads

\[{\cal C}_n \xi_{\ell mn}=\mbox{const.}\;\;\mbox{or}\;\;{\cal C'}_n \xi_{\ell mn}=\mbox{const.}\]

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:

(48)\[{\cal C}_n(r) W_{\ell mn} = 0 \;\;\mbox{at}\;\; r=r_i,r_o\;\;.\]

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

(49)\[{\cal C}_n(r) g_{\ell mn} = 0\;\;,\;\;{\cal C}'_n(r) g_{\ell mn} = 0 \;\;\mbox{and}\;\;{\cal C}_n(r) h_{\ell mn} = 0.\]

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:

\[V_{\ell m}(r) = r_{BC} V_{\ell m}^I \left(\dfrac{r_{BC}}{r}\right)^{\ell+1} + r_{BC} V_{\ell m}^E \left(\dfrac{r}{r_{BC}}\right)^{\ell}.\]

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

\[{\cal C}_n h_{\ell mn} = 0\]

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

\[\nabla_H \dfrac{\partial}{\partial r} g = -\nabla_H V^I\]

for the horizontal components and

\[\dfrac{\nabla_H^2}{r^2} g = \dfrac{\partial}{\partial r} V^I\]

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

\[{\cal C'}_n(r) g_{\ell mn} = V^I_{lm}\;\;\mbox{and}\;\; \dfrac{\ell(\ell+1)}{r^2} {\cal C}_n g_{\ell mn} = - \dfrac{\ell+1}{r} V^I_{lm}.\]

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

\[\left( {\cal C'}_n(r_o) + \frac{\ell}{r_o} {\cal C}_n(r_o) \right) g_{\ell mn} = 0\]

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

\[\left( {\cal C'}_n(r_i) - \frac{\ell+1}{r_i} {\cal C}_n(r_i) \right) g_{\ell mn} = 0.\]

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:

(50)\[\int d\Omega\; {Y_{\ell}^{m}}^\star\;\vec{e_r}\cdot\left[\vec{\nabla}\times \left(\vec{u^I}\times\vec{B^I}\right)\right] = - i\,\omega_I\,m\,\dfrac{\ell(\ell+1)}{r^2}\;g_{\ell m}^I(r)\; ,\]

(51)\[\int d\Omega\; {Y_{\ell}^{m}}^\star\;\vec{e_r}\cdot\left[\vec{\nabla}\times \vec{\nabla}\times\left(\vec{u^I}\times\vec{B^I}\right)\right] = - i\,\omega_I\,m\,\dfrac{\ell(\ell+1)}{r^2} \;h_{\ell m}^I(r)\;.\]

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:

(52)\[x(r)= \dfrac{r}{r_i}\;\;,\;\;-r_i\le r\le r_i\;\;.\]

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:

(53)\[g_{\ell m}^I(r) = \left(\dfrac{r}{r_i}\right)^{\ell+1} \,\sum_{i=0}^{M-1} g_{\ell m\,2i}^I\;{\cal C}_{2i}(r)\;\;.\]

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.