Observable-conditioned ΔC theory for Poisson PSF-fitting source detection (v30.6)

Rui Huang

Author B

Author C

Abstract

Poisson source-detection pipelines report a likelihood-based detection statistic for every candidate source — for example $\DETML$ in the XMM-Newton SAS — and use that number for catalog thresholding and sensitivity-map construction. Near the threshold, however, the relevant object is the sampling distribution of the underlying Cash-statistic likelihood-ratio improvement $\Delta C$, not merely the reported scalar statistic. We derive an observable-conditioned theory for the conditional mean of $\Delta C$ given the fitted source amplitude, the local background, the position-dependent PSF and the detection mask in the standard three-parameter single-source PSF fit, and use a Fisher projection to approximate the high-signal variance. The pipeline statistic is treated only as the final incomplete-gamma transform of $\Delta C$. We validate the framework on the XMM-Newton PN band-4 standalone emldetect chain. In self-consistent reference-fitter simulations spanning 17 detector positions and $\sim 8\times 10^5$ fits, the theoretical mean matches the simulated mean to within 0.07$\sigma$ in the high-signal regime, and the Fisher-projected variance reproduces the simulated variance to within ${\sim}5\%$. In the operational threshold band ($\DETML=6$–10), the asymptotic theory carries a known low-signal bias while the standalone emldetect adds an implementation-dependent residual. We absorb both into an empirical, position-dependent residual calibration, yielding a mean/median threshold fitted-amplitude calibration at $\DETML=6$ with an 11% uncertainty for the preferred 2D interpolation (17% for the nearest-calibration-position fallback), and a validated weak background dependence over $b=0.02$–$0.08\,\mathrm{counts\,pixel^{-1}}$. The resulting framework converts the $\DETML$ selection into a probabilistic mapping between fitted amplitude and detection statistic. It can be used backwards — to interpret traditional $\DETML$-selected catalogs together with a source-count prior — and forwards — to construct fitted-amplitude threshold-crossing cuts and percentile sensitivity maps; a production selection function requires additional validation of the dispersion and tail models.

Introduction

Source detection in Poisson-count photon imaging is usually done by fitting a parameterized source model on a known background and summarizing the result with a likelihood-ratio statistic. Under the standard Cash statistic formulation (Cash 1979), the natural statistical object is the Cash improvement $\Delta C = C_\mathrm{null} - C_\mathrm{best}$, where $C$ is proportional to twice the negative Poisson log-likelihood. Detection pipelines typically report a monotonic transform of $\Delta C$ as the “detection likelihood”; for example, in the XMM-Newton Science Analysis System (SAS), emldetect (Cruddace et al. 1988; Watson et al. 2009) fits a three-parameter source model — one free amplitude and two free position coordinates — to Poisson-count images and reports DET_ML, the incomplete-gamma transform of the underlying $\Delta C$. The reported detection-likelihood statistic is operationally important: it is used to decide which sources enter the catalog, to define detection thresholds and to build sensitivity maps.

Sensitivity maps and source-count methodology have been extensively developed and applied in major X-ray surveys (Cappelluti et al. 2009; Georgakakis et al. 2008; Laird et al. 2009; Luo et al. 2017, e.g.,; Mateos et al. 2008; Puccetti et al. 2009; Wang et al. 2016; Xue et al. 2011). The statistical treatment of Poisson source detection, limits and likelihood-ratio methods also has a long history in high-energy astrophysics (Broos et al. 2010, e.g.,; Gehrels 1986; Kashyap et al. 2010; Kraft et al. 1991; Protassov et al. 2002; Starck et al. 2002). The emldetect tool has been the standard for XMM-Newton serendipitous survey catalogs (Carrera et al. 2007; Rosen et al. 2016; Traulsen et al. 2019; Watson et al. 2009; Webb et al. 2020).

We emphasize at the outset that the $\DETML$ used in this paper is strictly an operational monotonic threshold scale inherited from the SAS convention. We do not interpret it as a calibrated frequentist false-alarm probability. The standard $\chi^2_3$ null distribution that supports the incomplete-gamma transform assumes internal regularity conditions (Chernoff 1954; Protassov et al. 2002; Self & Liang 1987) that are violated at the parameter-space boundary ($S=0$) and when the position parameters are unidentifiable under the null hypothesis. The applicability of that null distribution is itself contested: Stewart (2009) argue from Monte Carlo experiments that, under typical X-ray instrument conditions, the Cash statistic null distribution is quite close to the corresponding $\chi^2$ provided candidates with negative fitted amplitude are discarded from the sample, and on that basis reject Protassov et al. (2002)’s claim that likelihood-ratio statistics are unsuitable for source detection. More generally, the distribution of the Cash statistic itself departs from its high-count $\chi^2$ limit under low-count Poisson conditions (Bonamente 2020). This paper does not enter that dispute: our use of the $\DETML$–$\Delta C$ conversion is purely operational: it maps the detection threshold chosen by the survey analyst to the corresponding $\Delta C$ value, where the observable-conditioned theory applies.

The practically relevant threshold regime is low signal. For many XMM-Newton applications, the relevant sensitivity-map range is $\DETML$ of order 6–10. This is the regime where source counts are steep, Poisson fluctuations matter, and small changes in the detection threshold can change the inferred limiting flux. It is therefore also the regime where Eddington bias becomes important: random upward fluctuations are preferentially selected near the detection boundary, so the selected sample no longer represents the untruncated parent distribution (Eddington 1913; Wang 2004).

The core challenge is not the absence of a predictable mean. The expected $\Delta C$ at the fitted source has long been computable analytically and has been used in X-ray sensitivity-map construction as a single deterministic value. The direct origin of this idea is Stewart (2006, 2009): in a systematic comparison of the sensitivity of various source-detection methods under Poisson noise, Stewart (2009) give a closed-form approximation for $\langle U_\mathrm{Cash}\rangle$ in the single-parameter case where only the amplitude is fitted — after a sign conversion it is identical to the single-parameter limit of Eq. ([eq:mu3p]) in this paper (plus an empirical constant offset, to which we return in Section 4.5) — and explicitly propose it for Cash-type source-detection sensitivity-map construction. Inside SAS itself, from version 4.0 onward, esensmap carries an emldetect-style analytic Cash-statistic model, following the emldetect $\Delta C$ form of the 3XMM and 4XMM catalog pipelines (Rosen et al. 2016; Traulsen et al. 2019); outside SAS, the same expectation is evaluated at the detected (rather than injected) source flux to construct sensitivity maps for the New-ANGELS XMM legacy survey (Huang et al. 2025, appendix).

These prior works leave a gap in two respects. First, the conditioning variable they use is the model (injected) source amplitude, not the fitted amplitude actually reported by the detection pipeline. Second, and more fundamentally: although Stewart (2009) explicitly discuss the probability density $p(U;\alpha)$ of the detection statistic and weigh the mean versus the median for defining sensitivity, X-ray sensitivity-map practice ultimately used only the mean of that distribution; the width and shape of the distribution were never analytically modeled, validated by instrument position, or calibrated for a specific detection pipeline. Correspondingly, we find no precedent for using the conditional dispersion of $\Delta C$ in catalog-level threshold-crossing calculations. At the level of the moments of $\Delta C$, the closest work comes from spectral analysis rather than source detection: Kaastra (2017) give closed-form expressions for the expectation and variance of the C statistic under the null model, for goodness-of-fit testing.

It is, however, the distribution itself — not merely the mean — that controls the Eddington bias in the threshold regime and the precision of any downstream flux limit or source-count correction. We therefore work at the distribution level. Under standard regularity conditions, the conditional distribution of $\Delta C$ at the fitted source has a clean character — for a three-parameter single-source fit it is a non-central chi-square distribution with $\nu=3$ degrees of freedom, well approximated by a Gaussian in the high-signal limit — and both its mean and variance can be computed in closed form. This methodological lineage has been systematically developed in high-energy physics: based on the results of Wilks and Wald, Cowan et al. (2011) give the asymptotic distribution of the discovery test statistic $q_0$ under the signal hypothesis and introduce the Asimov data set as fluctuation-free reference data, thereby obtaining the expected sensitivity and its fluctuations analytically without Monte Carlo. The derivation in Section 3 of this paper is structurally an expansion of that framework to pixelized Poisson imaging: the Asimov mean gives the conditional expectation and the Fisher information gives the variance, the difference being that here one must explicitly handle the pixelized PSF, the fit mask, and the fluctuation directions absorbed by the free position parameters. The pipeline-reported detection-likelihood statistic, however, is a nonlinear incomplete-gamma transform of $\Delta C$ and is therefore not Gaussian; truncating the reported statistic near the threshold selectively truncates a non-Gaussian distribution, and that is precisely the regime where sensitivity maps and source-count corrections are most susceptible to Eddington bias.

In self-consistent reference-fitter simulations spanning 17 PN band-4 detector positions and including about $8 \times 10^5$ fits, the predicted conditional mean and Fisher-projected variance match the simulation results to within 0.07 $\sigma$ and about 5% respectively in the high-signal regime (Section 4). The working threshold band lies at lower fitted signal, where the asymptotic theory carries a residual low-signal bias and the standalone SAS emldetect introduces an implementation-dependent residual relative to the reference fit; both are absorbed into an empirical, position-dependent residual calibration (Section 5). The resulting sensitivity-limit equation (Section 6) reproduces the ground-truth limiting amplitude at $\mathtt{DET\_ML}=6$ to within 11% (preferred 2D interpolation) and 17% (nearest-neighbor fallback), and shows a weak background dependence over the tested background range $b=0.02\text{--}0.08$.

The observable-conditioning principle is set out in Section 2 and the three-parameter $\Delta C$ theory is derived in Section 3. Sections 4–6 cover reference-fitter validation, standalone emldetect threshold calibration, and the sensitivity-map recipe respectively. Section 7 lists the validated and unvalidated domains of this result.

Observable-Conditioned Residuals

The central statistical object is the conditional expectation of $\Delta C$ in observable space,

$$E[\Delta C | \widehat{S}, b, \mathrm{psf}, \mathrm{mask}],$$

where $\widehat{S}$ is the fitted source amplitude, $b$ is the local background in counts per pixel, $\mathrm{psf}$ is the pixelized PSF template used by the fit model, and $\mathrm{mask}$ is the image region included in the fit. Section 3 derives this expectation analytically; Section 5 turns it into an operational standalone emldetect threshold calibration, and Section 6 turns it into a sensitivity-map recipe. Figure 1 summarizes this logic.

Why the conditioning variable is S_hat, not S_true

In simulations, the injected source strength $S_\mathrm{true}$ is known and gives cleaner diagnostic plots than $\widehat{S}$, but it is not an available conditioning variable for detection theory: $S_\mathrm{true}$ is absent from every real image and every emldetect source table, and it is unavailable when constructing sensitivity maps. Without injecting information that does not exist in the use case, a theory conditioned on $S_\mathrm{true}$ cannot be applied to real detection or false-source recovery tests. Throughout this paper, $S_\mathrm{true}$ denotes only a simulation-generation label, while $\widehat{S}$ is the conditioning variable. Results conditioned on $S_\mathrm{true}$ test whether the model behaves sensibly at a known injected flux; results conditioned on $\widehat{S}$ test whether it is usable after the detection fit. Only the latter is relevant to standalone emldetect sensitivity-map work.

Residual decomposition

For each source or false-source trial, the theoretical prediction is evaluated at the observed fitted quantities:

$$\mu_{3p,i} = E_{\mathrm{theory}}[\Delta C | \widehat{S}_i, b_{i}, \mathrm{\mathrm{psf}}_i, \mathrm{\mathrm{mask}}_i],$$

The residual is:

$$R_{i} = \Delta C_{i} - \mu_{3p,i}$$

For a sample or bin:

$$\begin{aligned} \mathrm{mean}(\Delta C) &= \mathrm{mean}(R) + \mathrm{mean}(\mu_{3p}), \\ \mathrm{Var}(\Delta C) &= \mathrm{Var}(R) + \mathrm{Var}(\mu_{3p}) + 2 \mathrm{Cov}(R, \mu_{3p}) \end{aligned}$$

Only when $\mu_{3p}$ is approximately constant across the bin does the observed $\Delta C$ width reduce to $\mathrm{Var}(R)$; otherwise, mixed-$S$ or mixed-position samples introduce apparent variance effects that are not properties of the conditional theory. Single-$S$ residual bins show $\sigma^2_{3p}$ to be nearly correct on-axis and low at off-axis positions.

Thresholds are translated to DeltaC, not selected on observed DET_ML

DET_ML is a nonlinear monotonic tail-probability transform of $\Delta C$,

$$\mathrm{DET\_ML} = -\ln(\mathrm{gammaincc}(\nu/2, \Delta C/2)),$$

For a single-source three-parameter fit, $\nu=3$ and $\mathrm{shape}=1.5$. A Gaussian residual in $\Delta C$ is therefore no longer Gaussian in DET_ML, and selecting on the observed DET_ML means selecting on the very noisy quantity whose distribution we are measuring — it truncates the residual distribution and biases its mean upward near the threshold. Threshold ranges are therefore translated to $\Delta C$:

$$\begin{aligned} \mathrm{DET\_ML}&=6 \rightarrow \Delta C=14.338539, \\ \mathrm{DET\_ML}&=8 \rightarrow \Delta C=18.570072, \\ \mathrm{DET\_ML}&=10 \rightarrow \Delta C=22.755959, \end{aligned}$$

and used as the theory or predicted $\Delta C$ conditioning (for example $\mu_{3p,\rm ref}\in[\Delta C_\mathrm{ML6},\Delta C_\mathrm{ML10}]$), or as the $\widehat{S}$ bins in which the typical predicted $\Delta C$ falls in that interval. This makes the conditioning independent of the observed random fluctuations in $\Delta C_\mathrm{eml}$.

Three-Parameter DeltaC Theory

This section gives the closed-form derivation of the conditional mean and approximate variance of $\Delta C$ for the standard single-source PSF fit (one free amplitude and two free position coordinates).

Pixelized Poisson model

Consider a fit region containing pixels indexed by $i$. The observed counts are independent Poisson random variables,

$$Y_{i} \sim Poisson(\mu_i)$$

For a single source on a locally known background, the fitted source model is

$$\begin{aligned} \\mu_i(\theta) &= b_{i} + S p_{i}(x, y), \\ \theta &= (S, x, y) \end{aligned}$$

Here $b_i$ is the background expectation in pixel $i$, $S$ is the source amplitude in the same units as the image model, and $p_i(x,y)$ is the pixelized PSF template translated to position $(x,y)$. In the current PN band-4 validation runs, $b_i$ is usually a constant $b$ within the fit region, but the notation allows the background to vary per pixel.

The binary fit mask $M_i$ selects the pixels included in the fit:

$$\begin{aligned} M_{i} &= 1 \quad \text{included in the fit} \\ M_{i} &= 0 \quad \text{excluded from the fit} \end{aligned}$$

The mask represents the adopted fit footprint, for example the ‘ecut=15‘ pixel cut used in the standalone emldetect threshold calibration. A weighted-mask extension exists for fractional masks, but the main results of this paper use binary masks.

Cash statistic and DeltaC sign convention

Ignoring constants independent of the model parameters, the Cash statistic is

$$C(\theta) = 2 \sum_i M_{i} [ \\mu_i(\theta) - Y_{i} \log \\mu_i(\theta) ]$$

The null model is the background-only model,

$$\mu_{i,0} = b_{i}$$

The likelihood-ratio Cash improvement is

$$\Delta C = C_{\mathrm{null}} - C_{\mathrm{best}}$$

Under this sign convention, a better source fit gives a positive $\Delta C$. At a fixed fitted model $\mu_i$, the statistic can be written as

$$\Delta C = 2 \sum_i M_{i} [ b_{i} - \mu_i + Y_{i} \log(\mu_i / b_{i}) ]$$

This is the quantity later converted to the SAS ‘DET_ML‘. This paper never uses the catalog-column convention (in which the stored auxiliary quantity may represent $\Delta C/2$); the theoretical $\Delta C$ here is the full Cash difference $C_\mathrm{null} - C_\mathrm{best}$.

Conditional mean at fixed fitted result

The observable-conditioned theory evaluates the moments at the fitted result

$$\begin{aligned} \widehat{\theta} &= (\widehat{S}, \widehat{x}, \widehat{y}), \\ \mu_i &= b_{i} + \widehat{S} p_{i}(\widehat{x},\widehat{y}) \end{aligned}$$

The Asimov image of the fitted model is $Y_i = \mu_i$, the reference dataset that suppresses all statistical fluctuations at the fitted model (Cowan et al. 2011). Substituting it into the Cash improvement gives the conditional mean

$$\mu_{3p}(\widehat{S},b,\mathrm{psf},\mathrm{mask}) = 2\sum_i M_i\left[\mu_i\ln\left(\frac{\mu_i}{b_i}\right)-(\mu_i-b_i)\right] \label{eq:mu3p}$$

If the position parameters are frozen in this expression and only the amplitude degree of freedom is kept, it reduces to the sensitivity-mean expression given by Stewart (2009) for single-parameter Cash detection (whose form additionally contains an empirical constant offset, to which we return in Section 4.5). The generalizations of this paper are three: the conditioning variable is the fitted amplitude $\widehat{S}$ rather than the model amplitude, the fit mask $M_i$ is carried explicitly, and the projection term of the variance for the two free position parameters is supplied below.

For a constant background $b_i=b$ and a normalized PSF $p_i$, this is the expression implemented in deltac_moments_3param:

$$\begin{aligned} \mu_{3p} \\ &= 2 \sum_i M_{i} [ (b + \widehat{S} p_{i}) \log((b + \widehat{S} p_{i})/b) \\\ - \widehat{S} p_{i} ]\\ \end{aligned}$$

The mean depends on the pixelized PSF and the mask through the full pixel sum. It is not only a function of the Gaussian width or low-order PSF moments, although those summary quantities can serve as useful diagnostics.

Linearized score constraints

The best-fit parameters satisfy the score equations. For any fitted parameter $\theta_a \in (S,x,y)$,

$$dC/d\theta_a = 0$$

Linearizing around the fitted model, write

$$\begin{aligned} Y_{i} &= \mu_i + \eta_i, \\ E[\eta_i] &= 0, \\ \mathrm{Var}(\eta_i) &= \mu_i \end{aligned}$$

The first-order score constraints are

$$\sum_i M_{i} (d\mu_i/d\theta_a) \eta_i / \mu_i = 0$$

There are three such constraints:

$$\begin{aligned} d\mu_i/dS &= p_{i}, \\ d\mu_i/dx &= \widehat{S} dp_i/dx, \\ d\mu_i/dy &= \widehat{S} dp_i/dy \end{aligned}$$

These constraints are why the three-parameter fit is not equivalent to a simple unconstrained Poisson sum. Fluctuations along the fitted amplitude and position directions are absorbed by the fit and must be projected out when computing the residual variance.

DeltaC residual variance before projection

At the fixed fitted model, the first-order stochastic part of $\Delta C$ is

$$R_{\mathrm{linear}} = 2 \sum_i M_{i} \eta_i \log(\mu_i / b_{i})$$

If the fitted parameters were not constrained by the score equations, the variance of this linear residual would be

$$\begin{aligned} var_unconstrained \\ &= 4 \sum_i M_{i} \mu_i \log^2(\mu_i / b_{i}) \end{aligned}$$

This term is useful as a diagnostic, but it is not the final variance formula for the fitted likelihood ratio. It ignores the fact that $S$, $x$ and $y$ were chosen to maximize the likelihood.

Fisher projection for three fitted parameters

Define the Fisher information matrix of the three fitted parameters:

$$I_{ab} = \sum_i M_{i} (d\mu_i/d\theta_a)(d\mu_i/d\theta_b) / \mu_i$$

For $\theta=(S,x,y)$, its elements are

$$\begin{aligned} I_{SS} &= \sum_i M_{i} p_{i}^2 / \mu_i \\ I_{Sx} &= \sum_i M_{i} \widehat{S} p_{i} (dp_i/dx) / \mu_i \\ I_{Sy} &= \sum_i M_{i} \widehat{S} p_{i} (dp_i/dy) / \mu_i \\ I_{xx} &= \sum_i M_{i} \widehat{S}^2 (dp_i/dx)^2 / \mu_i \\ I_{xy} &= \sum_i M_{i} \widehat{S}^2 (dp_i/dx)(dp_i/dy) / \mu_i \\ I_{yy} &= \sum_i M_{i} \widehat{S}^2 (dp_i/dy)^2 / \mu_i \end{aligned}$$

The covariance between the $\Delta C$ residual and the score directions $\theta_a$ is encoded by

$$g_{a} = \sum_i M_{i} \log(\mu_i/b_{i}) (d\mu_i/d\theta_a)$$

Specifically,

$$\begin{aligned} g_{S} &= \sum_i M_{i} p_{i} \log(\mu_i/b_{i}) \\ g_{x} &= \sum_i M_{i} \widehat{S} (dp_i/dx) \log(\mu_i/b_{i}) \\ g_{y} &= \sum_i M_{i} \widehat{S} (dp_i/dy) \log(\mu_i/b_{i}) \end{aligned}$$

Under the Gaussian/Fisher approximation, conditioning on the fitted score constraints subtracts the projection of the residual onto the fitted-parameter subspace. This is the imaging realization of the Wald-type asymptotic argument (Cowan et al. 2011): the fluctuation directions of the fitted parameters are characterized by the Fisher information, and the components of the likelihood-ratio residual along those directions are absorbed by the fitting process. The resulting variance is

$$\sigma^2_{3p} = 4\left[\sum_i M_i\mu_i\ln^2\left(\frac{\mu_i}{b_i}\right) - g^T I^{-1}g\right] \label{eq:var3p}$$

This is the fitposition=True formula in lib/fit_mask_3param.py. If the position directions are removed and only $S$ is fitted, the formula reduces to the amplitude-only single-parameter expression.

Schur complement implementation

The code uses a block inverse to evaluate the three-parameter projection. Write

$$\begin{aligned} I &= [[I_{SS}, H], \\ [H^T, K]] \end{aligned}$$

where $H$ contains the cross-terms between $S$ and the positions, and $K$ is the $2\times2$ position block for $x$ and $y$. Then

$$\begin{aligned} g^T I^{-1} g \\ &= (g_{S} - H K^{-1} g_{pos})^2 / (I_{SS} - H K^{-1} H^T) \\\ + g_{pos}^T K^{-1} g_{pos} \end{aligned}$$

This is algebraically equivalent to inverting the full $3\times3$ Fisher matrix, but it makes the position projection explicit and is numerically convenient for handling large grids of $S$ and $b$ values.

Mask and weighted-mask conventions

For binary masks, the sum above runs only over pixels with $M_i=1$. This is the main convention used in the current PN band-4 standalone emldetect calibration.

For fractional weights $w_i$, the weighted statistic is

$$\Delta C_{w} = 2 \sum_i w_{i} [ b_{i} - \mu_i + Y_{i} \log(\mu_i/b_{i}) ]$$

The corresponding variance uses $w_i^2$ in both the Poisson variance and the score-covariance terms. Hence the weighted implementation in weighted_deltac_moments_3param uses

$$4 \sum_i w_{i}^2 \mu_i \log^2(\mu_i/b_{i})$$

and the weighted Fisher projection. When the weights are exactly 0 or 1, the weighted implementation reproduces the binary-mask formula. This paper does not claim to reproduce all internal SAS mask conventions; it uses the validated binary convention for the standalone emldetect calibration described here.

DET_ML convention

This paper uses the SAS emldetect natural-log convention:

$$\mathrm{DET\_ML} = -\ln(\mathrm{gammaincc}(\nu/2, \Delta C/2))$$

For a single-image fit with three fitted parameters, the chi-square degrees of freedom are

$$\begin{aligned} \nu &= 3, \\ \mathrm{shape} &= \nu/2 = 1.5 \end{aligned}$$

Therefore

$$\mathrm{DET\_ML} = -\ln(\mathrm{gammaincc}(1.5, \Delta C/2))$$

The second argument of gammaincc is $\Delta C/2$, not $\Delta C$. This log is the natural logarithm, not $\log_{10}$. The inverse transform is

$$\Delta C = 2\,\mathrm{gammainccinv}\left(1.5, \exp(-\mathrm{DET\_ML})\right) \label{eq:detml-inverse}$$

This is implemented in lib/detml.py with $\nu=3$.

The thresholds for the sensitivity-map band are

$$\begin{aligned} \mathrm{DET\_ML}&=6 \rightarrow \Delta C=14.338539 \\ \mathrm{DET\_ML}&=8 \rightarrow \Delta C=18.570072 \\ \mathrm{DET\_ML}&=10 \rightarrow \Delta C=22.755959 \end{aligned}$$

These values define the $\Delta C$ band of interest. They are not the observed DET_ML screening values used to measure the conditional distribution.

Validity range of the variance approximation

The variance expression $\sigma^2_{3p}$ is a local Gaussian/Fisher approximation around the fitted model. When the likelihood surface is close to quadratic and the fitted source is strong enough that position optimization does not absorb isolated Poisson noise features, the approximation is expected to be accurate. The corrected single-$S$ variance diagnostic is summarized in Figure 2.

Self-consistent reference-fitter validation

The predictions of Section 3 are validated on a self-consistent reference fitter in which source injection, source fitting and theoretical moment computation all use the same psfgen template PSF family. This test isolates the $\Delta C$ theory itself; the standalone SAS emldetect is handled separately in Section 5.

Simulation design

The reference validation uses I044 high-statistic psfgen template simulations:

instrument: PN
band: band 4
positions: 17
S_true grid: 8 values
b grid: 6 values
repetitions: 1000 per cell
total rows: 816,000
clean rows: 811,194

Each simulated source is fitted with a three-parameter Poisson likelihood fitter. For each row, the theoretical values are evaluated at the fitted source strength:

$$\begin{aligned} \mu_{3p} &= E[\Delta C | S_{\mathrm{fit}}, b, \mathrm{psf}, \mathrm{mask}], \\ \sigma^2_{3p} &= \mathrm{Var}[\Delta C | S_{\mathrm{fit}}, b, \mathrm{psf}, \mathrm{mask}] \end{aligned}$$

The diagnostic residual is

$$Z = (\Delta C_{\mathrm{obs}} - \mu_{3p}) / \sqrt{\sigma^2_{3p}}$$

This is an observable-conditioned validation because the prediction uses the fitted amplitude $S_\mathrm{fit}$, not the injected amplitude $S_\mathrm{true}$. Here $S_\mathrm{fit}$ is the name of the fitted amplitude in the simulation table; it is the same observable physical quantity $\widehat{S}$ used in the theory chapter.

High-signal (High-S) closure

For $S_\mathrm{fit} > 20$, the self-consistent reference fitter achieves closure at the level required by the paper’s theoretical claims:

Self-consistent reference fitter closure summary.

The stronger $S_\mathrm{fit} > 50$ subset is closer to zero mean and unit width. This is the expected asymptotic behavior of the Fisher/Gaussian approximation in the high-signal regime.

Figure 3 is the main theoretical proof figure. It shows that the cell-averaged $\Delta C_\mathrm{obs}$ follows $\mu_{3p}$, that the low-$S$ bias and the $S_\mathrm{fit} > 20$ closure band are clearly separated, that the high-$S$ standardized residuals approach the standard normal distribution, and that the closure is stable across detector positions. Table [tab:high-s-closure] gives the corresponding numerical summary.

Low-signal (Low-S) beyond-Fisher bias

The same validation shows that the theory does not close consistently at low fitted signal. The low-$S$ interval has a positive residual bias:

Low-signal beyond-Fisher bias in the self-consistent reference fitter.

This is the beyond-Fisher interval. Position optimization absorbs positive Poisson noise features, making the likelihood surface difficult to approximate by a local quadratic expansion. This bias is not a grid artifact; the v2/v10b audits identified it as a genuine limitation of the Gaussian/Fisher approximation at low fitted signal.

The m2 simulation correction model using $(S_\mathrm{fit},\,b,\,A_2,\,W_\mathrm{eff})$ — where $A_2$ is the second moment of the local PSF and $W_\mathrm{eff}$ is an effective PSF winginess parameter, both derived from position-dependent psfgen templates — reduces the leave-one-position-out cross-validation outlier residuals to about 0.095 $\sigma$. This can serve as a diagnostic tool for low-$S$ theory error. It is not an independent emldetect sensitivity-map recipe, since standalone emldetect has additional implementation residuals.

Figure 4 shows the beyond-Fisher low-$S$ bias and the m2 correction diagnostic.

Why the ML$\sim$10 interval is not covered by high-signal closure

The sensitivity-map band of interest is DET_ML around 6–10. Under the single-image three-parameter SAS convention,

$$\begin{aligned} \mathrm{DET\_ML}&=6 \rightarrow \Delta C=14.338539 \\ \mathrm{DET\_ML}&=8 \rightarrow \Delta C=18.570072 \\ \mathrm{DET\_ML}&=10 \rightarrow \Delta C=22.755959 \end{aligned}$$

For the current PN band-4 setup with $b=0.04$, this corresponds roughly to $\widehat{S}\sim 7.5$–$12.1$. This is below the high-$S$ validation range $S_\mathrm{fit} > 20$. Therefore high-$S$ closure alone cannot serve as the error budget for a $\mathtt{DET\_ML}\sim 10$ sensitivity map.

The threshold-band problem is different. It asks how precisely the source amplitude $S_\mathrm{lim}$ corresponding to the required $\Delta C_\mathrm{th}$ can be recovered after accounting for the low-$S$ theory bias and the standalone emldetect implementation effects. This operational calibration is the subject of Section 5.

Closed-form global corrections

The above low-signal bias can be described in compact closed form. In observable space, binning the reference fitter output by $\widehat{S}$ (using $\widehat{S}>4$ to constrain the low-signal end), across all 17 positions and 6 backgrounds, the mean residual and residual dispersion ratio can both be well described by a $(\widehat{S},b)$ expression. These global expressions empirically support the following interpretation: the residual beyond Fisher is a smooth, low-dimensional function of $(\widehat{S},b)$, with no significant off-axis dependence (Section 4.6). This smoothness motivates the empirical grid interpolation strategy deployed in Section 5: if the residual were chaotic or high-dimensional, sparse grid interpolation would fail. However, these specific closed-form global corrections are strictly diagnostic and are not directly used in the final operational calibration. The standalone emldetect calibration absorbs all low-signal bias and implementation-specific effects entirely into the empirical $R_\mathrm{eml}$ grid.

The empirical nature of these corrections is itself a constrained result, not merely a modeling convenience. Two direct validation tests (to be presented elsewhere) bound the simpler alternatives. First, a second-order expansion of $C_\mathrm{best}$ in the amplitude score (evaluated without fit coefficients) cannot reproduce the five coefficients of Eqs. ([eq:mean-corr])–([eq:sigma-corr]): its mean contribution approximates a constant position-rank gain rather than the observed $\widehat{S}^{-0.40}$ decay; its variance term cannot give the empirical $1/\widehat{S}$ dispersion deficit. This places the low-signal residual outside local-quadratic or Bartlett-type expansions, consistent with the position argmax mechanism separated in the reference fitter decomposition. That mechanism is structurally homologous to the look-elsewhere effect in high-energy physics (Gross & Vitells 2010): the observed $\Delta C$ is the supremum of the likelihood over a position search domain, and a Taylor expansion about an interior point cannot represent an argmax over a domain. Notably, Stewart (2009) likewise needs an empirical additive constant (about $+1$) in his single-parameter sensitivity-mean expression to match Monte Carlo results, and that constant degenerates after discarding negative-amplitude fits — qualitatively consistent with the low-signal additive mean shift we measure here, which can be viewed as an early sign of the same effect in the single-parameter case. The difference is that in our three-parameter production setup the shift is not constant but decays as $\widehat{S}^{-0.40}$ and depends on background. Second, a cross-PSF-transferability test using a Gaussian $\sigma=1.59$ pixel, a wide Gaussian $\sigma=2.79$ pixel, and a King $r_c=2.5$ pixel, $\alpha=1.5$ toy configuration (with the same $\mathtt{ecut}=15$ mask) shows that the correction coefficients are PSF-family-specific: the Gaussian toy mean law and the XMM-psfgen template production law differ at high significance in the most precise low $\widehat{S}$ bin, and neither $b A_\mathrm{psf}$ nor $b\,n_\mathrm{pix}$ allows merging the mean coefficients across PSF families. Therefore the low-signal correction must be calibrated per PSF template family. Within the XMM psfgen template family used in this paper, the 17-position off-axis stability shown below is precisely the validated operating domain.

Two further consequences of this PSF-family perspective are worth recording here. They build on a pilot study covering three EPIC cameras and two bands, using real psfgen templates prepared through a single template pipeline. First, once template preparation is unified, the per-instrument, per-band structure of the correction disappears: the five coefficients simplify to a tag-free form where three are shared constants and two are smooth univariate functions of scalar PSF features (effective PSF area and peak pixel fraction), and the model no longer retains any instrument or band label. The full presentation of this simplified tag-free correction law is deferred to a follow-up paper. Second, directly relevant to the application of any correction in this paper: the absolute value of the correction coefficients depends on the PSF template preparation convention (high-resolution sampling, recombination, cropping and normalization choices). Calibration and deployment must therefore use the same template pipeline; coefficients calibrated under one preparation convention must not be combined with templates prepared under another convention.

The conditional mean is corrected additively:

$$\mu_{3p}^\mathrm{corr}(\widehat{S},b) = \mu_{3p}(\widehat{S},b,\mathrm{psf},\mathrm{mask}) + \big(a_\mu + c_\mu\log_{10}b\big)\,\widehat{S}^{-p_\mu}, \label{eq:mean-corr}$$

where $a_\mu=0.59\pm0.12$, $c_\mu=-1.02\pm0.10$, $p_\mu=0.40\pm0.05$ (uncertainties from position-level bootstrap, see below), and the residual dispersion is corrected multiplicatively:

$$\sigma_C^\mathrm{corr}(\widehat{S},b) = \sqrt{\sigma^2_{3p}}\; \Big[\,1-\frac{a_\sigma + c_\sigma\log_{10}b}{\widehat{S}}\,\Big], \label{eq:sigma-corr}$$

where $a_\sigma=1.49\pm0.04$, $c_\sigma=0.39\pm0.03$. The dispersion correction has the form of a $1/\widehat{S}$ low-signal deficit that vanishes in the asymptotic limit and returns the Fisher-projected variance at high fitted signal. Eq. ([eq:sigma-corr]) reduces the bin-level scatter of the empirical-theory $\sigma_C$ ratio from $0.074$ (uncorrected) to $0.049$, with residuals consistent with the Monte Carlo uncertainty of the per-bin scatter estimate. At the highest background ($b=0.16\,\mathrm{counts\,pixel^{-1}}$) and lowest fitted signal there remains a mild residual where the empirical deficit is flatter in $\widehat{S}$ than the global form; this interval lies well below $\DETML=6$ and does not affect the validated operating range.

Off-axis independence of the corrections

The two corrections in Eqs. ([eq:mean-corr])–([eq:sigma-corr]) depend only on the fitted amplitude and background. Adding a linear off-axis term $d\,\theta$ to the scatter model barely changes the fit quality ($d=-5\times10^{-3}$ per arcmin; root-mean-square residual unchanged to three decimals), and the empirical $\sigma_C$ ratio stays flat at $0.92$–$0.93$ over the full $0$–$14'$ off-axis range. Thus the apparent "off-axis" variance shift is not an off-axis effect: it is the low $\widehat{S}$ deficit of Eq. ([eq:sigma-corr]), since those positions sample a larger fraction of low $\widehat{S}$ fits and so show up more strongly in off-axis rings. Within the validated domain, a single $(\widehat{S},b)$ correction suffices for all detector positions. This intra-family stability, contrasted with the cross-family dependence discussed above, delineates the validity boundary of Eqs. ([eq:mean-corr])–([eq:sigma-corr]).

The coefficients are stable under changes in the number of detector positions used in the fit. Position-level bootstrap (resampling the 17 positions with replacement, 400 realizations) recovers the central values cited above and provides the basis for the uncertainties in Eqs. ([eq:mean-corr])–([eq:sigma-corr]). Resampling random position subsets of fixed size $k$ keeps the central coefficients within the scatter and shrinks their spread as $k^{-1/2}$, showing that the positions independently sample the corrections with no systematic off-axis trend.

Limitations of high-signal asymptotic closure and the role of empirical calibration

The sensitivity-map band of interest is $\DETML$ around 6–10, corresponding to $\Delta C \approx 14.3$–$22.8$. For the current PN band-4 setup with $b=0.04\,\mathrm{counts\,pixel^{-1}}$, this corresponds roughly to $\widehat{S}\sim 7.5$–$12.1$. This is strictly below the high-signal validation range $\widehat{S} > 20$. We therefore explicitly state: the high-signal asymptotic closure of Section 4.2 does not directly validate the theory operationally in the threshold band. In this low-signal interval, the threshold-band residual $R_\mathrm{eml}$ is handled by pure empirical interpolation.

Physically, $R_\mathrm{eml}$ comes from three sources: (i) the beyond-Fisher bias of the $\Delta C$ theory itself near the $S=0$ boundary, where the non-negativity constraint violates Wilks’ theorem interior-point assumption; (ii) implementation-specific artifacts of standalone emldetect, including its internal optimizer, PSF sampling and numerical discretization details; (iii) the difference between the $\mathtt{ecut}=15$ fit footprint used in the theory and the mask actually applied by the pipeline. Although $R_\mathrm{eml}$ is in principle a function of $S$, it varies slowly enough over the narrow threshold band ($\widehat{S}\sim 7$–$12$) to be robustly approximated as a scalar offset per position. Paired-seed experiments comparing the self-consistent reference fitter and standalone emldetect show that independent $\Delta C$ tracks the reference fitter almost perfectly (correlation 0.996), with a small affine offset and an implementation jitter growing with off-axis angle; this supports treating the transfer from reference theory to independent pipeline as an empirical implementation layer rather than a change to the observable-conditioned form itself.

A standard Monte Carlo false-source-injection method (e.g. the esensmap simulation mode) would automatically capture all these effects: inject sources at each sky position and background level, run the full pipeline, and tabulate the detection fractions numerically. This method, which makes no theoretical approximation, offers excellent absolute accuracy. However, for full-sky mapping over the full range of background, off-axis angle and PSF geometry encountered in survey-scale sensitivity maps, it is computationally prohibitive. Large modern surveys therefore rely on extensive end-to-end simulations to characterize completeness, purity and selection thresholds, as in the eROSITA eFEDS source-detection calibration (Brunner et al. 2022; Liu et al. 2022).

The scientific value of the Section 3 theory is that it provides the correct functional form of the observable-conditioned $\Delta C$ expectation, since background, PSF geometry and mask vary continuously across the detector. By anchoring the computation to the deterministic theoretical mean $\mu_{3p}$, operational calibration only needs to empirically determine the slowly varying scalar residual $R_\mathrm{eml}$ on a sparse grid. This hybrid approach (analytic backbone plus empirical residual calibration) sacrifices the absolute accuracy of full Monte Carlo (the $S_\mathrm{lim}$ recovery error of 11–17% reported in Section 5) but avoids the extreme computational cost of dense false-source injection at every sky position. This operational calibration is the subject of Section 5.

Standalone emldetect threshold calibration

The standalone SAS emldetect uses its own source fitting, optimizer and PSF/model conventions; these implementation details generate position-dependent residuals relative to the reference theory derived in Section 3. Below we construct an empirical, position-dependent residual calibration that absorbs this bias for the threshold interval DET_ML=6–10, suitable for independent false-source sensitivity-map work.

Why standalone emldetect needs a separate residual

The raw theory prediction is

$$\mu_{3p}(S,b,\mathrm{psf},\mathrm{mask}) = E[\Delta C | S,b,\mathrm{psf},\mathrm{mask}]$$

If standalone emldetect matched this model exactly, the threshold limit could be found by solving

$$\mu_{3p}(S_{\mathrm{lim}},b,\mathrm{psf},\mathrm{mask}) = \Delta C_{\mathrm{th}}$$

In practice, this zero-residual formulation is only a borderline approximation within the threshold band. In the v13b validation, the raw recipe O0 has a maximum $S_\mathrm{lim}$ error of 19.9% at DET_ML=6, 15.6% at DET_ML=8, and 13.3% at DET_ML=10. The error is position-dependent, not a single scalar offset.

The failure of full observable-conditioning for standalone emldetect closure is shown in Figure 5. Figure 6 shows that increasing the detector footprint variants does not remove the off-axis bias.

Therefore, we define the standalone emldetect residual as:

$$R_{\mathrm{eml}} = \Delta C_{\mathrm{eml}} - \mu_{3p}(\widehat{S}_{\mathrm{eml}},b,\mathrm{psf},\mathrm{mask}) \label{eq:reml}$$

where $\Delta C_\mathrm{eml}$ is the standalone emldetect Cash improvement converted under the single-image $\nu=3$ convention, and $\widehat{S}_\mathrm{eml}$ is the fitted emldetect source amplitude. The calibrated threshold equation is:

$$\mu_{3p}(S_{\mathrm{lim}},b,\mathrm{psf},\mathrm{mask}) = \Delta C_{\mathrm{th}} - \widehat{R}_{\mathrm{eml}} \label{eq:sensitivity-limit}$$

This sign follows directly from the definition of the residual. If $\widehat{R}_\mathrm{eml}>0$, standalone emldetect produces a larger $\Delta C$ than the raw theory at the same fitted source strength, so the theoretical threshold required for the sensitivity limit is lower. If $\widehat{R}_\mathrm{eml}<0$, a higher source strength is required.

Threshold-band calibration sample

Calibration uses paired false-source simulations, where each realization is analyzed by both the reference fitter and the standalone SAS emldetect. The current production grid has 14 PN band-4 calibration positions; P3 is excluded due to a detector-edge anomaly.

The target threshold band is defined in $\Delta C$ space:

$$\begin{aligned} \mathrm{DET\_ML}&=6 \rightarrow \Delta C=14.338539 \\ \mathrm{DET\_ML}&=8 \rightarrow \Delta C=18.570072 \\ \mathrm{DET\_ML}&=10 \rightarrow \Delta C=22.755959 \end{aligned}$$

Rows are selected based on the reference theory prediction:

$$\mu_{3p,\mathrm{ref}} in [\Delta C_{\mathrm{ML6}}, \Delta C_{\mathrm{ML10}}]$$

This is very important. The threshold sample is not selected on the observed DET_ML$_\mathrm{eml}$, nor on the observed $\Delta C_\mathrm{eml}$, because these quantities carry the random fluctuations and implementation residuals we are trying to calibrate. The selection is based on the predicted/reference $\Delta C$ from the observable-conditioned theory.

The baseline calibration grid is built on:

$$b = 0.04 counts/pixel$$

The background shift toward $b=0.02$ and $b=0.08$ is tested separately in Section 5.6.

Recipe definitions

We compare four recipes.

O0: Raw theory.
Set $\widehat{R}_\mathrm{eml}=0$ and solve

$$\mu_{3p}(S_{\mathrm{lim}}) = \Delta C_{\mathrm{th}}$$

This tests whether the theory can be used without standalone emldetect calibration.

O4-global: Scalar residual correction.
Use a single mean $R_\mathrm{eml}$ for all positions. This tests whether the standalone emldetect residual can be treated as a global offset.

O4-lin: Position-dependent 2D linear interpolation.
_ is estimated by linear interpolation on the production coordinate sector_phi, restricted to leave-out positions inside the calibration-grid convex hull.

O4-nn: Nearest-neighbor fallback.
Estimate $\widehat{R}_\mathrm{eml}$ from the nearest calibration position in the same coordinate space. This is the fallback scheme for positions outside the O4-lin convex hull.

The production coordinates are:

$$(offaxis_arcmin, sector_phi)$$

sector_phi is a calibration binning coordinate with four quadrant codes (0, 90, 180, 270 degrees). It is not the physical detector azimuth. The physically defined DETX/DETY coordinates are retained as reference and future-work coordinates, but v15 shows they are not yet production-ready.

Leave-one-position-out validation

Validation uses leave-one-position-out cross-validation (LOPO): each position is held out in turn, the calibration is built from the remaining positions, and the predicted $S_\mathrm{lim}$ is compared with the position’s own calibrated truth. The reported error is

$$||S_{\mathrm{lim,pred}} - S_{\mathrm{lim,true}}| / S_{\mathrm{lim,true}}$$

The results in production coordinates are:

The main conclusions are:

• O0 is usable only as a coarse fallback. It is borderline at DET_ML=6.

• O4-global fails at DET_ML=6, so a single scalar correction is not sufficient.

• O4-lin is the preferred interpolation inside the hull.

• O4-nn is the full-position fallback scheme and passes the 20% limit grid in the current 14-position calibration set.

O0 passes the 20% limit grid at DET_ML=8 and DET_ML=10, but the most demanding edge in the sensitivity-map band is DET_ML=6. At this position O0 is only borderline while O4-global fails, motivating the position-dependent O4-lin/O4-nn residual calibration.

For O4-lin in sector_phi, the four leave-out positions outside the hull are P1, P8, P12 and P16. These positions need the O4-nn fallback.

Figure 7 and Table [tab:v13b-recipe] present this recipe comparison as the main result of this work.

Coordinate decision

Three coordinate choices were tested in v13b:

Therefore the production choice is sector_phi, with a mandatory warning:

$$\begin{aligned} sector_phi is sector-coded and should not be interpreted as physical detector \\ azimuth \end{aligned}$$

This choice is made for practical reasons. Although DETX/DETY is better grounded physically, the sparse grid does not give sufficient O4-lin coverage. v15 tested whether a denser DETX/DETY grid could fix this and it did not pass the current production limit grid; see Section 5.7.

Background transfer

The production grid is calibrated at $b=0.04$. To test whether the same grid can be used over a limited background range, v13b applies the $b=0.04$ grid to paired simulations at $b=0.02$ and $b=0.08$.

Compared with the LOPO interpolation error, the transfer error is very small. A rigorous statement suitable for the paper is:

R_eml has a weak background dependence between b=0.02–0.08, with a maximum tested S_limit transfer error of 5.7%.

Calling this calibration fully background-independent is incorrect. Figure 8 shows the same background-transfer validation pictorially.

Dense DETX/DETY test as a negative result

Calibrating directly on A coordinates on DETX/DETY would be more physically interpretable than sector_phi, but the sparse-grid result of Section 5.5 leaves only 6/14 leave-out positions inside the O4-lin hull. A denser DETX/DETY grid was therefore tested to ask whether spatial densification could recover the working coverage. The result is negative:

The underlying reasons are:

1. Hull topology: leave-one-position-out cross-validation leaves boundary positions outside the interpolation hull even when more interior points are added.

2. emldetect edge failure: detector-edge anchors fail and cannot safely extend the hull.

3. $S$-dependent $R_\mathrm{eml}$: a few problem positions show strong source-strength variation in the residual, violating the position-constant-residual approximation.

The v15 result does not supersede v13b. It shows that keeping sector_phi as the production coordinate and treating DETX/DETY as reference/future work is reasonable. Detailed dense-grid results are given in Table [tab:app-v15-dense-detx-dety].

Unified practical sensitivity-map recipe

The Section 3 theory and Section 5 residual calibration together produce a five-step procedure for evaluating the standalone emldetect sensitivity limit $S_\mathrm{lim}$ at a given detector position and detection threshold, valid within the PN band-4 single-image three-parameter domain established above.

Inputs and outputs

For each detector position, the required inputs are:

Target DET_ML threshold; per-pixel count background b; pixelized PSF template; binary fit mask / footprint; detector-position metadata; R_eml calibration grid

The output is:

$$S_{\mathrm{lim}}$$

where $S_\mathrm{lim}$ is the fitted source amplitude corresponding to the expected (i.e. Asimov) $\Delta C$ equal to the calibrated threshold after residual calibration. It is usually not the 50% completeness flux. The relation between this mean threshold limit and the completeness-percentage limit depends on the conditional $\Delta C$ distribution in the threshold band, which this paper does not characterize.

Output:

For the current v13b recipe, the validated domain is:

instrument: PN
band: band 4
fit: single-image, three-parameter source fit
nu: 3
background range tested: b=0.02--0.08 counts/pixel
threshold range: DET_ML=6--10
production coordinate: sector_phi

Step 1: Convert DET_ML to DeltaC

The SAS emldetect convention is

$$\mathrm{DET\_ML} = -\ln(\mathrm{gammaincc}(\nu/2, \Delta C/2))$$

For the single-image three-parameter case:

$$\begin{aligned} \nu &= 3 \\ \mathrm{shape} &= 1.5 \end{aligned}$$

Therefore the threshold inverse-transform formula is

$$\Delta C_{\mathrm{th}} = 2 \, \mathrm{gammainccinv}(1.5, \exp(-\mathrm{DET\_ML}))$$

The thresholds actually used are:

SAS single-image three-parameter threshold conversion used by the sensitivity-map recipe.

This conversion must use $\Delta C / 2$ as the second argument of gammaincc. It must use the natural-log convention. The old shape=2.5 transform is no longer used.

Step 2: Build the local theory function

At the detector position of interest, construct the local PSF and fit mask. The theory function is

$$\mu_{3p}(S,b,\mathrm{psf},\mathrm{mask}) = E[\Delta C | S,b,\mathrm{psf},\mathrm{mask}]$$

For a binary mask this is the pixel sum

$$\begin{aligned} \mu_{3p}(S,b,\mathrm{psf},\mathrm{mask}) \\ &= 2 \sum_i M_{i} [ (b + S p_{i}) \log((b + S p_{i})/b) - S p_{i} ]\\ \end{aligned}$$

The PSF should be normalized using the same convention as the calibration grid. In the current v13b recipe the mask is the 15-pixel fit footprint used in the paired standalone emldetect simulations. A weighted-mask variant exists in the code, but the production recipe described here uses the validated binary mask convention.

Step 3: Estimate standalone emldetect residual

The residual calibration grid is stored in the threshold band:

$$R_{\mathrm{eml}} = \Delta C_{\mathrm{eml}} - \mu_{3p}(\widehat{S}_{\mathrm{eml}},b,\mathrm{psf},\mathrm{mask})$$

The production interpolation coordinate is

$$(offaxis_arcmin, sector_phi)$$

sector_phi is a calibration block code, not the physical detector azimuth. For the current 14-position grid, O4-lin interpolation covers 10/14 leave-out positions. Positions outside the hull are

$$P1, P8, P12, P16$$

The interpolation rule is:

if query point lies inside the sector_phi convex hull:
    use O4-lin
else:
    use O4-nn fallback

This fallback scheme is required in current sparse-grid production use. Although O4-nn is not the preferred interpolator inside the hull, it is a validated full-position fallback scheme.

Step 4: Solve S_limit

The calibrated threshold equation is:

$$\mu_{3p}(S_{\mathrm{lim}},b,\mathrm{psf},\mathrm{mask}) = \Delta C_{\mathrm{th}} - \widehat{R}_{\mathrm{eml}}$$

The right-hand side is the effective $\Delta C$ target value in the theoretical reference model. Set:

$$\Delta C_{\mathrm{eff}} = \Delta C_{\mathrm{th}} - \widehat{R}_{\mathrm{eml}}$$

Then solve

$$f(S) = \mu_{3p}(S,b,\mathrm{psf},\mathrm{mask}) - \Delta C_{\mathrm{eff}} = 0$$

v13b validation used a bracketed root finder:

S_limit = brentq(f, 0.1, 500.0)

The approximate linear formula |dR| / dmu_dS should not be used as a production calculation method. This is the root cause of the deprecated v13 error (in which DETX/DETY nearest-neighbor errors were overestimated). This paper’s recipe uses the full nonlinear solve.

Step 5: Report uncertainty budget

The uncertainty report should state which interpolation branch was used and cite the validated error range.

For sector_phi at DET_ML=6:

Validated uncertainty budget for production coordinates at DET_ML=6.

For higher thresholds, the interpolation error decreases:

Interpolation errors of the v13b recipe at higher thresholds.

The recommended wording is:

Within the current sector_phi hull use O4-lin and quote the validated maximum S_limit error of 11.0% at ML=6. Outside the hull use O4-nn fallback and quote the validated maximum error of 17.2% at ML=6.

Distribution

The previous sections calibrated the behavior of the likelihood-ratio statistic in $\Delta C$ space. The calibrated distribution has two catalog-level uses. The first is interpretive: for traditional -selected catalogs — where the catalog may still be defined by an observed threshold $\Delta C>\Delta C_\mathrm{th}$ — the forward distribution provides the likelihood needed to infer the corresponding fitted-amplitude distribution. The second is constructive: catalog cuts can be defined directly in fitted-amplitude coordinates by requiring a selected percentile of the local $\Delta C$ distribution to reach the detection threshold. Both uses rely on the same forward kernel but apply it in opposite directions (Table [tab:two-catalog-uses]).

Here $S$ denotes the source-amplitude coordinate used by the observable-conditioned theory; throughout this paper we use the fitted-amplitude convention. It is not the intrinsic injected flux.

Figure 9 shows a visual version of the construction. The traditional threshold is a horizontal cut in $\Delta C$ space. Reading backwards from this cut gives the fitted-amplitude distribution of traditional threshold-selected catalogs once a source-count prior is specified. Reading forwards gives percentile crossing: the usual sensitivity map is mean crossing, while lower-percentile crossing defines a more conservative fitted-amplitude threshold-crossing cut. Standalone emldetect applications additionally include the pipeline residual calibration described in Section 5.

The forward kernel in $\Delta C$ space

For each detector position, the operational threshold is first converted to the corresponding single-image three-parameter $\Delta C$ threshold. For example, $\DETML=6$ corresponds to $\Delta C_\mathrm{th}=14.3$ under the $\nu=3$ convention. The distribution is then computed in $\Delta C$ space; remains only the final operational threshold scale.

In the calibrated standalone emldetect case, the mean of the forward kernel is

$$\Delta C_\mathrm{mean}(S,x,b) = \mu_{3p}(S,b,\mathrm{psf}_x,\mathrm{mask}_x) + \widehat{R}_\mathrm{eml}(x), \label{eq:mean-threshold-curve}$$

where $\widehat{R}_\mathrm{eml}$ is the empirical residual interpolated from the Section 5 calibration grid. Under the Gaussian residual approximation, the local kernel can be summarized as

$$p(\Delta C\mid S,x,b) \simeq \mathcal{N}\!\left[ \Delta C_\mathrm{mean}(S,x,b),\, \sigma_C^2(S,x,b) \right], \label{eq:deltac-forward-kernel}$$

The corresponding $p$-th percentile curve is

$$\Delta C_{p}(S,x,b) \simeq \Delta C_\mathrm{mean}(S,x,b) + z_p\,\sigma_C(S,x,b), \label{eq:deltac-percentile-curve}$$

where $z_p$ is the standard normal quantile. The percentile construction requires, besides the calibrated mean residual, the calibrated threshold-band dispersion model $\sigma_C$ in $\Delta C$ space. Thus $\Delta C_{50}\simeq\Delta C_\mathrm{mean}$, $\Delta C_{16}\simeq\Delta C_\mathrm{mean}-\sigma_C$, $\Delta C_{10}\simeq\Delta C_\mathrm{mean}-1.282\,\sigma_C$, $\Delta C_{2.3}\simeq\Delta C_\mathrm{mean}-2\,\sigma_C$.

In this proof-of-concept application we use the empirically corrected dispersion model from Section 4 as the distribution backbone. The Fisher-projected variance provides the analytic scale, but production use of the full kernel or low-percentile maps requires direct validation against standalone emldetect pass rates, since high-signal reference fitter closure alone cannot validate the –10 threshold band.

Application I: Interpreting traditional -selected catalogs

The first application retains the traditional catalog definition. Sources are selected by the observed detection statistic, e.g. >$T$, equivalently $\Delta C>\Delta C_\mathrm{th}$. The new element is that the catalog threshold can now be interpreted through the forward kernel of Eq. ([eq:deltac-forward-kernel]).

For a source with measured $\Delta C_\mathrm{obs}$, the fitted-amplitude posterior in catalog space is

$$p(S\mid \Delta C_\mathrm{obs},x,b) \propto p(\Delta C_\mathrm{obs}\mid S,x,b)\,\pi(S), \label{eq:posterior-single-source}$$

where $\pi(S)$ is the source-count prior expressed in fitted-amplitude coordinates. For a threshold-selected sample, a similar expression is

$$p(S\mid \Delta C>\Delta C_\mathrm{th},x,b) \propto P(\Delta C>\Delta C_\mathrm{th}\mid S,x,b)\,\pi(S), \label{eq:posterior-threshold-selected}$$

where

$$P(\Delta C>\Delta C_\mathrm{th}\mid S,x,b) = \int_{\Delta C_\mathrm{th}}^\infty p(\Delta C\mid S,x,b)\,d\Delta C . \label{eq:threshold-passing-probability}$$

Under the Gaussian residual approximation, this exceedance probability is the normal survival function evaluated in $\Delta C$ space. This path does not replace the traditional threshold. Rather, it adds the likelihood layer needed to interpret the fitted-amplitude distribution of -selected catalogs. This is where Eddington bias enters: the noisy threshold statistic must be combined with a source-count prior, especially when the underlying $dN/dS$ is steep (Eddington 1913; Georgakakis et al. 2008; Wang 2004). Therefore the specified value or threshold itself does not define a flux distribution; it defines a likelihood factor in $\Delta C$ space. To obtain $p(S|\Delta C_\mathrm{obs})$ or $p(S|\Delta C>\Delta C_\mathrm{th})$, the likelihood must be combined with a source-count prior, e.g. a power-law $dN/dS$ slope or index expressed in fitted-amplitude coordinates. Because this is an inverse problem, it necessarily depends on the population prior. If the required prior is the intrinsic $dN/dS_\mathrm{true}$ rather than the fitted catalog amplitude prior, an additional measurement model connecting $S_\mathrm{true}$ to the fitted amplitude $S$ is needed.

Accordingly, the corresponding sky coverage is probabilistic rather than the step-function area from a single limiting amplitude. At a given pixel $j$, a source with detected fitted amplitude $S$ has catalog-space density

$$q_j(S) \propto P(\Delta C>\Delta C_\mathrm{th}\mid S,x_j,b_j)\,\pi(S), \label{eq:pixel-detected-amplitude-density}$$

within the normalization range in fitted-amplitude coordinates. Thus the traditional threshold catalog assigns probabilities to detected amplitudes at each pixel rather than a single deterministic limiting amplitude. The corresponding effective area at fitted amplitude $S$ is

$$\Omega_\mathrm{th}(S) = \sum_j \Omega_j\, P(\Delta C>\Delta C_\mathrm{th}\mid S,x_j,b_j), \label{eq:traditional-probabilistic-sky-coverage}$$

where $j$ indexes sky or detector pixels. If the goal is the expected number of detected fitted amplitudes in interval $A$, the same probability-weighted area combined with the source-count prior gives

$$N_\mathrm{det}(A) \propto \int_A \Omega_\mathrm{th}(S)\,\pi(S)\,dS . \label{eq:traditional-detected-amplitude-distribution}$$

This is why Application I’s sky coverage follows the detected fitted-amplitude probability distribution, whereas the percentile-cut construction below starts from the selected fitted-amplitude threshold-exceedance probability.

Application II: Fitted-amplitude percentile cuts

The second application uses the same kernel forward. Instead of asking which fitted amplitudes are contained in a noisy -selected sample, it asks what fitted amplitude is needed to exceed the detection threshold in a specified fraction of realizations. This map-construction step fixes $S$ and evaluates $P(\Delta C>\Delta C_\mathrm{th}|S,x,b)$, so it does not require a $dN/dS$ slope or source-count index.

The usual sensitivity-map calculation is recovered as mean-curve crossing

$$\Delta C_\mathrm{mean}(S,x,b)=\Delta C_\mathrm{th}. \label{eq:mean-curve-crossing}$$

Under the symmetric Gaussian residual approximation, this is also the 50th percentile threshold-exceedance point. We avoid calling it the empirical 50% completeness flux because injection-recovery completeness is usually defined as $P(\mathrm{detected}|S_\mathrm{true})$ in injected-flux coordinates.

A more conservative fitted-amplitude cut is obtained by solving the threshold equation with a lower percentile curve,

$$\Delta C_p(S,x,b)=\Delta C_\mathrm{th}, \label{eq:percentile-threshold-crossing}$$

where $p<50$. Direction matters. Requiring 90% of realizations at a given fitted amplitude to exceed the threshold uses the lower 10th percentile, not the higher 90th percentile:

$$\Delta C_{10}(S,x,b)=\Delta C_\mathrm{th}.$$

Likewise, under the Gaussian approximation, using the $\mu-2\sigma$ lower envelope corresponds to about a 97.7% threshold-exceedance criterion. Table [tab:threshold-crossing-terms] summarizes the terms used in this paper.

Sky coverage for fitted-amplitude threshold-crossing cuts

Applying the same exceedance rule to each detector pixel yields a family of fitted-amplitude sensitivity maps: the mean-threshold map, the 84% threshold-exceedance map from the lower 16th percentile, the 90% threshold-exceedance map from the lower 10th percentile, and so on. Each map can be converted to a threshold-exceedance sky-coverage curve by counting the sky area where the local limiting amplitude is below the trial fitted amplitude,

$$\Omega_p(S) = \sum_j \Omega_j\, I\!\left[S_\mathrm{lim}^{(p)}(j) \le S\right],$$

where $j$ indexes sky or detector pixels and $p$ denotes the percentile curve used in Eq. ([eq:percentile-threshold-crossing]). This threshold-exceedance sky-coverage calculation is an instrument/catalog-space quantity in the same fitted-amplitude coordinates used by the observable-conditioned theory: it depends on local background, PSF, mask, threshold and residual calibration, but does not require an assumed $dN/dS$ slope. A source-count model is needed only for subsequent population inference.

Relation between the two uses and Eddington bias

Both applications come from the same conditional distribution $p(\Delta C|S,x,b)$ but answer different questions. Application I is backward: it infers the fitted-amplitude distribution from the observed value or -selected catalog. This is a natural setting for Eddington-bias modeling, since the asymmetry near the detection threshold is produced by the combination of the noisy threshold statistic and the steep source-count prior.

Application II is forward: it asks, from the fitted amplitude, for the probability of exceeding the threshold. This does not remove population-level Eddington bias from luminosity-function inference, but it does allow construction of an instrument catalog-space subsample with an explicitly stated threshold-exceedance probability. In this sense, the percentile sensitivity map is a catalog-construction tool, while the threshold-selected posterior is a catalog-interpretation tool.

Validation requirements and limitations

The Section 5 cross-validation validates the mean/median threshold fitted-amplitude mapping: at , the preferred bilinear interpolation reaches a maximum error of 11.0% within the calibration envelope, while the nearest-calibration-position fallback reaches 17.2

Thus the applications proposed here are a fast, distributed, catalog-oriented framework for threshold interpretation and fitted-amplitude threshold-exceedance sensitivity. Full injection-recovery simulation remains the absolute validation standard for production selection functions in complex survey fields.

The calibration grid, cross-validation tables, figure-generation scripts, configuration files and minimal scripts for reproducing the key validation results will be made available through a persistent Zenodo archive upon publication.

Scope and limitations

The proposed framework focuses on methodological progress: the conditional-mean formula, the Fisher-projected variance, and the observable-conditioning principle transfer algebraically to any Poisson PSF source-detection pipeline using a standard three-parameter fit; however, quantitative use requires pipeline-specific recalibration. Moreover, the quantitative calibration (specifically the position-dependent residual grid, leave-one cross-validation error budget and threshold-exceedance fitted-amplitude mapping) has been specifically validated using single-image false-source simulations in the XMM-Newton PN band-4 setup. Applying this method to different instruments, bands or detection pipelines will require re-deriving the pipeline-specific likelihood transform and fully recalibrating the implementation residuals with matched false-source grids.

Although the self-consistent reference fitter confirms the theoretical model structure in the high-signal region ($\widehat{S}>20$), these results cannot be blindly extrapolated to the threshold band ($\DETML\approx 6$–10). At lower fitted signal ($\widehat{S}\lesssim 12$), the reference fitter exhibits a positive beyond-Fisher bias. The independent calibration absorbs this bias together with other algorithmic effects into the empirical $R_\mathrm{eml}$ residual, providing operational predictions for local $\Delta C$ and distributions.

Moreover, we emphasize that our empirical calibration is not a comprehensive correction to the official XMM-Newton source catalogs (e.g. 4XMM). Catalog environments involve complexities such as complex background map variations, multi-observation mosaicking, vignetting weighting and source confusion effects, which are not fully captured in our single-image standalone simulations. Other current X-ray catalog pipelines handle related complexities by substantially different mechanisms, e.g. the eROSITA eSASS/eFEDS simulation-calibration detection chain (Brunner et al. 2022; Liu et al. 2022) and the Chandra Source Catalog Release 2 series with merged observations, local PSF modeling and Bayesian aperture photometry (Evans et al. 2024). These examples reinforce that the residual calibration developed here is pipeline-specific.

Regarding spatial interpolation, bilinear interpolation within the calibration envelope gives minimum error, while nearest-calibration-position assignment gives a validated full-position fallback. The adopted sector coordinates are calibration coordinates, not physical detector azimuths. Calibration retains a weak background dependence; transferring the $b=0.04\,\mathrm{counts\,pixel^{-1}}$ calibration to the $0.08\,\mathrm{counts\,pixel^{-1}}$ background level introduces up to 5.7% error in inferring sensitivity limits, indicating that extrapolation to substantially different backgrounds without recalibration should be avoided.

If no $R_\mathrm{eml}$ calibration grid is available, the uncalibrated theory can serve as an approximate estimate, with a maximum $S_\mathrm{lim}$ error of 19.9% at . This is a fallback, not a recommended calibration method. If the target setup differs from the validated domain, the calibration should be regarded as unvalidated. Examples include: MOS instruments; merged images or multi-sample fits; fixed-position fits with different degrees of freedom; backgrounds outside the tested range $b=0.02$–$0.08\,\mathrm{counts\,pixel^{-1}}$; PSF or mask conventions not matching the calibration grid; or official 4XMM catalog products.

Finally, Section 6 distinguishes two uses of the same forward kernel. The backward interpretation of traditional -selected catalogs is a population-inference problem requiring a source-count prior. The forward fitted-amplitude percentile map does not require such a prior at the map-construction stage, but source-count and Eddington-bias corrections still require an assumed or fitted source flux distribution, and where needed a measurement model connecting intrinsic flux to fitted flux.

Conclusions

The conditional mean improvement of the Cash statistic at the fitted source,

$$\mu_{3p}(\widehat{S}, b, \mathrm{psf}, \mathrm{mask}) = E[\Delta C | \widehat{S}, b, \mathrm{psf}, \mathrm{mask}],$$

Under the SAS single-image convention ($\nu=3$, $\mathrm{shape}=1.5$; $\DETML=-\ln Q(1.5,\,\Delta C/2)$, so $\DETML=6,\,8,\,10$ correspond to $\Delta C=14.3,\,18.6,\,22.8$), the standard three-parameter (one amplitude, two position) single-source PSF fit is derived in closed form. This mean itself is a generalization of existing analytic sensitivity-map work (Huang et al. 2025; Rosen et al. 2016; Stewart 2009; Traulsen et al. 2019), generalized in conditioning on the fitted amplitude rather than the model amplitude and explicitly carrying the fit mask. The essential difference of this framework from that line of work lies in the other two items: adding the Fisher-projected variance for the three fitted parameters, thereby giving the conditional dispersion of $\Delta C$ rather than only the conditional mean; and treating threshold selection as an exceedance event on a distribution rather than a single deterministic limit.

Self-consistent psfgen-template reference fitter simulations spanning 17 PN band-4 detector positions (about $8\times10^5$ fits) close the predicted mean to within 0.07$\sigma$ in the high-signal region ($\widehat{S}>20$) and close the Fisher-projected variance to within about 5%. The threshold band $\DETML=6$–10 lies at lower fitted signal (roughly $\widehat{S}\sim 7.5$–$12.1$ at $b=0.04\,\mathrm{counts\,pixel^{-1}}$), where the asymptotic low-signal bias and the implementation-dependent residual of standalone SAS emldetect require an empirical position-dependent residual calibration $R_\mathrm{eml}$. The validated mean/50th-percentile exceedance equation

$$\Delta C_\mathrm{mean}(S,b,\mathrm{psf},\mathrm{mask}) = \Delta C_\mathrm{th}$$

On the PN band-4 calibration grid at , it reproduces the reserved reference limiting amplitude with 11% accuracy (preferred bilinear interpolation, 10 of 14 hull positions retained) and 17% (nearest-calibration-position fallback, all 14 positions), with weak background dependence over $b=0.02$–$0.08\,\mathrm{counts\,pixel^{-1}}$ (maximum tested transfer error 5.7%). The adopted sector coordinates are calibration coordinates, not physical detector azimuths.

For catalog work, the calibration kernel can be used in both directions. Backward, it provides likelihood factors such as $p(\Delta C_\mathrm{obs}|S,x,b)$ or $P(\Delta C>\Delta C_\mathrm{th}|S,x,b)$ for combining with a source-count prior to interpret traditional -selected catalogs. Forward, it gives the local fitted-amplitude threshold-exceedance probability $P(\DETML>T|S,x,b)$. For example, the 90% threshold-exceedance fitted-amplitude limit is the source amplitude whose predicted $\Delta C$ distribution has its 10th percentile equal to the detection threshold. The sky area as a function of fitted amplitude is then obtained by counting detector pixels satisfying the selected threshold-exceedance condition.

The methodological core is independent of XMM-Newton: the observable-conditioned $\Delta C$ theory for three-parameter Poisson PSF source fitting, and threshold calibration in $\Delta C$ residual space rather than the pipeline-reported detection-statistic space. Extension to other Poisson PSF source-detection pipelines is a natural follow-up; the above XMM-Newton emldetect result is the first concrete validation and calibration application within this framework.

Data and code availability

The calibration grid, cross-validation tables, figure-generation scripts, SAS configuration files and minimal scripts for reproducing the key validation results (including exact commands, random seeds, psfgen parameters and emldetect invocation flags) have been deposited in a Zenodo archive. Private review links are available; the final DOI will be inserted upon acceptance.

The calibration grid is provided as machine-readable tables. Code is released under an open-source license.

Diagnostic role of injected source strength

Figure 10 is a diagnostic check of the conditioning variable. The simulation knows the injected source strength $S_\mathrm{true}$, but the theory and operational calibration condition on the fitted amplitude $\widehat{S}$. At fixed detected $\widehat{S}$, the normalized residual distributions from different $S_\mathrm{true}$ subsets overlap closely. The residual offset in the lowest $\widehat{S}$ bin is the same low-signal beyond-Fisher bias discussed in Section 4.2, not evidence for conditioning on $S_\mathrm{true}$.

Acknowledgments

We acknowledge support from the National Natural Science Foundation of China. This research used data obtained from the XMM-Newton satellite, an ESA science mission with instruments and contributions funded directly by ESA Member States and NASA.

References