Notes on the Numerical Methods for the PFSS Solver

This page was originally a Latex file that was converted to RST. The author is A. R. Yeates from the Department of Mathematical Sciences, Durham University, UK.

Basic Equations

The sunkit_magex.pfss code solves the for a magnetic field satisfying

\[\nabla\times\boldsymbol{B}=0,\qquad \nabla\cdot\boldsymbol{B}= 0\]

in a spherical shell \(1 \leq r \leq R_{ss}\), given boundary conditions

\[ \begin{align}\begin{aligned}&B_r(\theta,\phi) = g(\theta,\phi) \quad \textrm{on} \quad r=1\\&B_\theta=B_\phi=0 \quad \textrm{on} \quad r=R_{ss}.\end{aligned}\end{align} \]

The function \(g(\theta,\phi)\) is specified by the user.

Numerical Grid

Fig. 1 The numerical grid used in the solver.

The solver uses a rectilinear grid that is equally spaced in \(\rho\), \(s\), \(\phi\), where

\[\rho = \ln(r), \quad s=\cos\theta\]

in terms of spherical coordinates \((r,\theta\,\phi)\). The coordinate scale factors \(|\mathrm{d}\boldsymbol r/\mathrm{d}\rho|\), \(|\mathrm{d}\boldsymbol r/\mathrm{d}s|\), \(|\mathrm{d}\boldsymbol r/\mathrm{d}\phi|\) are

\[h_\rho = r = \mathrm{e}^\rho,\quad h_s = \frac{r}{\sin\theta} = \frac{\mathrm{e}^\rho}{\sqrt{1-s^2}}, \quad h_\phi = r\sin\theta = \mathrm{e}^\rho\sqrt{1-s^2}\]

The grid is illustrated in Fig. 1. Note that the longitudinal cell size goes to zero at the poles; these points are treated specially in the calculation of \(\boldsymbol{B}\). Note also that, since \(s\) is a decreasing function of \(\theta\), the components of a vector \(\boldsymbol{v}\) in \((\rho,s,\phi)\) are \(v_\rho = v_r\) but \(v_s = -v_\theta\).

We define the number of grid cells \(n_\rho\), \(n_s\), \(n_\phi\), with corresponding uniform spacings

\[\Delta_\rho= \frac{\ln(R_{ss})}{n_\rho}, \quad \Delta_s= \frac{2}{n_s}, \quad \Delta_\phi= \frac{2\pi}{n_\phi}\]

Note that the boundaries in \(s\) are at the poles \(s=\pm1\), at which points \(h_\phi\) is not defined. The solution is periodic in the longitudinal (\(\phi\)) direction.

Numerical method

Overall strategy

Rather than writing \(\boldsymbol{B}= \nabla\chi\) and solving \(\nabla^2\chi=0\), we write instead \(\boldsymbol{A}= \nabla\times\big(\psi \,\mathrm{e}_\rho\big)\). Then accounting for the unusual coordinates we get

\[\begin{split}\boldsymbol{B} &= \frac{1}{h_\rho h_sh_\phi} \begin{vmatrix} h_\rho e_\rho & h_\phi e_\phi & h_{s} e_{s} \\ \partial_\rho & \partial_\phi & \partial_s \\ 0 & \frac{h_\phi}{h_s}\partial_s\psi & -\frac{h_s}{h_\phi}\partial_\phi\psi \end{vmatrix}\\ &= -\frac{1}{h_sh_\phi}\left[\partial_s\left(\frac{h_\phi}{h_s}\partial_s\psi\right) + \partial_\phi\left(\frac{h_s}{h_\phi}\partial_\phi\psi\right) \right] e_\rho + \frac{1}{h_\rho h_\phi}\partial_\phi\partial_\rho\psi e_\phi + \frac{1}{h_\rho h_s}\partial_s\partial_\rho\psi e_{s}\\ &= -\Delta_\perp\psi e_\rho + \frac{1}{h_\phi}\partial_\phi\left(\frac{1}{h_\rho}\partial_\rho\psi\right)e_\phi + \frac{1}{h_s}\partial_s\left(\frac{1}{h_\rho}\partial_\rho\psi\right) e_{s}.\end{split}\]

This will take the curl-free form \(\boldsymbol{B}= \nabla\big(\tfrac1{h_\rho}\partial_\rho\psi\big)\) provided that

(1)\[\nabla^2_\perp\psi = -\frac{1}{h_\rho}\partial_\rho\left(\frac{1}{h_\rho}\partial_\rho\psi\right)\]

so our strategy is to solve Equation (1) for \(\psi\), then reconstruct \(\boldsymbol{A}\) and \(\boldsymbol{B}\). The reason for doing it this way is that it allows us to compute \(\boldsymbol{A}\) as well as \(\boldsymbol{B}\) (again, for legacy reasons).

Numerical solution

We follow the method described in [Ballegooijen], except that we modify the finite-difference discretization to suit our particular coordinates.

The discretization is chosen so that we will have \(\nabla\times\boldsymbol{B}=0\) to machine precision on a staggered grid, when the curl is taken using central differences. This property of essentially zero current density is required when using the PFSS solution to, e.g., initialize non-potential simulations. It would not typically be achieved by interpolating a spherical harmonic solution onto the numerical grid. However, we will see that the discrete solution effectively computes discrete approximations of the spherical harmonics, tailored to our particular difference formula.

In the following subsections, we describe the numerical solution in more detail.

Variables

Let the coordinate grid points be defined by

\[ \begin{align}\begin{aligned}&\rho^k = k\Delta_\rho, \qquad k=0,\ldots, n_\rho\\&s^j = j\Delta_s- 1, \qquad j=0,\ldots, n_s\\&\phi^i = i\Delta_\phi, \qquad i=0,\ldots, n_\phi\end{aligned}\end{align} \]

In the code the first two arrays are called rg and sg (that for pg is not required). There are also arrays rc, sc and pc corresponding to the cell centres, i.e. \(\rho^{k+1/2}\), \(s^{j+1/2}\) and \(\phi^{i+1/2}\).

To deal with the curvilinear coordinates, we define the edge lengths

\[ \begin{align}\begin{aligned}&L_\rho^{k+\frac{1}{2},j,i} = \int_{\rho^k}^{\rho^{k+1}} h_\rho\,\mathrm{d}\rho = \,\mathrm{e}^{\rho^{k+1}} - \,\mathrm{e}^{\rho^k}\\&L_s^{k,j+\frac{1}{2},i} = \int_{s^j}^{s^{j+1}} h_s\,\mathrm{d}s = \,\mathrm{e}^{\rho^k}\big(\arcsin(s^{j+1}) - \arcsin(s^j)\big)\\&L_\phi^{k,j,i+\frac{1}{2}} = \int_{\phi^i}^{\phi^{i+1}} h_\phi\,\mathrm{d}\phi = \,\mathrm{e}^{\rho^k}\sigma^j\Delta_\phi.\end{aligned}\end{align} \]

Here we used the fact that \(\Delta_\rho\), \(\Delta_s\) and \(\Delta_\phi\) are constant, and used the shorthand

\[\sigma^j := \sqrt{1 - (s^j)^2}.\]

Similarly we define the areas of the cell faces

\[ \begin{align}\begin{aligned}&S_\rho^{k,j+\frac{1}{2},i+\frac{1}{2}} = \int_{\phi^i}^{\phi^{i+1}}\int_{s^j}^{s^{j+1}} h_s h_\phi\,\mathrm{d}s\mathrm{d}\phi = \,\mathrm{e}^{2\rho^k}\Delta_s\Delta_\phi\\&S_s^{k+\frac{1}{2},j,i+\frac{1}{2}} = \int_{\rho^k}^{\rho^{k+1}}\int_{\phi^i}^{\phi^{i+1}} h_\rho h_\phi\,\mathrm{d}\phi\mathrm{d}\rho = \tfrac12\big(\,\mathrm{e}^{2\rho^{k+1}} - \,\mathrm{e}^{2\rho^{k}}\big)\sigma^j\,\Delta_\phi\\&S_\phi^{k+\frac{1}{2},j+\frac{1}{2},i} = \int_{\rho^k}^{\rho^{k+1}}\int_{s^j}^{s^{j+1}}h_\rho h_s\,\mathrm{d}s\mathrm{d}\rho = \tfrac12\big(\,\mathrm{e}^{2\rho^{k+1}}- \,\mathrm{e}^{2\rho^k}\big)\big(\arcsin(s^{j+1}) - \arcsin(s^j)\big)\end{aligned}\end{align} \]

In the code the face areas are stored in arrays Sbr, Sbs and Sbp (with only the dimensions required).

In the code the magnetic field \(\boldsymbol{B}\) is defined staggered on the face centres, so \(B_\rho^{k,j+\frac{1}{2},i+\frac{1}{2}}\), \(B_s^{k+\frac{1}{2},j,i+\frac{1}{2}}\), \(B_\phi^{k+\frac{1}{2},j+\frac{1}{2},i}\). These variables are called br, bs and bp.

The vector potential is located on the corresponding cell edges, so \(A_\rho^{k+\frac{1}{2},j,i}\), \(A_s^{k,j+\frac{1}{2},i+\frac{1}{2}}\), \(A_\phi^{k,j,i+\frac{1}{2}}\). In fact, these values are never required on their own, only multiplied by the corresponding edge lengths. So the variables alr, als and alp correspond to the products \(L_\rho A_\rho\), \(L_sA_s\) and \(L_\phi A_\phi\), respectively.

Finally, the potential \(\psi\) is located on the \(\rho\)-faces (like \(B_\rho\)), so \(\psi^{k,j+\frac{1}{2},i+\frac{1}{2}}\). It is stored in the variable psi.

Derivatives

Firstly, we have \(\boldsymbol{A}= \nabla\times\big(\psi\,\mathrm{e}_\rho\big)\). Numerically, this is approximated by

(2)\[A_s^{k,j+\frac{1}{2},i} = -\frac{\psi^{k,j+\frac{1}{2},i+\frac{1}{2}} - \psi^{k,j+\frac{1}{2},i-\frac{1}{2}}}{L_\phi^{k,j+\frac{1}{2},i}}, \qquad A_\phi^{k,j,i+\frac{1}{2}} = \frac{\psi^{k,j+\frac{1}{2},i+\frac{1}{2}} - \psi^{k,j-\frac{1}{2},i+\frac{1}{2}}}{L_s^{k,j,i+\frac{1}{2}}}\]

The magnetic field \(\boldsymbol{B}= \nabla\times\boldsymbol{A}\) is then approximated by

(3)\[(S_\rho B_\rho)^{k,j+\frac{1}{2},i+\frac{1}{2}} = (L_s A_s)^{k,j+\frac{1}{2},i+1} - (L_s A_s)^{k,j+\frac{1}{2},i} - (L_\phi A_\phi)^{k,j+1,i+\frac{1}{2}} + (L_\phi A_\phi)^{k,j,i+\frac{1}{2}},\]

(4)\[(S_s B_s)^{k+\frac{1}{2},j,i+\frac{1}{2}} = (L_\phi A_\phi)^{k+1,j,i+\frac{1}{2}} - (L_\phi A_\phi)^{k,j,i+\frac{1}{2}},\]

(5)\[(S_\phi B_\phi)^{k+\frac{1}{2},j+\frac{1}{2},i} = - (L_s A_s)^{k+1,j+\frac{1}{2},i} + (L_s A_s)^{k,j+\frac{1}{2},i}.\]

These formulae correspond to Stokes’ Theorem applied to the cell face. The condition \(\nabla\times\boldsymbol{B}=0\) may be expressed similarly as

(6)\[0 = (L_s B_s)^{k+\frac{1}{2},j,i-\frac{1}{2}} - (L_s B_s)^{k+\frac{1}{2},j,i+\frac{1}{2}} - (L_\phi B_\phi)^{k+\frac{1}{2},j+\frac{1}{2},i} + (L_\phi B_\phi)^{k+\frac{1}{2},j-\frac{1}{2},i}\]

(7)\[0 = (L_\phi B_\phi)^{k+\frac{1}{2},j+\frac{1}{2},i} - (L_\phi B_\phi)^{k-\frac{1}{2},j+\frac{1}{2},i} - (L_\rho B_\rho)^{k,j+\frac{1}{2},i+\frac{1}{2}} + (L_\rho B_\rho)^{k,j+\frac{1}{2},i-\frac{1}{2}}\]

(8)\[0 = (L_\rho B_\rho)^{k,j+\frac{1}{2},i+\frac{1}{2}} - (L_\rho B_\rho)^{k,j-\frac{1}{2},i+\frac{1}{2}} - (L_s B_s)^{k+\frac{1}{2},j,i+\frac{1}{2}} + (L_s B_s)^{k-\frac{1}{2},j,i+\frac{1}{2}}\]

Note that the factors \(L_\rho\), \(L_s\), \(L_\phi\) here are defined normal to the cell faces, not on the edges. But they have the same formulae.

In fact, condition (6) is automatically satisfied. This may be shown using equations (2) to (5), together with our formulae for \(L_s\), \(L_\phi\), \(S_s\) and \(S_\phi\).

Below, we will discretize (1) in such a way that conditions (7) and (8) are also satisfied exactly (up to rounding error).

Boundary conditions for \(\boldsymbol{B}\)

Boundary conditions are needed when br, bs, bp are averaged to the grid points for output. We use a layer of “ghost cells”, whose values are set by the following boundary conditions:

1. In \(\phi\), br and bs are simply periodic.

2. At the outer boundary \(\rho=\log(R_{ss})\), ghost values of bs and bp are set assuming constant gradient in \(\rho\).

3. At the inner boundary, \(\rho=0\), ghost values of bs and bp are set using equations (7) and (8) (effectively assuming zero horizontal current density).

4. At the polar boundaries, the ghost value of br is set to the polemost interior value from the opposite side of the grid. Similarly, bp is set to minus the polemost interior value from the opposite side of the grid. The values of bs at the poles are not meaningful as the cell faces have zero area. However, they are set to the average of the two neighboring interior values at that longitude (with the opposite one being reversed in sign).

Some of these conditions are chosen for compatibility with other codes, and are not necessarily the most straightforward option for a pure PFSS solver.

Discretization of Equation (1)

First, we approximate the two-dimensional Laplacian \(\nabla^2_\perp\psi\) by

\[ \begin{align}\begin{aligned}&\nabla^2_\perp\psi^{k,j+\frac{1}{2},i+\frac{1}{2}} = \frac{1}{S_\rho^{k,j+\frac{1}{2},i+\frac{1}{2}}}\left[ \frac{L_s^{k,j+\frac{1}{2},i+1}}{L_\phi^{k,j+\frac{1}{2},i+1}}\big(\psi^{k,j+\frac{1}{2},i+\frac{3}{2}} - \psi^{k,j+\frac{1}{2},i+\frac{1}{2}}\big) - \frac{L_s^{k,j+\frac{1}{2},i}}{L_\phi^{k,j+\frac{1}{2},i}}\big(\psi^{k,j+\frac{1}{2},i+\frac{1}{2}} - \psi^{k,j+\frac{1}{2},i-\frac{1}{2}}\big) \right.\\&\left. + \frac{L_\phi^{k,j+1,i+\frac{1}{2}}}{L_s^{k,j+1,i+\frac{1}{2}}}\big(\psi^{k,j+\frac{3}{2},i+\frac{1}{2}} - \psi^{k,j+\frac{1}{2},i+\frac{1}{2}}\big) - \frac{L_\phi^{k,j,i+\frac{1}{2}}}{L_s^{k,j,i+\frac{1}{2}}}\big(\psi^{k,j+\frac{1}{2},i+\frac{1}{2}} - \psi^{k,j-\frac{1}{2},i+\frac{1}{2}}\big) \right]\end{aligned}\end{align} \]

As shorthand we define the quantities

\[U^{j+\frac{1}{2}} = \left(\frac{L_s}{\Delta_s\Delta_\phi L_\phi}\right)^{j+\frac{1}{2}}, \qquad V^j = \left(\frac{L_\phi}{\Delta_s\Delta_\phi L_s}\right)^j,\]

noting that these both depend on \(j\) only. In the code these are called Uc and Vg. Then we can write our discretization as

\[ \begin{align}\begin{aligned}\nabla^2_\perp\psi^{k,j+\frac{1}{2},i+\frac{1}{2}} = \frac{1}{\,\mathrm{e}^{2\rho^k}}\Big[U^{j+\frac{1}{2}}\big(\psi^{k,j+\frac{1}{2},i+\frac{3}{2}} - \psi^{k,j+\frac{1}{2},i-\frac{1}{2}} \big) + V^{j+1}\psi^{k,j+\frac{3}{2},i+\frac{1}{2}} + V^{j}\psi^{k,j-\frac{1}{2},i+\frac{1}{2}}\\- \Big(2U^{j+\frac{1}{2}} + V^{j+1} + V^{j}\Big)\psi^{k,j+\frac{1}{2},i+\frac{1}{2}}\Big].\end{aligned}\end{align} \]

This is the left-hand side of (1).

To discretize the right-hand side of (1), we use the approximation

(9)\[-\frac{1}{h_\rho}\partial_\rho\left(\frac{1}{h_\rho}\partial_\rho\psi\right)^{k,j+\frac{1}{2},i+\frac{1}{2}} = -\frac{c(\Delta_\rho)}{L_\rho^{k,j+\frac{1}{2},i+\frac{1}{2}}}\left(\frac{\psi^{k+1,j+\frac{1}{2},i+\frac{1}{2}} - \psi^{k,j+\frac{1}{2},i+\frac{1}{2}}}{L_\rho^{k+\frac{1}{2},j+\frac{1}{2},i+\frac{1}{2}}} - \frac{\psi^{k,j+\frac{1}{2},i+\frac{1}{2}} - \psi^{k-1,j+\frac{1}{2},i+\frac{1}{2}}}{L_\rho^{k-\frac{1}{2},j+\frac{1}{2},i+\frac{1}{2}}} \right),\]

where

\[c(\Delta_\rho) = \frac{2\,\mathrm{e}^{\Delta_\rho/2}}{\,\mathrm{e}^{\Delta_\rho} + 1} = \mathrm{sech}\left(\frac{\Delta_\rho}{2}\right).\]

Combining this with (9), we arrive at

(10)\[ \begin{align}\begin{aligned}U^{j+\frac{1}{2}}\big(\psi^{k,j+\frac{1}{2},i+\frac{3}{2}} - \psi^{k,j+\frac{1}{2},i-\frac{1}{2}} \big) + V^{j+1}\psi^{k,j+\frac{3}{2},i+\frac{1}{2}} + V^{j}\psi^{k,j-\frac{1}{2},i+\frac{1}{2}} - \Big(2U^{j+\frac{1}{2}} + V^{j+1} + V^{j}\Big)\psi^{k,j+\frac{1}{2},i+\frac{1}{2}}\\= -\frac{c(\Delta_\rho)\,\mathrm{e}^{2\rho^k}}{L_\rho^{k,j+\frac{1}{2},i+\frac{1}{2}}}\left(\frac{\psi^{k+1,j+\frac{1}{2},i+\frac{1}{2}} - \psi^{k,j+\frac{1}{2},i+\frac{1}{2}}}{L_\rho^{k+\frac{1}{2},j+\frac{1}{2},i+\frac{1}{2}}} - \frac{\psi^{k,j+\frac{1}{2},i+\frac{1}{2}} - \psi^{k-1,j+\frac{1}{2},i+\frac{1}{2}}}{L_\rho^{k-\frac{1}{2},j+\frac{1}{2},i+\frac{1}{2}}} \right).\end{aligned}\end{align} \]

The reader may verify algebraically that conditions (7) and (8) follow if this finite-difference equation is satisfied.

Method of solution

Equation (10), together with the appropriate boundary conditions, yields a large (but sparse) system of \(n_\rho n_sn_\phi\times n_\rho n_sn_\phi\) linear equations to solve. Fortunately, following [Ballegooijen], we can reduce this to a series of symmetric tridiagonal eigenvalue problems, if we look for eigenfunctions of the form

(11)\[\psi^{k,j+\frac{1}{2},i+\frac{1}{2}} = f^kQ_{lm}^{j+\frac{1}{2}}\,\,\mathrm{e}^{2\pi I mi/n_\phi}.\]

Here the \(k\) in \(f^k\) is a power, not an index, and \(I\) is the square root of \(-1\) (since we already used \(i\) and \(j\) for indices). This reduction will enable very efficient solution of the linear system.

Substituting (11) in Equation (10) gives

\[ \begin{align}\begin{aligned}-V^{j}Q^{j-\frac{1}{2}}_{lm} + \left(V^{j} + V^{j+1}+ 4U^{j+\frac{1}{2}}\sin^2\left(\tfrac{\pi m}{n_\phi}\right) \right)Q^{j+\frac{1}{2}}_{lm} - V^{j+1}Q^{j+\frac{3}{2}}_{lm}\\= \frac{c(\Delta_\rho)\mathrm{e}^{2\rho^k}}{L_\rho^{k,j+\frac{1}{2},i+\frac{1}{2}}}\left(\frac{f - 1}{L_\rho^{k+\frac{1}{2},j+\frac{1}{2},i+\frac{1}{2}}} - \frac{1 - f^{-1}}{L_\rho^{k-\frac{1}{2},j+\frac{1}{2},i+\frac{1}{2}}} \right)Q_{lm}^{j+\frac{1}{2}}.\end{aligned}\end{align} \]

The right-hand side can be simplified since the dependence on \(k\) cancels out. This leaves the tridiagonal eigenvalue problem

\[-V^{j}Q^{j-\frac{1}{2}}_{lm} + \left(V^{j} + V^{j+1}+ 4U^{j+\frac{1}{2}}\sin^2\left(\tfrac{\pi m}{n_\phi}\right) \right)Q^{j+\frac{1}{2}}_{lm} - V^{j+1}Q^{j+\frac{3}{2}}_{lm} = \lambda_{lm}Q_{lm}^{j+\frac{1}{2}},\]

where \(f\) is determined from the eigenvalues \(\lambda_{lm}\) by solving the quadratic relation

\[\lambda_{lm} = \frac{c(\Delta_\rho)}{\mathrm{e}^{\Delta_\rho/2} - \mathrm{e}^{-\Delta_\rho/2}} \left(\frac{1-f^{-1}}{1-\mathrm{e}^{-\Delta_\rho}} - \frac{f-1}{\mathrm{e}^{\Delta_\rho} - 1}\right).\]

This may be rearranged into the form

\[f^2 - \left[1 + \mathrm{e}^{\Delta_\rho} + \mathrm{sech}\left(\frac{\Delta_\rho}{2}\right)\lambda_{lm}(\mathrm{e}^{\Delta_\rho}-1)(\mathrm{e}^{\Delta_\rho/2} - \mathrm{e}^{-\Delta_\rho/2}) \right]f + \mathrm{e}^{\Delta_\rho} = 0,\]

with two solutions for each \(l\), \(m\) given by

\[f_{lm}^+, f_{lm}^- = F_{lm} \pm \sqrt{F_{lm}^2 - \mathrm{e}^{\Delta_\rho}}, \quad \textrm{where} \quad F_{lm} = \tfrac12\Big[1 + \mathrm{e}^{\Delta_\rho} + \lambda_{lm}(\mathrm{e}^{\Delta_\rho}-1)\sinh(\Delta_\rho) \Big].\]

In the code, the eigenvalues are called lam and the corresponding matrix of eigenvectors is Q. The solutions \(f_{lm}^+\) and \(f_{lm}^-\) are called ffp and ffm respectively.

The solution may then be written as a sum of these two sets of radial eigenfunctions:

\[\psi^{k,j+\frac{1}{2},i+\frac{1}{2}} = \sum_{l=0}^{n_s-1}\sum_{m=0}^{n_\phi-1}\Big[c_{lm}(f_{lm}^+)^k + d_{lm}(f_{lm}^-)^k) \Big] Q_{lm}^{j+\frac{1}{2}}\,\mathrm{e}^{2\pi I mi/n_\phi}. \label{eqn:psisum}\]

The coefficients \(c_{lm}\) and \(d_{lm}\) are then determined by the radial boundary conditions:

1. At the inner boundary \(\rho = 0\), where \(k=0\), we want \(B_\rho = -\nabla^2_\perp\psi\) to match our given distribution \(g^{j+\frac{1}{2},i+\frac{1}{2}}\). We have

   \[B_\rho^{k,j+\frac{1}{2},i+\frac{1}{2}} = \sum_{l=0}^{n_s-1}\sum_{m=0}^{n_\phi-1}\frac{\lambda_{lm}}{\mathrm{e}^{2\rho^k}}\Big[c_{lm}(f_{lm}^+)^k + d_{lm}(f_{lm}^-)^k) \Big] Q_{lm}^{j+\frac{1}{2}}\mathrm{e}^{2\pi I mi/n_\phi},\]

   so

   \[g^{j+\frac{1}{2},i+\frac{1}{2}} = \sum_{l=0}^{n_s-1}\sum_{m=0}^{n_\phi-1}\frac{\lambda_{lm}}{\mathrm{e}^{2\rho^0}}\Big[c_{lm} + d_{lm}\Big] Q_{lm}^{j+\frac{1}{2}}\mathrm{e}^{2\pi I mi/n_\phi}.\]

   We take the discrete Fourier transform of \(g^{j+\frac{1}{2},i+\frac{1}{2}}\) in \(i\), so that (noting \(\,\mathrm{e}^{\rho^0}=1\)),

   \[\sum_{l=0}^{n_s-1}\sum_{m=0}^{n_\phi-1}\lambda_{lm} \Big[c_{lm} + d_{lm}\Big] Q_{lm}^{j+\frac{1}{2}}\,\mathrm{e}^{2\pi I mi/n_\phi} = \sum_{m=0}^{n_\phi-1}b_m^{j+\frac{1}{2}}\,\mathrm{e}^{2\pi I mi/n_\phi}.\]

   Then the orthonormality of \(Q_{lm}^{j+\frac{1}{2}}\) for different \(l\) allows us to determine

   (12)\[ c_{lm} + d_{lm} = \frac{1}{\lambda_{lm}}\sum_{j=0}^{n_s-1}b_m^{j+\frac{1}{2}}Q_{lm}^{j+\frac{1}{2}}.\]

2. At the source (outer) surface \(\rho=\ln(R_{ss})\), where \(k=n_\rho\), there are two options.

   • Radial field - We impose \(\partial_\rho\psi = 0\), in the form \(\psi^{n_\rho,j+\frac{1}{2},i+\frac{1}{2}}=\psi^{n_\rho-1,j+\frac{1}{2},i+\frac{1}{2}}\), which gives

   (13)\[ \frac{d_{lm}}{c_{lm}} = \frac{(f_{lm}^+)^{n_\rho} - (f_{lm}^+)^{n_\rho-1}}{(f_{lm}^-)^{n_\rho-1} - (f_{lm}^-)^{n_\rho}}.\]

   (Numerically it is better to compute this ratio the other way up, to prevent overflow.)

   • Imposed \(B_r\) - In this case the boundary condition is treated similarly to the inner boundary. We require

   \[\hat{g}^{j+\frac{1}{2},i+\frac{1}{2}} = \sum_{l=0}^{n_s-1}\sum_{m=0}^{n_\phi-1}\frac{\lambda_{lm}}{\mathrm{e}^{2\rho^{n_\rho}}}\Big[c_{lm}(f_{lm}^+)^{n_\rho} + d_{lm}(f_{lm}^-)^{n_\rho}\Big] Q_{lm}^{j+\frac{1}{2}}\mathrm{e}^{2\pi I mi/n_\phi},\]

   so we may again take the discrete Fourier transform to end up with

   (14)\[ c_{lm}(f_{lm}^+)^{n_\rho} + d_{lm}(f_{lm}^-)^{n_\rho} = \frac{\mathrm{e}^{2\rho^{n_\rho}}}{\lambda_{lm}}\sum_{j=0}^{n_s-1}\hat{b}_m^{j+\frac{1}{2}}Q_{lm}^{j+\frac{1}{2}},\]

Solving (13) simultaneously with either (12) or (14) gives \(c_{lm}\) and \(d_{lm}\). These are called clm and dlm in the code.

Remark: as we increase the grid resolution, the eigenfunctions \(Q_{lm}^{j+\frac{1}{2}}\), which are functions of \(\theta\), should converge to the corresponding associated Legendre polynomials \(P_l^m(\cos\theta)\), up to normalization. This is illustrated in Fig. 2.

Fig. 2 Comparison of \(P_l^m(\cos\theta)\) (coloured lines) with the discrete eigenfunctions \(Q_{lm}^{j+\frac{1}{2}}\) (black dots), for \(m=6\) and \(l=0,\ldots,4\), at resolution \(n_s=60\) and \(n_\phi=120\).

References

[Ballegooijen] (1,2)

Mean Field Model for the Formation of Filament Channels on the Sun, ApJ, 539, 983-994, 2000.