SRG/ART-XC Galactic Bulge deep survey. I. Maximum likelihood source detection algorithm for X-ray surveys

ABSTRACT

We describe an X-ray source detection method entirely based on the maximum likelihood analysis, in application to observations with the ART-XC telescope onboard the Spectrum Roentgen Gamma observatory. The method optimally combines the data taken at different conditions, a situation commonly found in scanning surveys or mosaic observations with a telescope with a significant off-axis PSF distortion. The method can be naturally extended to include additional information from the X-ray photon energies, detector grades, etc. The likelihood-based source detection naturally results in an optimal use of available information for the sources detection and stable and uniform definition of detection thresholds under different observing conditions (PSF, background level). This greatly simplifies the statistical calibration of the survey needed to, e.g. obtain the $\log N-\log S$ distribution of detected sources or their luminosity function. The method can be applied to the data from any imaging X-ray telescope.

1 INTRODUCTION

The Mikhail Pavlinsky ART-XC telescope is one of the two X-ray telescopes onboard the Spectrum Roentgen Gamma observatory Sunyaev et al. (2021). The peak of its effective area is reached in 6–11 keV energy band. Its primary science objective is detection of the Galactic and extragalactic source populations in the hard band (Pavlinsky et al. 2019). The most frequently used mode of ART-XC operations is wide-field surveys done with scanning observations. In this paper, we describe the source detection algorithm developed specifically for ART-XC surveys, but which also is widely applicable to any imaging X-ray telescope.

The detailed description of the ART-XC telescope is given in Pavlinsky et al. (2021). Here, we briefly list its characteristics pertaining to our analysis. ART-XC consists of seven co-aligned grazing incidence Wolter I (Wolter 1952) telescopes. The full operating energy range is $4-35$ keV, but the telescope is most sensitive in the $4-12$ keV band. The field of view (FoV thereafter) is approximately circular with a radius of 18${}^{\prime }$⁠. The telescope uses CdTe semiconductor double-sided strip detectors with a spatial resolution of 45${}^{\prime \prime }$ (Levin et al. 2014). The particle-induced background is a dominating component and its distribution on the detector is nearly flat. It is highly stable due to SRG operating in the Earth–Sun second Lagrange point, except for periods of Solar flares (Pavlinsky et al. 2021).

The main limiting factor for the ART-XC source detection is a relatively wide point spread function (PSF thereafter) and a relatively high particle-induced background. The net effect of these two factor is that the source photons are spread over a large area and the signal should be detected against a significant level of Poisson noise. An additional complication for classical source detection methods is a strong variation of the ART-XC response (the PSF, mirror vignetting, background) within the field of view. Examples of the ART-XC PSF measured during ground tests (Krivonos et al. 2017) are shown in Fig. 1. The PSF size and shape strongly changes as the source is moved away from the optical axis. It is circular with a $\sim 30$${}^{\prime \prime }$ HPD in the central 9${}^{\prime }$ but then widens to $\sim 2{}^{\prime }$ and becomes strongly elongated towards the FoV edge (see fig. 13 in Pavlinsky et al. 2021, for more information). During the scans, the source spends $\approx$25 per cent of the time in the central FoV zone with a compact PSF. Approximately half of the source photons are collected in this zone. The other half of the source photons are collected when it is in the outer regions of the FoV with the poor PSF. Each source is therefore observed under very heterogeneous conditions, and there is no clear boundary one can use to simply filter out observing periods with poor data quality.

The optimal source detection should take all these changes into account instead of operating on a single, combined image as is typically done for pointed observations (see, e.g. for a recent review Masias et al. 2012). Our detection method does just that by combining the likelihoods computed individually at the time of registration of each photon, as well as by using additional information from the photon energy and detector grades. The basic information we use to identify sources is the map of the likelihood ratio test, which for each sky location provides a measure of probability that the observed photon distribution can be attributed entirely to the background. A high value of the likelihood ratio at a given sky location corresponds a low probability for the ‘background-only’ hypothesis, and hence implies an existence of a source. In many settings, this approach represents a statistically optimal test for detection or rejection of sources (Neyman & Pearson 1933).

The fully likelihood-based source detection approach has indeed been historically used for analysing data from the $\gamma$-ray observatories, which are characterized by a complex detector response and where one is often forced to work at lower levels of statistical significance (Mattox et al. 1996; Hartman et al. 1999; Abdo et al. 2010). For the imaging X-ray telescopes, the maximum-likelihood analysis has been applied in a more limited manner. It has been used for modelling of detected sources in the analysis pipelines, e.g. for ROSAT (Cruddace, Hasinger & Schmitt 1988), Swift(Evans et al. 2020), XMM–Newton (Gabriel et al. 2004), and the Chandra Source Catalogue (Evans et al. 2010). However, in all prior implementations, the maximum-likelihood fitting was the second stage of a two-step procedure in which the candidate sources were first detected by much cruder methods. For example, the Swift source catalogue (Evans et al. 2020) uses the sliding box detection while the Chandra Source Catalogue approach is equivalent to using an approximately Gaussian detection filter. The prime reason for using these simplified methods is a high computational cost of performing the likelihood-based analysis at all locations. We solved this issue by heavily optimizing the process of maximizing the likelihood function. As a result, our detection technique uses the likelihood-based approach at the crucial step of initial source detection. This naturally leads to other advantages such as the statistically optimal combination of data taken at non-uniform conditions (background, PSF, etc.) or even a combination of data from different telescopes.

A conceptually similar algorithm has been described previously by Ofek & Zackay (2018) and Lynx Team (2019), and applied to obtain the catalogue of a point sources from the first year of the ART-XC all-sky survey (Pavlinsky et al. 2022). This algorithm is based on the derivation of a ‘matched filter’ whose cross-correlation with the data provides a close approximation of the likelihood ratio test for sources near the pre-defined detection limit. In terms of completeness and purity of detection, this matched filter approach performs nearly identically to the full likelihood ratio method (see further discussion in Appendix  B). The critical improvement introduced in our approach is a full likelihood analysis at each sky location that not only established existence of a source but also returns its best-fitting flux. As a result, we automatically measure the fluxes of detected sources. Our method also leads to a stable and uniform definition of detection thresholds and estimation of the number of false detections.

While the bulk of our work is generic and applicable to the X-ray data from any imaging X-ray telescope (Chandra, XMM–Newton, Swift, NuSTAR, eROSITA), there are several ART-XC specific issues that we had to address. They include a procedure for estimating the local backgrounds and reducing the illumination from bright off-axis sources during the survey scans. Our solutions to these problems are instructive and are given below.

The paper is organized as follows. In Section 2, we describe the detection method. The mathematical basics are given in Section 2.1; the iterative maximum likelihood fit and flux measurement procedure are described in Section 2.2; our treatment of the detection thresholds in terms of the likelihood ratio is described in Section 2.3, estimation of the positional and flux count rate uncertainty is described in Section 2.5. In Section 3.1, we discuss how the illumination from bright sources is suppressed. The background estimation procedure for the ART-XC scans is detailed in Section 3.2.

2 SOURCE DETECTION

The ART-XC detector background is very stable, excluding periods of Solar flares (Pavlinsky et al. 2021) but relatively high (⁠$\sim 0.4$ cts s${}^{-1}$ per the 1020 arcmin${}^2$ field of view in 4–12 keV energy band in each telescope). The expected sensitivity of the typical ART-XC deep survey programs (a few ksec total exposure) is a $\mathrm{few}\times {10}^{-13}$ erg s cm${}^{-2}$ (Pavlinsky et al. 2019). The background contribution in the PSF beam is non-negligible for sources in this flux range and type of exposures: $\sim 42$ background versus $\sim 20$ sources counts detected within the 1 arcmin PSF beam for 3 ks. For the medium-exposure survey used as an example in this paper below, the sensitivity is $8\times {10}^{-13}$ erg s${}^{-1}$ cm${}^{-2}$⁠, and $N_{\mathrm{src}}=12$$N_{\mathrm{bg}}=11$ for a 800 s exposure. In the typical shallow regions of the all-sky survey, $t_{\exp }=100$ s, $N_{\mathrm{src}}=6$$N_{\mathrm{bg}}=1.4$⁠, and a typical sensitivity is $4\times {10}^{-12}$ erg s${}^{-1}$ cm${}^{-2}$⁠. The sensitivities are computed for the same detection threshold, which results in increasing of the number of the false sources 1.5 times between 100 s and 800 s effective exposure and 1.86 times between 100 s and 3 ks. Moreover, the majority of the ART-XC observations are done in the ‘scanning’ mode, which results in each source being exposed roughly uniformly within the field of view. This creates a wide mix of significantly varying PSF (see Fig. 1), and also the source-to-background ratios due to non-uniform vignetting. Under these conditions, only a highly optimal source detection method allows us to address the main scientific objectives of ART-XC of surveying large regions in the sky in the hard X-ray band (Pavlinsky et al. 2021).

As we discussed above, the likelihood-based analysis is the theoretically most optimal method for detecting faint sources. It is also the optimal approach for combining non-uniform data sets (Lehmann & Romano 2005), such as exposures of the same source with different PSF in our case. Therefore, we base our method on the likelihood function analysis.

The expected source density at a flux limit of ${10}^{-13}$ erg s cm${}^{-2}$ in the 4–12 keV band is 5.8 deg${}^{-2}$ in the extragalactic sky (Gilli, Comastri & Hasinger 2007) and $\approx 15$ deg${}^{-2}$ in the Galactic Ridge (Hands et al. 2004). This translates to 0.003 and 0.01 sources per the PSF area, respectively. Therefore, ART-XC is unaffected by source confusion except for a small area around the Galactic Centre, and its typical detection regime is that of finding isolated point sources. We provide a rigorous treatment of the likelihood-based source detection for such a case. In Section 2.1 below, we start the description of our method with the discussion of the likelihood calculation for a source in a given sky position.

2.1 Likelihood ratio test for existence of a source

The basic observable for an X-ray telescope is a set of events with a known arrival time, position on the detector, and energy. For a source with a given flux and sky position, we have a model for the count rate at this detector location and arrival time (see discussion near equation 3 below). Following the approach of Cash (1979), we can write the Poisson probability to observe a set of events in a single pixel, given a model, as follows:¹

where the observation is split into very short intervals $\delta t$ such that at most one event is detected in each interval; μ$(=\mu (t)$⁠) is the model count rate; and $n=1$ if the event is detected and 0 if not. With this fine binning, equation (1) can be rewritten as

where the product is over the detected set of events. For a set of independent pixels, the probability is a product of probabilities given by equation (2). It can be simplified as follows

The model count rate μ is, generally, a sum of the source (r) and background (b) models: $\mu =r+b$⁠, and so equation (3) can be written as

Let us now consider how the source model, r, can be factored (the background model is discussed in Section 3.2 below). Assuming that the source flux is constant over the observing period, we can write $r=Rs$⁠, where R is the count rate for a source located at the optical axis, and s can be written as

where $\mathrm{PSF}(\mathrm{\Delta }x,\mathrm{\Delta }y,\xi ,\eta ,E)$ is the position- and energy-dependent PSF, i.e. the probability of photon scattering to location $(\mathrm{\Delta }x,\mathrm{\Delta }y)$ relative to the nominal source location within the field of view, $(\xi ,\eta )$⁠; $V(\xi ,\eta ,E)$ is the vignetting factor at the nominal source location; and the integral is over one pixel, necessary if the pixel and PSF sizes are comparable as is the case for ART-XC. The time dependence in equation (4) is encoded in equation (5) through the time dependence of $\xi$ and $\eta$ that describes how the source moves through the field of view during scans.

From the probability of the observed data given the model (equation 4, 5), we can now construct the classical likelihood ratio test, broadly used for models comparison, for rejecting a case that the observed photons are entirely produced by the background (Wilks 1938):

where e is the exposure-averaged vignetting correction (cf. equation 5)

The application of the likelihood ratio test for simple nested hypothesis test (⁠$R=0$ versus $R\ne 0$ in our case) requires that the null values of the tested parameter may not be on the boundary of the set of its possible values. In many situations this may represent a serious problem (see extensive discussions in Protassov et al. 2002). However, this concern is not applicable in our case because in the presence of non-zero background solutions with negative R are mathematically permissible. Our iteration scheme described below does enforce the non-negative values of R, but this does not affect the basis for applicability of the likelihood ratio test.

The summation in equation (6) is formally done over all photons, but in practice only the photons sufficiently close to the given location contribute significantly. For photons with large $\mathrm{\Delta }x$⁠, $\mathrm{\Delta }y$⁠, $s\to 0$⁠, and $\ln ((Rs+b)/b)\to 0$⁠, also. Therefore we set the ART-XC PSF model to zero outside of the $5^{\prime }$ radius. This aperture contains 95 per cent of the total source flux near the FoV edge, and nearly 100 per cent of that on-axis. The photons that enter the likelihood calculations are then also selected within $5^{\prime }$ of the given position. This area contains $0.04(f/{10}^{-13}\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}{)}^{-1.5}$ sources on average (Gilli et al. 2007), hence the assumption of a single source implicitly used in our derivation is valid.

The value of R that maximizes the likelihood ratio in equation (6) can be found, as usual, by solving the equation

After taking derivatives, equation (8) can be rewritten as follows

Many existing source detection packages (see e.g. xmm sas Gabriel et al. 2004) implement an equivalent of the procedure described above, for a subset of trial source locations, identified by other means – usually, via simplified source detection methods such as sliding box or wavelet transform (e.g. Voges et al. 1999; Brunner et al. 2022, check also a more comprehensive background subtraction technique which can advance such initial analysis Ehlert et al. 2022). Clearly, this first step leads to an overall sub-optimal detection, although in the regime of very high background levels the matched filter becomes identical to the likelihood ratio test (Helstrom 1968; Vio & Andreani 2018). The prime reason for using these simplified method is the absence of analytical solutions to equation (9) and a high computational cost of solving equation (9) numerically because in general it has no analytic solution. The source rate solution can be straightforwardly performed with a family of gradient-free root finding solvers (Rios & Sahinidis 2013), which usually require more then 10–20 iterations to achieve a per cent level of convergence. Note that each evaluation of equation (9) involves a substantial amount of vector algebra, and if one aims to find the solution for R at each location covered by a data set, the computational costs accumulate. Below we describe a faster and highly accurate method that enables the implementation of the likelihood-based detection scheme for maximizing the likelihood on a desktop-class computer.

2.2 Iterative calculation of the likelihood map

We start with modifying equation (9) as follows:

This equation can be used as a basis for the iterative scheme,

where $R^{\prime }$ and $R^{\prime \prime }$ are the previous and new iterations of R, respectively. Let us show that this iterative scheme does indeed converge. Let R is the exact solution of equation 9 (and, equivalently, 11), $\delta R^{\prime }=R^{\prime }-R$ and $\delta R^{\prime \prime }=R^{\prime \prime }-R$ are the previous and new deviations of iterations from that value. Then (cf. equations 11 and 12),

Now note that per equation (11), $(eR{)}^{-1}\sum (R{s}_i)/(R{s}_i+b_i)=1$ and therefore $\delta R^{\prime \prime }<\delta R^{\prime }$ if $R^{\prime }>0$ and thus the iterations always converge. In practice, $<20$ iterations are needed to obtain a solution for R with a 1 per cent accuracy. We also checked that the convergence is largely independent on the initial value of $R^{\prime }$⁠.

Implementation of this method opens us a possibility of computing the map of the likelihood ratio over the entire observed region (all sky in the ART-XC case). Local maxima of a sufficient height in such a map correspond to real sources. Moreover, the value of R obtained by solving equation (11) is, automatically, the maximum likelihood estimate of the source countrate.

2.3 Detection thresholds in terms of $\mathbf{\Delta }\ln L$

Real sources correspond to local maxima in the likelihood ratio map. We need to set a reasonably high threshold to avoid accepting a large number of false sources, i.e. the source is detected if

This definition of the detection threshold leads to a robust and stable estimation of the number of false detections that can be done nearly from the first principles. Indeed, according to the Wilks theorem (Wilks 1938), the probability distribution of the quantity $2\mathrm{\Delta }\ln L$ should approach to the distribution ${\chi }_1^2$ for the case when there really is no source in the sky (⁠$R=0$⁠). This implies that that the distribution of $\mathrm{\Delta }\ln L$ should be independent of the background level. A direct Monte Carlo simulation shown in Fig. 2 a demonstrates that indeed the distribution of $\mathrm{\Delta }\ln L$ does not depend on the background and follows the ${\chi }_1^2$ law.

Note, however, that to characterize the probability of the false detection we need the distribution of the peak heights of the local maxima, not all $\mathrm{\Delta }\ln L$ values across the field. The $\mathrm{\Delta }\ln L$ map is obtained, effectively, from cross-correlation of the data and a kernel $\sim \ln (s/b+1)$ (Ofek & Zackay 2018, cf. also equation 6), which has a width approximately equal to that of the PSF. This makes the $\mathrm{\Delta }\ln L$ map smooth and its values not independent for close separations. The distribution of the local maxima peak heights in this case are most easily obtained through Monte Carlo simulations. The simulation results for a representative ART-XC survey are shown in Fig. 2b. Note that because the distribution of $\mathrm{\Delta }\ln L$ is background-independent, we can expect that the distribution of peaks is also not very sensitive to the background level, and this is indeed the case (cf. solid black and dash-dotted red lines in Fig. 2b).

Vio & Andreani (2018) discuss the difficulty of analytical computation of the peak height distribution function in the noise field with non-Gaussian distribution, and outline an approximate approach based on the order statistics method (Hogg, McKean & Craig 2013). In the limit of high peak values (or small probabilities of chance fluctuations), this approach gives the false alarm probability that equals the complimentary cumulative distribution function (CCDF) of the noise field, $\mathrm{\Phi }(z)$⁠, times the effective number of independent pixels in the image (cf. equation 17 in Vio & Andreani 2018), where z is the peak height. In our case, $z=\mathrm{\Delta }\ln L$ and $\mathrm{\Phi }(z)$ is the CCDF for the ${\chi }_1^2$ distribution, $\mathrm{\Phi }(z)=1-\mathrm{\Gamma }(2z,1/2)$⁠. This approach based on the order statistics is approximate, but we find that it can be modified as follows to very accurately describe the distribution of false detection observed in the Monte Carlo simulations.

Let us consider a grid of pre-defined locations on the $\mathrm{\Delta }\ln L$ map, separated by a distance larger than the size of the effective filter (⁠$\approx$the PSF size). The distribution of the maximal value on this grid can be precisely obtained with the order statistics approach. However, the true peak in the smoothed image will almost never coincide with the nearest grid point, and hence its value will be somewhat higher than the maximum on the grid. The true peak value in the image is then the maximum on the grid plus an additional random variable whose distribution reflects the random offset between the true peak and grid points, and how fast the $\mathrm{\Delta }\ln L$ map falls off from the points of its local maxima. The resulting distribution is a convolution of the probability density function given by the order statistics approach with a kernel skewed towards positive values. For a broad class of the CCDFs, this operation is approximately equivalent to simply shifting CCDF to the right by a constant log-offset, or, equivalently using $\mathrm{\Phi }(kz)$ with $k\lesssim 1$⁠. Indeed, we find that in our case, the function $1-\mathrm{\Gamma }(0.89\times 2\mathrm{\Delta }\ln L,1/2)$ provides an excellent description of the distribution of false detections in our Monte Carlo simulations (dashed line in Fig. 2b). We will explore this technique further in a follow-up paper (Semena & Vikhlinin, in preparation).

The stability with respect to the background is a very useful property as it allows us to greatly simplify the statistical calibration of the detection process in a survey. Estimation of the number of false detection still requires a Monte Carlo simulation, but the simulations are needed for different PSF conditions and not for different levels of the background. In the ART-XC survey case, there are multiple scans of any point in the sky and the resulting effective PSF is uniform. Therefore, a single set of simulations is required for statistical calibration of the survey.

The exact value of the $\mathrm{\Delta }\ln L$ threshold is determined from a tradeoff between the desired low number of false detections and sensitivity. For ART-XC surveys of a few tens of square degrees, a good choice is $\mathrm{\Delta }\ln L_{\mathrm{l}\mathrm{i}\mathrm{m}}=11.4$⁠, resulting in 0.025 false sources per square degree (Fig. 2). Note that for applications such as measuring the $\log N-\log S$ distribution of detected sources or their luminosity function, we need detection thresholds defined as a function of flux. Fortunately, a relationship between flux and $\mathrm{\Delta }\ln L$ is well defined (Fig. 3). We defer further discussion to a future work where we will present the $\log N-\log S$ of ART-XC sources.

2.4 Additional information from source spectra and grades

So far, we have considered a case of source detection in a single (broad) energy band. In principle, there is additional information that can be used to improve the detection sensitivity. For example, the source spectra are typically very distinct from that of the background (see Fig. 4). This difference can be incorporated into the likelihood model if one is aiming to detect sources of a specific spectral type. In addition, the background and real source events typically have a different detector grade distribution (see Section  A and description of the ART-XC detectors in Levin et al. 2014). Grades can significantly help with rejecting the background in X-ray detectors (see, e.g. Bailey et al. 1983; Lumb & Holland 1988). For ART-XC case, the relative frequency of different grades from the background and from real photons from the onboard calibration source are shown in the right panel of Fig. 4. Events with $\mathrm{grade}>8$ are dominated by the background and they are by default excluded from the analysis. However, there are non-negligible differences between the source and background grade distribution for $\mathrm{grade}\le 8$⁠, and this can also be incorporated in the likelihood model. Specifically, we can replace s and b in equations (5) and (4) with

E is the registered photon energy, G is the event grade, $f(E)$ is a normalized source spectrum convolved with the telescope response, $p_{\mathrm{S}}(G,E)$ is the detector grade distribution for ‘sky’ events, $V(\xi ,\eta ,E)$ and $\mathrm{PSF}(\mathrm{\Delta }x,\mathrm{\Delta }y,\alpha ,\beta ,E)$ are the energy dependent vignetting and PSF, respectively. For the background model,

where $f_{\mathrm{B}}(E)$ is the normalized background spectrum, $p_{\mathrm{B}}(G,E)$ is the detector grade distribution for background events, and B is the overall background intensity. For a source with a Crab-like spectrum, using the spectral and grade information for source detection leads to an increase of $\mathrm{\Delta }\ln L\approx 0.6$⁠, which is equivalent to an additional $\approx 1\sigma$ signal being added to the source detection process.

2.5 Source parameter uncertainties

The value of R that corresponds to the local and the maximum’s location are automatically the maximum likelihood estimation of the source flux and position, respectively. A small additional analysis is required to derive uncertainties on the source parameters, such as flux and position. They can be obtained, e.g. with a Markov chain sampling technique (we use an implementation described in Salvatier, Wiecki & Fonnesbeck 2015). This is a natural choice in our case because the entire detection approach is based on computing the likelihood (probabilities) given the source parameters (equation 4). Alternatively, one can use the Cash (1979) approach based on varying the source parameters around the best-fitting value until $\ln L$ decreases by a specific value; for a example a change of $\ln L$ by $-0.5$ corresponds to the boundaries of a single-parameter 68 per cent confidence interval. We note that based on this approach, the source position uncertainties can be obtained directly from the likelihood ratio map we use for source detection: a region around a local maximum within $-1.15$ of the peak value² corresponds to the 68 per cent CL uncertainty of the source position. We have verified via bootstrap approach that this method does result in the correct positional uncertainties.

3 ART-XC SPECIFIC ISSUES

Detection method outlined about is generic and can be applied to any X-ray imaging telescope. There are always however instrument-specific issues. In this section, we describe a treatment of two such issues in the ART-XC case.

3.1 Illumination from bright sources

The ART-XC PSF has strong wings, resulting in ‘haloes’ around bright sources extending to a distance of up to 1 degree. These haloes are formed due to single reflections in the mirror system of the telescope (see Nariai 1988; Buitrago-Casas et al. 2017; Pavlinsky et al. 2019; Buitrago-Casas et al. 2020). They are notoriously hard to model (but see e.g. Raimondi & Spiga 2015).

For a fixed position of the source relative to the optical axis of the mirror module, singly reflected photons occupy a limited area in the detector plane. During the scan, the affected areas move and change morphology with a speed different from that of the scan (see top panels in Fig. 5). At any given moment in time we can predict where the most affected areas are and mask those areas from further analysis. This approach strongly suppresses the effect of scattering haloes on detection of other source and at the same time does not lead to a complete loss of data at any sky position.

Based on observations of bright sources during the calibration phase, we created a library of illumination masks from bright point sources for a set of offsets from the telescope optical axis, ranging from 5 to 85 arcmin. The masks were created such that the illumination is suppressed by a factor of 15 at distances in a range of $15-50{}^{\prime }$ from the bright source (see Fig. 5). Such filtering results in a $\sim 50\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}$ exposure loss in a $\sim {30}^{\prime }$ region around bright sources (Fig. 5).

3.2 Background model

The ART-XC background consists of three main components: the detector particle-induced background, diffuse radiation from the sky, and the residual halo illumination from bright sources (Section 3.1). Except for the regions of bright diffuse emission in the Galactic Centre region and the Galactic ridge, by far the dominant background component is the particle-induced one (⁠$\sim 99.5\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}$ of the total). It is very stable (except for short episodes of solar flares, see Pavlinsky et al. 2022) and had median value of $2\times {10}^{-4}$ s${}^{-1}$ per pixel in 4–12 keV energy band during the first two years of operation. Large background data sets were accumulated during ART-XC scanning observations in the source-free regions. With a small renormalization using the observed background intensity above $\sim 40$ keV,³ these data sets can be used to accurately predict the particle-induced background at all locations on the detector and throughout the useful energy band. The renormalization accuracy is better than 1 per cent for time intervals of $1$ ksec, which is much shorter than the typical time scale of the secular changes in the particle-induced background.

The second component we need to model is the astrophysical diffuse emission. In practice, it is significant only in the Galactic Centre and Ridge regions. These are also the regions densely populated with bright point sources. Despite the application of the illumination suppression algorithm (see the 3.1 paragraph), a certain fraction of their extended halo remains. A typical angular scale of this residual emission is $\sim 1{}^{\circ }$⁠. We do not distinguish these residual haloes from the truly diffuse Galactic emission or from emission from a large number of point sources well below the ART-XC detection threshold (cf. Revnivtsev et al. 2009). The diffuse sky background model for ART-XC is derived as follows.

First, we subtract the particle-induced background model from the observed photon map. This map is flat-fielded using the exposure map including mirror vignetting. Then we remove small-scale features associated with detectable sources using a wavelet decomposition method with an implementation similar to that described by Vikhlinin et al. (1998). With wavelet decomposition we remove all detectable structures on spatial scales of 20, 40, 80, and 160${}^{\prime \prime }$ are discarded, i.e. comparable to the width of the ART-XC PSF. The residual map is smoothed with a Gaussian filter with a width of 10${}^{\prime }$⁠. This smoothing scale was chosen because it sufficiently suppresses the photon counting statics noise, while at the same time remains small compared to the typical scale of the Galactic background. The resulting map gives an estimate of the diffuse background count rate from a given region in the sky, which is used – together with the particle-induced background map – in the likelihood ratio calculations.

4 SUMMARY AND CONCLUSIONS

Putting all components of our source detection pipeline together, we show in Fig. 6 the application of our algorithm to the data in the ART-XC ‘L20’ field (a $5{}^{\circ }\times 4{}^{\circ }$ field in the Galactic Plane region centred at $l=20{}^{\circ }$⁠, $b=0{}^{\circ }$⁠, also check Semena et al. 2023 (Paper II), for the analysis of Galactic Centre field, populated with several bright X-ray sources). The top panel shows the total photon image correlated with the field of view averaged PSF. This is the type of data used in the traditional source detection pipelines. In the middle, we demonstrate the photon image with the background model subtracted, correlated with the average PSF and divided by the exposure map. This map shows close resemblance to the likelihood map despite not being exactly applicable for the purpose of the source detection. Notice significant exposure and background non-uniformities across the field in both cases. As a result, the surface brightness peaks at the location of detected sources (marked by red circles) are not very prominent. Obviously, a confident source detection using such data is possible only at a very high level of detection threshold, and there will be significant variations of the threshold value across the field.

The situation is drastically better in the $\mathrm{\Delta }\ln L$ map (Section 2.1), in the sense that the fluctuations in the source-free regions in this map are uniform (as expected, Section 2.3) and have a much lower amplitude compared to the peaks corresponding to sources. Twenty five sources are clearly detected above a threshold of $\mathrm{\Delta }\ln L>11.4$⁠. Less than one false detection is expected in 20 deg${}^2$ for this choice of threshold (Section 2.3), so most if not all of these detections are real sources (Karasev et al. 2023). For comparison, the PSF convolution algorithm produces 110 false detection for the same flux limit. Alternatively, if we set the threshold for the PSF convolution such that there is only one false detection, then only 14 of 25 sources are detected. This comparison shows that for the ART-XC data, we are far from the high-background regime where the likelihood-based approach and PSF convolution become equivalent. The improvements in sensitivity of the source detection using the $\mathrm{\Delta }\ln L$ compared with the PSF-smoothed image are due to the maximally optimal use of available information (Section 2.1). Primarily, the main improvement is due to using separately the data with the sharp PSF on-axis and those with the degraded PSF off-axis. Additional small improvements in sensitivity are due to using photon energies and grades (Section 2.4).

To conclude, we presented an X-ray source detection method based on the likelihood analysis and demonstrated the performance using the ART-XC data. Our method results in uniform detection conditions even using very non-uniform data (PSF, background level). This stability greatly reduces the effort required for further statistical calibration needed, e.g. to obtain the $\log N-\log S$ or the source luminosity function. We will present this analysis in application to ART-XC in a future paper. In the next paper in this series (Paper II), we present a catalogue of sources detected using this method in a deep Galactic Bulge survey with ART-XC.

ACKNOWLEDGEMENTS

This work is based on observations with Mikhail Pavlinsky ART-XC telescope, hard X-ray instrument onboard the SRG observatory. The SRG observatory was created by Roskosmos (the Lavochkin Association and its subcontractors) in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI) in the framework of the Russian Federal Space Program, with the participation of Germany. The ART-XC team thanks the Russian Space Agency, Russian Academy of Sciences and State Corporation Rosatom for the support of the ART-XC telescope design and development. The science data are downlinked via the Deep Space Network Antennae in Bear Lakes, Ussurijsk, and Baykonur, funded by Roskosmos. Alexander Vikhlinin’s work on this project was partially supported by the Chandra X-ray Observatory grant AR1-22013X.

DATA AVAILABILITY

At this moment SRG/ART-XC data and corresponding data analysis software have a private status. However, we plan to provide public access to the ART-XC scientific archive in the future.

Footnotes

REFERENCES

APPENDIX A: ART-XC GRADES

The ART-XC possess semiconductor double-sided strip detectors, which have their own specifics in events grades determination. Each event is stored in the telemetry with six amplitudes from three top and three bottom strips centred at event. The grade is determined by the comparison of these amplitudes to the specified signal thresholds. This provides us with sixteen independent grades, which could be classified as single (grade = 0), double (i.e. with two adjacent stripes in one or both surfaces, 0 < grade < 8) or multiple (grade > 8). Multiple events are nearly impossible to produce with X-ray photon, therefore they are dominated by cosmic ray hits.

APPENDIX B: ON THE STATISTICAL POWER OF THE LIKELIHOOD RATIO TEST FOR SOURCE DETECTION

The likelihood ratio test (LRT) we use in this work is often assumed to be the most optimal method for source detection. The theoretical basis for this assumption is the Neyman–Pearson Lemma (Neyman & Pearson 1933) which proves that for simple hypothesis testing (⁠$R=0$ versus $R=R_1$ in our case) the likelihood ratio approach does indeed provide the uniformly most powerful (UMP) test. The Karlin–Rubin theorem (Karlin & Rubin 1956) proves an extension of the UMP nature of the LRT for non-simple hypotheses (e.g. $R=0$ versus $R>0$ in our case) in certain cases, e.g. when the likelihood is in fact a monotonic function of a scalar combination of observed values. We also note that the original Neyman & Pearson (1933) paper considered an extension of simple hypothesis comparison to the set of alternatives in the case where the critical regions of observables are bounded by the single family of non-intersecting multidimensional surfaces. In our case, the source detection in the images is a non-simple hypothesis test, the conditions for Karlin–Rubin theorem are not met, and we were unable to prove whether or not the critical regions considered by Neyman & Pearson (1933) are bounded by a single family of surfaces. The question, therefore, is in order on how close is the likelihood ratio approach implemented in this paper to the most powerful test for source detection? In other words, how close to the theoretically possible maximum is the probability of detecting real sources given a desired rate of false detections?

We have performed a study to this effect using Monte Carlo simulations. The starting point is the notion that the optimal filter detection method proposed by Ofek & Zackay (2018) does correspond by construction to the condition of the Neyman–Pearson Lemma for detecting sources with a specified flux (the flux for which the optimal filter is constructed). The optimal filter is then the UMP test for detecting such sources. Therefore a comparison between the optimal filter and our method gives us an idea for how close we are to theoretical maximal performance.

To implement this analysis, we simulated a large number of images containing only background or background plus a source of a given flux. The background level and the source PSF model were chosen to be typical for the ART-XC surveys with 1000 s vignetted exposure. These simulated data were processed with our pipeline and also cross-correlated with the optimal filter contrasted for this flux. The thresholds in both methods were set such that rate of false detections was identical (this rate was close to that expected for $\mathrm{\Delta }\ln L=7$⁠). For these threshold levels, we then compared the completeness for detecting real sources by the optimal filter and our LRT-based method. We found the completeness of the two methods to be identical within 0.15 per cent (which in turn was within the statistical uncertainty of the simulation). We thus conclude that the LRT for the specific case of the image source detection does indeed approach the ‘optimality’ in terms of Neyman & Pearson (1933) lemma.