A Data-driven Approach to X-ray Spectral Fitting: Quasi-Deconvolution

X-ray spectral fitting of astronomical sources requires convolving the intrinsic spectrum or model with the instrumental response. Standard forward modeling techniques have proven success in recovering the underlying physical parameters in moderate to high signal-to-noise regimes; however, they struggle to achieve the same level of accuracy in low signal-to-noise regimes. Additionally, the use of machine learning techniques on X-ray spectra requires access to the intrinsic spectrum. Therefore, the measured spectrum must be effectively deconvolved from the instrumental response. In this note, we explore numerical methods for inverting the matrix equation describing X-ray spectral convolution. We demonstrate that traditional methods are insufficient to recover the intrinsic X-ray spectrum and argue that a novel approach is required.


INTRODUCTION
The intrinsic emission spectrum of an X-ray source, F (E), is defined by the underlying physical emission processes over a continuum of photon energies, E (e.g. Kahn & Blissett 1980;Weisskopf 1999;Weisskopf et al. 2002). The spectrum has units of photons m −2 s −1 keV −1 . We relate it to the observed spectrum, S(E ), through a convolution with the instrumental response, R(E , E): Note that E denotes the discrete photon energies captured by the detector. Fitting over a broadband spectrum in this manner is common practice in X-ray astronomy due to the limited resolution of the detectors. The instrumental response matrix for the Chandra X-ray Observatory is captured in two parts: the redistribution matrix file (rmf) and the ancillary response file (arf) which are encapsulated in R(E , E). The rmf contains the mapping from the continuous energy space to the detector position space. Analogously the arf contains the Corresponding author: Carter Rhea carterrhea@astro.umontreal.ca effective area of the detector as well as its quantum efficiency as a function of time-averaged energy. Although we focus on Chandra response matrices, they are ubiquitous in X-ray astronomy; thus, this problem extends to other existing and future X-ray observatories. The standard method for determining the true spectrum requires the following prescription: choose an appropriate parametric, physically-derived model to explain the intrinsic emission and fit the model using equation 1. The fit is generally optimized by reducing the chi-squared statistic (e.g. Arnaud 1996). Alternatively, it is possible to deconvolve equation 1 and directly solve for F (E). However, the illconditioning of the response matrix makes this method unstable and thus infrequently used (e.g. Blissett & Cruise 1979). Matrices are considered ill-conditioned when their rows are not linearly independent of one another. Certain applications, such as extracting spectral parameters using machine learning techniques (e.g. Rhea et al. 2021b), require the intrinsic spectrum rather than the observed spectrum. Therefore, the response matrix must be taken into account in order to isolate F (E) from equation 1. Since the response matrix greatly affects the observed spectrum and changes significantly across the Chandra field-of-view, the handling of it is crucial to proper analysis. Thus, we investigate several arXiv:2105.09470v1 [astro-ph.IM] 13 May 2021 numerical methods to deconvolve the intrinsic spectrum from the response matrix.

METHODOLOGY
In the following section, we will describe the derivation of (and reasoning behind) the matrix formulation of equation 1 and methods for solving the matrix equation.

Deriving the Matrix Formulation
Despite the simplicity of equation 1, a direct convolution of the instrumental response function and model spectrum poses several issues. Due to the finite spectral resolution, the rows of the response matrix are not independent. We confine the integral to the energies covered by the detector. Additionally, since the sampling of the detector space E is discrete, we can rewrite equation 1 as a matrix equation (Kaastra & Bleeker 2016): We have replaced the instrumental response function by its matrix counterpart, R ij . In this formulation, S i represents the observed photon count rate in units of counts s −1 for a given detector energy bin. F j is the model spectrum flux in units of counts m −2 s −1 in emitted energy bin j. For simplicity, we will write equation 2 in the following form:

Solution Methods
Several methods for standard matrix equation solutions exist (Ax = b; ref); however, our application poses an additional constraints: the response matrix is illcondition (condition number >> 100; e.g. Wilkinson 1972). The most straight-forward solution is to use a Moore-Penrose pseudo-inverse (e.g. Penrose 1955) to invert the response matrix. In doing so we can directly calculate the intrinsic spectrum with a single matrix multiplication: where R † is the pseudo-inverse. The pseudo-inverse allows for the computation of an unique inverse matrix for non-square systems. Although this method is computationally efficient and direct, it suffers spurious oscillations owing to the ill-conditioning of R (e.g. Varah 1973). In order to diminish the effects, we instead solve the normal equations: Doing so has the added benefited of inverting a square matrix. Since RR T is still ill-conditioned, we again use the Moore-Penrose pseudo inverse: Numerous methods exists to solve ill-posed problems such as that described by the Fredholm Integral Equation of the First Kind illustrated in X-ray spectral analysis (equation 1; e.g. Hansen 1992). We explore two methods: preconditioning (e.g. Estatico 2002) andregularization (e.g. Neumaier 1998). We apply standard preconditioning using the normal equations; however, the regularization technique is more involved.

Regularization
Matrix regularization is a family of algorithms designed to overcome ill-conditioned matrices by imposing a strict condition such as smoothness on a least-squared solution. We apply Tikhonov regularization which augments the standard least-squares formalism by a Lagrangian mutliplier (e.g. Calvetti & Reichel 2004): where λ > 0 is the regularization parameter, D is the regularization matrix, and · L2 is the L2 (or Euclidean) norm (Horn & Johnson 2012). A standard choice for the regularization matrix is the identity matrix, I. We optimize the value of λ using the standard L-curve analysis ( Kindermann & Raik 2019). The result of minimizing equation 7 is the following expression for the intrinsic spectrum: We must create a set of mock X-ray spectra in order to test the feasibility of solving the matrix equation (equation 2) by inversion or a least-squares method.

Creation of Data
We use the sherpa (v4.13) tool fake pha in order to create mock spectra with an approximate signal-to-noise ratio of 20. We take two rmf and arf files from different regions on ACIS-I3 from the observation 7253 (we note that the choice in observation is arbitrary and used to demonstrate the spatial changes in the convolved response matrix; the choice of chip is also arbitrary). We tested several chips and several different observations covering a range of Chandra observation cycles. We report no difference in our results. We test two intrinsic emission types: a simple 1-dimensional powerlaw (powerlaw) and an absorbed thermal plasma emission model (phabs*apec). The powerlaw index parameter is set to −0.5, the column density, n H is set to 10 20 cm −2 . Although these values are arbitrary, we tested several values with no change in the results. Two spectra are created for each emission type; they differ only in the arf and rmf used in their creation.

RESULTS & DISCUSSION
When applied to the powerlaw model, the full spectral unfolding using Tikhonov regularization paired with a normalized preconditioner recovers the intrinsic spectra correctly up to several percent (<3% errors). However, when applied to a physically motivated thermal model (such as MEKAL or APEC), the algorithm only successfully recovers the underlying continuum spectrum. Unfortunately, the method fails to fully capture the prominent emission lines (such as Fe K-α). Recovering the precise shape and peak amplitude of these emission lines is crucial for subsequent calculations of temperature and metallicity. Therefore, this method is ill-suited to solve the inverse problem posed in equation 1. The authors are currently exploring the use of recurrent neural networks to solve the equation.