© ESO 2020

1. Introduction

The Sun’s coronal magnetic field configuration is an important component in understanding the physics of the heliosphere and the solar dynamo. Due to the vast extent of the heliosphere, a direct measurement of its global structure is not possible. Therefore, computational modeling tools are employed to approximate its structure. Modern high-accuracy algorithms, like full magnetohydrodynamic (MHD) solvers, compute a vast set of physical phenomena and require significant computing power. These models rely on additional modeling assumptions that are difficult to verify and on boundary conditions that have to be modeled because they cannot be observed directly.

Simpler models exist that produce less precise results, but that can be computed orders of magnitude more quickly. For some scientific efforts like long-running data evaluation tasks, for example characterizing solar wind streams, rapid computation of the heliospheric magnetic field over long periods is crucial, while the overall accuracy can be lower without affecting the large-scale results of these studies.

One of the earliest models that was used to model the solar coronal magnetic field is the potential field source surface (PFSS) model (Altschuler & Newkirk 1969; Schatten et al. 1969). A brief description of this model is presented in Sect. 2.2. The PFSS model returns an analytic expression that allows us to predict the magnetic field configuration of the Sun between the photosphere and a virtual spherical surface, called the source surface, at a height of a few solar radii in the corona. The boundary condition for the PFSS model is that all magnetic field lines have to be oriented radially at the source surface, which is in accordance with observations farther out in the heliosphere.

An improved version of this model is the current sheet source surface (CSSS) model (Zhao & Hoeksema 1995), which adds a second virtual sphere, the cusp surface, between the photosphere and source surface. Above this surface, all magnetic field lines are required to be open, but do not have to be oriented radially. This allows more freedom when modeling the magnetic field.

Slightly more sophisticated models build upon the force-free approach which neglects external forces on a restricted domain near the solar surface. The most general group of these models are the nonlinear force-free models (Aly 1989; Wiegelmann 2008). Simplifying this physical model by assuming a current-free domain in addition to it being force-free leads to the PFSS model that allows the extrapolation of the coronal magnetic field using only the line-of-sight component of the photospheric magnetic field configuration (also called a synoptic magnetogram). The photospheric magnetograms can be obtained from direct observations and the procedure has been applied for many decades. While the PFSS model itself makes more simplifying assumptions than the full MHD approach, it relies less on modeled parameters, which can potentially introduce false or inaccurate assumptions to the model.

The assumption that the source surface is spherical simplifies the mathematical framework as well as the computations significantly. Without hints as to exactly what the source surface looks like, as well as a lack of computing power back when this model was created, the only reasonable starting point was to assume a spherical source surface. Also, the very low resolution of available magnetic observations of the photosphere at that time made more accurate model assumptions meaningless.

Since then, several suggestions have been made that question the merit of the spherical source surface. Schulz et al. (1978), Schulz (1997) have developed an algorithm to alter the shape of the source surface. They proposed surfaces of constant magnetic flux (isogauss) to function as the source surface, and presented a model with a surface that has greater heights above the poles compared to the equator. In this model the free parameter of the PFSS model (the source surface height R_ss) is substituted for the constant magnetic flux B₀ of the isogauss surface, which is determined by a hypothetical solar internal dipole and the constraint of magnetic field lines to be normal to this surface. The strategy behind this approach is to better match field lines to the predictions of an MHD implementation that produces field lines that are not quite radially oriented at the height of the spherical source surface.

Expanding on this idea Levine et al. (1982) proposed a nonspherical source surface that is determined by three free parameters, the mean height and two parameters defining the shape. Levine et al. (1982) also deviated from the constraint of magnetic field lines that are strictly perpendicular to the source surface, and compared the orientation of computed field lines to total solar eclipse observations.

Riley et al. (2006) computed isosurfaces of B_r/|B| = 0.97 by employing a magnetohydrodynamic solver and found shapes that resemble prolate spheroids with indentations at the poles. While these isosurfaces do not constitute the source surface, they illustrate the deviation of the magnetic field configuration from spherical symmetry.

The spherical source surface is the simplest assumption and it allows fast computation of the solar magnetic field. However, observations show that the physical conditions of the outer solar coronal plasma depend on longitude as well as latitude, and that polar regions exhibit a more radial orientation of the magnetic field than at the equator (see, e.g., McComas et al. 1998). An oblate elliptical source surface creates a magnetic field that is consistent with this observation, as is discussed in Sect. 2.4.

In the past decades the increase in computing power and space-bound magnetic field observations, have given rise to more intricate models. MHD solvers have allowed a wider array of physical phenomena to be incorporated into the derivation of the solar magnetic field. While potentially more accurate, these models require a larger set of input parameters and model assumptions that are difficult to verify with current observations. The simplicity of the PFSS model allows calculating the rough structure of the solar magnetic field comparably quickly. It is therefore desirable to improve this simple model without adding the complexity of a full MHD approach.

In this work we suggest modifying the classical PFSS model to incorporate elliptical source surfaces. This increases the number of free parameters from one (the source surface height R_ss) to two (adding the ellipticity or deviation from sphericity of the source surface A). In contrast to the approaches mentioned above, we chose to implement a finite difference solver rather than an analytical solver. This allows easier adjustments of the utilized source surface and gives us the possibility to increase solution accuracy by adding more grid points to regions of strong magnetic gradients. We hope that this minor adjustment aids in the modeling of the large-scale structure of the solar magnetic field more accurately and gives insight into the reliability of the PFSS model for different stages of the solar activity cycle.

Evaluating the validity and accuracy of the various existing models is a field of study of its own. Since the phenomena in question cannot be recreated in their entirety in a laboratory environment, the scientific community is forced to browse sparsely available spacecraft data to look for indicators of correctness or lack thereof. In this work we focus on presenting alterations to the classical PFSS model and the computational results produced by them. In Sect. 4 we discuss options of evaluating the results produced by the different PFSS models.

2. Methods

The PFSS model is often used to investigate the link between the interplanetary medium and the solar surface. Because of its simplicity it has strengths and weaknesses, which we investigate here by implementing three different versions of the PFSS. In Sect. 2.1 we give a brief summary of the mathematical framework underlying all the implementations. The first implementation recreates the framework developed by Zhao & Hoeksema (1993) so that subsequent alterations to the model can be compared to a version that has been thoroughly evaluated. This classical implementation is briefly described in Sect. 2.2. We created an implementation of this approach and call it the spherical harmonic coefficient (SHC) version. The second version is a numerical finite differences solver that solves Laplace’s Equation (Eq. (1)) on several grid points throughout the computational domain of the PFSS model and is described in Sect. 2.3. The third version is an alteration of the second; it allows for an elliptical source surface and is described in Sect. 2.4. A summary of our treatment of input data as well as the parameters for the computational solver are presented in Sect. 2.5.

2.1. The PFSS model

Altschuler & Newkirk (1969) and Schatten et al. (1969) independently proposed a magnetostatic model for the solar corona. They partition the domain of interest into three regions (Fig. 1).

Fig. 1. Three regions underlying the PFSS model: region 1 is the inside of the Sun, region 2 is the computational domain between the photosphere and the source surface, and region 3 is the interplanetary space where the solar wind flows radially outward. Figure adapted from Schatten et al. (1969).

Employing several instruments the magnetic field can be measured at the photosphere, which is the boundary separating regions 1 and 2. Region 2 is the computational domain of the PFSS model. The boundary between regions 2 and 3 is called the source surface. Above the source surface, the solar wind dominates the magnetic field configuration (region 3). At least during the quiet times of the solar activity cycle, and aside from violent eruptions, the Sun’s photosphere displays features that persist for several Carrington rotations. Therefore, in a first approach the lower region of the solar corona (region 2) can be assumed to be electrostatic, or $\frac{\partial \mathit{E}}{\partial t}=\mathbf{0}$. Furthermore, it is assumed that due to the sharp decrease in particle density above the photosphere and with a smaller decrease in magnetic field strength, the electric current density can be neglected to some point, or j = 0. Ampére’s law then states

$$\begin{array}{c}\hfill \mathrm{\nabla }\times \mathit{B}={\mu }_0(\mathit{j}+{ϵ}_0\frac{\mathit{\partial }\mathit{E}}{\partial t})=\mathbf{0},\end{array}$$

where B is the magnetic flux density, E is the electric field, j is the electric current density and μ₀ and ϵ₀ are the permeability and permittivity of free space, respectively.

A curl-free vector field can be described as the gradient of a scalar potential (∇ × ∇f = 0 for twice the continuously differentiable f, the curl of a gradient vanishes everywhere). Therefore, we write B = −∇Ψ, where Ψ is the scalar magnetic potential. Gauss’s law then states

$$\begin{array}{c}\hfill \mathrm{\nabla }·\mathit{B}=-\mathrm{\nabla }·\mathrm{\nabla }\mathrm{\Psi }=-\mathrm{\Delta }\mathrm{\Psi }=0.\end{array}$$(1)

Equation (1) can be integrated within region 2, given two boundary conditions, one at the photosphere and one at the source surface. At some height the solar wind carries the magnetic field outward, where it is said to be frozen in. The solar wind is advected outward radially, so at the source region of this flow, the magnetic field lines have to be aligned radially as well. Therefore the upper boundary condition is given by the restriction of the magnetic field lines to be perpendicular to the source surface (Neumann boundary condition). For the lower boundary, the photosphere, the magnetic field configuration is known and is supplied to the algorithms by synoptic (line of sight) magnetograms; therefore, a Dirichlet boundary condition is applied.

For the lower boundary condition, two main approaches are employed regularly by the scientific community, called the radial approach and the line-of-sight approach (Wang & Sheeley 1992; Altschuler & Newkirk 1969). Due to its simplicity and the widespread availability of data for comparison, we employ the radial approach for this study.

2.2. Spherical harmonic coefficient implementation

The mathematical framework for the approach to model the solar coronal magnetic field using spherical harmonic coefficients and associated Legendre polynomials can be found in Altschuler & Newkirk (1969) and Chapman & Bartels (1940). All that is necessary to recreate the magnetic field is an implementation of the associated Legendre polynomials $P_l^m$ as well as the computed harmonic coefficients g_lm and h_lm (for explanations of these symbols, see Chapman & Bartels 1940). The magnetic field configuration at any point in region 2 of the PFSS model can then be acquired by evaluating an analytic expression (see, e.g., Eqs. (8)–(10) in Altschuler & Newkirk 1969).

Today the harmonic coefficients are computed by several groups, including the John M. Wilcox Solar Observatory at Stanford University (Zhao & Hoeksema 1993) and published on their website (Hoeksema 2020). We used synoptic magnetograms from the Wilcox Solar Observatory (Duvall et al. 1977) and compared the results from the Stanford PFSS implementation with our own. Within floating-point precision, our SHC solver produces the same coefficients as published by Stanford.

2.3. Grid approach

Instead of fitting the observed photospheric magnetograms to spherical harmonic functions, we have developed a numeric solver that works on a three-dimensional grid and employs finite differences. This allows us to deform the grid in other implementations to incorporate elliptical source surfaces, while also giving us a tool to compare our results with the implementations of other groups.

The computational grid stretches from the photosphere (at r = R_⊙) to the source surface in the radial direction, from the northern boundary supplied by the magnetogram to the southern boundary in the meridional direction and around the sphere in the zonal direction without boundaries. Grid spacing is equidistant in zonal direction and follows a sine-latitude distribution in meridional direction, as do the underlying synoptic magnetograms. In the radial direction spacing between grid points increases from the photosphere to the source surface geometrically (see Sect. 2.5 for more details).

The solution method is an explicit time-stepping algorithm that solves Laplace’s Eq. (1) at each grid point for the potential field Ψ. This is done by evaluating the analytic expression utilizing finite differences and solving for the magnetic potential Ψ_ijk, where i, j, and k denote the radial, meridional, and zonal position in the numerical grid, respectively.

In the initial state, the potential field is set to zero at each grid point except the lowest grid shell which is derived from the synoptic line-of-sight magnetograms. The potential at the uppermost grid shell (at the source surface, r = R_ss) is kept constant throughout the iteration process, thus implementing the radial boundary condition at the source surface. The potential at each grid point is stored and compared to that obtained in the following time step. Let ${\mathrm{\Psi }}_{\mathit{ijk}}^t$ denote the magnetic potential at position i, j, k in time step t, and ${\mathrm{\Psi }}_{\mathit{ijk}}^{t-1}$ the potential computed in the previous time step. The algorithm terminates if the maximum relative deviation to the previous time step at all grid points drops below a specified accuracy threshold p, or

$$\begin{array}{c}\hfill {\displaystyle \underset{\mathit{ijk}}{max}}\left(e_{\mathit{ijk}}\right)={\displaystyle \underset{\mathit{ijk}}{max}}\left(\frac{|{\mathrm{\Psi }}_{\mathit{ijk}}^t-{\mathrm{\Psi }}_{\mathit{ijk}}^{t-1}|}{|{\mathrm{\Psi }}_{\mathit{ijk}}^t|}\right)<p.\end{array}$$(2)

The numerical grid solver produces nearly the same polarity configuration at the source surface as the SHC implementation. Closer to the photosphere the numerical solver can resolve a finer structure than the SHC solver, which is discussed in Sect. 3.

2.4. Grid approach with an elliptic source surface

Ideally, an accurate algorithm for the computation of the solar magnetic field would accommodate arbitrary heights of the source surface at every longitudinal and latitudinal position. We propose a small step in this direction by implementing an ellipsoidal source surface controlled by a single parameter A. An oblate source surface in the PFSS model creates a magnetic field that reaches radial orientation closer above the poles compared to equatorial regions, as can be seen in the renderings in Figs. 2–5. Let R_ss be the source surface height of the classical PFSS model. The actual source surface height r_ss in our model is given by

$$\begin{array}{cccc}& \mathrm{oblate}\mathrm{case}:\hfill & \hfill r_{\mathrm{ss},\mathrm{equator}}& =A·R_{\mathrm{ss}}\hfill \\ \hfill & r_{\mathrm{ss},\mathrm{poles}}\hfill & \hfill =R_{\mathrm{ss}}\\ \hfill & \mathrm{prolate}\mathrm{case}:\hfill & \hfill r_{\mathrm{ss},\mathrm{equator}}& =R_{\mathrm{ss}}\hfill \\ \hfill & r_{\mathrm{ss},\mathrm{poles}}\hfill & \hfill =A·R_{\mathrm{ss}}.\end{array}$$

Fig. 2. Magnetic field polarity configuration at different heights. Results shown for a prolate ellipsoidal source surface with ellipticity 2.0 (left column), for the classical spherical source surface (center), and for an oblate ellipsoidal surface again with ellipticity 2.0 (right). All figures were created using our grid solver. The source surface height for all models is 2.5 R⊙ (minor half-axis in the ellipsoidal cases). Depicted is Carrington rotation 2066. Data for the lower boundary was obtained from MDI onboard SOHO. First row: three-dimensional rendering of a few magnetic field lines for each model as well as a cut through the height levels depicted below. Rows 2 to 5: magnetic field polarity configuration at height levels 1.0R⊙, 1.5R⊙, 2.0R⊙, and 2.5R⊙ (minor half-axis/radius) in a synoptic Carrington format. The height levels correspond to ellipsoids at distances of 0%, 33%, 67%, and 100% between the photosphere and source surface. Red magnetic field lines and pixels are directed inward, blue lines are directed outward, and cyan indicates closed field structures. The projection of Earth on the height levels is drawn as a black ascending line near the solar equator. Black ellipses, circles, and squares are inserted to highlight differences in the regimes of interest (see text for details).

As the Sun still needs to be approximated as a sphere, we cannot simply employ a homogeneous elliptical grid. Hence, we created a grid that exhibits spherical symmetry at the photosphere while incrementally deforming higher grid shells to the ellipsoidal shape. A brief explanation of the grid is summarized in Appendix A.

The elliptical implementation works in the same manner as the spherical one by repeatedly solving Laplace’s equation at all grid points until the accuracy threshold p is surpassed at all grid points (see Eq. (2)).

Magnetic field lines are oriented perpendicular to the source surface, which defines an isopotential. In the spherical case, this means field lines are already oriented radially. To obtain the radial orientation of the magnetic field lines in the elliptical model, a spherical surface can be positioned above the source surface, and some form of interpolation technique, for example employing splines, can be applied to “bend” the magnetic field lines into the radial orientation.

2.5. Parameters for the grid solver and treatment of input magnetograms

For all model evaluations we employ a source surface at a heliocentric height of 2.5R_⊙ (minor half-axis in the elliptical version). As input data, we used synoptic magnetograms from the Wilcox Solar Observatory for testing purposes and comparisons with data products published by Hoeksema (2020). Because the resolution of these maps is low, no image processing needs to be applied. For high-resolution grid tests and for the plots presented in this work we employed synoptic magnetograms from the Michelson Doppler Imager (MDI) on board the Solar and Heliospheric Observatory (SOHO; Scherrer et al. 1995). The high-resolution synoptic magnetic maps produced by MDI are scaled down to 87 × 175 pixels using a Lancosz filter. In a second step these magnetograms are corrected for the monopole offset which is introduced by small changes of the photospheric magnetic field during the data acquisition period of about 27 days and for data gaps near the poles.

A well-known difficulty of numerical models is the trade-off between high accuracy and computational demand. After systematically testing and comparing results of several distributions and densities of grid points we chose our numerical grid to have 35 × 87 × 175 grid points in radial, meridional, and zonal directions. Radial spacing increases geometrically from the photosphere up to the source surface with a geometric factor q ≈ 3.3%. The radial position of grid points on a height level i (i ∈ {1, …, N_r − 1}) in the spherical grid is r_i = r_i − 1 + (r₁ − R_⊙)⋅q^i − 1, where N_r = 35 is the number of height levels in the radial direction and r₁ − R_⊙ ≈ 18.7 Mm. For the elliptic case, this height dependence is distorted according to the stretching treatment presented in Appendix A. This geometric increase in radial grid point spacing allows higher accuracy of the solver near the photosphere where magnetic gradients are strongest, while reducing the computational footprint in the higher regions where a lower accuracy is sufficient for acceptable results.

Polar boundaries are introduced by computing virtual values at the poles as being the average value of all northernmost and southernmost grid points at that specific height. The virtual polar points serve as the outermost neighbor for all adjacent grid points that make up its average. This decreases computation time as information of the solution is allowed to directly travel over the poles rather than only in the zonal direction.

As termination criterion (Eq. (2)) we used p = 0.01 which means the value change from one time step to the next is below 1% at all grid points. Magnetic field line tracking is done by employing an adaptive Runge-Kutta-Fehlberg method of fourth order (RKF45), tri-linearly interpolating between computational grid points.

The spherical PFSS model is a special case of the elliptical PFSS model with ellipticity A = 1.0. The elliptic solver gives the same values at all grid points within floating-point rounding accuracy for this special case compared to the purely spherical solver.

Our implementation has been written in C/C++ and employs the CUDA framework (version 9.2) by NVIDIA. The time-stepping solution process is quite basic and can be sped up significantly by employing more sophisticated solution processes. However, the computation time for one Carrington rotation with the resolution and accuracy presented here takes between 10 and 20 minutes on the NVIDIA GTX Titan which was released in 2013. For the studies presented here, the computation time is sufficiently fast.

3. Results

Without answering the question of which parameter for the ellipticity gives the most realistic coronal magnetic field configuration (if any), we illustrate here the qualitative differences between the spherical and elliptical PFSS models. We chose Carrington rotation 2066 during the minimum between solar activity cycles 23 and 24 in early 2008 to illustrate the differences the model parameters incur. There are only a few and weak CMEs registered for this Carrington rotation (Yashiro 2020), and a few coronal holes and active regions (Barra et al. 2009).

Figures 2 and 3 depict the solar coronal magnetic field configuration. The figure consists of a three-dimensional rendering of the magnetic field configuration as seen from the vernal equinox in the first row. The other rows show the magnetic field polarities at different heights between the photosphere and the source surface. For illustrative purposes only, we chose a strong ellipticity (A = 2) of the source surface which we assume is higher than a realistic configuration would exhibit.

Fig. 3. Same as in Fig. 2, but the height levels are spheres. The heights are again 1.0R⊙, 1.5R⊙, 2.0R⊙, and 2.5R⊙. The dashed line in the first row depicts the source surface, which here is not the same as the uppermost height level examined. The uppermost levels exhibit open field structures in the ellipsoidal cases because the sphere of 2.5R⊙ is partially within the source surface and touches it at the equator (prolate case) or the poles (oblate case).

The three-dimensional renderings consist of two magnetic field line mappings from the source surface down to the photosphere and vice versa. At the specific height (photosphere, source surface), an equidistant two-dimensional grid is spanned in sin(latitude) and longitude with −14.5/15.0 ≤ sin (latitude) ≤ +14.5/15.0 and 0 ≤ longitude  <   2π. Each of these grid points constitutes the starting point of a magnetic field line, which is traced throughout the computational domain. Blue magnetic field lines have a positive sign pointing outwards, while red lines have a negative sign pointing inwards. Closed field lines originating and ending on the photosphere are in cyan. In the renderings, 15 × 30 magnetic field lines are illustrated starting on the source surface and the same number starting on the photosphere.

Similarly, rows 2 to 5 are cuts of magnetic field line mappings originating at intermediate heights but with a higher resolution of 200 × 400 field lines, corresponding to 200 × 400 pixels. The format is similar to synoptic Carrington maps. The uppermost and lowermost pixels again are positioned at sin (latitude_max) = +14.5/15.0 and sin (latitude_min) =−14.5/15.0, respectively. Each pixel represents the polarity of the magnetic field line at that pixel center position with the same color scheme as the field lines in the renderings. Due to the sine-latitude spacing pixels near the poles cover a larger area than pixels at the equator, effectively reducing resolution at high latitudes.

In Fig. 2 the intermediate maps correspond to ellipsoidal surfaces which grow in ellipticity, as does the underlying computational grid when approaching the source surface from below. Figure 3 shows intermediate height maps originating on spheres between the two computational boundaries. These represent height levels of 1.0R_⊙, 1.5R_⊙, 2.0R_⊙, and 2.5R_⊙ heliocentric distance. We note that the middle column for the spherical model is the same in both figures.

Figure 2 considers global differences among the three presented models. Each height level depicts the magnetic polarity structure at the same relative distance between the boundaries, which correspond to different absolute heliocentric heights for each column. Figure 3 helps to visualize local differences. Each pixel represents the same physical position in all three columns. While most structures are present in every model, their general appearance can vary greatly. In Figs. 2 and 3 the pronounced current sheet warp (rectangles) at longitude 270°, for example, is less pronounced in the prolate case and extends to lower latitudes in the oblate case compared to the spherical reference model. Another example is the open positive structure slightly south of the equator at longitude 180° (ellipses and circles) which is present in the spherical and prolate models, but missing in the oblate version. Also, closed field structures extend higher in the computational domain in the oblate case, as can be easily seen in the third row of Fig. 2.

Another characteristic that can be computed by only considering magnetic field data is the flux-tube expansion factor $f_s=B_{\odot }{R}_{\odot }^2/B_{\text{ss}}{r}_{\text{ss}}^2$ (Wang & Sheeley 1990), where B_⊙ is the magnetic flux density at the photosphere, B_ss is the magnetic flux density at the source surface, R_⊙ is the heliocentric distance of the photosphere, and r_ss is the heliocentric distance of the source surface.

The flux-tube expansion factor is inversely related to solar wind speed and offers a means of comparing model prediction to in situ spacecraft data (ibid). In Figs. 4 and 5, we illustrate the expansion factor for the three PFSS models. The format is similar to Figs. 2 and 3, but here each pixel is color-coded with the expansion factor for that specific field line. The oblate PFSS model exhibits higher expansion factors at lower latitudes than the other two models. This suggests lower solar wind speeds in the oblate model, which should be verifiable by analyzing spacecraft data.

Fig. 4. Height levels show the expansion factor. Each pixel again corresponds to one field line, and the color represents its expansion factor between photosphere and source surface.

Fig. 5. Same as Fig. 4, but for height levels as spheres.

The expansion factor also illustrates some differences in resolution between the classic implementation employing spherical harmonic coefficients and our grid solver. Figure 6 depicts the expansion factor height levels for the SHC implementation of orders 9 and 20 as well as our implementation. Unsurprisingly, the middle column, which corresponds to a maximum principle order 20, exhibits more detail than the left column, which corresponds to the classical approach of order 9. Our grid implementation in the right column offers even more detail at the cost of longer computation times. While the classical SHC implementation of order 9 takes less than 1 minute using a single thread on the CPU (Intel Xeon E5-1650), our grid solver takes about 15 minutes employing the massively parallel architecture of a GPU (GTX Titan).

Fig. 6. Comparison of expansion height maps for the classical PFSS implementation using spherical harmonic coefficients with order 9 in the left column, with order 20 in the middle column, and our grid solver in the right column. Again Carrington rotation 2066 is depicted. All models employ a spherical source surface at 2.5R⊙ heliocentric distance. The height levels in rows 1 to 4 are spheres at heliocentric distances of 1.0R⊙, 1.5R⊙, 2.0R⊙, and 2.5R⊙.

4. Discussion and conclusions

The alteration to the PFSS model to incorporate elliptical source surfaces allows us to tweak the magnetic field computations without needing to employ the full set of MHD assumptions and computational complexity. Only one parameter has been added that needs to be determined and evaluated. In this regard, it is comparable with other improvements of the PFSS model like CSSS, which includes two additional parameters, namely the height of the cusp surface and a length scale of the assumed horizontal electric currents in the corona (Zhao & Hoeksema 1995).

The CSSS model still relies on a spherical symmetry which is a strict constraint for a physical system that we know deviates from this type of symmetry. It is known that the “true” source surface, if it even exists, is most probably not elliptical either (see, e.g., Cohen 2015; Schulz et al. 1978; Riley et al. 2006; Panasenco et al. 2020). It remains to be seen whether the elliptical PFSS approximates the true magnetic field configuration better or worse than the other improvements that have been developed by the community in the past decades. It might also be useful to combine the elliptical source surface with these models, thereby merging possible advantages either model has over the original spherical PFSS model.

The finite difference solver allows more complex shapes to act as source surfaces. For example, the polar indentations found by Riley et al. (2006) may be modeled using the same techniques presented here.

Phenomenologically, there are a few differences between the three grid models presented in this work. In the oblate elliptical model, magnetic field lines tend to bend towards the poles, while in the prolate model they bend away from the poles (see the renderings in Figs. 2–5). Hence, particles traveling along these field lines will have slightly different trajectories. In the oblate case, closed magnetic field structures extend higher relative to the source surface compared to the other two models.

One way to evaluate PFSS models is to compare in situ spacecraft measurements of the heliospheric magnetic field with predictions of the PFSS model. After measuring the solar wind speed at the spacecraft, the likely footpoint of the solar wind plasma package on the source surface can be computed by tracing the Parker spiral. The magnetic polarity at the footpoint, as computed by the PFSS implementation, can then be compared to the magnetic field measured at the spacecraft. While there are phenomena above the source surface that might alter the magnetic field orientation as well as inaccuracies in tracing the correct Parker spiral, for example due to solar wind acceleration processes, a better magnetic model of the corona is expected to increase the agreement of spacecraft data with predicted model data.

Going one step further it is possible to track the source surface footpoint of the spacecraft down to the photosphere along the magnetic field lines computed by the model. The photospheric footpoint positions can then be compared to the appropriate positions in extreme ultraviolet (EUV) images, which originate close to the photosphere. It is generally assumed that open magnetic flux structures map to darker regions in EUV images (Huang et al. 2019). A better magnetic field model should therefore map more often to dark regions in EUV maps than do worse models.

Evaluating predictions of the PFSS models utilizing spacecraft data is complicated due to the distances between spacecraft (typically at heliocentric distances of about 1AU) and the computational domain (below a few R_⊙). Connecting spacecraft positions to the source surface introduces errors increasing with distance. Fortunately, two recent missions might help this endeavor. The Parker Solar Probe (Fox et al. 2016) is already collecting data and Solar Orbiter (Müller & Marsden 2013) was launched in February 2020. The Parker Solar Probe is scheduled to reach a perihelion heliocentric distance of less than 10R_⊙, while the Solar Orbiter will go as close as 60R_⊙. Solar Orbiter will have a higher ecliptic inclination allowing for measurements closer to the solar poles.

Having instruments close to the computational domain of the PFSS model will help to evaluate the alterations made in this work significantly. Data from both Solar Orbiter as well as the Parker Solar Probe are expected to be available during the 2020s.

Badman et al. (2020) already employed a spherical PFSS model and found that a lower source surface between 1.3R_⊙ and 1.5R_⊙ matched the observations better than the traditional source surface height of 2.5R_⊙. They also pointed out that this exceptionally low source surface might compensate for washed out small-scale structures due to traditional modeling parameters, which the model alterations presented here might help to alleviate.

Acknowledgments

Wilcox Solar Observatory data used in this study was obtained via the web site http://wso.stanford.edu at 2019:10:31_09:49:01 PDT courtesy of J.T. Hoeksema. This work uses data from the Michelson Doppler Imager (MDI) onboard the Solar and Heliospheric Observatory (SOHO). SOHO is a project of international cooperation between ESA and NASA to study the Sun, from its deep core to the outer corona, and the solar wind.

References

Appendix A: Mathematical framework of the elliptical grid solver

The framework of the employed grid is not trivial, in particular the transition from the spherical configuration at the base to elliptical at the source surface. Here we briefly describe the framework employed for this work. The derivation follows the procedure laid out by Piercey (2007) for non-orthogonal curvilinear coordinate systems.

To simplify the mathematical description, we distinguish between two domains. The first domain is a basic spherical coordinate system on which the algorithm performs its operations. This domain is called the computational domain. By stretching this domain along one or two axes we obtain the physical domain, which corresponds to real-world coordinates that we are interested in. By employing this two-domain approach we circumvent the need for complicated derivations of the mathematical framework by hiding it in a very simple transformation from the computational to the physical domain. Stretching the computational grid along one axis (z) produces a prolate ellipsoidal source surface. Stretching the grid along two axes (x and y) produces an oblate ellipsoid. In the following, we concentrate on the oblate case. The prolate version is obtained in an analogous manner.

Figure A.1 shows cuts through the two domains. Coordinates in the computational domain are denoted by bars above the symbols ($\overline{x}$, $\overline{\mathit{y}}$, $\overline{z}$, $\overline{r}$, $\overline{\theta }$, $\overline{\varphi }$), whereas the same symbols without bars are used for the physical domain (x, y, z, r, θ, ϕ). They are in two groups: x, $\overline{x}$, y, $\overline{\mathit{y}}$, z and $\overline{z}$ are the well-known cartesian coordinates and r, $\overline{r}$θ, $\overline{\theta }$, ϕ, and $\overline{\varphi }$ are the spherical coordinates in their respective domains. The stretching of the computational grid is performed by an analytic stretching function along the $\overline{x}$- and $\overline{\mathit{y}}$-axes according to

$$\begin{array}{cc}\hfill x& =a\overline{x}=a_A(\overline{r})\overline{r}sin\overline{\theta }cos\overline{\varphi }\hfill \\ \hfill y& =a\overline{y}=a_A(\overline{r})\overline{r}sin\overline{\theta }sin\overline{\varphi }\hfill \\ \hfill z& =\overline{z}=\overline{r}cos\overline{\theta }\hfill \end{array}$$

Fig. A.1. Cut through the $\overline{x}\overline{\mathit{y}}/x\mathit{y}$-plane of the computational and physical domain for an oblate ellipsoidal source surface. The number of grid points shown is considerably reduced to improve clarity. The presented grid has an ellipticity of A = 2 at the source surface, which results in a source surface height of 2.5R⊙ over the poles and of 5.0R⊙ over the equator.

As the stretching function, we chose

$$\begin{array}{c}\hfill a_A(\overline{r})=1+\frac{A-1}{{\overline{r_u}}^2-{\overline{r_l}}^2}({\overline{r}}^2-{\overline{r_l}}^2)=2\alpha +a_s{\overline{r}}^2,\end{array}$$

where A is the ellipticity parameter at the source surface; $\overline{r_u}$ and $\overline{r_l}$ are the radial positions of the upper and lower computational boundaries, respectively; $\alpha =(1-a_s{\overline{r_l}}^2)/2$; and $a_s=(A-1)/({\overline{r_u}}^2-{\overline{r_l}}^2)$. We chose the squared dependence on radial distance in the computational domain to decrease the rate of change in the lower region where the algorithm requires a higher computational accuracy compared to the outer region near the source surface.

With these relations between computational and physical domain the gradient basis vectors of the physical coordinate system are

$$\begin{array}{cc}\hfill {\mathit{g}}^r& =\frac{1}{a+r·{sin}^2\theta ·\frac{\partial a}{\partial r}}\left(\begin{array}{c}sin\theta cos\varphi \\ sin\theta sin\varphi \\ acos\theta \end{array}\right)\hfill \\ \hfill {\mathit{g}}^{\theta }& =\frac{1}{r(a+r·{sin}^2\theta ·\frac{\partial a}{\partial r})}\left(\begin{array}{c}cos\theta cos\varphi \\ cos\theta sin\varphi \\ -sin\theta (a+r·\frac{\partial a}{\partial r})\end{array}\right)\hfill \\ \hfill {\mathit{g}}^{\varphi }& =\frac{1}{a·r·sin\theta }\left(\begin{array}{c}-sin\varphi \\ cos\varphi \\ 0\end{array}\right).\hfill \end{array}$$

With these vectors and a twice continuously differentiable scalar function Ψ, the Laplace operator in general curvilinear coordinates can be expressed as

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\mathrm{\Psi }=\frac{1}{\sqrt{g}}{\displaystyle \sum _i\sum _j\frac{\partial }{\partial q^i}\left(\sqrt{g}{g}^{\mathit{ij}}\frac{\partial \mathrm{\Psi }}{\partial q^j}\right),}\end{array}$$

where g = det(G) is the determinant of the metric coefficient matrix G with entries g^ij = gⁱ ⋅ g^j, and i, j ∈ {1,2,3} corresponding to the coordinates r, θ, and ϕ, respectively.

All Figures

	Fig. 1. Three regions underlying the PFSS model: region 1 is the inside of the Sun, region 2 is the computational domain between the photosphere and the source surface, and region 3 is the interplanetary space where the solar wind flows radially outward. Figure adapted from Schatten et al. (1969).
In the text

Fig. 2. Magnetic field polarity configuration at different heights. Results shown for a prolate ellipsoidal source surface with ellipticity 2.0 (left column), for the classical spherical source surface (center), and for an oblate ellipsoidal surface again with ellipticity 2.0 (right). All figures were created using our grid solver. The source surface height for all models is 2.5 R⊙ (minor half-axis in the ellipsoidal cases). Depicted is Carrington rotation 2066. Data for the lower boundary was obtained from MDI onboard SOHO. First row: three-dimensional rendering of a few magnetic field lines for each model as well as a cut through the height levels depicted below. Rows 2 to 5: magnetic field polarity configuration at height levels 1.0R⊙, 1.5R⊙, 2.0R⊙, and 2.5R⊙ (minor half-axis/radius) in a synoptic Carrington format. The height levels correspond to ellipsoids at distances of 0%, 33%, 67%, and 100% between the photosphere and source surface. Red magnetic field lines and pixels are directed inward, blue lines are directed outward, and cyan indicates closed field structures. The projection of Earth on the height levels is drawn as a black ascending line near the solar equator. Black ellipses, circles, and squares are inserted to highlight differences in the regimes of interest (see text for details).

In the text

Fig. 3. Same as in Fig. 2, but the height levels are spheres. The heights are again 1.0R⊙, 1.5R⊙, 2.0R⊙, and 2.5R⊙. The dashed line in the first row depicts the source surface, which here is not the same as the uppermost height level examined. The uppermost levels exhibit open field structures in the ellipsoidal cases because the sphere of 2.5R⊙ is partially within the source surface and touches it at the equator (prolate case) or the poles (oblate case).

In the text

	Fig. 4. Height levels show the expansion factor. Each pixel again corresponds to one field line, and the color represents its expansion factor between photosphere and source surface.
In the text

	Fig. 5. Same as Fig. 4, but for height levels as spheres.
In the text

Fig. 6. Comparison of expansion height maps for the classical PFSS implementation using spherical harmonic coefficients with order 9 in the left column, with order 20 in the middle column, and our grid solver in the right column. Again Carrington rotation 2066 is depicted. All models employ a spherical source surface at 2.5R⊙ heliocentric distance. The height levels in rows 1 to 4 are spheres at heliocentric distances of 1.0R⊙, 1.5R⊙, 2.0R⊙, and 2.5R⊙.

In the text

Fig. A.1. Cut through the $\overline{x}\overline{\mathit{y}}/x\mathit{y}$-plane of the computational and physical domain for an oblate ellipsoidal source surface. The number of grid points shown is considerably reduced to improve clarity. The presented grid has an ellipticity of A = 2 at the source surface, which results in a source surface height of 2.5R⊙ over the poles and of 5.0R⊙ over the equator.

In the text