Richardson-Lucy deconvolution with a spatially Variant point-spread function of Chandra: Supernova Remnant Cassiopeia A as an Example

Richardson-Lucy (RL) deconvolution is one of the classical methods widely used in X-ray astronomy and other areas. Amid recent progress in image processing, RL deconvolution still leaves much room for improvement under a realistic situations. One direction is to include the positional dependence of a point-spread function (PSF), so-called RL deconvolution with a spatially variant PSF (RL$_{\rm{sv}}$). Another is the method of estimating a reliable number of iterations and their associated uncertainties. We developed a practical method that incorporates the RL$_{\rm{sv}}$ algorithm and the estimation of uncertainties. As a typical example of bright and high-resolution images, the Chandra X-ray image of the supernova remnant Cassiopeia~A was used in this paper. RL$_{\rm{sv}}$ deconvolution enables us to uncover the smeared features in the forward/backward shocks and jet-like structures. We constructed a method to predict the appropriate number of iterations by using statistical fluctuation of the observed images. Furthermore, the uncertainties were estimated by error propagation from the last iteration, which was phenomenologically tested with the observed data. Thus, our method is a practically efficient framework to evaluate the time evolution of the remnants and their fine structures embedded in high-resolution X-ray images.


INTRODUCTION
Imaging analysis is critically important for studying diffuse celestial sources. X-ray astronomy, starting with the first space application of X-ray CCD in ASCA (Burke et al. 1994), has delivered detailed images of various celestial objects; e.g., supernova remnants such as SN1006 (Bamba et al. 2003) and Cassiopeia A (hereafter Cas A; Hwang et al. 2004), and galaxy clusters such as A2142 (Markevitch et al. 2000). Since the X-ray mirrors used in Chandra are the largest and most precisely built, exceeding the angular resolution of Chandra is considered to be challenging. Therefore, enhancing the technique of imaging analysis has been an essential direction to utilize the highest spatial resolution and data accumulated over decades.
X-rays are collected primarily by total reflection from the surface of an X-ray mirror, therefore the response function for the distribution of focused X-rays, called the point-spread function (PSF), is nearly energy independent. On the condition that a PSF is independent of incoming photon energy and the position of the focal plane, a reverse calculation of a convolution of PSF, so-called image deconvolution (see the review on deconvolution in astronomy by Starck et al. (2002)), is highly simplified. There are various deconvolution methods proposed by assuming that a PSF is constant during the deconvolution process, such as the deconvolution of Suzaku XIS (Sugizaki et al. 2009). One of the latest examples is the image restoration algorithm Expectation via Markov chain Monte Carlo (Esch et al. 2004). It is applied to the double active galactic nuclei in NGC 6240 (e.g., Fabbiano et al. 2020;Paggi et al. 2022), succeeding in finely resolving the two cores. Similarly, a classical method, Richardson-Lucy (RL) deconvolution proposed by Richardson (1972) and Lucy (1974), is often used (e.g., Grefenstette et al. 2015;Thimmappa et al. 2020;Sobolenko et al. 2022).
The choice of method depends on the trade-off between accuracy and computational cost. Relaxing the condition that a PSF is positional-independent and/or energy-independent, the deconvolution methods increase the complexity of the calculation. RL deconvolution is one of the simplified methods but still has room for improvement in practical situations. In gamma-ray astronomy, a PSF can change by one order of magnitude with energy and incident angle; it is calculated for each event, e.g., RL algorithm optimized for Fermi-LAT and EGRET (Tajima et al. 2007). In contrast, as the number of photons is much larger in X-ray astronomy, event-by-event reconstruction is less practical; image-based reconstruction thus can be the first choice in X-rays. However, there are few studies on extending the RL method, especially their application to diffuse sources obtained by Chandra. We therefore explored its applicability to the Chandra data and considered the associated systematic errors.
In this paper, we implement RL deconvolution with a spatially variant PSF (RL sv ) algorithm, assuming it to be used for Chandra images. Section 2 describes the principle of the RL sv method. One of the technical difficulties is reducing computational cost in calculating PSFs. This is solved by decimating the sampling interval of PSFs, while a side-effect is discussed in Section 5.1. Section 3 presents an example of its application to a diffuse source observed by Chandra. We apply the RL sv method to the supernova remnant of Cas A as an example, because Cas A is bright and extended over the entire field of view of the ACIS detector, which would be the best target for the first application. The remnant is intensively studied because of its unique structure and evolution, e.g., the velocities and thickness of shocked filaments (Patnaude & Fesen 2009;Sato et al. 2018;Tsuchioka et al. 2022), where the method can contribute to advancing our understanding of the phenomena. In Section 4, we propose a reliable number of stop iterations and uncertainties of the method. We develop the method to estimate the number of convergent iterations by generating fluctuations due to statistical errors during iteration. Furthermore, the uncertainty on the RL sv -deconvolved image is estimated by using the law of error propagation (e.g., Ku 1966). As a result, filaments and ambiguous structures of Cas A are deconvolved to be sharper with some knowledge of the statistical uncertainties.

RL Deconvolution
The RL algorithm iteratively estimates a true image from an observed image using Bayesian inference. It generally assumes that the PSF does not change with a position in the image. The RL algorithm is expressed by where i and j are mapping the image in the sky, and k is mapping the image on the detector. The indices of the summation run through all the pixels. W (r) is the restored image after r iterations, and H is the observed image on the ACIS detector. P jk is the probability that a photon emitted in sky W bin j is measured in data space H bin k, or P (H k |W j ).

RL with a Spatially Variant PSF
Previous Chandra image deconvolution approaches (e.g., Thimmappa et al. 2020;Sobolenko et al. 2022) used a simplified approximation for the P jk values, i.e., they used the same PSF for each j bin. Here we assume that the PSF changes as a function of the off-axis angle and the roll angle. As a consequence, the Chandra RL algorithm is extended. The formula for RL sv is obtained by rewriting Equation (1) as (2) P jjk refers to a PSF at a position of j (first index) which returns a probability that an event emitted at W j (second index) is observed at H k (third index), or P j (H k |W j ). Computational cost and memory requirements need to be minimized for calculating the thirdorder tensor of the PSF, which is a distinctive feature of the RL sv algorithm. When H is corrected for slight differences among pixels in effective area and exposures, its normalization can be chosen arbitrarily. Here we use where N k is the detector count image, and A k is the Hadamard product of effective area and exposure time.

RL sv Deconvolution with Total Variation Regularization
There are regularization techniques to enhance the RL method, which are also readily available for the RL sv method. Among these, total variation (TV) regularization (Rudin et al. 1992) is effective in handling statistical errors, which is used in the RL method (Dey et al. 2006). The formula for RL sv with the regulariza- tion is obtained by rewriting Equation (2) as (3) The only difference from the RL sv algorithm of Equation (2) i . In this paper, we utilize the parameter of λ TV = 0.002, as proposed by Dey et al. (2006).

Comparison to other methods
Deconvolution methods require an understanding of their applicability to a practical condition, as well as optimization of computation cost and accuracy (for features of various methods see Naik & Sahu (2013)). The RL method is well studied and has been used to incorporate regularization (e.g., van Kempen & van Vliet 2000; Dey et al. 2006;Yuan et al. 2008;Yongpan et al. 2010) and the recent trend of deep learning (e.g., Agarwal et al. 2020). For the Chandra users, RL deconvolution with a single PSF is frequently used because the method is already implemented as arestore in the Chandra Interactive Analysis of Observations (CIAO; Fruscione et al. 2006), Chandra's standard data processing package.
Compared to other methods, the RL method forces the deconvolved image of each iteration to be nonnegative, and its integral value is conserved. Additionally, the method converges to the maximum likelihood solution for a Poisson noise distribution (Shepp & Vardi 1982), which is suitable for Chandra images with noise from counting statistics. Depending on the application, it is less prone to ringing artifacts than inverse PSFbased methods (e.g., Sekko et al. 1999;Neelamani et al. 2004); see the results of the comparison by Dalitz et al. (2015). According to White (1994), it is robust against small errors in the PSF. Because Cas A is a bright and diffuse X-ray source with a moderately large apparent diameter, it is an ideal target to demonstrate the RL sv method. It has been observed by Chandra almost every year since 1999. The Chandra data of Cas A used in this paper are listed in Table 1: ACIS-S observation of 2004 in Obs. ID=4636, 4637, 4639, and 5319. The image size is 1489×1488 pixels, or 743 ′′ × 742 ′′ given a unit pixel of 0. ′′ 492. Data processing and analysis were performed using CIAO version 4.13. The data were reprocessed from the level 1 event files by chandra repro. Since the roll angle and optical axis of the four observations are almost the same (the maximum difference of the optical axis location is about 4 unit pixels), all the events were merged into one by merge obs. The total exposure time was 428.29 ks.

Generating the PSF of Chandra
The Chandra telescope system consists of four pairs of nested reflecting surfaces, configured in the Wolter type I geometry. The high energy response is achieved by coating the mirrors with iridium. It has attained the highest angular resolution of 0. ′′ 492 among existing X-ray telescopes. Its mirror of Chandra has been extensively calibrated on the ground and in orbit (Jerius et al. 2000). The Chandra PSF is positional-dependent, mainly due to aberrations. The RL sv method includes the position dependence of the PSF, which is useful for highly extended X-ray sources. Because creating a PSF for each position is computationally expensive, it is decimated at some intervals. For reference, creating a PSF takes several seconds, depending on the computational environment and desired accuracy. The sampling interval of PSFs was chosen to be 35 × 35 pixels (total of 43 × 43 = 1849). The interval was determined empirically by trying several different ranges. In general, the PSFs simulated from each observation should be merged for a precise calculation. Here, the PSF of the Obs. ID of 4636 was used as a representative since its sampling of the PSF is decimated. The PSFs at the lattice points were generated by CIAO's simulate psf using the Model of AXAF Response to X-rays (Wise et al. 1997;Davis et al. 2012) at a monochromatic energy of 2.3 keV. They were applied to the observed image with energies from 0.5 to 7.0 keV. Figure 1 shows all the PSFs sampled every 35×35 pixels. The optical axis is located at the northeast in the image, where the spread of the PSF is minimum. As the position is away from the optical axis, the tail of the PSF increases with its gradual shift of the elliptical axis. Although it is a trade-off with photon statistics, it is effective to run the RL sv method with the optimal monoenergetic PSF for each of the multiple energy decompositions (see Figure 6). This is because, at shorter wavelengths, the effect of diffuse reflection due to the roughness of the mirror surface is not negligible (Jerius et al. 2000).

3.3.
Results of the RL sv method Figure 2(a) is an observed ∼400 ks image using the energy range from 0.5 to 7.0 keV, as explained in Section 3.1. We applied the RL sv method to the image with a sampling interval of each PSF of 35 × 35 pixels. The number of iterations is 200. Note that the choice of the iteration number is discussed in Section 4.1. The result of the RL sv method is presented in Figure 2(b). The unit of flux and its range in Figure 2(b) is the same as in Figure 2(a). The overall structures in the RL sv image become more vivid than the original ones. The images around the off-axis are significantly improved compared to those around the optical axis.
To make the differences more precise, we present magnified images of the original image in Figures 2(a-1, 2, 3) and the RL sv -image in Figures 2(b-1, 2, 3). The three regions represent a sharp filament in the northeast, complicated filaments in the north, and a slightly diffuse area in the south. The filamentary structures in the northeast and north become sharper in the RL sv -image. We will quantify the filament width in detail and discuss the systematic uncertainties associated with the method in Section 4.2.

Assessment of the reasonable number of iterations
We considered a way of assessing an appropriate number of iterations, which is one of the issues with the RL method. This is because the method has the property of excessive amplification of noise as the number of iterations increases. We propose a method to suppress convergence using statistical errors during the iteration. The formula of the method is written as The only difference from the RL sv algorithm, Equation (2), is the G(N k )/A k term. N (counts) is the map of detector counts. A (photons cm −2 s −1 ) is the Hadamard product of effective area and exposure time. G(N k ) is a random number generator following a Poisson distribution with a count in the kth pixel of N k ; i.e., G(N k )/A k is a flux in units of photons cm −2 s −1 . The reason for normalizing G(N ) by dividing it with A is to account for the slight variations in effective area and exposure time among pixels. The performance of the RL sv algorithm using Equation (4) is compared to that using Equation (2). The convergence is evaluated by using the mean squared error (MSE) of the two images: one step before and after iterations. Figure 3 shows the history of the MSE during the iteration. The curve obtained by the RL sv algorithm, Equation (4) saturates at a certain level, while the other continues to decrease. The saturation level is caused by the injection of Poisson fluctuation at each step, which is considered as an indicator of stopping. In Figure 3, the iteration number of ∼30 seems appropriate.

Assessment of image blurredness
We then designed a simplified method for evaluating a certain amount of confidence. The RL sv -deconvolved image should have a similar amount of fluctuation accompanying the observation image. The principle of the method is to propagate errors of the converged RL sv image into the next step. Our choice of using the last step for the error propagation is just to simplify the task. Here, each error in the observed image is considered statistically independent. Assuming only uncertainties on H k , using the law of error propagation, the image uncertainty can be expressed as where W ′ is the image of the next iteration number of any estimated true image of W , and √ N k is the statistical error of N k .
We compared the off-axis RL image with the error of Equation (5)  outer regions. These radial profiles were created using dmextract in CIAO. In Figure 4(d), the framed regions from Figures 4(a)-(c) are color-coded as blue, orange, and black, respectively. The error bars of the radial profiles in Figure 4(d) correspond to statistical errors, represented by blue and black, and the result obtained by applying Equation (5) to the 199th iteration of the RL sv image, is indicated by orange. From Figure 4(d), the profile of the off-axis RL sv image agreed with that of the on-axis image within the statistical errors. This method gives a guideline for a certain level of confidence associated with the RL sv method.

Enhancement Technique and Possibilities
In this section, further enhancements to the RL sv method are discussed. The first is to reduce the loss of down-sampling PSFs. The positional dependence of the PSF does not contain high-frequency components, so the decimation of PSF sampling should work to some extent. For a small image such as a core plus jet structure in an active galactic nucleus, keeping a high sam- pling rate of the PSFs might be possible. However, for a largely extended source such as a supernova remnant or galaxy cluster, to minimize the sampling rate is critically important for practical use. The higher the decimation, the more emphasized the boundary of the segment. Taking Cas A as an example, the edges of specific segments clearly appear when the sampling interval is 35 × 35 pixels. To smooth out the edges, we propose that the PSFs' boundaries be randomly selected from nearby PSFs (see more details in the Appendix). Second, this method can be developed by incorporating several regularization methods. We implemented an RL sv method incorporating the TV regularization expressed in Equation (3). Finally, the RL sv method is naturally applied to color images. By decomposing observed images into several colors (or energy bands) and generating PSFs for an appropriate energy in each band, an energy-dependent RL sv method can be realized.
We compare the RL sv method with and without the TV regularization. We use the same image as in Section 3.1. The PSF sampling is 35 × 35 pixels and the number of iterations is 200. The PSFs' boundaries are randomly selected following the Appendix. Figure 5 presents the enlarged eastern image after applying the RL sv method to the entire region. Figures 5(a) and (b) show the results of the RL sv method of Equation (2) and the regularization version of Equation (3), respectively. The TV regularization preserves the sharp structure to remain and smoothes out statistical errors. In this way, regularization can be added to the RL sv method.
We implement the RL sv method including these enhancements and adapt it to the Cas A observational data described in Section 3.1. Figure 6(a) is the observed image of Section 3.3 divided into three energies in RGB: 0.2-1.2 keV (red), 1.2-2.0 keV (green), and 2.0-7.0 keV (blue). Cas A is dominated by the thermal radiation in 2.15×10 -8 9.00×10 -8 3.07×10 -7 9.90×10 -7 3.16×10 -6 ⩽4 keV and the nonthermal radiation in ⩾4 keV. We applied the RL sv method with the TV regularization Equation (3) to each energy image of Figure 6(a) using the appropriate energy of the PSF (red is 0.92 keV, green is 1.56 keV, and blue is 3.8 keV) based on the official CIAO page. 1 The sampling interval of PSFs is 35×35 pixels. PSF is randomly selected at the sampling boundaries according to the Appendix. The number of iterations is 30, according to Section 4.2. The result of the RL sv method is presented in Figure 6(b). The energy dependence in Figures 6(b-1) and (b-2) are clearly visible by this method.

Constraint on the Uncertainty
We propose two complementary ways to obtain a guideline on the stop condition of the RL sv method. One is to obtain a minimum of residuals by inserting statistical uncertainties into each update during the iteration process. This gives a rough estimate of the limit of the iterations. The other is to include the statistical uncertainties in the last step of the iteration. By combining the two methods, it is possible to derive uncertainties using the errors obtained by the latter method at the optimal iteration number estimated by the former. This is a quick and convenient way to derive the systematic uncertainties associated with the RL sv method.
The systematic uncertainty in the former method is a way of defining the optimal number of iterations. One easy way is to use the same level of residuals as shown in Section 4.1. It is intrinsically difficult to distinguish the signal of the celestial objects from the statistical noise. This difficulty needs to be overcome by comparing the deconvolved images without errors and the errorestimated images around the optimal iteration number recommended by the former method. This method is based on a compromise between the computational cost and the simplicity of use while keeping a reasonable statistical error.

CONCLUSION
We have improved the processing capability of RL deconvolution by incorporating the positional dependency of the Chandra PSF. The RL sv method is applied to the entire region of Cas A with an estimation of its limit and errors, which are based on the phenomenological method for evaluating a reasonable number of iterations and uncertainties. It shows that the features of shock waves and jets are sharper than those measured in the original image, with a certain amount of knowledge of the associated errors. The RL sv -deconvolved profile of the off-axis image at the southeastern filament became shaper and agreed with that of the on-axis observation within the statistical errors. This method is useful for a detailed diagnosis of other extended X-ray sources obtained by Chandra.
The code used in this paper is available at doi:10.5281/zenodo.8020557.
We would like to thank the anonymous referee for helpful comments and feedback on this paper. This research has made use of data obtained from the Chandra Data Archive and the Chandra Source Catalog, and software provided by the Chandra X-ray Center (CXC) in the application packages CIAO. This work was supported by JSPS KAKENHI grant Nos. 20K20527, 22H01272, and 20H01941.

APPENDIX BOUNDARIES OF THE PSF
Decimating the sampling number of PSFs is an effective approach to minimize computational cost. However, this technique can introduce side-effects at the boundary of segments when switching between PSFs. The variation in shape between neighboring PSFs, caused by the sampling interval, leads to discontinuities in the deconvolution process. To mitigate this issue, a simple countermeasure is to randomly select adjacent PSFs at their boundaries, which helps to smooth out the discontinuities. The weights of the probabilities for selecting PSFs are illustrated in Figure 7. The presence and severity of artifacts depend on factors such as the dissimilarity in shape between neighboring PSFs, statistical characteristics of the observed images, and other relevant factors. Therefore, the presented technique serves as an example, and the problem is optimized by the range of pixels to be randomized.  Figure 7. Illustration of the 9 × 9 pixels around the intersection of the PSF switchover, overlaid with the probability weights on selecting PSFs. The quadrants are named A, B, C, and D for illustrative purposes. The weights of the probabilities are either 1/3 or 2/3 at the two boundaries, while they are 1/9, 2/9, or 4/9 at the corners. The result of applying the selection rule in Figure 7 is shown in Figure 8. We compared the RL sv method for the observed data in the eastern region in Section 3.1 with a PSF sampling of 35×35 pixels and 200 iterations.  (c), where the white lines are used to clarify the border lines. Figures 8(b) and (c) are the RL svdeconvolved images with and without using the randomization of PSFs, respectively. To illustrate the differ-ences, we presented magnified images of Figures 8(b) and (c) as Figures 8(b-1) and (c-1), respectively. It appears that the discontinuity at the boundaries of the PSFs, indicated by the green arrows, is smeared out to some extent.