OBTAINING POTENTIAL FIELD SOLUTIONS WITH SPHERICAL HARMONICS AND FINITE DIFFERENCES

Magnetograms provide the radial magnetic field on the visible surface of the Sun. The actual measurement is for the line-of-sight (LOS) component of the magnetic field, which is then transformed into the radial component assuming an (approximately) radial field near the solar surface. As the Sun rotates, the individual magnetograms can be combined into a synoptic magnetogram that covers the whole spherical surface. Synoptic magnetograms are provided by many observatories, including Wilcox Solar Observatory (WSO), the Michelson Doppler Imager (MDI) instrument on the Solar and Heliospheric Observatory (SOHO), the Global Oscillation Network Group (GONG), Solar Dynamics Observatory (SDO), and the Synoptic Optical Long-term Investigations of the Sun (SOLIS) observatory. Today's magnetograms contain hundreds to thousands of pixels along each coordinate direction. These magnetograms can be used to extrapolate the magnetic field into the solar corona.

The simplest model (Schatten et al. 1969) assumes a current-free, in other words potential, magnetic field that matches the radial field of the magnetogram on the surface, while it satisfies a simple boundary condition at the outer boundary at some radial distance R. The outer boundary condition is usually taken at R = 2.5 R_s (solar radii), and a purely radial field is assumed at this "source surface." Mathematically the problem is the following: given the magnetogram data that define the radial component of the magnetic field as M(θ, ϕ) at r = 1 R_s, find the scalar potential Φ so that

$$\begin{equation} \nabla \cdot (\nabla \Phi) = 0 \end{equation} \tag{ 1 }$$

$$\begin{equation} \left.\frac{\partial \Phi }{\partial r}\right|_{r=1} = M(\theta,\phi) \end{equation} \tag{ 2 }$$

$$\begin{equation} \left.\Phi \right|_{r=R} = 0. \end{equation} \tag{ 3 }$$

Here, θ ∈ [0, π] and ϕ ∈ [0, 2π] are the co-latitude and longitude coordinates, respectively. Once the solution is found, the potential field solution is obtained as

$$\begin{equation} \mathbf {B}= \nabla \Phi, \end{equation} \tag{ 4 }$$

and it will trivially satisfy both the divergence-free and the current-free properties

$$\begin{eqnarray} \nabla \cdot \mathbf {B}= \nabla \cdot (\nabla \Phi) &=& 0 \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} \nabla \times \mathbf {B}= \nabla \times (\nabla \Phi) &=& 0. \end{eqnarray} \tag{ 6 }$$

We note that the current is only zero inside the domain. If the solution is continued out to r > R with a purely radial magnetic field B_r(r > R) = (R/r)²B_r(R), there will be a finite current at r ⩾ R; on the other hand, the divergence will be zero for all r > 1.

The potential field solution is often obtained with a spherical harmonics expansion (Altschuler et al. 1977). Here we briefly summarize the procedure in its simplest possible form. The base functions φ_nm are the spherical harmonic functions Y_nm multiplied with an appropriate linear combination of the corresponding radial functions rⁿ and r^{−n − 1} so that the boundary condition φ_nm(R, θ, ϕ) = 0 is satisfied:

$$\begin{equation} \varphi _{nm}(r,\theta,\phi) = (r^n - R^{2n+1}r^{-n-1}) Y_{nm}(\theta,\phi). \end{equation} \tag{ 7 }$$

The indexes n ⩾ 0 and m (|m| ⩽ n) are the integer degree and order of the spherical harmonic function, respectively. The functions φ_nm are solutions of the Laplace equation (1), satisfy the boundary condition at r = R, and they form an orthogonal base in the θ, ϕ coordinates. The magnetic potential solution can be approximated as a linear combination of the base functions

$$\begin{equation} \Phi (r,\theta,\phi) = \sum _{n=1}^N \sum _{m=-n}^{n} f_{nm} \varphi _{nm}(r,\theta,\phi), \end{equation} \tag{ 8 }$$

where N is the highest degree considered in the expansion and the n = m = 0 harmonic is not included as it corresponds to the monopole term. The coefficients f_nm can be determined by taking the radial derivative of Equation (8) and equating it with the magnetogram radial field at r = 1:

$$\begin{eqnarray} M(\theta,\phi) &=& \sum _{n=1}^N \sum _{m=-n}^{n} f_{nm} \left.\frac{\partial \varphi _{nm}}{\partial r}\right|_{r=1} \nonumber\\ & =& \sum _{n=1}^N [n + (n+1)R^{2n+1}]\sum _{m=-n}^{n} f_{nm} Y_{nm}(\theta,\phi). \end{eqnarray} \tag{ 9 }$$

Exploiting the orthogonality of the base functions, we can take the inner product with Y_nm to determine f_nm as

$$\begin{eqnarray} f_{nm} &=& \frac{1}{4\pi [n + (n+1)R^{2n+1}]} \int _0^\pi d\theta \sin \theta\nonumber\\ && \times \int _0^{2\pi } d\phi M(\theta,\phi) Y_{nm}(\theta,\phi), \end{eqnarray} \tag{ 10 }$$

where the 1/(4π) coefficient results from the normalization of the spherical harmonics. An alternative approach to obtaining the harmonic coefficients f_nm is to employ a (least-squares) fitting procedure in Equation (9). This is much more expensive than evaluating the integral in Equation (10), but it can be more robust if the magnetogram does not cover (well) the whole surface of the Sun.

Using the spherical harmonic coefficients, the potential can be determined on an arbitrary grid using Equation (8) and the magnetic field can be obtained with finite differences. Alternatively, one can calculate the gradient of the base functions analytically and obtain the magnetic field as

$$\begin{equation} \mathbf {B}(r,\theta,\phi) = \sum _{n=1}^N \sum _{m=-n}^{n} f_{nm} \nabla \varphi _{nm}(r,\theta,\phi), \end{equation} \tag{ 11 }$$

for 1 ⩽ r ⩽ R. Spherical harmonics provide a computationally efficient and very elegant way of solving the Laplace equation on a spherical shell. However, one needs to be cautious of how the integral in Equation (10) is evaluated, especially when a large number of harmonics are used in the series expansion.

We will use the GONG synoptic magnetogram for Carrington Rotation 2077 (CR2077, from 2008 November 20 to December 17) as an example to demonstrate the problem. The magnetogram contains the radial field on a 180 × 360 latitude–longitude grid on the solar surface. The grid spacing is uniform in cos θ (or sine of the latitude) and in longitude ϕ. Figure 1 shows the radial field.

Zoom In Zoom Out Reset image size

Section 2 discusses the naive and more sophisticated ways of obtaining the potential field solution with spherical harmonics. Section 3 describes an alternative approach using an iterative finite difference. The various methods are compared in Section 4, where we also demonstrate the ringing effect that can arise in the spherical harmonics solution, and we draw our conclusions.

To turn the analytic prescription given in the introduction into a scheme that works with real magnetograms, one has to pick the maximum degree N and evaluate the integrals in Equation (10) for each pair of n and m up to the highest order. The resulting f_nm coefficients can be used to construct the three-dimensional (3D) potential magnetic field solution at any given point using Equation (11).

2.1. Naive Spherical Harmonics Approach

The simplest approximation to Equation (10) is a discrete integral using the original magnetogram data:

$$\begin{eqnarray} f_{nm} &=& \frac{1}{4\pi [n + (n+1)R^{2n+1}]} \nonumber\\ && \sum _{i=1}^{N_{\theta }}\sum _{j=1}^{N_{\phi }} (\Delta \cos \theta)_i (\Delta \phi)_j M_{i,j} Y_{nm}(\theta _i,\phi _j), \end{eqnarray} \tag{ 12 }$$

where M_{i, j} is the radial field in a pixel of the N_θ by N_ϕ sized magnetogram. The pixel is centered at the (θ_i, ϕ_j) coordinates, and the area of the pixel is given by (Δcos θ)_i(Δϕ)_j.

Unfortunately, the uniform cos θ mesh used by most of the magnetograms is not optimal to evaluate the integral in Equation (10). In fact this procedure will only work with maximum order N that is much less than N_θ. Figure 2 shows the P_{90, 0} associated Legendre polynomial discretized in different ways. The red curve shows the discretization on 180 grid points that are uniform in cos θ. Clearly, the red curve is a very poor representation near the poles, where cos θ = ±1. This is important, because the amplitude of the Legendre polynomial is actually largest near the poles. This means that the orthonormal property is not satisfied in the discrete sense, and the coefficients obtained with Equation (12) are very inaccurate. The Legendre polynomial can be represented much better on a uniform θ grid (shown by the green curve in Figure 2), as we will discuss below.

Zoom In Zoom Out Reset image size

A clear signal of this problem is that the amplitudes of the higher-order spherical harmonics are not getting smaller with increasing indexes n and m, i.e., the harmonic expansion is not converging. The black line in Figure 3 shows the amplitudes S_n = ∑_mf²_nm which oscillates wildly for n > 60 for this 360 × 180 magnetogram. The oscillations are almost exclusively due to the f_n0 coefficients; the m ≠ 0 coefficients are well behaved. This plot can be directly compared with Figure 15 in Altschuler et al. (1977), where the power spectrum is more-or-less exponentially decaying. We believe that the reason is that these authors used a least-square fitting to the LOS magnetic field instead of calculating the spherical harmonics from the radial field as shown above. While the two methods are identical analytically (assuming that the LOS and radial fields correspond to the same solution), the use of least-squares fitting mitigates the lack of orthogonality among the discretized Legendre functions, while the naive approach described above relies heavily on the orthogonality property.

Zoom In Zoom Out Reset image size

Given the non-converging series expansion, the resulting potential field will be very inaccurate in the polar regions and will have essentially random values depending on the number of spherical harmonics used. This is demonstrated by Figure 4 that shows the radial magnetic field reconstructed with various numbers of harmonics using the original magnetogram grid. One would expect the radial component of the potential field to reproduce the magnetogram shown in Figure 1. Instead, we find that the solution deviates strongly in the polar region if the harmonics expansion is continued above N = 60. For N = 60 (or lower) the solution looks reasonable, but is strongly smoothed due to the insufficient number of harmonics. This is most obvious around the active regions in the top panel of Figure 4.

Zoom In Zoom Out Reset image size

We note that these numerical errors are not related or comparable to the observational uncertainties of the magnetograms, which are usually also quite large in the polar regions. The observational uncertainties are essentially unavoidable but are within some well understood range. On the other hand, these numerical artifacts are definitely avoidable while the errors are essentially unbounded if one uses too many harmonics.

2.2. Spherical Harmonics with Remeshed Magnetogram

One can get much more accurate results if the magnetogram is remeshed to a grid that is uniform in the co-latitude θ, has an odd number of nodes, and contains the two poles θ = 0 and θ = π. In fact this is the standard grid used in spherical harmonics transform (e.g., Suda & Takami 2001) and it is often referred to as using the Chebyshev nodes, since the uniform θ grid points correspond to the Chebyshev nodes in the original cos θ coordinate which is the argument of the Legendre polynomials. Figure 2 shows that the Legendre polynomial is much better represented on the uniform θ grid than on the uniform cos θ grid. Remeshing the magnetogram introduces some new adjustable parameters into the procedure: the number of grid cells N'_θ on the new mesh, and the interpolation procedure.

If the remeshing is done with the same number of grid points as in the original magnetogram grid, the latitudinal cell size at the equator will be a factor of π/2 larger than that in the uniform cos θ grid. On the other hand, the uniform θ grid will contain many more points than the original in the polar regions, so the interpolation procedure may create some unwanted artifacts. To maintain the resolution of the original data around the equator, we set N'_θ to (π/2)N_θ rounded to an odd integer. For the remeshing we chose a simple linear interpolation procedure, which works satisfactorily, but one could certainly use higher-order interpolation procedures such as splines. Before doing the interpolation, we add extra grid cells corresponding to the north and south poles of the magnetogram grid, and the values at these two extra cells are set as the average of the pixels around the poles:

$$\begin{eqnarray} M_{0} &=& \frac{1}{N_{\phi }}\sum _{j=1}^{N_{\phi }} M_{1,j} \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} M_{N_{\theta }+1} &=& \frac{1}{N_{\phi }}\sum _{j=1}^{N_{\phi }} M_{N_{\theta },j}. \end{eqnarray} \tag{ 14 }$$

The co-latitude coordinates of the uniform θ mesh are

$$\begin{equation} \theta ^{\prime }_{i^{\prime }}=\pi \frac{i^{\prime }-1}{N^{\prime }_{\theta }-1} \end{equation} \tag{ 15 }$$

for i' = 1...N'_θ. We use simple linear interpolation from the extended magnetogram mesh to the uniform θ mesh:

$$\begin{equation} M^{\prime }_{i^{\prime },j} = \alpha M_{i,j} + (1-\alpha) M_{i^{\prime }+1,j}, \end{equation} \tag{ 16 }$$

where the index i is determined so that $\theta _i \le \theta ^{\prime }_{i^{\prime }} \le \theta _{i+1}$ and

$$\begin{equation} \alpha = \frac{\theta _{i+1} - \theta ^{\prime }_{i^{\prime }}}{\theta _{i+1} - \theta _i}. \end{equation} \tag{ 17 }$$

Finally, the spherical harmonics coefficients are determined with the integral approximated as

$$\begin{eqnarray} f_{nm} &=& \frac{1}{4\pi [n + (n+1)R^{2n+1}]}\nonumber\\ && \times \sum _{i=1}^{N^{\prime }_{\theta }} \sum _{j=1}^{N_{\phi }} \epsilon _i w_i(\Delta \phi)_j M^{\prime }_{i,j} Y_{nm}(\theta ^{\prime }_i,\phi _j), \end{eqnarray} \tag{ 18 }$$

where $\epsilon _1=\epsilon _{N^{\prime }_{\theta }}=1/2$ and _i = 1 for all other indexes. The w_i coefficients are the Clenshaw–Curtis weights (Clenshaw & Curtis 1960; Potts et al. 1998) defined as

$$\begin{equation} w_i = \frac{1}{H} \sum _{k=0}^{H} \epsilon ^{\prime }_k \frac{-2}{4k^2-1}\cos (2\pi k\theta ^{\prime }_i), \end{equation} \tag{ 19 }$$

where H = (N'_θ − 1)/2 and '₀ = '_H = 1/2 and '_k = 1 for all other indexes. We note that for N'_θ order of 10 or more

$$\begin{equation} w_i \approx \frac{{|}\cos \theta ^{\prime }_{i_+} - \cos \theta ^{\prime }_{i_-}{|}}{2}, \end{equation} \tag{ 20 }$$

where i₊ = min (i + 1, N'_θ) and i₋ = max (i − 1, 1) are the indexes of the neighboring cells, or the cell itself at the poles.

Using the proper grid allows us to use a larger number of spherical harmonics. N is limited only by the N ⩽ 2N'_θ/3 and N ⩽ N_ϕ/3 alias-free conditions (Suda & Takami 2001). For our example 180 × 360 magnetogram, we remesh it to a 283 × 360 uniform θ grid, and we can obtain accurate solutions up to N = 120 deg spherical harmonics. Figure 5 shows the potential field solution obtained with the remeshed magnetogram grid with N = 120. Compared with the naive method, the solution is much more reasonable in the polar regions. There is some smoothing when compared to the magnetogram shown in Figure 1, most obvious near the active regions.

Zoom In Zoom Out Reset image size

While the remeshing is definitely a big improvement over using the original magnetogram, it would be nice to be able to use the original magnetogram data without remeshing and interpolation. The following section shows that this can be easily achieved with a finite difference solver.

The Laplace equation (1) with the boundary conditions (2) and (3) can be solved quite easily with an iterative finite difference method. The advantage of finite differences compared with spherical harmonics is that the boundary data given by the magnetogram directly affect the solution only locally, while the spherical harmonics are global functions and their amplitudes depend on all of the magnetogram data. If the magnetogram contains large discontinuities, we expect the finite difference scheme to be better behaved.

The finite difference method has advantages if the solution is to be used in a finite difference code on the same grid, because one can guarantee zero divergence and curl for the magnetic field in the finite difference sense. The solution obtained with the spherical harmonics has zero divergence and curl analytically, but not on the finite difference grid, which may severely underresolve the high-order spherical harmonic functions in some regions (see Figure 2).

The finite difference method was applied to the solar potential magnetic field problem as early as 1976 (Adams & Pneuman 1976), but the method was limited by the computational resources available at the time. Solving a 3D Laplace equation on today's computers is an almost trivial problem. We implemented the new Finite Difference Iterative Potential-field Solver (FDIPS) code in Fortran 90. The serial version does not require any external libraries, while the parallel version uses the Message Passing Interface (MPI) library for communication. FDIPS can solve the Laplace equation on a 150 × 180 × 360 spherical grid to high accuracy on a single processor in less than an hour. The parallel code can solve the same problem in less than 5 minutes on 16 processors.

We briefly describe the algorithm in FDIPS. We use a staggered spherical grid: the magnetic field is discretized on cell faces while the potential is discretized at the cell centers. We use one layer of ghost cells to apply the boundary conditions so the cell centers are located at r_i, θ_j, ϕ_k with i = 0, ..., N_r + 1, j = 0, ..., N_θ + 1, and k = 0, ..., N_ϕ + 1. The θ_j and ϕ_k coordinates of the real cells are given by the magnetogram, while the ghost cell coordinates are given by θ₀ = −θ₁, $\theta _{N_{\theta }+1} = 2\pi - \theta _{N_{\theta }}$, $\phi _0 = \phi _{N_{\phi }}$, and $\phi _{N_{\phi }+1} = \phi _1$. We allow for a non-uniform radial grid extending from r = 1 to R, but for the sake of simplicity in this paper a uniform radial grid is used with r_i = 1 + (i − 1/2)Δr with Δr = (R − r)/N_r.

The radial magnetic field components are located at the radial cell interfaces at (r_{i + 1/2}, θ_j, ϕ_k), where r_{i + 1/2} = (r_i + r_{i + 1})/2 for i = 0, ..., N_r, j = 1, ..., N_θ, and k = 1, ..., N_ϕ. Similarly, the latitudinal components are at r_i, θ_{j + 1/2}, ϕ_k with cos θ_{j + 1/2} = (cos θ_j + cos θ_{j + 1})/2, and j = 0, ..., N_θ. Note that the interface is taken half-way in the cos θ coordinate and not in θ, because this makes the cells equal in area when the magnetogram is given on a uniform cos θ grid. Finally, the longitudinal field components are located at (r_i, θ_j, ϕ_{k + 1/2}), where ϕ_{k + 1/2} = (ϕ_k + ϕ_{k + 1})/2 for k = 0, ..., N_ϕ.

The staggered discretization keeps the stencil of the Laplace operator compact and makes the boundary conditions relatively simple. The magnetic field is obtained as a discrete gradient of Φ:

$$\begin{eqnarray} B_{r,i+1/2,j,k} &=& \frac{\Phi _{i+1,j,k}-\Phi _{i,j,k}}{\Delta r}\nonumber \\ B_{\theta,i,j+1/2,k} &=& \frac{\sin \theta _{j+1/2}(\Phi _{i,j+1,k}-\Phi _{i,j,k})}{r_i \Delta \cos \theta } \\ B_{\phi,i,j,k+1/2} &=& \frac{\Phi _{i,j,k+1}-\Phi _{i,j,k}}{r_i \sin \theta _j \Delta \phi }.\nonumber \end{eqnarray} \tag{ 21 }$$

Note the sin θ/Δcos θ factor in the θ derivative for the uniform cos θ grid. For the uniform θ grid this is replaced with 1/Δθ. The divergence of the magnetic field, i.e., the Laplace of Φ, is obtained as

$$\begin{eqnarray} 0 = (\nabla ^2 \Phi)_{i,j,k} &=& \frac{r^2_{i+1/2}B_{r,i+1/2,j,k} - r^2_{i-1/2,j,k}B_{r,i-1/2,j,k} }{r^2_i \Delta r} \nonumber \\ &&{+}\frac{\sin \theta _{j+1/2}B_{\theta,i,j+1/2,k} - \sin \theta _{j-1/2}B_{\theta,i,j-1/2,k} }{r_i \Delta \cos \theta } \nonumber\;\; \\ &&+\frac{ B_{\phi,i,j,k+1/2} - B_{\phi,i,j,k-1/2} }{r_i\sin \theta _j\Delta \phi }. \end{eqnarray} \tag{ 22 }$$

Again 1/Δcos θ is used for the uniform cos θ grid, while on the uniform θ grid this is replaced with 1/(sin θΔθ).

The magnetogram boundary condition is applied by setting the inner ghost cell as

$$\begin{equation} \Phi _{0,j,k} = \Phi _{1,j,k} - \Delta r M^{\prime }_{j,k}, \end{equation} \tag{ 23 }$$

where $M^{\prime }_{j,k} = M_{j,k}-\bar{M}$ is the magnetogram with the average field (i.e., the monopole due to observation errors)

$$\begin{equation} \bar{M} = \frac{1}{4\pi }\sum _{j,k} (\Delta \cos \theta)_j (\Delta \phi)_k M_{j,k} \end{equation} \tag{ 24 }$$

removed. The zero potential at $r_{N_r+1/2}=R$ is enforced by setting the ghost cell as

$$\begin{equation} \Phi _{N_r+1,j,k} = -\Phi _{N_r,j,k}. \end{equation} \tag{ 25 }$$

The boundary conditions at the poles are a bit tricky. Cells (i, 1, k) and (i, 1, k') are on opposite sides of the north pole if k' = mod(k − 1 + N_ϕ/2, N_ϕ) + 1. Therefore, the ghost cells in the θ direction are set as $\Phi _{i,0,k}=\Phi _{i,1,k^{\prime }}$ and $\Phi _{i,N_{\theta }+1,k}=\Phi _{i,N_{\theta },k^{\prime }}$. We note here that N_ϕ is assumed to be an even number. The periodic boundaries in the ϕ direction are simple: $\Phi _{i,j,0}=\Phi _{i,j,N_{\phi }}$ and $\Phi _{i,j,N_{\phi }+1}=\Phi _{i,j,1}$.

We need to find Φ_{i, j, k} that satisfies the discrete Laplace equation (22) with the boundary conditions applied via the ghost cells. The initial guess is Φ = 0, which provides a non-zero residual because of the inhomogeneous boundary condition at the inner boundary applied by Equation (23). We use this residual with a negative sign as the right-hand side of the Poisson equation (∇²Φ)_{i, j, k} = R_{i, j, k}, and use Φ_{0, j, k} = Φ_{1, j, k} as the new homogeneous inner boundary condition instead of Equation (23). We use the Krylov-type iterative method BiCGSTAB (van der Vorst 1992) to find the solution. The linear system is preconditioned with an Incomplete Lower–Upper decomposition (ILU) preconditioner to speed up the convergence. We use ILU(0) with no fill-in compared to the original matrix structure, so the preconditioner is a diagonal matrix, but its elements depend on all elements of the original matrix. We have implemented a serial as well as a parallel version of the algorithm. In the parallel version, the preconditioner is applied separately in each sub-domain. FDIPS finds an accurate (down to 10⁻¹⁰ relative error) solution on a 180³ grid in less than 1000 iterations. Even running serially, this takes less than an hour on today's computers.

Once the solution is found in terms of the discrete potential Φ_{i, j, k}, we apply the original boundary conditions including Equation (23) and calculate the magnetic field with Equation (21). The divergence of the magnetic field will be zero to the accuracy of the Poisson solver. The curl of the magnetic field will be zero in a finite difference sense simply because it is constructed as the discrete gradient of the potential. The boundary condition at the inner boundary is also satisfied exactly: B_{r, i = 1/2, j, k} = M'_{j, k}. Averaging rB_ϕ and rB_θ to the $r_{N_r+1/2} = R$ position at the outer boundary also gives exactly zero tangential fields due to Equations (25) and (21). Depending on the application, we may interpolate the potential magnetic field onto a co-located grid, or use it on the original staggered grid.

Figure 6 shows the solution of the magnetic field obtained with the finite difference solver FDIPS on a 150 × 360 × 180 grid. Since we use the same uniform cos θ grid as the magnetogram, the radial field obtained is identical with the magnetogram at r = 1. The tangential components agree well with the remeshed spherical harmonics solution shown in Figure 5. It took 1166 iterations to get a solution with a relative accuracy of 10⁻¹⁰. The run time was almost exactly one hour on a 2.66 GHz Intel CPU.

Zoom In Zoom Out Reset image size

We have discussed various ways to obtain the potential field solution based on solar magnetograms. While spherical harmonics provide an efficient and elegant method, there are some subtle restrictions that require attention. If one wants to use many spherical harmonics (the same order as the number of magnetogram pixels in the co-latitude direction), the magnetogram data on the N_θ × N_ϕ grid have to be remeshed onto a uniform θ grid with N'_θ × N_ϕ points, N'_θ must be an odd number, and the new grid must include both poles. After the remeshing, the maximum degree of harmonics N is limited only by the anti-alias limit to min (2N'_θ/3, N_ϕ/3). We used a simple linear interpolation for the remeshing.

The remeshing can be avoided by the use of a 3D finite difference scheme. One can use the original magnetogram grid, and the only freedom is in choosing the radial discretization. The finite difference scheme provides a solution that is fully compatible with the boundary conditions, and the solution has zero divergence and curl in the finite difference sense.

Figure 7 compares the solutions obtained with the three methods along the radial direction for a fixed latitude θ = 765 and longitude ϕ = 30°. The spherical harmonics series were truncated at N = 120 for both the naive and remeshed methods. The naive spherical harmonics algorithm gives incorrect results close to the solar surface where the high-order harmonics dominate. This is most obvious for the radial component, which is given at r = 1 by the magnetogram, and it is exactly reproduced by the finite difference scheme. The latitudinal component at r = 1 is also very different from the values given by the remeshed harmonics and the finite differences. The latter two methods agree reasonably well with each other. For radial distances above r = 1.05, all three methods agree quite well.

Zoom In Zoom Out Reset image size

So far we restricted our example to a GONG magnetogram taken at the solar minimum. If one uses an MDI magnetogram during solar maximum, the largest magnetic fields are much stronger (order of 1000 G) and the resolution of the magnetogram is much finer (order of 1000 pixels). The finer magnetogram resolution allows the use of a larger number of harmonics, even when using the original magnetogram grid (naive approach). But the strong and sharp gradients in the magnetogram will bring out another problem with the spherical harmonics approach: the ringing effect. The ringing is due to the so-called Gibbs phenomenon: the step-function-like magnetogram data result in high-amplitude high-order harmonics in Fourier space. The ringing effect and other artifacts are discussed in great detail by Tran (2009).

Figure 8 demonstrates this effect on the 3600 × 1080 resolution MDI magnetogram for Carrington Rotation 2029 (from 2005 April 21 to May 18), with the maximum radial field strength around ±3000 G. The remeshed harmonics method with N = 90 is compared with the finite difference method on a 150 × 360 × 180 grid (the magnetogram data are coarsened to a 360 × 180 grid). In the spherical harmonics solution, the ringing is very clearly visible around the active regions, both in the radial and latitudinal components. The finite difference scheme, on the other hand, shows no sign of ringing in either component. This is obvious for the radial component, which simply coincides with the coarsened magnetogram, but for the latitudinal component it is due to the fact that the finite difference solution of the Laplace equation does not suffer from ringing artifacts even for discontinuous boundary data. For the spherical harmonics approach, the ringing becomes weaker with increased number of harmonics, but is still quite apparent even for N = 180 (not shown). The results of the remeshed and naive harmonics methods are essentially the same up to N = 180, i.e., the ringing is not due to the remeshing of the magnetogram.

Zoom In Zoom Out Reset image size

In terms of computational efficiency, a good implementation of the spherical harmonics scheme is much faster than the finite difference scheme. In fact, it may be more costly to construct the potential field solution on a 3D grid from the spherical harmonics coefficients than to obtain the coefficients themselves. Our Fortran 90 code can obtain the spherical coefficients up to N = 60, 90, and 120 deg in 1, 1.8, and 3.3 s, respectively, while the reconstruction of the solution on the 151 × 361 × 180 grid takes 5, 12, and 20 minutes, respectively. All timings were done on a single 2.66 GHz Intel CPU. The reconstruction cost can be improved by running the code in parallel, and/or truncating the series in parts of the grid where the higher-order harmonics have a negligible contribution. We also note that going beyond about N = 360 harmonics becomes fairly complicated (Potts et al. 1998).

The computational cost of the finite difference scheme scales with the number of grid cells and the number of iterations. The number of iterations is fairly constant for multigrid-type methods, but for the Krylov sub-space schemes it grows with the problem size, although slower than linearly. The finite difference scheme can be speeded up by parallelizing the code, which is fairly straightforward for the Krylov sub-space schemes. Since we limit the ILU preconditioning to operate independently on the sub-domains of each processor, the preconditioner becomes less efficient as the number of processors increases, which results in an increase in the number of iterations. To minimize this effect, the parallel FDIPS code splits the grid in the θ and ϕ directions only, so the sub-domains in each processor contain the full radial extent of the grid. Our experiments confirmed that by using this domain decomposition, the number of iterations indeed does not depend much on the number of processors. Our largest test so far involves a 450 × 540 × 1200 grid with 30 times more cells than the 150 × 180 × 360 grids discussed in most of this paper. For the largest problem we need about 8500 iterations to reach the 10⁻¹⁰ relative accuracy, a factor of nine increase relative to the smaller problem discussed here. Using 108 CPUs, the solution is obtained in about 5.3 h.

Despite the various limitations, for some applications the spherical harmonics approach may still be preferred, for example if the solution is needed to obtain a spherical power spectrum of the solar magnetic field. If the solution is to be used in a finite difference code, the finite difference solution is probably preferable. We are using the FDIPS code to generate the potential field solution as the initial field for our solar corona model (van der Holst et al. 2010).

This paper attempts to call the attention of astrophysicists and solar physicists to the limitations and potential pitfalls of using the spherical harmonics approach to obtain a potential field solution. The spherical harmonics representation of the potential field solutions is available from several synoptic magnetogram providers, although the details of the method used to obtain the spherical harmonics are not always clear. A spherical harmonics based PFSS package implemented in IDL is available as part of the Solar-Soft library (http://www.lmsal.com/solarsoft). This package uses the magnetogram remeshing technique either onto the Chebyshev (uniform θ) or the Legendre co-location points.

We are not aware of any publicly available code that uses finite differences to solve this particular problem. To allow other researchers to use and compare the two approaches, we make our finite difference code FDIPS publicly available at the http://csem.engin.umich.edu/fdips/ Web site.