An atlas of coronal electron density at 5Rs II: A spherical harmonic method for density reconstruction

This is the second of a series of three papers that present a methodology with the aim of creating a set of maps of the coronal density over a period of many years. This paper describes a method for reconstructing the coronal electron density based on spherical harmonics. By assuming a radial structure to the corona at the height of interest, line-of-sight integrations can be made individually on each harmonic basis prior to determining coefficients, i.e. the computationally-expensive integrations are calculated only once during initialization. This approach reduces the problem to finding the set of coefficients which best match the observed brightness using a regularized least-squares approach, and is very efficient. The method is demonstrated on synthetic data created from both a simple and an intricate coronal density model. The quality of reconstruction is found to be reasonable in the presence of noise and large gaps in the data. The method is applied to both LASCO C2 and STEREO COR2 coronagraph observations from 2009/03/20, and the results from both spacecraft compared. Future work will apply the method to large datasets.


Introduction
Reliable maps of the coronal density are important for linking various solar wind structures to the low solar atmosphere, for studies of the coronal response to the solar cycle, and for space weather applications: either as an inner boundary conditions for solar wind models, or for direct ballistic extrapolation into interplanetary space. Estimates of the coronal electron density can be made through inversion of coronal visible light observations. This has been achieved using several methods of varying complexity during eclipses, or by coronagraphs, for several decades. The introduction of Morgan (2015) gives a summary of the field, including discussion of the difficulties involved and examples of applications. A comprehensive review is given by Aschwanden (2011). This paper presents a new inversion method based on spherical harmonics for the extended inner solar corona, valid for regions where the large-scale structure is close to radial. Spherical harmonics as a basis for 3D reconstruction is used in some branches of medicine and geophysics (e.g. Merrill et al. 1996;Arridge & Schotland 2009;Levis et al. 2015, and references within). The method is described in section 2, and is tested on a simple set of synthetic data in section 3. A more complicated set of synthetic data is discussed in section 4. An approach to regularizing the higher-order spherical harmonics is presented in section 5. A discussion of datagaps, noise and temporal changes is given in section 6. Application to observations are demonstrated in section 7. Conclusions are in section 8. The appendix presents an alternative method to calculate the spherical harmonic coefficients based on iteration rather than least-squares.

Outline
For a spherical surface at a constant height r = r 0 , the coronal density ρ at Carrington longitude φ and latitude θ may be approximated by a spherical harmonic basis, where the c i are coefficients and S i are the real-valued spherical harmonics, with the i index n sph = (L + 1) 2 . By increasing the order L to large values, any sufficiently continuous density structure can be well approximated by equation 1.
If a radial coronal density structure is assumed above the height of interest, the profile f (r ≥ r 0 ) of density with height can be described by a simple function. For example, considering mass flux conservation for a spherically-expanding corona under acceleration for heights at around 5R , with α = 2.2. Thus the coronal density can be described by ρ(φ, θ, r) = ρ(φ, θ, r 0 )f (r), r ≥ r 0 For a volume segmented into discrete voxels, the observed K-coronal (electron) brightness B k is the line-of-sight summation of the product of density and a factor g which contains known constants, Thomson scattering coefficients and the length of each line-of-sight segment through each voxel (see for example section 2.1 of Quémerais & Lamy (2002), and references within): where the j index labels voxels lying along the line of sight, thus S ij is the value of the spherical harmonic at order level i and voxel j.
Each spherical harmonic S ij may be summed independently of the other harmonics along the line of sight to give the brightness contribution resulting from each harmonic.
Defining A i : the total brightness is given by This describes a linear relationship between the contribution from each spherical harmonic density distribution and the observed brightness. For the purpose of finding an unknown density distribution from observed brightness, a reconstruction space with prescribed S ij , f (r j ) and g j is created. The line of sight summations of equation 5 are calculated, and the problem is reduced to finding the coefficients c i -thus the line-of-sight integrations are made only once, leading to high efficiency. Given a large number of observations (n obs n sph ), the system is overdetermined and can be solved using least squares. The ability to perform the line of sight summations independently for each spherical harmonic is based on the assumption of a radially-structured corona at heights above the height of interest, and a uniform profile to the decrease in density with height (e.g. equation 2). The assumption of a radial corona is reasonable at r =5R , and the approximation of an assumed radial density profile is discussed later.

Application
Consider a set of observed coronal images recording brightness B k , taken over an extended time period (e.g. half a solar rotation, ∼2 weeks). Circular samples of data at constant distance from Sun center, at a height at which the coronal structure is deemed close to radial (e.g. 5R ), are extracted over the time period, giving b, which records B k as a function of position angle and time. For each member of b, a geometrical line-of-sight is defined through the corona, extending to large heights behind and in front of the point of closest approach to the Sun (similar to the description in the following section for the creation of synthetic observations). A set of S ij , g j , and f (r) are prepared (with the unknown f (r) set according to equation 2). The line-of-sight summation of equation 5 is then implemented. This gives a set A i , one for each spherical harmonic, each of size n obs .
Assuming a normal distribution to observational errors, the problem is reduced to solving min c |b − Ac| 2 , with matrix A of size n sph × n obs , b of size n obs and c the coefficients of size n sph . The least-squares solution to equation 7 is For numerical stability, before applying equation 8, A and b are divided by the mean of the absolute values of A (both contain very small numbers).

A simple test
Synthetic observations are made from a known density distribution. For this example, a spherical distribution of density at height 5R is created using equation 1, with L = 11 (n sph = 144). The c i are created from a set of random numbers in the range −1 to 1, divided by weight l + m + 1, so that higher-order components are reduced in amplitude.
The distribution is then scaled between a minimum at a typical value for electron density in a coronal hole (Doyle et al. 1999), and a maximum within a streamer (Gibson et al. 2003).
This distribution is shown in figure 1a. This will be the target density distribution against which the method is tested. The distribution is simple in the sense that it is based directly on spherical harmonics -it is not similar to a true coronal density distribution, yet it serves as an initial test of the method.  et al. (1999) to fix the minimum density at each height. Similarly, the formulation of Gibson et al. (2003) can be used to set the maximum density at each height. The g i are then calculated, and the resulting emission summed along each line of sight. The 'observed' K-coronal brightness b, as a function of position angle and time, is shown in figure 2.
An important choice in reconstructing the density is the choice of L, or the maximum number of orders. For the sake of this first simple test, this is set at L = 11, to match the order of the input distribution. Solving equation 7 takes a few seconds on a 2.8GHz Intel Core i7 desktop computer with 16Gb memory. The reconstructed density map is shown in figure 1b. The percentage difference between target (ρ t ) and reconstructed (ρ r ) density is shown in figure 1c. The mean absolute percentage deviation is 3.8%, whilst the distribution correlation C over the sphere, given by is 99.8% (theρ are means). Figure  The algorithm is close to giving a perfect reconstruction for this simple test case. This is perhaps not surprising given that the test density is based directly on spherical harmonics, and that the information on the number of orders (L = 11) has been used for the solution.
Note that the original density distribution used to create the synthetic observations has a density decrease with height based on the formulation of Doyle et al. (1999) and Gibson et al. (2003). This gives a decrease with height which is proportional to the relative density of each point, but which does not follow the spherically uniform decrease of equation 2. For the reconstruction, the true decrease is assumed unknown, and equation 2 is used. It is obvious from the success of the reconstruction that this leads to only a minor error.
The Appendix describes an alternative method for finding the coefficients c, based on the properties of the spherical harmonics and iteration. The alternative method performs well in the case where the target density is directly based on spherical harmonics. In general, and for the rest of this work, it is not used since its performance degrades (in both accuracy and efficiency) in comparison to the least-squares method on more complicated density distributions. It is included in the Appendix since it is an interesting approach and may prove useful in other contexts.

A more realistic test
In this section, a complicated, narrowly-peaked, density distribution is used to test the reconstruction method. In contrast to the previous simple test, the density distribution is not based directly on a spherical harmonic basis, and therefore the distribution cannot be exactly fitted by a limited order of spherical harmonics, and the number of orders required in the reconstruction cannot be determined beforehand. This distribution is where ρ 1 , ρ 2 are summed spherical harmonic series with weighted random coefficients (as in the simple case of the preceding section), with L = 11 and M = 9, and with ρ 1 scaled -10between 0 and 1. The exponential term forms ridges centered on where the ρ 2 function passes through zero, and these ridges can be made narrow by setting ω to a small value.
The ρ 1 term introduces variability to the value of both the ridges and the background. This initial density distribution is scaled to appropriate coronal values of density in a similar way to the simple case above. The resulting density distribution is shown in figure 5a. Through the exponential function, this distribution has extended, narrow and intricate structures, and is more similar to the expected form of the true coronal density distribution, being distributed along narrow sheets along polarity inversion regions and pseudostreamers (e.g. Despite the decent structural correlation in density distribution, and the almost identical match between model and observed brightness, the reconstruction suffers from high-frequency longitudinal oscillations, leading to large inaccuracy near the equator and regions of low density (including a small negative region). These oscillations are caused by large spherical harmonic coefficient values at high frequencies as the data is overfitted. Figure 7 shows the optimal density that can be achieved using a 25 th order spherical -11harmonic basis. The coefficients are calculated directly from integrating the product of each spherical harmonic basis with the true input density over the spherical shell by

Morgan
Steep jumps in density cause high-frequency oscillations (Gibbs oscillations), which can be seen in figure 7b. These are minor compared to the large-amplitude errors in the least-squares tomographical reconstruction.
The tendency of the reconstruction to contain negative densities near high-density regions is a problem which plagues coronal tomography. This test shows that it is a problem which arises not solely due to rapid temporal changes in the streamer belt or due to contamination by coronal mass ejections (this test data has zero noise and no temporal changes). It is a problem intrinsic to the observations -of convolution of linear lines of sight through an extended spherical structure, and related to missing information at heights below the height of interest r 0 for any single observation. Even for tomography at heights below 5R , this problem is unavoidable at the limit of the instrument field of view. The problem of extreme oscillations in reconstructed density is worse near the equator: for a given observation, the LOS integrations at the equator pass through only a limited range of longitude and through only a very small range of latitude. At the poles, the LOS observations pass through the whole polar corona, near to the axis of rotation, giving a more stable reconstruction.
The results of this section show that some form of regularization is required to impose smoothness on the reconstruction and to avoid negative densities.

Regularisation of the higher-order harmonics
Other coronal tomography methods impose a condition on the spatial smoothness of the reconstruction (e.g. Frazin 2000) to avoid unphysical high-frequency components. A similar and necessary extension of the spherical harmonic approach is given here. It is desirable to increase the highest order of the spherical harmonics in order to reconstruct the density structure at the finest possible resolution, yet this leads to greater instability of the highest orders. Coronal tomography methods achieve stability by imposing a weighted penalty term for lack of spatial smoothness in the reconstructed density -thus the optimal reconstruction is given by a compromise between the best fit to the data and the spatial smoothness of the reconstruction (regularization).
The noise σ at each position angle and time bin is estimated from the original pre-binned data by isolating the highest-frequency spatial and temporal component. To where λ is a regularisation factor that sets the balance between fitting the data and imposing a priori constraints on the solution. w is a square matrix, with diagonal elements i = 0, 1, ..., n sph − 1 given by and non-diagonal elements are zero (the l and m are the spherical harmonic longitudinal and latitudinal order). w takes the place of the more commonly-used identity matrix so that the regularisation has a larger direct impact on higher frequency harmonics.
In previous work on regularization in coronal tomography, the commonly-used positivity constraint on the density selects values of λ where density is everywhere zero or positive. From our own tests on this approach, this gives an overly-smooth solution -that is, for all small values of λ the positivity constraint is not satisfied, and only at large values does the density become everywhere positive. A different approach is taken here. Our fitting routine finds an optimal solution using two parameters. One is λ (the smoothing parameter), and the other is a minimum density threshold ρ . The main steps in this approach are: 1. Values λ k , with index k = 0, 1, ..., n k − 1 are set by a logarithmic increment between the minimum entry of the diagonal of the co-variance matrix A σ A σ divided by 10, and the maximum entry multiplied by 2. Typically we set n k = 25.

A minimum density is estimated from the observed brightness values through a
spherically-symmetric inversion of the 2 nd percentile minimum of brightness. Values of ρ j , with index j = 0, 1, ..., n j − 1 are set between the minimum density divided by 5, and the minimum density multiplied by 2. Typically we set n j = 20.
3. For each value of λ k , an initial solution is given by equation 12. This solution gives an initial density distribution on a longitude-latitude map at the coronal height of interest (e.g. 5R ).
4. For each value of ρ j the initial reconstruction solution at the current λ k is thresholded to a minimum value of ρ j . A new set of spherical harmonic coefficients are calculated directly from this thresholded density map via equation 11. These adjusted coefficients c k,j are used to give a measure of goodness-of-fit to data for the current value of λ and ρ by: Thus a 2D array χ k,j is gained that maps the goodness of fit as a function of λ and ρ . The final task is to define an optimal point within this array. Figure 8 shows χ k,j for the complicated density distribution, calculated over a grid of 25 λ and 20 ρ points.
As expected, χ increases with increasing λ -a smoother density reconstruction gives a poorer fit to data. χ also increases with increasing ρ , since the reconstructed density is 6. Missing data, noise and rapid temporal changes Figure 10a shows the brightness test data degraded through the addition of random normally-distributed noise at 5% of the mean signal level. Regularized tomography applied to this noisy dataset gives the density of figure 11a. The reconstructed density has a mean absolute fractional deviation of 12.1% from the target, with a structural correlation of 95%.
The brightness values are fitted with a mean absolute deviation of 4.3%. The solution has a minimum density of ρ = 9.96 × 10 3 cm −3 , and λ = 1.75 × 10 3 .
The largest reconstruction errors are near the equator, where high-density regions are underestimated, and low-density regions overestimated -that is, the reconstruction gives density which is too smooth over longitude compared to the sharply-defined structures and large gradients of the true density. This is an important point to remember when interpreting tomography results applied to real data -the equatorial regions are the most important regions in the context of space weather studies, yet is where the reconstruction errors are greatest.
All coronagraphs suffer from occasional datagaps, with the potential to seriously degrade tomographical reconstructions. Figure 10b shows a half-solar-rotation set of noisy synthetic observations with 4 missing days of data (around one-third missing) split into 3 gaps of 2 days, 1 day and 1 day. The reconstructed density for this data is shown in figure 11b. It deviates from the target density by 14.1%, with a spatial correlation of 94%.
The reconstructed and observed brightness deviate by 4.3%. The solution has a minimum density of ρ = 9.9 × 10 3 cm −3 , and λ = 2.19 × 10 3 . Thus the spherical harmonic basis provides stability in the presence of even quite substantial datagaps.
The most detrimental noise in coronagraph data is perhaps not a normal distribution, but rather isolated pixels or groups of pixels of spurious high/low values caused by, for example, sporadic bursts of energetic particles which can seriously deteriorate some images, or the passage of bright planets. The weighted fitting can help reduce the impact of these on the results. More importantly, rapid changes in brightness and structure caused by CMEs have a large detrimental effect on reconstruction. Paper I introduces several processing steps to reduce these problems. In particular, the dynamic separation technique (DST) reduces the effect of CMEs, and also results in a smoother signal with reduced salt-and-pepper noise. Observations which are seriously degraded (possibly due to bursts of energetic particles), can be identified and discarded, as described in Paper I. Occasionally, telemetry or read errors can lead to missing blocks of data within an image. Discarding bad images, or missing data blocks, will result in short datagaps, which seems acceptable for the spherical harmonic method as shown above.
Lastly, coronal structure must change, either slowly, or rapidly, and may reconfigure very rapidly during the passage of large CMEs. Time-dependent coronal tomography (based on regularisation methods) has been successfully applied by Vibert et al. (2016). In principle, the spherical harmonic approach can be extended to include time-dependency, with the coefficients becoming functions of time. Initial experiments with a time-dependent density model shows that this is a very challenging task -particularly if a step-change in density is needed to account for rapid changes. Further development is necessary, reserved for a future publication.

Application to observations
This section applies the tomography to observations made by the LASCO C2 and the STEREO SECCHI COR2 A coronagraphs for a half-Carrington rotation period centered on 2009/03/20 12:00. At this time, the STEREO A spacecraft is separated by 60 • from SOHO.
The data are processed and calibrated according to the method of Paper I. The height of interest is set at 5.5R , and the data rebinned into a position-angle and time array with 180 position angle bins, 200 time steps. The data array for LASCO C2 is shown in figure   12a, and for COR2 A in figure 12c. The data binning can be set at higher resolution, at the expense of computational time. The binning here allows reconstructions to be made in approximately 10 minutes.
The choice of period, and height, is to allow convenient comparison with figure 5 of Frazin et al. (2010). The density reconstruction for LASCO C2 is shown in figure   13a, and for COR2 A in figure 13b. The LASCO C2 data is fitted with a mean absolute deviation of 10.6%, with a smoothing parameter of λ = 6.2 × 10 4 and minimum density ρ min = 1.4 × 10 3 cm −3 . For COR2 A the values are 7.0%, λ = 5.1 × 10 4 and ρ min = 6.5 × 10 3 cm −3 . The mean absolute fractional difference between the two reconstructed densities is 38%, with a spatial correlation of 81%. Comparing with figure   5 of Frazin et al. (2010), these density maps are smoother, and have maximum densities at around half the values of Frazin et al. (2010). Currently there is no other empirical verification for density maps such as these. From figure 13, COR2 A seems to give a better reconstruction, in that the streamer belt is narrow, and is fitting the data more closely.
Comparison with future in situ measurements of the coronal density by the Parker Solar Probe will be invaluable for coronal tomography.

Conclusions and future work
For heights where the coronal structure can be well-approximated as radial with an uniform density decrease with increasing height (i.e. the extended inner corona), a model of the density based on spherical harmonics leads to a very efficient and stable method for reconstruction. This is demonstrated for a simple and complex model coronal density distribution. The method is robust to large datagaps of several days. Without regularisation, the smoothness of the reconstructed density is dictated by the highest order of the spherical harmonic basis. However, the true coronal density is likely to have steep gradients between regions of low and high density, or very narrow regions of high density, and a high order is required to approximate these. To counteract this problem, we provide a method for regularised solutions where the smoothness of the reconstructed density, and a minimum density threshold, is taken into consideration.
The application of this method to a large dataset will be presented in the third paper of this series. Other future work involves finding a robust time-dependent extension to the spherical harmonic approach, where the harmonic coefficients can change as a function of time. We also aim to experiment with other approaches similar to spherical harmonics that have proved useful in geophysics, including wavelet-based spherical functions (Chambodut et al. 2005). We anticipate these may prove useful for the non-radial corona, in particular for extreme ultraviolet (EUV) observations of the low corona.
Spherical harmonics are a simple yet powerful basis for inversion of coronal density, and should be a consideration for other coronal applications such as EUV diagnostics in the low corona, or 3D reconstructions of the coronal magnetic field with future spectropolarimetric instruments.

Appendix
This appendix describes an iterative procedure to find the spherical harmonic coefficients c i . For this procedure, the observed data b and the A i (see equation 5 of section 2.1) are first normalized to achieve numerical stability -both are very small numbers (b and |A i | on the order of 10 −10 and 10 −16 respectively). b is normalized to a mean of zero and unity standard deviation by whereb is the mean and σ b is the standard deviation. The A i are normalized by the mean of their absolute value (calculated over all orders): Starting with an initial estimate of coefficients (labelled with a prime, c i , since they are operating on normalized arrays) all set to zero, the following iterative algorithm, with iteration counter k cycling through equations 3-5, converges towards a solution: where σ b mod is the standard deviation of b mod and λ ( 1) is a parameter that controls the rate of convergence. At values too large, the process does not converge. This becomes important as the number of spherical harmonic orders becomes high. λ = 1 n sph gives good results for the examples in this work (where n sph is the number of spherical harmonics).
The iterations continue until k reaches a set value, or the convergence rate drops below a set threshold. Note that equation 5 is not strictly necessary, it is included to greatly increase the rate of convergence.
After convergence is reached, the c i are scaled to account for the normalizations of equations 1 and 2, to give solution c i : where σ b and σ b mod are the standard deviations of the observed and modelled brightness.
Finally, the mean density which should be included in the zeroth-order DC component, c 0 , is estimated directly from the observed brightness by: C is a correction factor based on the curtailing of the line of sight to a limited range. Due to the curtailing, the summation in the denominator is too small, leading to an overestimate of the mean density by a few percent. This correction is easily quantified by calculating n los j=1 g j f (r j ) for a single case of a very long line of sight (where the emission essentially drops to zero at large heights), and comparing the same value for the curtailed line of sight.
This gives the correction factor C directly.
To fit any function on a sphere to a set of spherical harmonics, the coefficient of a spherical harmonic at a given order can be found by integrating the product of the function and the spherical harmonic over the sphere (see equation 11). In this case, where the spherical harmonics are multiplied by geometrical and other factors and integrated over extended lines of sight, the iterative algorithm of equations 3-5 in essence implements a similar approach. For the simple test case of section 3, this iterative method gains a more accurate reconstruction than the least-squares approach, with a mean absolute fractional deviation of 2% between the reconstructed and target densities. For the more complicated cases, it loses accuracy compared to the least-squares approach, and with increasing number of spherical harmonic orders, it becomes considerably less efficient.