Multi-component Decomposition of Astronomical Spectra by Compressed Sensing

The signal measured by an astronomical spectrometer may be due to radiation from a multi-component mixture of plasmas with a range of physical properties (e.g. temperature, Doppler velocity). Confusion between multiple components may be exacerbated if the spectrometer sensor is illuminated by overlapping spectra dispersed from different slits, with each slit being exposed to radiation from a different portion of an extended astrophysical object. We use a compressed sensing method to robustly retrieve the different components. This method can be adopted for a variety of spectrometer configurations, including single-slit, multi-slit (e.g., the proposed MUlti-slit Solar Explorer mission; MUSE) and slot spectrometers (which produce overlappograms).


Introduction
From solar flares to quasars, spectrometers are used to investigate a wide variety of astrophysical phenomena. To obtain measurements of the physical conditions of the emitting material (e.g., extreme UV emission lines of many-times ionized Fe atoms) inevitably requires a forward model with the following components: 1. P-A physics-based model of the radiative process operating in the astrophysical object of interest (e.g., emission, absorption, scattering, and gravitational redshift); 2. O-An optical model of the telescope (including the point-spread function (PSF), instrumental spectral broadening, and if the telescope is ground-based, atmospheric seeing); and 3. D-A detector model capturing the properties of the sensing system (e.g., nonlinearity, dark current, gain patterns, and sources of noise).
The goal of spectroscopic measurements is to provide observational constraints of the physical properties of the system. For the sake of discussion, suppose one has perfect knowledge of the optical system and the detector. Even then, a physics model is still required for the most basic of spectroscopic measurements, such as the Doppler shift of a spectral line. For instance, consider spectroscopic measurements of extreme UV (EUV) emission lines from solar coronal plasma in the optically thin regime in the absence of scattering. In order to extract the Doppler shift of the line, the local enhancement (emission line) or deficit (absorption line) in the detected spectrum must first be associated with a known spectral line from a certain atomic species. This provides a reference rest wavelength of the line against which a Doppler shift (in wavelength and in velocity) can be measured. Accounting for, or in the absence of, gravitational redshift, the Doppler shift informs us about the motion of plasma along the line of sight (LOS), v LOS . Given an atomic model(e.g., CHIANTI, Dere et al. 1997;Young & Landi 2009;Landi et al. 2013) we can also attribute the emission line to plasma in a certain range of temperatures. Assuming thermal equilibrium conditions (e.g., thermal collisional excitation rates balanced by spontaneous radiative deexcitation), the atomic model also provides the thermal width σ th of the line. In the general case, the spectrum measured at the detector may be due to a heterogeneous mixture of plasmas at different temperatures and Doppler velocities. For example, measurement of a spectral line with an observed width σ obs >σ th suggests the emitting plasma has multiple components moving at different LOS velocities (which, depending on the physics model, may be interpreted as a sign of turbulence; Brooks & Warren 2016;van Ballegooijen et al. 2017;Polito et al. 2018). The detection of multiple spectral lines associated with different temperatures suggests multiple thermal components in the emitting plasma. The presence of multiple components contributing to a single spectrum on the detector may be due to spatial inhomogeneities along the LOS or within the plane-ofsky area visible to the slit (or multiple slits in the case of a multi-slit instrument).
The aim of this work is to describe a method for decomposition of the spectra into constituent components of the emitting spectra by techniques of compressed sensing. The driver behind this work is to address the complexities introduced by the use of multi-slit instrumentation in solar physics (e.g., the Multi-slit Solar Explorer or MUSE, proposed as a NASA Heliophysics mission), which promises to greatly increase the field of view and/or drastically improve the temporal cadence of spectroscopic data. In this paper, the method is demonstrated in the context of MUSE, for which the goal is to account for any effect of blends on the primary lines of interest, rather than to use the decomposition results directly. However, the method described is very general and can be applied to a host of other astrophysical problems. The article is structured as follows. Section 2 gives a mathematical description of the problem for the general case. Section 3 briefly discusses the parameter space we explore. Section 4 describes in some detail the compressed sensing approach, while Section 5 discusses some examples of application of the method to solar spectral data. Finally, in Section 6 we discuss the presented method and its results.

Problem Statement
Consider a multi-slit sensing system with a set of parallel 0, , 1 m { } , each of width w, and a common detector at the focal plane (see Figure 1). Assume the regular slit spacing , where λ is the wavelength of observation and f/A is the f-ratio of the light beam at the slit ( f and A being the effective focal length and aperture diameter, respectively, of the telescope that feeds the spectrograph). The detector has a 2D array of pixels indexed [i, j], and we assume that i corresponds precisely to the spectral dimension and j corresponds precisely to the spatial dimension along the slit, as can be accomplished through geometric corrections when instrument alignment is not perfect. Variations along i and j are therefore separable, and we need only to consider variations in i for the decomposition. The measured spectrogram y[i] (intensity at the ith pixel) can have contributions from photons originating from multiple slits.
is the spectrogram due solely to radiation from slit m. The emitting material seen by slit m has a density function Ψ over some parameter space of physical properties (e.g., temperature, density, and Doppler velocity) and radiates a spectrum denoted I λ =P(Ψ). The operator O m acting on P(Ψ) represents how the radiation is processed by the optical system, including how light through an individual slit is dispersed and focused onto the detector to form a spectrum. For two adjacent identical slits that happen to observe the same packet of material residing in the astrophysical object of interest with physical property Ψ would lead to the same spectrum O P ( (Ψ)), Figure 1. Schematic of a multi-slit spectral sensing system: multi-wavelength radiation I λ emitted by an extended astrophysical object (e.g., the coronal plasma of a solar active region) is focused by a telescope onto a system of M parallel slits S that transmits light sampling a picket fence subset of the overall field of view. A disperser optical element (e.g., a diffraction grating), or system of elements, then disperses the transmitted radiation as a spectrum and focuses it on the detector consisting of an array of N pixels. The measured spectrogram is the superposition of spectra originating from all slits, and includes contributions from the detector gain G and noise e. Accurate interpretation of the spectrogram requires correct identification of the spectral lines (the same lines represented by the same color) and their slit origin, a problem addressed in Section 4. except Y + O P m 1 ( ( )) would be translated from Y + O P m 1 ( ( )) by some pixel offset Δ i, namely For an ideal detector, the operator D gives an identity mapping (i.e., it preserves the spectrograph exactly). A noisy, linear detector can be described as , where G is the gain and e is the stochastic (and perhaps systematic) noise introduced by the detector.
In general, different slits will be exposed to incident radiation from different parts of the astrophysical scenery. So the challenge is to decompose the net measured spectrum y[i] into constituent components by identifying 1. the slit(s) that contributed to the net spectrogram, and 2. the physical properties of the radiating material along the column of integration in the LOS seen by each slit.
Though the above description of the problem applies to both radiation from optically thick and optically thin plasmas, the rest of this paper will be concerned with the simple case of optically thin plasmas. This allows us to express the physical model P as a linear operator over the parameter space density function Ψ. The aim of measuring the physical properties of emitting material is equivalent to finding the function Ψ such that Equation (1) is satisfied. The following sections (in particular Section 5) will provide some concrete examples.

Differential Emission Measure (DEM)
A common concept encountered in EUV and X-ray observations of solar plasma is the DEM(see, e.g., Boerner et al. 2012), defined by the relation where y[i] is the detected spectrogram (after dark subtraction, flatfielding, correction for nonlinearities, etc.), K i (T) is the temperature response function of the ith spectral channel, and DEM is the electron density squared, integrated along the LOS, contained in a temperature bin of width dT. K i (T) encapsulates assumptions about the radiative properties of the emitting plasma (i.e., P) and the optical properties of the system (i.e., O), such as PSF, effective area, etc.
The aim of the DEM inversion problem in previous work was to recover DEM(T) given a set of measurements y[i] (see Cheung et al. 2015, and references therein). In that context, DEM(T) is the parameter space density function Ψ, which spans over the temperature range [T 0 , T 1 ], but did not include the dimensions of Doppler velocity, as we assumed that the observing system in question(e.g., the EUV channels of the Atmospheric Imaging Assembly (AIA), see Boerner et al. 2012;Lemen et al. 2012) does not have sufficient spectral resolution to resolve Doppler shifts. Furthermore, the physical model of radiation assumes the emissivity of EUV spectral lines is only a function of T (though as shown by Martínez-Sykora et al. 2011;Testa et al. 2012, there is also a slight dependence on plasma density).

Velocity DEM (VDEM)
As an extension of the DEM inversion problem, consider the situation where the sensing system has sufficient spectral resolution and sampling to be sensitive to Doppler shifts. In this case the parameter space density function Ψ(T, v) spans temperature and Doppler velocity space. Solving the VDEM problem means we seek to quantify how much plasma (in terms of emission measure n e 2 ) is at a certain temperature T moving with Doppler velocity v. In other words, find Ψ(T, v) such that the following relation holds: In the case of multiple slits, the rhs of Equation (4) also includes a sum over m (the index for the M slits) and K i,m may be slit-dependent. In Section 4, we describe how a compressed sensing technique is used for the class of problems similar to Equation (4).

Compressed Sensing Solution Approach
Suppose Ψ0 is a density function over an n-dimen- [ ], and f j,0 and f j,1 are the lower and upper bounds of the jth parameter. The operator P acts on f to generate (and propagate) radiation to the sensing system. The operator O takes this incident radiation arriving at the sensing system (e.g., a telescope and its optical system) and produces the M-tuple y (each component of y is the spectrogram measured in a pixel of the detector). P and O,and the deterministic part of D (e.g., gain) are assumed to be known.
The solution strategy begins with generating response functions. Consider each dimension of parameter space is discretized into a finite number of points N j . For each point f in parameter space, compute the detector response f r according to Equation (1). The set of all response vectors are used to generate the response matrix 1 . Parameter estimation (i.e., measuring the physical properties of the radiating plasma) is then equivalent to solving the following linear system for x is an N-tuple of coefficients for the response functions. In other words, we seek to express the measured spectrogram y as a linear superposition of response functions f r . For a multi-dimensional parameter space, this linear system may be underdetermined. There are a variety of compressed sensing schemes that can be used to tackle such types of problems. For instance, for DEM reconstruction from narrowband EUV data taken by the AIA (Boerner et al. 2012;Lemen et al. 2012) on board NASA's Solar Dynamics Observatory(SDO; Pesnell et al. 2012), Cheung et al. (2015) presented a validated inversion scheme based on basis pursuit (Chen et al. 2001). That particular problem has M=6 (six EUV channels used) and N≈20 (number of T log bins). For much larger problem sizes, we found the lasso method (Tibshirani 1996) implemented in the Python scikit learn module (Pedregosa et al. 2011) to give reliable results. Lasso seeks a solution for Equation (5) by finding: In other words, # x is the argument x which minimizes the objective function in the square brackets. By adding an L1 penalty term, Lasso promotes sparsity in the solution.α is a hyperparameter (i.e., a parameter that is not fitted by the algorithm) used to control the level of sparsity. A larger value of α tends to yield solutions that have smaller L1 norm (more sparse).

Single-slit Spectrometer: Hinode/EIS
Perhaps the most common type of astronomical spectrometer instruments are those with a single slit,i.e.,M=1. An example is the EUV Imaging Spectrometer(EIS; Culhane et al. 2007) on board the Hinode mission (Kosugi et al. 2007). Hinode/EIS is sensitive to emission lines of many-times ionized metallic species (e.g., Fe, O, and Ni) found in solar coronal plasmas at temperatures betweenT log K 5.0 andT log K 7.3.
In general, an EIS spectrogram consists of multiple emission lines corresponding to different plasma components of various temperatures and Doppler velocities. Sample spectrograms are shown in Figure  To invert the spectrograms for VDEMs, we follow the procedure outlined in Section 4. Consider a unit of emission measure (e.g., for observations of solar plasmas at 1 au, an emission measure of = EM 10 0 25 cm −5 or above is generally detectable given a sufficiently bright emission line and sufficient effective area and exposure time). Consider a VDEM distribution Ψ(T, v)=EM 0 δ(T 0 , v 0 ), i.e., an isothermal plasma at temperature T 0 at a single LOS velocity v 0 . Using the physical model P (CHIANTI) and instrumental response model O ( ), compute the response vector r T v , 0 0 ( ) for this VDEM distribution, and repeat for other values of T 0 and v 0 in VDEM space. This allows us to construct the response matrix R used for VDEM inversions. The parameter space has a velocity range of ±400 km s −1 with a velocity bin size of 10 km s −1 and the temperature ranges from Although there are imperfections in the recovered VDEM, they reproduce the salient features in the ground truth VDEMs. Over a range of temperatures, the inversion correctly reproduces the spread of emission measure in Doppler velocity space. For instance, the range of Doppler velocities with significant emission measure is much narrower for We also tested our VDEM inversion method on actual Hinode/EIS observations. Even though in this case the ground truth is unknown, we have confirmed that the results we obtain are in good agreement with findings of previous detailed studies. We analyzed EIS spectral data of AR 11190, observed on 2011 April 15 01:17:19. This data set has been analyzed by Warren et al. (2012; their selected region number 13). We used for the inversion the spectral lines listed in Table 1, which provide a good temperature coverage fromT log K 5.8 [ ] to ∼6.7. In order to include the Ca XVII line, we also included the strongest known blends (from O V and Fe XI; see Warren et al. 2012 and reference therein). The response function is constructed for a velocity range of±200 km s −1 (with a velocity bin size of 5 km s −1 ) and a temperature range ), by using CHIANTI v.8 emissivity functions (assuming CHIANTI ionization equilibrium, and coronal abundances). We calibrate the Hinode/EIS spectra with standard routines available in SolarSoft to remove the CCD dark current, cosmic-ray strikes on the CCD, and take into account hot, warm, and dusty pixels. In addition, the radiometric calibration is applied to convert the data from photon events to physical units, and we correct for the CCD offset, and for the orbital variations of wavelength calibration. Residual shifts in wavelength calibration are corrected by fitting an Fe VIII line (in the short-wavelength channel) and a Si VII line (in the longwavelength channel), in a patch of quiet Sun in the field of view, and assuming these velocities to be zero(see, e.g., Young et al. 2012).
The results of our VDEM inversions applied to these EIS data are shown in Figure 3. In particular we show two sets of panels (one for "box 1" on the left, and one for "box 2" on the right; see  [ ] bin). We also note that Warren et al. (2012) included further constraints to the high temperature end of the DEM by using the SDO/AIA 94 Å observations (estimating their Fe XVIII contribution). Such constraints were not applied in our current study.

Multi-slit Spectrometer
The MUSE is a science mission proposed for NASA's Heliophysics program (Tarbell & De Pontieu 2017). Like Hinode/EIS, it measures atomic EUV lines emitted by coronal plasma. However, by using 37 spectrally dispersing slits (i.e., M = 37), it allows spatial rasters of solar active regions up to two orders of magnitude faster than EIS or any other existing or planned spectrometers. This design allows the MUSE instrument to capture many more solar eruptions and flares, and, for the first time, to capture them with sufficient spatio-temporal resolution to reveal the dynamic evolution of the active corona. The extremely rapid (12 s cadence), subarcsecond (0 4) resolution rasters (170×170 arcsec 2 ) with broad temperature coverage, accompanied by large FOV context imaging in several EUV lines (Fe XII 195 Å and He II 304 Å) will allow MUSE to address its top-level science goals: (1) to determine which mechanisms drive coronal heating and the solar wind, (2) to understand the genesis and evolution of the unstable solar atmosphere, and (3) to investigate fundamental physical processes in the solar corona. The three spectral passbands are dominated by spectral lines with wavelengths around 108 Å (Fe XIX, Fe XXI), 171 Å (Fe IX), and 284 Å (Fe XV). These lines are formed around T log K [ ]=7.0, 7.1, 5.7, and 6.4, respectively. Because the passbands are spectrally wider than the (wavelength) separation between neighboring slits, the multi-slit design can, in principle, lead to overlap of spectral information from neighboring slits. This is minimized by: (1) the selection of band passes to study bright, well-isolated lines as primary diagnostics, (2) the selection of a slit spacing that minimizes possible blends from other slits. This typically limits multi-slit confusion to regions in which the primary lines are not bright, or where the plasma has unusual emission measure distributions (e.g., a predominance of very cool plasma, e.g., in coronal loop fans, which can lead to contamination by secondary lines). Our spectral decomposition code has been shown to be very effective in disambiguating the multi-slit confusion even for these difficult conditions (as shown below, and in more detail, in a follow-up paper focusing on MUSE applications).
To satisfy the science requirements for MUSE, it is not necessary to accurately determine a VDEM distribution for each slit S: Y T v S , , ( ). Instead this VDEMS distribution is only used as an intermediate step to disambiguate any multi-slit confusion of the primary lines, e.g., by calculating, for each slit, the primary lines and secondary lines from the VDEMS distribution. Similar to the example in Section 5.1, we use the physical model P (CHIANTI) and instrumental response model (O, which takes into account the position of the 37 slits in the spectrogram) to compute the response vector r T v S , , 0 0 0 ( )for this VDEMS distribution, and repeat for other values of T 0 and v 0 in VDEM space and the 37 slits S 0 . This allows us to construct the response matrix R used for VDEMS inversions. The matrix R has a velocity range of±400 km s −1 with a velocity bin size of 10 km s −1 , a temperature range from with a bin size of 0.2, and is calculated for all 37 slits. For the detector model D, we added photon noise based on a Poisson distribution, with exposure times of 1.5 s and the MUSE effective area.
Using the same snapshot of the 3D radiative MHD simulated solar flare from the previous section, we synthesized the optically thin MUSE spectrum using the CHIANTI database for all three passbands (108, 171, and 284 Å, i.e., N=600). The synthetic MUSE observations take into account the instrument response for the 37 slits, the spectral and spatial resolution, and all spectral lines from the CHIANTI database with wavelengths that could fall on the detector from any slit. Specifically, spectra were generated using all lines from CHIANTI v8.0.7, with updated data for Fe VII. We use the spectral decomposition code on the synthetic MUSE data from all three passbands to determine the VDEMS. As mentioned, this inverted VDEMS is not the end goal, but instead can be used to reconstruct the dominant spectral lines without contribution from other nondominant spectral lines and removing the confusion from adjacent slits. The example in Figure 4 compares the synthetic Fe IX171 Å intensity map from the 3D radiative MHD flare simulation (ground truth, panel (A)) with a disambiguated map from the inverted VDEMS (based on synthetic data that includes photon noise, panel (B)). The correlation between the ground truth and intensities derived from the inversions are very good (both for the cases without photon noise, panel (C), and with photon noise, panel (D)), illustrating that the spectral decomposition code accurately disambiguates the MUSE data, even for the challenging scene presented by a flare (in which secondary lines and strong Doppler shifts can potentially lead to multi-slit confusion). Note that this is only shown as an illustration of how the decomposition code can disambiguate multi-slit spectra: in the concept of operations of MUSE, the disambiguation code would be used to identify locations of multi-slit confusion so that users can isolate signals from secondary lines or a neighboring slit before analyzing the primary lines. For the application of this method to MUSE data, a more detailed description of the various trade-offs and optimal choices for velocity bins, temperature bins, and inversion parameters, as well as the dependence on abundance, density, and noise, will be provided in a follow-up paper (J. Martínez-Sykora et al. 2019, in preparation).

Slot and Slitless Spectrometers
A special case is when the slit spacing d is equal to the slit width w (in terms of angular coverage over the plane of the sky). This allows us to treat slitless spectrometers and spectrometers with slot modes within this multi-slit framework. An example of a slitless spectrometer is the S082A instrument on Skylab (Tousey et al. 1977). The Coronal Diagnostic Spectrometer on board the Solar and Heliospheric Observatory and the Hinode/EIS instrument (Culhane et al. 2007) are spectrographs with slot modes. The validation of this method to decomposing spectra from these types of instruments like the proposed COronal Spectroscopic Imager in the EUV (COSIE) is detailed in a companion paper (Winebarger et al. 2018).

Discussion
In this paper, we outlined a general approach to performing the decomposition of spectrograms from astronomical spectrographs with a variety of configurations (e.g., single-slit, multi-slit, and slot mode). The decomposition of spectrograms is not only helpful for removing possible ambiguities (e.g., from which slit this detector signal originated), but such decomposition is also useful for estimating the physical parameters of interest. In this paper and in the companion paper (Winebarger et al. 2018, on decomposing overlappograms from slot spectrometers), we demonstrate application of the technique for interpreting spectrograms from different instruments observing the solar corona.
The applications considered in this paper were for EUV radiation emanating from optically thin plasmas. However, the method is also useful for the interpretation and inversion of spectra from optically thick plasmas. Suppose the spectrometer is designed to detect an atomic absorption line formed in the (partially) optically thick photospheres or chromospheres of astrophysical objects. The emergent spectrum from the optically thick material depends on a number of physical parameters describing the atmosphere including the number density of the absorbing species, the ambient temperature, the slope of the source function S ν (τ) as a function of optical depth τ, the local magnetic field (if magneto-optical effects like Zeeman splitting are important), and more (especially if the line does not form in local thermodynamic equilibrium). Nevertheless, given a physics model including the desired effects, one can construct a library of emergence spectra over parameter space, and fold them through the optical model O to compute the response matrix R. Seeking a solution along the lines of Equation (6) remains an attempt to express the measured spectrogram as the linear superposition of spectra from the library of atmospheric models. However, since the operator P is not linear for optically thick radiation, the components of the solution vector x cannot be interpreted as a density of emitting material (in the case of optically thin EUV radiation, this density is the emission measure) along a single ray. Instead, the linearity should be attributed to the optical characteristics of the instrument (i.e., the operator O).
Within the framework presented here, a single component inversion is equivalent to seeking to describe the spectrogram y with a single response vector f r ( ). Unless one is certain the object is spatially resolved given the telescope PSF, there is perhaps no compelling justification (other than computational expediency) for a single component inversion. For instance, a spatially extended PSF would lead to the addition of signal associated with photons from different parts of an astrophysical object.
The multi-slit configuration of an instrument like MUSE can be thought of as having a spatially extended (over multiple slits) PSF. But even in the absence of multiple slits, the PSF of a telescope can be sufficiently extended that variations of the physical parameters of interest in the plane of the sky are not resolved. In such cases, the inversion of an observed spectrogram with a single atmospheric model (i.e., single set of physical parameters) may yield systematic errors. It is then perhaps more appropriate to fit the spectrogram with a linear combination of atmospheric models. Whether such an approach yields superior results over single component inversion remains to be tested and validated. The metrics of interest would depend on the specific use case and are outside the scope of this paper.
In this article, we outlined a novel framework for multicomponent decomposition of astronomical spectra. We have illustrated application of the method to solar observing instruments, but it can also be used for the interpretation of spectra from other astrophysical sources.
The MUSE team acknowledges support from NASA contract 80GSFC18C0012 to LMSAL. M.C.M.C., B.D.P., J.M.S., and P.T. acknowledge support by NASA's Heliophysics Grand Challenges Research grant Physics and Diagnostics of the Drivers of Solar Eruptions (NNX14AI14G to LMSAL). P.A. acknowledges funding from his STFC Ernest Rutherford  Figure 2 (panel (A)) is nicely reproduced by the synthetic Fe IX171 Å intensity (panel (B)) calculated after applying the spectral decomposition code on synthetic multi-slit MUSE observations including photon noise. The good correlation is illustrated by the lower panels, which show 2D histograms of the ground truth (vertical axis) and inverted intensities (horizontal axis) without and with photon noise for an exposure time of 1.5 s in panels (C) and (D), respectively. Vertical error bars show, for a given inverted Fe IX 171 Å intensity, the range of true intensities within the central 95th percentile range. The discrepancy in the reconstruction is consistent as would be expected from Poisson noise statistics (here for an exposure of 1.5 s). Hence the reconstructed intensities are comparable regardless of whether or not noise is added.