SPICE PSF Correction: General Framework and Capability Demonstration

To begin, let us suppose SPICE observes some static (meaning it does not change over the time of observation) cube on the sky, with spectral radiance $L(x_{s},y_{s},\lambda_{s})$ as a function of sky longitude and latitude (which we will write as $x_{s},y_{s}$) and wavelength $\lambda_{s}$. These are projected onto the detector plane by the PSF, $P(x_{d},y_{d},\lambda_{d},x_{s},y_{s},\lambda_{s})$, which is a function of both sky coordinates ($x_{s},y_{s},\lambda_{s}$) and detector plane coordinates ($x_{d},y_{d},\lambda_{d}$). In reality, $x_{d}$ and $\lambda_{d}$ are time-multiplexed onto the same surface (although not necessarily the same axis) by the grating, slit, and scanning mechanism, but we will treat them as independent for clarity of demonstration. The flux density $E(x_{d},y_{d},\lambda_{d})$ on the detector plane observed for this cube is the integral of the cube against the PSF over all sky angles:

Where we have made the small angle approximation and placed the Sun at the equator of the coordinate system which allows us to set the $\cos y_{s}$ weight in the integral to one. The sensors divide the detector plane into pixels which we will index by $i,j,$ and $k$, each of which have weight functions (typically non-overlapping ‘bins’) $\theta_{ijk}(x_{d},y_{d},\lambda_{d})$. The fluxes $\Phi_{ijk}$ into each of these pixels are the integrals of the detector plane flux densities against each of these weight functions:

At this point we point out what we call the overall response function of the $\{ijk\}$th pixel of the instrument,

The name of the game is to find the spectral radiance, $L$, which is an unknown continuous function and therefore has an infinite number of degrees of freedom. To limit the number of degrees of freedom we can without loss of generality (in practice if not in principal) define $L$ in terms of a linear combination of a set of appropriately chosen basis functions, $B_{lmn}(x_{s},y_{s},\lambda_{s})$:

where $c_{lmn}$ are the coefficients of the linear combination. The basis function set ought to be linearly independent (i.e., no basis function can be expressed as a linear combination of the others). It is useful for the present purpose if each one is spatially compact and does not overlap much (or at all) with the other basis functions. They should also cover the sky plane at least as well as the pixels cover the detector plane. A set of pixel like ‘bins’ in $\{x_{s},y_{s},\lambda_{s}\}$ with equivalent (or denser) spacing to the pixels will suffice and is the choice we use in this paper. The indices $\{l,m,n\}$ are intended to reflect the identification of the bins with the coordinate axes. In terms of this, the pixel fluxes are

If we recognize an array of response coupling terms,

(the consituents of this are all known and computable, being composed of the basis functions and the in principal known instrument response functions) we are left with the following linear equation relating the unknowns ($c_{lmn}$) to the knowns ($\phi_{ijk}$):

We have used three indices each for the coefficients and the observables, but this is simply a notational convenience that allows us to identify each of the coordinate axes with its own index. Mathematically, we can ‘flatten’ all these indices down to a pair of indices (one for the coefficients, and one for the observables) by relabeling according to the following one-to-one correspondence:

Our linear equation then becomes

Which we recognize as a familiar matrix equation, with the array of coupling terms becoming a forward response matrix. Any N-dimensional linear discrete or integral transform can be similarly treated, and nonlinear ones may be amenable to a linearization and iteration approach, so this can be a very powerful and general method for dealing with these problems.

At this point, a pair of lower dimensional examples may be illustrative. Consider the spectrum for a single pixel, i.e., a function only of (for example) wavelength. This is equivalent to the DEM problem as presented by Plowman & Caspi (2020), the only difference being that instead of the instrument temperature response functions found there, we have a combination of the wavelength PSF and pixel bin functions:

The forward matrix is

This is illustrated graphically in Figures 4 and 5. The only difference between them is that in Figure 4, the individual steps are depicted, whereas Figure 5 combines all of these steps into the forward matrix – This is possible because all of the individual steps are linear operations.

Figures 6 and 7 are of the same format and show how the same treatment can be scaled up to two dimensions (spatial, in this example, although it makes no difference in the mathematical treatment). Essentially, the additional dimension in the source and the detector are simply multiplexed down to the one ‘input’ (conventionally the columns) and one ‘output’ dimension of the forward matrix. Mathematically, the only difference is the addition of a multiplexing (e.g. by ‘flattening’ in numpy or reforming to a 1D vector in IDL) step at the beginning and, if desired, a demultiplexing step at the end (e.g., by using reform in IDL or reshape in numpy) – But these are trivial 1-to-1 operations. This additional multiplexing step is shown in the top left and bottom center panels of Figure 7, which otherwise is identical in format to the one-dimensional figure (5). The output resolution of these examples differs from the input, because in general the resolution of solar sources differs than that of our detectors, to demonstrate that this framework can accomodate such resolution differences.

The SPICE forward problem for the effective ($y-\lambda$) PSF is equivalent to this 2D correction, while the full forward problem follows the same trend from 1D to 2D to 3D – the same multiplexing and demultiplexing operations suffice to convert the problem into matrix form. The practical implementation of this can be conceptually thorny, but we have developed a general N-Dimensional, coordinate aware framework in Python for computing these matrices. This will be described in more detail in Section 3.4, and the framework is also available for download (see Section 3.4).

It is clear from the examples provided above that the forward matrix is very large, even for images of modest size. However, only detector elements within PSF range of the source elements have non-zero entries, and the PSF is much smaller than the overall image. Therefore, this matrix is sparse. A wide variety of specialized sparse matrix algorithms exist to work with such matrices and they drastically reduce the real-world memory and processing requirements of working with these matrices. Standard examples can be found in, for instance, Press et al. (2002). Additionally, solving linear problems based on sparse matrices is well developed field in computer science, and by casting the problem in this form we can avail ourselves of the tools that have been developed. Our framework for computing the forward matrices is built on these algorithms, specifically making use of the scipy.sparse sparse matrix package.

Regularization typically takes the form of an additional term to minimize which is added to $\chi^{2}$. In this derivation, we will take this to be of the form

In other words, instead of minimizing $\chi^{2}$ by itself, we seek to minimize the merit function

Where $\mathbf{X}$ is a matrix specifying the regularization. This can be as simple as the identity matrix, in which case the regularization makes the solver prefer solutions with a smaller $l^{2}$ norm (i.e., ${\bf c}\cdot{\bf c}$). Exactly how much smaller depends on the size of the $\lambda$ parameter, which tunes the strength of the regularization. In Plowman & Caspi (2020), the operator was composed of the derivatives (WRT temperature, in that case) of the basis functions integrated against each of the other basis functions. This results in a regularization operator that makes the solver prefer solutions with the smaller derivatives – specifically the log of the derivative in that case.

The addition of the regularization term changes the matrix equation we wish to solve from Equation 14 to

In other words, the matrix that needs to be inverted simply has the regularization matrix times the tunable paramater ($\lambda X$) added to it.

For this class of problem, $\lambda$ is often chosen to make each the residuals of the fit to be some number $\chi_{0}$ (typically equal to the number of data points, or ‘reduced’ $\chi^{2}$ of order unity). Noting that the residuals are ${\bf b}-A\cdot{\bf c}$, and with a slight rearrangment of Equation 19, this requirement is equivalent to

We want to solve this for $\lambda$. But this is just one parameter and can’t satisfy each component of this vector equation. Instead, minimize the RMS of the LHS of this equation over lambda. The result is

But we don’t know ${\bf c}$, either. We therefore use an initial guess for ${\bf c}$ to determining our $\lambda$. We’ll need one later on when we introduce the positivity constraint, anyway.

It turns out that a simple initial guess can be formed using $F^{T}\cdot{\bf d}$. An analogy can be made to upscaling an image by an integer factor. More generally, the results of this procedure are roughly a model which has been subjected to the source response function twice. It will be much smoother, but features will be in the right place. This is desirable for regularization and it is more conservative and therefore avoids overfitting. It does result in an over-smoothed solution as well, however, the amount of which depends on how much sharper the source is than the results of the instrument resolution degradation. To mitigate this, we apply a final correction factor to $\lambda$. The SPICE sources are significantly sharper than the SPICE PSF, so the $\lambda$ used is about 0.1.

The amplitude needs to be rescaled (since $F$ is not normalized). We choose this scaling factor to be the one that minimizes the $\chi^{2}$ of this initial guess compared with the data, resulting in the following initial guess for ${\bf c}$:

The framework as currently constructed uses regularization based the $l^{2}$ norm. The derivative-based regularization used in Plowman & Caspi (2020) could also be done here, but it would require constructing gradient operator matrices with respect to image position and wavelength, which is doable but too involved to include in the present publication. It will be implemented in future work, however, which will open up a variety of possibilities for including more physical constraints on inversions using this framework, as well as some other intrinsic benefits.

One condition that the sources of these data must satisfy is positivity: spectral radiances cannot be negative. This is a nonlinear constraint, which can’t be directly captured by a linear formalism. Instead, we use the same nonlinear mapping on the coefficients as in Plowman & Caspi (2020). We express $c_{j}$ in terms of a second set of parameters, $s_{j}$, such that

We then choose the mapping function, $g$, to be one whose range is limited to the non-negative numbers. Two such examples are $e^{s_{j}}$ and $s_{j}^{2}$; we use the former (exponential) mapping here, as it seems to produce better convergence properties. We then linearize and iterate beginning from an initial guess. The results of this are described in more detail in Plowman & Caspi (2020), and we will not repeat all of the steps of the derivation here although we have changed the notation somewhat in an attempt at better clarity. The one difference is that in Plowman & Caspi (2020) only the case where the regularization acts on the $s_{j}$, rather than the $c_{j}$ was considered. If the regularization acts on the $c_{j}$ instead, there are additional terms (shown below).

The result is the following matrix equation

and

These equations suffice to ensure positivity, which in addition to being required by the physics makes the solutions more well-posed in many situations. It does require an iteration, but the convergence is generally good provided the mapping functions are monotonic (this makes the overall optimization problem convex, so the figure of merit has a single global minimum).

Since these matrices are very large, they require sparse solvers. These generally rely on some sort of iterative process in order to avoid dealing with the fully populated matrix. The code for applying these corrections is developed in Python, so we have been using the solver algorithms included in the sparse package of SciPy (Virtanen et al., 2020a). Specifically the stabilized biconjugate gradient (bicgstab van der Vorst, 1992) and lgmres (Baker et al., 2005) algorithms. Both of these algorithms have a long heritage – they are variants of algorithms mentioned in the 1992 edition of Press et al. (2002), namely the Biconjugate Gradient and Generalized Minimal Residual methods. We have tested both and found their convergence and performance to be comparable and more than acceptable for the problem at hand.

For the work shown next, we use a ‘Boilerplate’ model of the PSF which consists of a two component, 2D Gaussian where each component can be raised by an additional exponent to weight it more toward its wings or its core. We found that setting this exponent to $\sim 1.5$ made for a closer match between the IRIS line profiles (with our PSF added) and the SPICE line profiles, although the effect on the corrected Doppler shifts was very small. The mathematical expression for these Gaussians is:

where $\Delta\mathbf{r}$ is the difference between the coordinate vector of the source and the coordinate vector on the image plane, e.g., $\{x_{d}-x_{s},y_{d}-y_{s},\lambda_{d}-\lambda_{s}\}$ in equation 1. $\Sigma$ is the covariance matrix of the Gaussian, and $\gamma$ is the exponent (an exponent of one is a standard Gaussian). We specify this in terms of a set of initial axis lengths (the uncertainties of the principal axes without rotation) and rotation angles, in which case one begins with a diagonal covariance matrix that has the principal axis uncertainties on the diagonals and then rotates the matrix to its final position with a succession of rotation matrice(s) – one in the 2D case, 3 in the 3D case – Our implementation includes an example implementation of this PSF that works in both 2D and 3D.

Last, we show correction of a second part of the coordinated IRIS and SPICE observing campaign. This region is more quiet sun, with a couple of isolated bright features (The mean and brightest features are both about a factor of 4 dimmer in this region than in the first one). Removal of artifacts from these features is another good test of the correction of the PSF. A context image for these observations is shown in Figure 13, while Doppler images for the region are shown in Figure 14. The correction removes almost all the artifacts and results in a Doppler image similar to the IRIS data, although somewhat less so than before. A small relatively bright region (circled in blue) in the upper left of the region appears to still have some strong Doppler signals that are not present in IRIS, although they are reduced in magnitude. However, this region has a highly complex spectral structure and varies rapidly spatially and temporally in IRIS, so a close match between SPICE and IRIS isn’t expected here. The corrected ‘synthetic’ SPICE (IRIS degraded with noise added) do not show this artifact, so it would not seem to be an issue with the correction method itself. Nevertheless, it appears as if the ‘corrected’ SPICE data may not be fully corrected in that location. Other differences can be attributed to the lower signal in the region and possibly also to more dissimilarities between the region as seen by IRIS vs. SPICE (the intensities appear more different than in the previous region). We have also had to increase the rotation angle of the PSF from $15^{\circ}$ to $20^{\circ}$ to best remove the artifacts; the rotation angle does appear to change across the field (we tried a variety of PSF rotation angles for both region, and the ones which best reduced the Doppler artifacts are shown). The correction can account for per-pixel and/or field-wide variation in the PSF, however; it is just a (perhaps not so simple) matter of quantifying that variation.