Image Deconvolution To Resolve Astronomical X-ray Sources in Close Proximity: The NuSTAR Images of SXP 15.3 and SXP 305

1. Introduction

1.1. Astronomical Point Sources and X-Ray Telescopes

Astronomical systems as observed by modern telescopes are divided into two types: point-like sources and extended sources. The former are unresolved and present images whose sizes and shapes are dictated purely by the combined point spread function (PSF) of the telescope optics and the detector; whereas extended sources and crowded fields (containing many closely-spaced point sources) present images that represent the true brightness patterns in the sky convolved with the imaging PSF [1]. If a reliable estimate of the PSF can be made, then in principle, mathematical deconvolution can reduce the instrumental blurring and provide a sharper image.

Most observations from space-based X-ray telescopes are used to study distant galactic and extragalactic point sources. Examples of such point sources are X-ray binaries [2], novae [3], and active galactic nuclei [4], while extended sources are found in supernova remnants [5], galaxy clusters [6], star-forming regions [7], and galactic halos [8]. The goal is always to obtain a more detailed image of the X-ray source irrespective of its nature. In extended sources, the finer the image, the more detailed the structure of the source. In the case of point sources, higher-resolution images can be crucial in distinguishing between two closely-spaced sources (actually their 2D projections in the sky).

Among presently operating X-ray observatories, Chandra [9] has the best spatial resolution and can distinguish even between sources projected to 0.5″ apart; whereas NuSTAR[10] has a broader energy coverage (viz. 3-79 keV versus 0.5-8 keV for Chandra), but also (suffers from) a broader PSF (a half power diameter (HPD) of 58″). The NuSTAR PSF also deteriorates considerably for off-axis sources (observed far out from the image center). This can create a grave issue in distinguishing between point sources lying within the full radius of the PSF when there are no simultaneous observations from the Chandra or the XMM-Newton [11] telescope. Image deconvolution can then help resolve the problem to a certain extent. For instance, Sakai et al. [12] used a Richarson-Lucy image deconvolution technique on a high-resolution Chandra image of the supernova remnant Cas A to bring out more detailed features in the low-flux regions. A non-iterative deconvolution technique has also been used on Suzaku-XIS [13] images to enhance the resolution from its 110″ HPD to $\sim {20}^{\prime \prime }$ [14].

1.2. Image Deconvolution

Image deconvolution is a mathematical technique used to improve the resolution or clarity of images by reversing the effects of blurring and distortion. These distortions arise from the instrument’s PSF, atmospheric effects, or other turbulence during observation. Deconvolution aims to reconstruct a more accurate version of the original image by using prior knowledge of the PSF and other known properties of the telescope. There are several methods for image deconvolution, each with specific advantages and limitations. Iterative methods, such as the Richardson-Lucy algorithm [15,16], apply successive approximations to recover the image, and they are particularly effective in dealing with noisy data. Fourier-based methods [17] work in the frequency domain to perform deconvolution efficiently, but they can be more sensitive to noise. Wavelet-based deconvolution [18] is often used for images with structures and edges, allowing for better edge detection and detail recovery. More recently, machine learning techniques [19] have become very common; they use deep learning models trained on known PSFs and image pairs, and they can achieve results comparable to or better than those from deconvolution methods; but these new techniques are heavily dependent on the availability of training data.

Applications of image deconvolution have appeared in diverse fields. In astronomy, deconvolution is used to sharpen images of distant stars, galaxies, and various other celestial objects captured by ground-based and space-based observatories. In medical imaging, deconvolution improves the quality of microscopy and radiology scans, aiding in more precise diagnoses. It is also widely used in remote sensing, surveillance, and increasingly in consumer photography to enhance image quality and detail.

1.3. High-Mass X-Ray Binary Pulsars

X-ray pulsars are rapidly spinning neutron stars (NSs) that emit X-rays modulated by their spin periods. These pulsars are found in NS X-ray binary (XRB) systems with high-mass or low-mass companion stars from which the NS accretes material. The material falls onto the star’s magnetic poles, forming bright spots (`hotspots’) that emit X-ray photons. As the NS rotates, the X-ray beams sweep across space, appearing as regular pulses to any observer. The pulse period tells us how fast the neutron star is spinning, and changes in this period over time reveal details about the interaction of the NS with its companion star [20,21].

High-mass X-ray binaries (HMXBs) are systems in which the NS companion is a massive ($>8M_{\odot }$) hot star, usually of spectral type O or B. A special type of HMXBs is the subset of Be X-ray binaries in which the companion is a Be star. Be stars are a class of B-type stars that show Hydrogen emission lines (e.g., H$\alpha$) in their optical spectra (hence the `e’ in the class name). They form large decretion disks around them due to material shed from their rapidly rotating surfaces enhanced by processes thought to be related to stellar pulsations [22,23]. Decretion disks play a crucial role in driving the observed X-ray behavior in Be-XRBs. When the NS passes through or interacts with the disk material, then the process of accretion occurs leading to powerful X-ray outbursts [24].

1.4. The Magellanic Pulsars SXP 15.3 and SXP 305

The Small Magellanic Cloud (SMC) is a nearby dwarf galaxy that is home to more than 100 HMXBs [25], about 70 of which are X-ray pulsars, and all but one are Be-XRBs [21]. The low metallicity of the SMC makes it a great environment for forming massive stars and studying their life cycles [26]. With a plethora of known HMXBs, the SMC offers unique opportunities to understand how these systems evolve.

In this work, we address the difficult problem of resolving point sources in NuSTAR X-ray images, especially when there is no other major distinguishable spectral or timing feature, and the observer has to rely entirely on the imaging capability of the telescope. The particular case that we explore concerns the closely-spaced pair of SXP 15.3 and SXP 305 [27,28]. These two pulsars have widely different spin periods, and their X-ray outbursts allow us to study the extreme conditions that develop near their accreting NSs.

SXP 15.3 [29,30,31] and SXP 305 [28] are both Be-HMXB pulsars [21,26], and they are separated by only 7″ in the sky (see Figure 1 of Monageng et al. [28]). Lazzarini et al. [27] determined the pulse period of SXP 305 for the first time using a NuSTAR observation (also used in this work), but the NuSTAR image did not resolve the pair as the projected separation of the two sources is much smaller than the NuSTAR HPD of 58″.

2. Methods

2.1. General Framework for Image Deconvolution

Image deconvolution is a computational technique used to reconstruct the `true’ image $I_{\mathrm{true}}(x,y)$ that has been blurred by the instrument’s PSF and the superimposed noise. The observed image $I_{\mathrm{obs}}(x,y)$ is expressed as

$$I_{\mathrm{obs}}(x,y)=I_{\mathrm{true}}(x,y)\ast PSF(x,y)+N(x,y),$$

(1)

where ∗ denotes convolution, and $N(x,y)$ represents the noise at each image point $(x,y)$. The goal of deconvolution is to estimate $I_{\mathrm{true}}(x,y)$ by reversing the effects of the PSF while accounting for noise.

Deconvolution is an ill-posed problem because noise in the observed image can be amplified during the reduction process. To address this issue, regularization techniques and additional constraints are applied. Common deconvolution approaches include:

• Direct Methods: These include Fourier-based deconvolution and Wiener filtering, which are computationally efficient but very sensitive to the noise content.
• Iterative Methods: Methods such as the Richardson-Lucy algorithm that refine the estimate during successive iterations and that are very effective against Poisson noise.
• Bayesian Methods: These incorporate prior knowledge about the image, such as smoothness or sparsity, to guide the cleaning algorithm.
• Machine Learning Approaches: Neural networks trained on simulated data can directly predict the `true’ image from the observed outcome [19].

The general update rule for iterative deconvolution is written as

$$I_{\mathrm{true}}^{(k+1)}(x,y)=I_{\mathrm{true}}^{\left(k\right)}(x,y)+\lambda ·\mathcal{R}(I_{\mathrm{true}}^{\left(k\right)},I_{\mathrm{obs}},PSF),$$

(2)

where $I_{\mathrm{true}}^{\left(k\right)}(x,y)$ is the estimate after k iterations, $\lambda$ is a step size, and $\mathcal{R}(\cdots )$ is a residual term based on the convolution mismatch. Iterations continue for as long as the residual is minimized and an adopted convergence criterion is met.

2.2. The Richardson-Lucy Algorithm

The Richardson-Lucy (RL) algorithm [15,16] is an iterative deconvolution technique optimized for low-flux data, such as many X-ray observations of SMC sources. It is based on maximum likelihood estimation under Poisson statistics. The algorithm updates the image estimate by

$$I_{\mathrm{true}}^{(k+1)}(x,y)=I_{\mathrm{true}}^{\left(k\right)}(x,y)·{\displaystyle \frac{I_{\mathrm{obs}}(x,y)\ast PS{F}^T(x,y)}{I_{\mathrm{true}}^{\left(k\right)}(x,y)\ast PSF(x,y)}},$$

(3)

where $PS{F}^T(x,y)$ is the transpose (or flipped) version of the PSF. The numerator represents the cross-correlation of the observed image with the flipped PSF, while the denominator represents the deblurred estimate of the true image.

The execution steps for the RL algorithm are:

• Initialization: Start with an initial guess for the true image, typically a uniform image or a given $I_{\mathrm{obs}}(x,y)$.
• Update Rule: Correct the estimate iteratively based on the ratio of observed data to the simulated blurred estimate.
• Stopping Criterion: Stop the iteration when the changes in $I_{\mathrm{true}}(x,y)$ are negligible or when the residual error is sufficiently small.

The RL algorithm assumes an accurately known PSF and works quite well in images with superimposed Poisson noise. Various extensions of this algorithm can also handle spatially varying PSFs. This work uses a Python `scikit’ image restoration package called skimage.restoration.richardson_lucy (URL: https://scikit-image.org/, accessed on 26 October 2024).

2.3. NuSTAR Data Analysis

NuSTAR (Harrison et al. [10]) is the first orbiting telescope capable of focusing very high-energy X-rays. It was launched on June 13, 2012 and operates in the 3–79 keV energy range, surpassing the  10 keV limit of earlier focussing X-ray telescopes. NuSTAR comprises two co-aligned hard X-ray grazing incidence telescopes (Wolter I design) that focus X-rays onto two separate solid-state focal-plane detectors with a focal length of approximately 10 meters. The two focal-plane modules (FPMA and FPMB) are designed to be identical, enabling their combined data to enhance the signal.

The archival data from both FPMA and FPMB modules were processed using the NuSTAR Data Analysis Software (NuSTARDAS v2.1.4) with calibration files from CALDB (v20240325). This software is included in the HEASoft (v6.34) package distributed by the High Energy Astrophysics Science Archive Research Center (HEASARC). The nupipeline command was employed to convert level-1 event files into level-2 event files, applying standard data cleaning procedures. High-level products, including barycenter-corrected source event files, skymap images, light curves, and spectra, were generated using the nuproducts command. For this study, skymap images, created in the broad energy band of 3-79 keV, were utilized.

2.4. NuSTAR PSF

The PSF of NuSTAR describes how the telescope and detector system responds to point X-ray sources, effectively characterizing the blurring introduced by the instrument’s optical components. It is effectively determined by its focusing optics, which consist of multilayer-coated, grazing-incidence mirrors designed to reflect high-energy X-rays. The on-axis PSF has a full width at half maximum of approximately 18″ and an HPD of about 58″ [10].

However, the PSF varies with off-axis angle due to optical aberrations inherent in the telescope’s design. As sources appear offset from the optical axis (the image center), the PSF broadens and becomes more asymmetric, resulting in decreased spatial resolution and increased overlap between nearby sources. Additionally, the PSF exhibits a slight energy dependence because of reflectivity variations in the multilayer mirror coatings. Accurately modeling the NuSTAR PSF—including its dependence on energy and off-axis angle—is crucial for image deconvolution techniques aiming to resolve sources in close proximity, such as SXP 15.3 and SXP 305 [27,28]. In this study, we have utilized the calibrated PSF models stored in CALDB, and we have accounted for the specific offsets and energy ranges that pertain to the particular FPMA and FPMB data sets.

Figure 1. Central 100×100 pixel regions of the NuSTAR PSF at different off-axis angles ranging from 0′ to 7′. Each panel represents the PSF corresponding to a specific off-axis angle, highlighting the variation in spatial distribution due to the telescope’s optics. The color scale indicates relative intensity, and the images are shown in a logarithmic scale for a better view of the PSF structure.

Figure 1. Central 100×100 pixel regions of the NuSTAR PSF at different off-axis angles ranging from 0′ to 7′. Each panel represents the PSF corresponding to a specific off-axis angle, highlighting the variation in spatial distribution due to the telescope’s optics. The color scale indicates relative intensity, and the images are shown in a logarithmic scale for a better view of the PSF structure.

3. Results

3.1. Deconvolution of Simulated NuSTAR Images

The NuSTAR PSF was used to simulate two nearby point sources. The sources were placed at various distances and the deconvolution algorithm was tested for efficiency. An example is shown in Figure 2, where a 4′ PSF was used for both sources separated by 10″.

3.2. Deconvolution of the SXP 15.3 – SXP 305 Pair

Figure 3 demonstrates that deconvolution of the SXP 15.3/305 image enhances the visibility of fine structures, particularly in the region around the sources. The cropped image isolates the target and reduces background noise during computations. Then, image rotation helps to align the raw cropped image with the adopted input PSF image. After rotation, the alignment of features remains consistent, preserving thus the various morphological details. Finally the RL iteration procedure sharpens the image and reveals details that were previously blurred due to the instrument’s PSF.

In the resulting image (right panel in Figure 3), the maximum X-ray luminosities are not the same in the two flux-enhanced regions marked by boxes. This is never the case in NuSTAR observations of single sources, so this feature signifies the presence of two X-ray emitters in the cleaned image. The sources are however not entirely resolved owing to their small separation by about $7^{\prime \prime }$ for which the performance of the RL algorithm begins to degrade somewhat (simulations show strong performance for separations of $\ge {10}^{\prime \prime }$). Nevertheless, we believe that the brighter region corresponds to SXP 305 in outburst for which pulsations of 305 seconds were measured in this observation. Before this observation, only 15.3-second pulsations had been detected from this region, when SXP 15.3 was bright and SXP 305 was presumably in quiescence.

A comparison between the original, rotated, and deconvolved images of Figure 3 shows the effectiveness of th RL method in improving spatial resolution. The deconvolved image (right panel) can then be processed by other source detection algorithms (e.g., Gaussian fitting of the signal) to pinpoint the precise positions of the two sources that were unresolved in previous observations [27,28].

4. Discussion

We have demonstrated that utilizing advanced image deconvolution techniques, specifically the RL algorithm tailored to the NuSTAR PSF, we can significantly enhance the spatial resolution of NuSTAR X-ray images. By choosing the correct model of the telescope’s PSF and by incorporating it into the deconvolution procedure, it is possible to resolve closely-spaced X-ray sources that are indistinguishable by standard imaging techniques. Our application of the RL iterative method to the case of the HMXB pulsars SXP 15.3 and SXP 305 (right panel in Figure 3) illustrates this approach.

Observations of these nearly overlapping sources with NuSTAR have been challenged by the instrument’s broad and off-axis-dependent PSF, leading to blended images in which pulsations from one or the other source are detected at different epochs without the two sources ever being spatially resolved. By employing a refined PSF model accounting for energy dependence and off-axis variations inherent in the telescope’s optics, we have increased the effectiveness of the RL deconvolution algorithm.

The ability to resolve such X-ray sources has significant implications for astrophysical research. For HMXBs like SXP 15.3 and SXP 305, accurately determining their individual fluxes, spectral properties, and temporal behaviors is essential for understanding their accretion processes, magnetic fields, and evolutionary states [28]. In crowded pulsar fields (as in the SMC), where multiple X-ray sources may be located in close angular proximity, enhanced resolution facilitates a more precise characterization of the source population and their individual contributions to the X-ray luminosity of the entire region.

Furthermore, the improved detection capabilities enabled by refined PSF modeling and deconvolution can also aid the studies of transient phenomena. Detecting outbursts, pulsations, or other temporal variations from faint sources requires not only temporal resolution but also the ability to spatially isolate the sources in order to avoid contamination from nearby objects. By enhancing the spatial resolution, we improve the signal-to-noise ratio for individual sources, enabling thus the possibility of more sensitive timing analyses.

Despite the advancements in this field of research, certain challenges and limitations remain: (a) The RL deconvolution procedure is sensitive to inaccuracies in the adopted PSF model; any mismatch between the modeled PSF and the actual instrumental response can introduce artifacts or erroneous features in the cleaned image. This necessitates ongoing efforts to validate and update the existing PSF models using observations of bright isolated point sources across the detector field. (b) Noise amplification is an inherent risk in deconvolution algorithms, particularly in low-flux scenarios typical of faint (galactic and extragalactic) X-ray sources. Regularization techniques and careful selection of iteration parameters are required to mitigate this effect.

5. Conclusions

Our study reveals the significant impact that PSF modeling and RL deconvolution can have on the analysis of NuSTAR observations of X-ray point sources by enhancing the spatial resolution and enabling the detection of multiple sources in close proximity (≳ 10″). In a case study of pulsars SXP 15.3 and SXP 305, we have tested these tools in an investigation of a faint crowded X-ray field in the sky. This approach, not only improves our understanding of specific astrophysical objects, but also contributes to the broader capability of X-ray astronomy to explore the very high-energy universe with greater precision and clarity.

Future work should explore the integration of spatially varying PSFs directly into deconvolution algorithms, allowing thus for the continuous and seamless incorporation of off-axis signal variations. Bayesian deconvolution and machine learning approaches may offer alternative ways to handle complex PSF attributes and noise characteristics. Furthermore, applying the same methodology to a broader set of observations carried out at different wavelengths and targeting more crowded fields of high source density can help evaluate the effectiveness of the procedure and identify telescope-specific requirements.