11institutetext: IRAP, CNRS, UPS, CNES, 14 avenue Edouard Belin, 31400
Toulouse, France
11email<EMAIL_ADDRESS>22institutetext: Max Planck Institute for
Extraterrestrial Physics, Giessenbachstrasse 1, 85748 Garching, Germany
33institutetext: Dr. Karl-Remeis-Sternwarte and ECAP, Friedrich-Alexander-
Universität Erlangen-Nürnberg, Sternwartstr. 7, 96049 Bamberg, Germany
44institutetext: Universität Innsbruck, Institut für Astro- und
Teilchenphysik, Technikerstr. 25/8, 6020 Innsbruck, Austria 55institutetext:
Argelander-Institut für Astronomie (AIfA), Universität Bonn, Auf dem Hügel 71,
53121 Bonn, Germany

# The SRG/eROSITA All-Sky Survey

X-ray selection function models for the eRASS1 galaxy cluster cosmology

N. Clerc, email<EMAIL_ADDRESS>J. Comparat 22 R. Seppi 22 E.
Artis 22 Y. E. Bahar 22 F. Balzer 22 E. Bulbul 22 T. Dauser 33 C. Garrel 22 V.
Ghirardini 22 S. Grandis 44 C. Kirsch 33 M. Kluge 22 A. Liu 22 F. Pacaud 55 M.
E. Ramos-Ceja 22 T. H. Reiprich 55 J. Sanders 22 J. Wilms 33 X. Zhang 22

(Received -; accepted -)

###### Abstract

Aims. Characterising galaxy cluster populations from catalog of sources
selected in astronomical surveys requires knowledge of sample incompleteness,
known as selection function. The first All-Sky Survey (eRASS1) by eROSITA
onboard Spectrum Roentgen Gamma (SRG) has enabled the collection of large
samples of galaxy clusters detected in the soft X-ray band over the Western
Galactic hemisphere. The driving goal consists in constraining cosmological
parameters, which puts stringent requirements on accuracy, flexibility and
explainability of the selection function models.

Methods. We use a large set of mock observations of the eRASS1 survey and we
process simulated data identically to the real eRASS1 events. We match
detected sources to simulated clusters and we associate detections to
intrinsic cluster properties. We train a series of models to build selection
functions depending only on observable surface brightness data. We develop a
second series of models relying on global cluster characteristics such as
X-ray luminosity, flux, and expected instrumental count-rate as well as on
morphological properties. We validate our models using our simulations and we
rank them according to selected performance metrics. We validate the models
with datasets of clusters detected in X-rays and via the Sunyaev-Zeldovich
effect. We present the complete Bayesian population modelling framework
developed for this purpose.

Results. Our results reveal the surface brightness characteristics most
relevant to cluster selection in the eRASS1 sample, in particular the
ambiguous role of central surface brightness at the scale of the instrument
resolution. We have produced a series of user-friendly selection function
models and demonstrated their validity and their limitations. Our selection
function for bright sources reproduces well the catalog matches with external
datasets. We discuss potential inconsistencies in the selection models at low
signal-to-noise revealed by comparison with a deep X-ray sample acquired by
eROSITA during its performance verification phase.

Conclusions. Detailed modelling of the eRASS1 galaxy cluster selection
function is made possible by reformulating selection into a classification
problem. Our models are used in the first eRASS1 cosmological analysis and in
sample studies of eRASS1 cluster and groups. These models are crucial for
science with eROSITA cluster samples and our new methods pave the way for
further investigation of faint cluster selection effects.

###### Key Words.:

X-rays: galaxies: clusters – Surveys – Catalogs – Methods: statistical

##  1 Introduction

Galaxy clusters sit at the most massive nodes of the cosmic web. They form
last in the cosmic evolution by accreting groups and smaller structures. Their
distribution is sensitive to the underlying cosmological model and for this
reason they are recognised as a key cosmological probe (Clerc & Finoguenov,
2023). Being rare objects, their census requires surveys spanning a large
fraction of the sky with sensitive instrumentation. The diffuse X-ray emitting
gas, filling the Megaparsec-wide intracluster medium (ICM), signposts
virialized halos and enables their discovery up to large distances. In the
context of X-ray cluster studies, accurate knowledge of a selection function
has been required for cosmological abundance studies (e.g. Böhringer & Chon,
2015; Mantz et al., 2015; Finoguenov et al., 2020; Garrel et al., 2022),
cluster clustering studies (e.g. Marulli et al., 2018; Lindholm et al., 2021),
scaling relation studies (e.g. Pacaud et al., 2006; Mantz, 2019; Bahar et al.,
2022), investigation of extreme objects (e.g. Hoyle et al., 2012).

The eROSITA instrument onboard the Spectrum Roentgen Gamma Mission
(SRG/eROSITA, Predehl et al., 2021) has surveyed the entire sky during its
first six months of operations (Merloni et al., 2024), collecting enough
photons to discover several thousands of galaxy clusters in the Western
Galactic hemisphere (Bulbul et al., in press). Those extended sources are
identified in multi-band optical surveys and their redshift is measured, hence
providing a distance to us observers. The first catalogs of clusters
discovered in the first eROSITA All-Sky Survey (eRASS1) Western Galactic
hemisphere (Bulbul et al., in press; Kluge et al., submitted) are the support
of detailed individual cluster studies (e.g. Veronica et al., 2023; Liu et
al., 2023) and population studies revealing for the first time the large-scale
evolution of clusters and groups up to $z\simeq 1$ and beyond. Among them the
cosmological analyses particularly stand out, being a driver of the eROSITA
mission design and of the construction of cluster catalogs. The first
cosmological results based on cluster number counts are presented in
Ghirardini et al. (submitted); Artis et al. (in prep.). In this article, we
explain how the selection function model supporting these results is
constructed. In fact, virtually any cluster population study based on the
published eRASS1 catalogs needs to account for incompleteness at some stage.
It is indeed inefficient, if not ill-defined, to require complete samples for
population analyses (Rix et al., 2021); modelling incompleteness is in general
more fruitful and in fact, rather inevitable. The current blooming of large
astronomical surveys fosters development of accurate selection models with
powerful statistical methods. The recent Gaia survey is an illustrative
instance where selection functions for the parent catalog (Boubert & Everall,
2020) and for its subsets (Rix et al., 2021; Boubert & Everall, 2021) are
modelled using detailed accounting for the magnitudes, position and parallaxes
(and uncertainties) of the stars.

This paper is organised as follows. We first present detailed motivation for
this work in Sect. 2, where we also introduce our notations. We present the
twin eRASS1 simulations in Sect. 3, in particular we discuss our procedure
matching detections and simulated clusters. We describe a first class of
selection function models relying uniquely on X-ray counts profiles in Sect.
4, highlighting some salient features of cluster selection in the observable
space. In Sect. 5 we introduce a second class of models relying on
intermediate variables. We show the outcome of a validation of our
intermediate models using external catalogs (X-ray and millimeter-band
detections) in Sect. 6. We discuss our results in Sect. 7 and summarize our
findings in Sect. 8.

Unless stated differently, we assume a flat $\Lambda$CDM cosmological model
with parameters from Planck Collaboration et al. (2016). We use $\ln$ for
natural logarithm, base 10 logarithm writes $\log_{10}$. The notation
$\mathcal{N}\left(\mu\mathrel{};\mathrel{}\sigma\right)$ indicates a normal
distribution centred on $\mu$ with standard deviation $\sigma$. We use symbol
$\arcsec$ for arcsecond units.

##  2 Motivation

Let us follow an illustrative example based on cluster count cosmology
studies, although this demonstration may be generalized to other kinds of
population analyses. Most cosmology analyses involve a numerical likelihood:

$L({\rm model})\propto p\left(\rm data\mathrel{}\middle|\mathrel{}\rm
model\right)$

This quantity represents the probability that a set of observed data derives
from a certain model. One assumes the model to be dependent on a set of
parameters $\@vec{\Theta}_{\rm model}$ fully describing it (see Table 1):

$L(\@vec{\Theta}_{\rm model})\propto p\left(\rm
data\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm model}\right)$ (1)

Posterior distributions $p\left(\@vec{\Theta}_{\rm
model}\mathrel{}\middle|\mathrel{}\rm data\right)$ are final products of
cosmology analyses, obtained by application of the Bayes theorem involving the
likelihood (Eq. 1). Most often the required integrals are computed by repeated
samplings of the parameter space, and by formulating as many queries to the
likelihood function value $L$. A practical consequence is that likelihood
evaluation must be computationally efficient.

In experiments involving counts of galaxy clusters – which include cluster
abundance and cluster clustering studies – objects are grouped in bins of one
or several measured quantities. The bins are drawn in a predefined parameter
space $\@vec{\Theta}_{\rm bin}$: measured mass, redshift, flux, pairwise
distances, etc. The likelihood functions combine observed numbers
$\hat{N}_{j}$ and model-predicted numbers $n_{j}$ by means of a statistics
($j$ is an index for the bins):

$L\equiv
L\left(\left\\{\hat{N}_{j}\right\\},\left\\{n_{j}\right\\},\dotsc\right)$

A famous example is the Poisson log-likelihood (see also Cash, 1979), adequate
for shot-noise dominated bins:

$\ln L=\sum_{j}\hat{N}_{j}\ln(n_{j})-\sum_{j}n_{j}\ \ (+{\rm constant})$

The size of the bins is part of the likelihood design and it is usual practice
to consider infinitesimally small bins containing $\hat{N}_{j}=0$ or $1$
object. In such a case, $n_{j}$ transforms into
$\mathrm{d}n/\mathrm{d}\@vec{\Theta}_{\rm bin}$ and the sum signs into
continuous integrals in the expression above. In presence of additional
variance terms, this likelihood must be modified. For instance, Garrel et al.
(2022) present a likelihood accounting for sample variance (e.g. Hu &
Kravtsov, 2003). This component is an important source of variance when
dealing with large numbers of objects (relatively to Poisson noise that
affects small samples).

A model should then predict the value of $n_{j}$. In the ideal case of
unrestricted computational power, excellent fidelity may be achieved by means
of end-to-end models of the eRASS1 sky, and by retrieving the number $n_{j}$
from the simulated catalogs. However it is prohibitive to run one new end-to-
end simulation each time a value is queried by the likelihood function (Eq.
1). The cosmological analysis process necessarily resorts to approximations,
interpolations and emulators. For instance, a parameterized halo mass function
(HMF) model would replace series of full N-body simulations; intra-cluster
medium (ICM) scaling relations replace detailed hydrodynamical simulations;
while a selection function model replaces mock image generation and
processing.

We write $p\left(I\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm sel}\right)$
the probability of selecting an object in the sample given the value of
parameters $\@vec{\Theta}_{\rm sel}$. These values must be a prediction of the
model. The expression for the modelled counts $n_{j}$ may write formally:

$n_{j}=\int\mathrm{d}\@vec{\Theta}_{\rm bin}\mathrm{d}\@vec{\Theta}_{\rm
sel}\mathrm{d}\@vec{\Theta}_{\rm model}\mathbb{1}(\@vec{\Theta}_{\rm bin}\in
j)\\\ \times p\left(I\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm
sel}\right)p\left(\@vec{\Theta}_{\rm bin},\@vec{\Theta}_{\rm
sel}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm
model}\right)\frac{\mathrm{d}n}{\mathrm{d}\@vec{\Theta}_{\rm model}}$ (2)

In the latter expression we have introduced
$\mathrm{d}n/\mathrm{d}{\@vec{\Theta}_{\rm model}}$, the predicted number
distribution of galaxy clusters under a certain model assumption. The
indicator function $\mathbb{1}$ ensures only those objects with
$\@vec{\Theta}_{\rm bin}$ taking values in bin indexed by $j$ are counted. The
three sets of parameters may share part or all of their members.

An important preparatory task in the cosmological analysis consists in
carefully selecting $\@vec{\Theta}_{\rm bin}$ and $\@vec{\Theta}_{\rm sel}$ so
the model is accurate and the modelling effort is well balanced between
$p\left(I\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm sel}\right)$ and
$p\left(\@vec{\Theta}_{\rm bin},\@vec{\Theta}_{\rm
sel}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm model}\right)$. Depending
on how close $\@vec{\Theta}_{\rm sel}$ is from purely observable quantities,
selection function models involve a varying amount of astrophysical modelling.
It is also important to recognize that some aspects of selection may result
from human intervention and thus they are difficult to model with analytic
formulae. A certain degree of empiricism may then be introduced in building
selection models. We will discuss several options in this paper and compare
their advantages and shortcomings.

Table 1: Glossary of symbols and conventions used throughout this paper.

Symbol |  Signification
---|---
$\@vec{\Theta}_{\rm model}$ |  Variables governing the population of clusters, e.g. cosmological parameters, or $(M_{500},z)$
$\@vec{\Theta}_{\rm sel}$ |  Variables describing the selection of clusters, e.g. count-rate
$\@vec{\Theta}_{\rm meas}$ |  Variables describing measurement associated to clusters, e.g. measured flux
$I_{\rm{catalog}}$ |  Booleanaa$a$Boolean random variables values may be mapped to integers 0 and 1, following usual practice.  random variable indicating the presence of a cluster in a catalog, e.g. $I=1$
$I_{\rm{eRASS1}}$ |  A source appears in the X-ray source detection list, it is detected
$I_{\rm{main}}$ |  A source is detected and classified as extended in the eRASS1 primary cluster catalog ($\mathcal{L}_{\rm ext}>3$)
$I_{\rm{cosmo}}$ |  A source is detected and classified as extended in the eRASS1 cosmology catalog ($\mathcal{L}_{\rm ext}>6$)
$L_{X}$ |  True cluster rest-frame 0.5-2 keV luminosity integrated in a cylinder of radius $R_{500}$
$z$ |  True cosmological redshift of a cluster (no peculiar motion)
$f_{X}$ |  True, noiseless absorbed 0.5-2 keV flux of a cluster within an aperture of radius $R_{500}$
$CR$ |  True, noiseless unabsorbed 0.2-2.3 keV count-rate of a cluster within $R_{500}$ aperture
$T_{\rm exp}$ |  Local exposure time in the eRASS1 survey
$N_{H}$ |  Hydrogen column density along the line-of-sight direction
bkg |  Local background surface brightness in the eRASS1 survey in the 0.3-2.3 keV band
$\mathcal{H}$ |  Abridged notation for $(N_{H},T_{\rm exp},{\rm bkg})$
111

##  3 The eROSITA galaxy cluster survey and its simulated twin

###  3.1 The eRASS1 galaxy cluster samples

The primary eRASS1 cluster catalog in the Western Galactic hemisphere
(hereafter eRASS1-primary or eRASS1-main) is described in detail in Bulbul et
al. (in press) and Kluge et al. (submitted). This catalog builds upon the
eRASS1 X-ray source detection list in the soft band ($0.2-2.3$ keV) presented
in Merloni et al. (2024). The source detection is achieved with the eROSITA
Science Analysis Software (eSASS in version 010), and it incorporates two
essential steps in its final stages (see Brunner et al., 2022, for further
information). First, sources are detected by searching counts significantly
above the background level; those sources are individually characterized using
a point-spread function (PSF) fitting algorithm. The detection likelihood
(DET_LIKE_0 in the parent catalog, hereafter shortened to $\mathcal{L}_{\rm
det}$) and the extent likelihood (EXT_LIKE, hereafter $\mathcal{L}_{\rm ext}$)
are among the most relevant parameters for selecting extended sources. In
particular, the eRASS1-primary catalog is compiled based on a low threshold of
$\mathcal{L}_{\rm ext}>3$ to maximize completeness. A rigorous cleaning
process is applied to the X-ray source catalog to ensure clean and
uncontaminated images are suitable for extended source identification. A
series of additional cleaning steps remove spurious split detections; we will
assume their impact on completeness is null, as suggested by the number of
secondary matches found in the digital twin matching the fraction of split
sources in real data (Seppi et al., 2022, and next section).

The description of the eRASS1 cosmology sample (hereafter eRASS1-cosmo) is
detailed in Bulbul et al. (in press). It is compiled with a more conservative
cut on the measured parameter $\mathcal{L}_{\rm ext}>6$ compared to the
eRASS1-primary sample to maximize the purity. Exclusively, galaxy clusters
with measured photometric redshifts between 0.1 and 0.8 are selected, keeping
only areas of the sky where the photometry can be uniformly applied in
redshift measurements in the Legacy DR10-south area (Kluge et al., submitted).
The cosmology sample comprises 5 259 identified clusters with low
contamination levels below 5%.

This paper focuses only on X-ray selection function models. However, both the
eRASS1-primary and eRASS1-cosmo catalogs are derived from extensive optical
follow-up designed for identification and redshift measurement. This
identification may act as an additional selection filter impacting
completeness and purity estimates and should be accounted for in science
analyses with these catalogs (see Ghirardini et al., submitted). Most of the
discussion in the paper neglects this effect, and we refer to Kluge et al.
(submitted) for a detailed description of the optical selection function.

###  3.2 The eRASS1 digital twin

Understanding selection effects requires mock observations reproducing as many
of the characteristics of the actual data as possible. We briefly summarize
here the main features of the eRASS1 digital twin depicted in Seppi et al.
(2022). The parent halo catalog originates from UNIT1i dark-matter simulations
(Chuang et al., 2019). A full sky light cone is created by replicating shells
of the individual snapshots of the simulation (Comparat et al., 2020). The
large, albeit finite, box size prevents simulating the very nearby, massive
halos that constitute a tiny fraction of the sources in eRASS1 and require a
separate treatment. X-ray sources are associated to the light cone using
models for the emission of AGN (Comparat et al., 2019) and galaxy clusters and
groups (Comparat et al., 2020). The cluster model is generated from a set of
real clusters by accounting for the covariances between the surface brightness
profile, halo mass, temperature, and redshift. The baryon profiles are painted
onto halos based on observations; therefore, the predicted emissivity profiles
are representative of observations present in the literature. Low-mass groups
$M_{500c}<5\times 10^{13}M_{\odot}$ require a flux rescaling correction (see
appendix of Seppi et al., 2022). It replaces simulated objects along the
observed stellar mass–X-ray luminosity relation of groups and clusters
(Anderson et al., 2015; Comparat et al., 2022, Zhang Y. et al. in
preparation). A simulated half-sky twin contains more than $10^{6}$ simulated
halos up to $z\sim 1.5$. Besides AGN and clusters, foreground X-ray emitting
stars and galactic foregrounds also have their share of simulated photons.

The SIXTE software (Dauser et al., 2019) generates mock event lists and images
of these photons, as would be seen by eROSITA in the eRASS1 survey. In
particular, instrumental characteristics and scanning laws enable reproducing
the response of the telescopes to the X-ray simulated sky. In order to examine
variability due to stochastic Poisson noise, several tens of simulations are
produced (using the same X-ray parent population). Processing of event files
in tiles of size $3.6\times 3.6\deg^{2}$ took place in a similar way as for
real data. In particular the eSASS (version eSASSusers_201009) source finding
algorithm runs over each tile and delivers a list of detections with
associated measurements $\mathcal{L}_{\rm det}$ and $\mathcal{L}_{\rm ext}$.

All clusters simulated in a twin mock are associated to a set of properties,
such as mass $M_{500}$, true redshift $z$, X-ray luminosity $L_{X}$, flux
$f_{X}$ or count-rate $CR$ (all within $R_{500}$, see Table 1). These ‘true’
properties will serve in establishing selection models in Sect. 5. However, we
consider important at this stage to recall that such properties are ‘labels’
associated to simulated sources. The link between those labels and the mock
events depends obviously on the models imprinted in the twin simulation.

###  3.3 Matching input and detected sources

The next step consists in attributing a binary flag to each simulated halo
stating whether it is selected in the cluster sample. This operation relies on
matching the input catalog to entries in the detection list. We explore and
compare two different procedures. The first method, hereafter ‘photon-based’,
takes advantage of the SIXTE simulator as it individually tracks simulated CCD
events back to the source that emitted the photon. The second method,
hereafter ‘position-based’, uses positions and sizes of sources on sky, aided
with prior knowledge of the flux distribution of selected sources.

####  3.3.1 Photon-based matching

The photon-based matching is in fact our baseline method. We refer the reader
to Seppi et al. (2022) for a comprehensive explanation of the procedure. The
input object contributing most photons to the counts making a detection is
identified as the best matching counterpart. If a detection is split into
several sources by the detection software, again preference is given to the
source comprising most photons originating from the identified counterparts.
This enables assigning a unique match (ID_Uniq) to a given detection and to a
given simulated object. This is the definition we take for matches in the
following. Sources with a unique match are flagged as selected ($I_{\rm{\rm
main}}=1$). We note that this procedure involves sources of all types:
simulated AGN, clusters, stars; and detected point-like and extended sources.

####  3.3.2 Positional matching

The positional matching technique is a two-way match between the input and
detection catalogs. It was employed in the context of establishing a selection
function for the eROSITA Final Equatorial-Depth Survey (eFEDS, see Liu et al.,
2022). This technique does not assume a fixed cross-matching radius; instead
it takes into account the sizes of sources. Each simulated halo is listed with
a central coordinate and an angular extent (we take 10% of the virial radius).
Each extended detection is reported with a coordinate, a 1-$\sigma$ error
circle, and an angular extent (we take the core-radius of the best-fit
$\beta$-model). We combine source extents with the source positional
uncertainty; effectively it is similar to spreading location of the source
over a small region. We use the NWAY algorithm (Salvato et al., 2018), since
it is well-suited to cross-matching catalogs in presence of positional
uncertainties. We first take the input catalog as reference and look for
matches in the detection catalog. For each reference source, the algorithm
returns a list of all matches located within a large buffer region of radius 3
arcmin. Each of the matches is assigned a probability $p_{i}$ that it is
indeed a valid association, based on its proximity and positional
uncertainties. A second value $p_{\rm any}$ is computed, representing the
probability for the input source to have at least one counterpart among the
detections. The probabilities $p_{i}$ and $p_{\rm any}$ account for chance
associations, through the source density estimated over a degree-scale region.
Sorting matches by the value of $p_{i}$ provides a ranking of most likely
counterparts. More often than not our catalogs enclose complex configurations
with multiple input sources projected along neighboring line-of-sights, and
detections split in multiple sources. To deal with these effect we update the
probabilities with a prior distribution of the flux of selected sources, as
enabled by the NWAY formalism. The exact shape of this prior is unimportant in
this context, and we simply take it from the simpler pre-launch selection
functions derived in Clerc et al. (2018). As a consequence, NWAY will
preferentially up-weight the probability of brighter sources. Converting
probabilities into binary flags requires setting thresholds on $p_{i}$ and
$p_{\rm any}$. Our experiments showed that the result is rather insensitive to
their exact values. For this first pass (input catalog as reference) we kept
only pairs with $p_{\rm any}>0.9$ and $p_{i}>0.1$. Pairs with the highest
$p_{i}$ among all associations are called primary matches. We reiterate the
above procedure after exchanging the role of the input and detection catalogs.
For this second pass (detection catalog as reference) we do not input any
prior and we chose $p_{\rm any}>0.5$ and $p_{i}>0.1$ to select valid pairs and
primary matches. Finally, matches that are primary in both runs were selected
as solid matches. For each input source, it is flagged as a selected cluster
($I_{\rm{\rm main}}=1$) if it belongs to a solid match. Such a double matching
procedure is more conservative than a single-pass procedure; as a primary
matched source will not be available for a second pair.

####  3.3.3 Comparison of matching techniques

Figure 1: Comparing the outcome of two matching procedures relating simulated
halos (black plain histogram) to detected sources (photon-based matching as a
blue dashed line, position-based matching as a thin orange line). A simulated
halo is flagged as selected ($I_{\rm{\rm main}}=1$) whenever it is (solidly)
matched to an extended source in the detection catalog. The details of the
matching algorithms impact the flux distribution of detected sources,
especially at flux below $\sim 10^{-14}$ erg s-1 cm-2, with deviations up to a
factor 2 in certain bins.

Figure 1 illustrates the comparison between both matching procedures. They
agree well on a wide range of fluxes. Deviations are visible at low flux, i.e.
below a few $10^{-14}$ erg s-1 cm-2. One possible cause (not unique) stems
from the presence of a bright AGN near a relatively dim cluster. An example is
shown in Fig. 2. In photon-based matching, even if there is an extended
detection nearby, it is matched to the input simulated AGN due to the large
amount of photons (and events) it deposits in the vicinity of the extended
detection. In the position-based procedure, the extended source is matched to
the input cluster, since it does not include input AGN sources in the
analysis. The problem of characterising a source as detected in such case is
ill-posed. Another subtle difference in the treatment of the catalogs plays a
role, due to the definition of an extended detection. The position-based
matching takes as input the list of sources with extent likelihood greater
than a predefined threshold (6 for the cosmology sample). The photon-based
matching takes as input the entire list of sources, associates a detection in
the list and only then flags as selected those halos with a detection having
extent likelihood above threshold. Reversing the order of operations has an
impact in crowded regions with multiple split sources, or multiple halos along
neighouring line-of-sights.

Figure 2: Example case of blending, where a bright AGN drives the detection
and the presence of a faint extended source favors its classification as
extended. This figure is a cut-out of a simulated count image, smoothed with a
Gaussian. The red 2- and 5-arcmin circles are centred on a faint simulated
cluster (filled triangle symbol) with flux $\sim 5\times 10^{-14}$ erg s-1
cm-2. Detected sources appear with blue squares. A detection close to the
cluster centre is classified as extended (magenta circle). It also coincides
with a bright simulated AGN (orange ‘x’). The position-based matching
algorithm does consider the simulated cluster as selected ($I_{\rm{\rm
main}}=1$), contrary to the photon-based matching procedure. Two other faint
simulated clusters are shown with open triangles, they have no impact on the
detection.

The comparison exercise described in this section serves in bracketing the
uncertainty associated to the matching procedure. In practice, model
uncertainties in the low count regime are expected to overcome such
variations. The photon-based matching is taken as a baseline for the rest of
the paper.

###  3.4 Training and test samples

Each realization of the twin eRASS1 sky comprises about $10^{6}$ halos, less
than 4% being detected ($I_{\rm{\rm eRASS1}}=1$) and less than 1% being
selected ($I_{\rm{\rm main}}=1$) according to the matching algorithm. The
sample is split into two parts: two thirds are saved for training a model, the
other third is left untouched and will serve to test the validity of the
model. Splitting is performed after shuffling the list of halos at random.

We also create supersets of the simulated sky by concatenating training sets
of eleven realizations, and randomly shuffling their content. Although the
realizations are not independent stricto sensu (they share the exact same
population of objects), this procedure helps in reducing the impact of photon
noise. Each super-training set is $11\times 2/3\simeq 7$ times larger than the
actual eRASS1 sky.

##  4 Selection models with surface brightness profiles

Equation 2 involves three important factors that require modelling. The task
of modelling the selection function consists in building
$p\left(I\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm sel}\right)$, this is
the main interest of this paper. The choice of $\@vec{\Theta}_{\rm sel}$
usually results from a compromise between precision of the selection function
and complexity of the model $p\left(\@vec{\Theta}_{\rm bin},\@vec{\Theta}_{\rm
sel}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm model}\right)$.

A natural choice for the selection parameters $\@vec{\Theta}_{\rm sel}$ is the
collection of pixel values that form the images of clusters. Then
$p\left(I\mathrel{}\middle|\mathrel{}\rm image\right)$ is obtained by feeding
eSASS with an image and applying thresholds on the values ($\mathcal{L}_{\rm
ext}$, EXT) returned by the algorithm; because eSASS is deterministic,
$p\left(I\mathrel{}\middle|\mathrel{}\rm image\right)$ takes either value 0 or
1. The selection function model is thus very precise – it is actually the most
precise model since it exactly reproduces the actual data processing. However,
analysing a typical $10\arcmin\times 10\arcmin$ image with eSASS takes about
three seconds with standard CPUs. Computing a likelihood (Eq. 1) for about 12
000 clusters would then amount to about three CPU$\cdot$hours, which is
unrealistically long for a cosmological analysis. To this cost must be added
the time to compute a model cluster image given $\@vec{\Theta}_{\rm model}$.

Realizing that galaxy cluster images are almost circularly symmetric, we may
accept to lose precision and to reduce the complexity of modelling a cluster.
Let us consider the number of photon hits (counts) deposited by sources onto
the detectors, split in several sky annular apertures around the centre of a
putative dark matter halo. We emphasize that the number of counts should be a
prediction of the model $p\left(\rm
counts\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm model}\right)$. It is not
the outcome of a measurement process – that is available for detected sources
only. It is not the scope of this paper to discuss how such a generative model
is constructed, we simply assume it exists. For constructing our selection
model we rely on the twin simulations, that have kept a record of the origin
of each count deposited in the image.

###  4.1 Illustration with a single feature: the $90\arcsec$ cluster counts

Figure 3: Distribution of clusters in the training set as a function of
$\Sigma$ (units counts arcsecond-2), the average surface brightness in the
$90\arcsec$ radius around their centre. The x-axis is rescaled with $c=2\times
10^{-4}$ arcsecond-2. _Top panel:_ histogram of all simulated clusters (blue)
and histogram of the subset of those found as extended by eSASS (orange).
_Bottom panel:_ dots indicate the ratio of the histograms (empirical selection
rate). Error bars are the 68% confidence range estimated according to App. C.
The vertical dashed line indicates the transition $\hat{\eta}=1.94$ in the
logistic model described by Eqs. 3 and 4 – corresponding to $N\simeq 30$
counts.

Let us initiate our procedure by involving only one feature, namely the number
of cluster counts received in a circular aperture around the (RA, Dec) centre
of the simulated halo. The aperture radius is set to $R=90\arcsec$. We
transform the counts (integer values $N$) with:

$\eta=\ln\left(1+\frac{\Sigma}{c}\right)$ (3)

We set $\Sigma=N/(\pi R^{2})$ and $c=0.0002$ counts arcsecond-2. In our
training set, $N$ ranges from 0 to 3194 counts within $90\arcsec$ aperture.
Fig. 3 displays the histogram of $\eta$ associated to all halos in the
training set, and for those flagged as selected ($I_{\rm{\rm main}}=1$). The
ratio of the two histograms provides an empirical estimate of the probability
of detecting a halo given the number of counts deposited on the detectors. The
empirical probability follows a characteristic ”S-shaped” curve. We perform
logistic regression to fit a model to the probability of detecting a cluster
as extended ($\mathcal{L}_{\rm ext}>3$) given the value $\eta$. We use the
scikit-learn implementation (Pedregosa et al., 2011) with the Broyden-
Fletcher-Goldfarb-Shanno optimizer and fit two coefficients, namely the
intercept $w_{0}$ and the slope $w_{\eta}$ of the linear model
$f(\eta)=w_{0}+w_{\eta}\eta$. The model probabilities are such that

$p\left(I_{\rm{main}}\mathrel{}\middle|\mathrel{}\eta\right)=\left(1+e^{-f(\eta)}\right)^{-1}.$
(4)

The cost function to minimize during the fitting procedure is the log-loss
$\mathcal{C}(w_{0},w_{\eta})$, its expression is shown in Eq. 22.

Thanks to the one-to-one relation between $\eta$ and $N$ (Eq. 3), we thus
obtain a model $p\left(I_{\rm{main}}\mathrel{}\middle|\mathrel{}N\right)$. We
can interpret the values of the best-fit coefficients as follows: $w_{\eta}$
governs the sharpness of the transition in the S-shaped curve,
$\hat{\eta}=(-w_{0}/w_{\eta})$ is the value at which this transition occurs.
In our experiment, we find $w_{\eta}=4.59\pm 0.04$ and $w_{0}=-8.92\pm 0.06$,
leading to $\hat{\eta}=1.94$. This corresponds to a transition of the
”S-shaped” curve taking place at $N\simeq 30$ cluster counts in the
$90\arcsec$ circular region. Positiveness of $w_{\eta}$ reflects the fact that
detection probability increases with number of counts $N$. The computation of
coefficient uncertainties is detailed in App. A.

###  4.2 A model using light profiles from all sources

The model presented so far involves only one feature, it is too simplistic to
explain subtle variations in the probability of selecting a given cluster. The
model improves by involving more of the quantities available from the twin
simulations, namely the number of detector hits deposited by each of the five
sky source components: the cluster of interest, Active Galactic Nuclei (AGN),
stars, back- and foreground and other clusters. Each of these contribute to
the count image and deposit photons into 7 annular regions around the central
coordinate of the cluster. Radial boundaries (units arcseconds) are $0-20$,
$20-40$, $40-60$, $60-90$, $90-120$, $120-150$ and $150-180$. For each
simulated cluster, we thus construct a $5\times 7=35$-element feature vector
<EMAIL_ADDRESS>and perform logistic regression,
involving a total of 36 coefficients $\\{w_{j}\\}_{j=0..35}$ (one for each
feature, plus the intercept). Introducing $\Sigma_{j}$ the average surface
brightness in an annulus, we write:

$f(\eta)=w_{0}+\sum_{j=1}^{35}w_{j}\eta_{j}{\ \rm,with\ \
}\eta_{j}=\ln\left(1+\frac{\Sigma_{j}}{c}\right)$ (5)

One may interpret the value of the coefficients $w_{j}$ as the sensitivity of
the detection rate $p\equiv
p\left(I_{\rm{main}}\mathrel{}\middle|\mathrel{}\@vec{\eta}\right)$ to a small
variation of the surface brightness of one component in one annulus. Indeed we
have:

$\left.\frac{\partial p}{\partial\eta_{j}}\right|_{\eta_{k}}=p(1-p)w_{j}$

where the derivative is performed at constant $\eta_{k}$ ($k\neq j$). It
yields:

$\left.\frac{\partial
p}{\partial\Sigma_{j}}\right|_{\Sigma_{k}}=\frac{p(1-p)}{\Sigma_{j}+c}w_{j}$
(6)

In other terms, everything else maintained fixed, a small relative variation
$\epsilon=\Delta\Sigma_{j}/\Sigma_{j}$ of the surface brightness of one
component in one given annulus leads to a variation of the detection
probability proportional to $w_{j}\epsilon$. In particular, positive
coefficients relate to those features that marginally contribute to selecting
a cluster, and conversely for negative coefficients. Fig. 4 conveniently
represents the value of the 35 coefficients $w_{j}$ $(j\neq 0)$ split into
each component and annulus; the last coefficient not shown on this figure is
$w_{0}=-8.82\pm 0.08$.

Figure 4: Representation of the 35 coefficients $w_{j}$ of a logistic
regression model $p\left(I_{\rm{main}}\mathrel{}\middle|\mathrel{}\rm
counts\right)$, trained to predict whether a cluster is selected in the
primary cluster sample. The 35 features are surface brightness in 7 radial
annuli, associated to the 5 components indicated in legend (counts from the
galaxy cluster of interest, from neighbouring AGN, fore- and background
counts, counts from foreground stars, counts from other neighbouring galaxy
clusters.) High absolute value of a coefficient indicates high importance of
the associated feature (Eq. 6).

Keeping in mind that $w_{j}$ are proportional to the marginal increase of the
detection probability, comparing together values of coefficients provides an
approximate way of assessing the ‘importance’ of a feature in promoting
detectability of a cluster. For instance, it appears that a small increase of
surface brightness in the $20-40\arcsec$ annulus from a given cluster has the
strongest impact on the detectability of this cluster. The fact that it is
twice as important as the $0-20\arcsec$ surface brightness is not surprising.
A detected source ($I_{\rm{\rm eRASS1}}=1$) must be classified as extended in
order to be selected ($I_{\rm{\rm main}}=1$); from this perspective, it is
much more profitable to increase the counts beyond the eRASS point-spread
function radius ($\sim 30\arcsec$). This argument may appear more clearly from
Fig. 32, which shows the coefficients $w_{j}$ of a model trained to predict
the detection (not the selection) of a cluster. Clearly in this case the
central $20\arcsec$ surface brightness is the most relevant feature
contributing to detectability of a cluster.

The effect of star-emitted photons is barely constrained, mostly due to their
paucity in the simulation. AGN photons within one arcmin of the centre tend to
marginally increase the detection rate, while those located beyond 1 arcmin
tend to decrease the detection rate. This is because the source detection
algorithm would put a mask and remove photons useful for cluster detection. An
excess of photons from instrumental and astrophysical fore- and background
within $1.5\arcmin$ of a cluster tend to increase its probability of being
selected; while at larger distances they decrease its probability by a large
factor. The impact of neighbouring clusters is similar, although less
pronounced.

###  4.3 Internal validation of the selection models

The previous interpretation of the $w_{j}$ coefficients in terms of marginal
probabilities should not hide it is only relevant to a specific model, here a
logistic model, that is not a perfect classifier. We assess performance of
this model on the left-apart test sample. Fig. 5 represents the so-called
‘precision-recall’ curve. This curve is obtained by scanning values of a
threshold set to convert probabilities
$p\left(I_{\rm{main}}\mathrel{}\middle|\mathrel{}\rm counts\right)$ into a
binary classification $I_{\rm{main}}$. Once this threshold is set, one
retrieves the number of true positives (TP), false positives (FP) and false
negatives (FN) by comparison with the actual detection flag in the test
sample. The higher the recall222Recall (also known as sensitivity) and
precision are respectively equivalent to completeness and purity, defined with
respect to the true classification in the test sample. However we prefer using
the former terminology in order to avoid possible confusion with the usual
measures of completeness and purity of a certain sample, defined with respect
to the entire population of clusters. $TP/(TP+FN)$ and the precision
$TP/(TP+FP)$, the better performance of a model. For any value of the
threshold, the model using 35 surface brightness indicators performs
significantly better than the model using only the 1.5-arcmin cluster counts
as input. We have also displayed the curve obtained from a logistic model
using two features, namely the 1.5-arcmin cluster and background counts. Its
performance lies between the simple and the sophisticated ($7\times 5$)
models. This agrees with common sense, in that adding more detailed
information provides more precise models.

A second test is shown on Fig. 6. This evaluation is performed by binning test
clusters by their value of
$p\left(I_{\rm{main}}\mathrel{}\middle|\mathrel{}\rm counts\right)$ returned
by the model. The actual fraction of objects selected in the sample is
reported on the vertical axis. The 35-feature model produces the plain
envelope, in satisfactory agreement with the one-to-one line. This result
indicates reliable probabilities delivered by the model in a statistical
sense. The probabilities from the model using only one cluster count indicator
(round dots and dashed envelope) are less reliable. This is reflected in Fig.
3 where the ”S”-shaped curve does not reach unity around values of $\eta\simeq
2.5-4$ (corresponding to $N\simeq 60-270$ counts in the $90\arcsec$ aperture).

Figure 5: Performance of three different logistic models
$p\left(I_{\rm{main}}\mathrel{}\middle|\mathrel{}\rm counts\right)$ to predict
a cluster as selected given the values of surface brightness in radial bins.
The plain black line is obtained with a model using surface brightness values
in 7 annuli for each of the 5 sky components. The dot-dashed line is for a
model taking as input both cluster and background counts within $90\arcsec$.
The dashed grey line stands for a model using as input only the $90\arcsec$
cluster counts. Each precision-recall curve results from model evaluations on
the test sample. It is obtained by varying the threshold over which a cluster
should be considered as selected and counting the number of true positives
(TP), false negatives (FN) and false positives (FP). Figure 6: Reliablity of
two different logistic models
$p\left(I_{\rm{main}}\mathrel{}\middle|\mathrel{}\rm counts\right)$ to predict
the probability of a cluster being selected given the values of surface
brightness in radial bins. Each series of points is obtained with the test
sample, by comparing probabilities predicted by the model (horizontal axis)
with the actual fraction of selected objects (vertical axis) in each bin of
probability. The envelopes materialise the 68% uncertainties (App. C). An
ideal model would align along the 1:1 curve (thin dotted line). Black squares
and plain lines correspond to the model taking $7\times 5$ parameters as
input, grey dots and dashed lines for the model using only cluster counts
within $90\arcsec$ to make its prediction.

Despite its apparent simplicity, a linear logistic model using circularly
averaged surface brightness as training features does output probabilities
that are representative of the actual selection probabilities. It is therefore
well-suited to selection function problems involving populations of clusters.
However, an average precision of $77\%$ (Fig. 6) indicates that it is not
well-suited to the classification on an individual basis. This floor
performance is due to using circularly averaged profiles, hence neglecting
two-dimensional effects (e.g. blending, masking, gradients) in the selection
process. Furthermore, logistic regression is a simple linear model. An
interesting development would involve more sophisticated models such as neural
networks (e.g. building upon the multi-layer perceptron). Their use may be
able to add the right level of non-linear feedback and complexity to
accurately reproduce the selection process, at the expense of explainability
and simplicity.

A selection function model based on surface brightness as above is attractive
due to its use of observable features. In the context of cosmological studies
(as idealized in Eq. 2), it is necessary to construct a second, independent
model able to predict the radial surface brightness profile of five
components: clusters, AGN, stars, background and neighbouring clusters. Such a
complex model $p\left(\rm counts\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm
model}\right)$ represents a computational bottleneck.

##  5 Selection models with intermediate variables

We may instead choose $\@vec{\Theta}_{\rm sel}$ such as to reduce this
complexity, at the cost of a less precise selection model. We consider
variables (i.e. selection model features) such as cluster mass, redshift,
luminosity, size, central emissivity, etc. We also consider exposure time,
background level, galactic hydrogen column density. For some of these
variables (e.g. flux, luminosity) it is more convenient to manipulate their
logarithmic values.

Having fixed the set of selection variables $\@vec{\Theta}_{\rm sel}$, we
build a model assuming that the detection probability varies slowly over the
range of interest of these (rescaled) features. However, we do not make prior
assumption on the exact shape, nor on monotonicity of the function. A kernel-
based model is well-suited to describe this problem, as it reflects the fact
that similarly-looking clusters should have similar detection probabilities.
We use Gaussian Process classification (GPC) to build our model, using the
SVGP implementation in the GPy package (GPy, since 2012). Details on the
procedure and algorithm are discussed in App. D. The very nature of GPC allows
to issue not only an estimate of
$p\left(I\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm sel}\right)$, but also
a range of statistical uncertainty. Our implementation does not provide
extrapolation properties and a model returns a default value ($p=0.5$) when
$\@vec{\Theta}_{\rm sel}$ is outside the domain where the training set lives.
For the purpose of cosmological analysis it is not deemed an issue, as long as
care is taken in handling the integration bounds in Eq. 2.

In the following we highlight three models relevant to this paper. In
particular we demonstrate the interest of a model based on count-rate $CR$,
which was selected for the eRASS1 cosmological analysis (Ghirardini et al.,
submitted). In App. E we describe several other models, notably those
including cluster morphology in their input parameters.

###  5.1 Mass and redshift

Setting $\@vec{\Theta}_{\rm sel}=\\{M_{500c},z\\}$ allows to gain insights on
the selection process in the fundamental parameter space of cluster cosmology.
Fig. 7 is a visualization of a model adjusted to predict the selection of a
cluster with extent likelihood above 3. The probability of selecting a cluster
located at $z=0.3$ monotonically increases with mass, with a transition taking
place at $M_{500c}\simeq 5\times 10^{14}$ $M_{\odot}$.

Figure 7: Representation of the model
$p\left(I_{\rm{main}}\mathrel{}\middle|\mathrel{}M_{500c},z\right)$ predicting
a cluster to be detected and selected with an extent likelihood above 3. The
explanatory variables (features) are labels attached to simulated clusters in
the twin simulation, standing for galaxy cluster $M_{500c}$ mass and
cosmological redshift. The bottom-left panel represents contours of equal
detection probability, labelled in steps of 0.1. Both one-dimensional curves
(top-left and bottom-right panels) are slices through the function displayed
in the bottom-left corner, at fixed $z=0.3$ and $M_{500c}=3\times 10^{14}$
$M_{\odot}$ (indicated with dotted lines).

Despite its formal interest, this model is not of practical use in context of
cosmology studies. The two features entering this model (mass and redshift)
are mere labels attached to clusters in the twin simulation. The model is
therefore valid only under assumption of physical models implemented in the
simulation. In particular, it is only valid for the set of scaling relations
and cosmological parameters (Planck Collaboration et al., 2016) imprinted in
the simulation.

###  5.2 Luminosity and redshift

Setting $\@vec{\Theta}_{\rm sel}=\\{L_{X},z\\}$ advantageously removes the
strong dependence of the model upon the $M_{500}\rightarrow L_{X}$ relation
that is imprinted in the twin simulation. It is therefore up to the
cosmologist to model the relation
$p\left(L_{X}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm model}\right)$
with (most likely) parameterized scaling laws. A visual representation of the
model is shown in Fig. 8. As expected, clusters brighter and closer to us are
more likely to be selected.

Figure 8: Representation of the model
$p\left(I_{\rm{main}}\mathrel{}\middle|\mathrel{}L_{X},z\right)$ predicting a
cluster to be detected and selected with an extent likelihood above 3. The
explanatory variables (features) are labels attached to simulated clusters in
the twin simulation, standing for galaxy cluster luminosity $L_{X}$ measured
in in the 0.5–2 keV energy band at the cluster rest-frame (units erg s-1,
rescaled in $\log_{10}$) and for cosmological redshift $z$. Other details are
similar as in Fig. 7.

Although useful, this model cannot account for variations of the detection
probability as a function of sky position: galactic absorption, background and
exposure time values all have an impact on cluster detection. Offering these
degrees of freedom in the selection function allows one to detach from the
assumptions imprinted in the twin simulation. We have trained a second model
with those three features, using the GPC formalism. Fig. 9 is a partial
visualisation of the model output, all following quantities but exposure time
$T_{\rm exp}$ being fixed at $z=0.3$, $N_{H}=3\times 10^{20}$ cm-2, background
brightness $5.2\times 10^{-15}$ erg s-1 cm -2 arcmin-2 in the 0.3–2.3 keV
energy band, $L_{X}=10^{43}$ or $10^{44}$ erg s-1. The dependence of the model
output on the exposure time appears clearly on this figure and thus the model
is more precise in delivering a selection probability for a given cluster.

Figure 9: Representation of the models
$p\left(I_{\rm{main}}\mathrel{}\middle|\mathrel{}L_{X},z,N_{H},T_{\rm
exp},{\rm bkg}\right)$ and
$p\left(I_{\rm{cosmo}}\mathrel{}\middle|\mathrel{}L_{X},z,N_{H},T_{\rm
exp},{\rm bkg}\right)$ for fixed values of cluster redshift $z$, two different
values for the luminosity $L_{X}$, and local galactic absorption $N_{H}$, a
nominal background level and a range of exposure times $T_{\rm exp}$ (along
the x-axis). The shaded regions indicate the approximate 68% confidence range
output of the model.

###  5.3 Unabsorbed count rate and redshift

In the course of the development of the eROSITA cosmology pipeline, trade-off
discussions led to setting $\@vec{\Theta}_{\rm sel}$ as a five-component
feature, namely: the cluster redshift $z$, sky position-dependent quantities
$\mathcal{H}\equiv$ ($N_{H}$, $T_{\rm exp}$ and background) and the cluster
$CR$. Here $CR$ represents the unabsorbed, 0.2–2.3 keV survey-average count-
rate collected within a $R_{500}$ aperture and unaffected by shot noise. It
derives from the formula:

$CR=\frac{L_{X}}{4\pi d_{L}(z)^{2}ECF(z,k_{B}T,N_{H}=0)}\ \ \ {\rm(units\
s}^{-1}{\rm)}$ (7)

In this expression, $ECF$ is the Energy Conversion Factor from the $0.5-2$ keV
band in the cluster rest-frame to the $0.2-2.3$ keV band in the observer
frame. It is computed for a APEC emission spectrum at temperature $k_{B}T$,
redshift $z$, unabsorbed (hence $N_{H}=0$) and folded through the survey-
averaged instrumental response of the 7 telescopes. This quantity is readily
available for each simulated cluster in the twin simulation, since we know its
redshift and temperature. The luminosity-distance $d_{L}$ is computed using
the same cosmology as for the twin simulation. Fig. 10 represents a few slices
through the multi-dimensional model
$p\left(I_{\rm{cosmo}}\mathrel{}\middle|\mathrel{}CR,z,N_{H},T_{\rm exp},{\rm
bkg}\right)$, fixing the value of $N_{H}$ and background level as previously,
for two different redshift values and a range of count-rate ($CR$) values. The
dependence on count-rate shows the usual ”S”-shaped curve and larger exposure
times imply higher sensitivities. The slight dependence on redshift is less
intuitive: all other quantities being fixed, a source at higher redshift is
more compact, thus its surface brightness is concentrated over a smaller area,
making its detection comparatively easier.

Figure 10: Representation of a model
$p\left(I_{\rm{cosmo}}\mathrel{}\middle|\mathrel{}CR,z,N_{H},T_{\rm exp},{\rm
bkg}\right)$ for two values of cluster redshift $z=0.2$ and $z=0.5$, fixed
local galactic absorption $N_{H}$, a nominal background level, two values of
exposure times $T_{\rm exp}$ (units s) and a range of count-rate values
(x-axis). The shaded regions indicate the approximate 68% confidence range
output of the model.

The selection model obtained this way is quasi-independent of the cosmology
imprinted in the twin simulation. It is also capable of reflecting the
variations of selection depth over the sky. The cluster model is relatively
simple (one cluster is represented by one count-rate and one redshift), which
makes the computation of
$p\left(CR,z\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm model}\right)$ fast
and easy in likelihood (Eq. 2). However this simplicity comes at a cost, since
all variations of the selection due to e.g., cluster morphology, are
marginalized over, relying on the distribution of morphologies in our twin
simulation.

###  5.4 Internal validation of the selection models

We now quantify the absolute performance of models relying on intermediate
variables $\@vec{\Theta}_{\rm sel}$. By doing so we will also assess the
relative performance between models involving various parameters as input. We
use our test sample to make such tests and compare the predicted model
outcomes to the actual detection flags in the sample. Similarly as in Sect. 4
we will highlight two performance tests: the precision-recall curve and the
calibration curve.

Fig. 11 compares five of our models predicting the presence of a cluster in
the eRASS1-main sample. They differ from each other by their input parameters
(features $\@vec{\Theta}_{\rm sel}$). The model involving only mass and
redshift has quite a low precision and recall (average precision of $34\%$).
Predicting the detectability of a cluster requires indeed more detailed
characterisation. Changing to a model taking luminosity and redshift notably
improves the overall performance, since luminosity is more tightly linked to
photon counts than mass. However, it still has an average precision of $56\%$,
unsatisfactory enough for detailed analyses. Including sky position-dependent
parameters $\mathcal{H}=(N_{H},T_{\rm exp},{\rm bkg})$ dramatically boosts the
classification performance (average precision of $67\%$). As expected, adding
information on the local background, exposure time and foreground absorption
increases the capability of the model. An additional performance gain is
obtained by adding a morphological parameter, here the $EM_{0}$ parameter (as
presented in App. E). The average precision increases to $71\%$, closer to
that obtained with the $7\times 5$ surface brightness model depicted in Fig.
5. Finally, a model involving (absorbed) flux, apparent $R_{500}$ radius and
$\mathcal{H}$ shows good performance with an average precision reaching
$68\%$.

Figure 11: Precision-recall curves obtained with selection models trained to
predict the selection of a cluster in the primary sample. The input parameters
vary from one model to the other, as indicated in legend. A model involving
luminosity, redshift, morphological features and local sky information
performs best (yellow dot-dashed curve).

A similar assessment is provided in Fig. 12, now for models predicting the
presence of a cluster in the eRASS1 cosmology sample. The baseline model used
for the eRASS1 cosmological analysis relies on $CR$ (count-rate) and redshift,
as well as on sky local information. Its good performance (average precision
70%) is comparable to a model involving luminosity and redshift. Again, adding
morphological information (here via $EM_{0}$, see App. E) provides a
noticeable enhancement of the model precision (average precision $74\%$).

Figure 12: Similar as Fig. 11, comparing three models predicting the presence
of a cluster in the cosmology. The model represented with the black curve
corresponds to the baseline selection model used in the analysis of Ghirardini
et al. (submitted).

We now assess the reliability of the probabilistic output of the models. In
Fig. 13 we compare the same five models for the primary catalog selection
function. All curves lie very close to the one-to-one line: the probability
outputs reflect well the actual fraction of positive detections in the test
catalog. The result for the cosmology sample classifiers appear in Fig. 14 and
they all are very well aligned along the one-to-one line.

Figure 13: Evaluating the reliability of five models predicting the
probability of a cluster to appear in the eRASS1 primary cluster catalog
(x-axis). The y-axis represents the actual fraction of selected objects in the
test sample. The bottom panel shows the residuals (difference) of each curve
with respect to the one-to-one line. Shaded regions indicate 68% confidence
intervals in each bin of probability. Figure 14: Similar as Fig. 13, comparing
three models predicting the presence of a cluster in the cosmology sample.

The tests depicted above illustrate good prediction performance of the models,
in the sense that probability outputs correctly reflect the actual selection
probabilities. This performance indicator is valid for the ensemble of
clusters in the simulated eRASS1 sky. In the course of our model production
efforts, we have monitored several other indicators in order to assess the
performance and reliability of the models. Among them were: the Receiver-
Operator Curve (ROC) comparing recall and fall-out rates $FP/(FP+TN)$; the
Detection-Error Tradeoff Curve (DET) comparing fall-out and miss rates
$FN/(FN+TP)$ and the correlation coefficient suggested by Matthews (1975). All
together, these performance indicators helped in selecting the best models to
be used in the cosmological analysis.

##  6 External validation of the selection function

We have described several models and we validated their performance on a
simulated test set. This test set is produced in a very similar manner as the
training sample, splitting the original twin simulation set in two parts.
Additional tests are required using new samples independent from the twin
simulation. We aim to test eRASS1 selection models with i) a sample of
clusters selected from a much deeper dataset (eFEDS) and ii) two SZ-selected
sample selected from South Pole Telescope surveys of different depths:
SPTpol-100d and SPT-Early Cluster Sample (SPT-ECS). Selection biases affecting
these samples must be taken into account; moreover we have to link their
measured properties to eRASS1 selection variables. We present our formalism in
the most generic manner, then we apply it to the specific cases of eFEDS and
SPT-SZ.

###  6.1 General presentation of the formalism

Given two surveys, $R$ taken as a reference and the other $T$ to be tested, we
wish to validate our selection function models for both $R$ and $T$. We do so
by comparing the populations selected by each of the survey. In the optimal
case where $R$ is much more complete (e.g. deeper) than $T$, the problem does
not require precise knowledge of the selection process leading to catalog $R$.
In general both catalogs comprise a biased selection of the true underlying
population of clusters, the latter we take as prior in the following analysis.

Extending the formalism laid out in previous work (e.g. Grandis et al., 2020)
we calculate for each object $i$ in the reference catalog an associated
probability
$P_{i}=p\left(I_{\rm{T}}\mathrel{}\middle|\mathrel{}I_{\rm{R}},\@vec{\Theta}_{\rm
meas}\right)$. There, $I_{\rm{T}}$ and $I_{\rm{R}}$ are random boolean
variables indicating whether an object is listed in catalog $T$ and $R$
respectively. $\@vec{\Theta}_{\rm meas}$ is a collection of measured
properties attached to the cluster: its measured flux, redshift, luminosity,
extent, etc. are likely candidates. Using Bayes’ rule we obtain:

$P_{i}=p\left(I_{\rm{T}}^{i}\mathrel{}\middle|\mathrel{}I_{\rm{R}}^{i},\@vec{\Theta}_{\rm
meas}\right)=\frac{p\left(I_{\rm{T}}^{i},I_{\rm{R}}^{i},\@vec{\Theta}_{\rm
meas}^{i}\right)}{p\left(I_{\rm{R}}^{i},\@vec{\Theta}_{\rm
meas}^{i}\right)}=\frac{N_{i}}{D_{i}}$ (8)

where:

$N_{i}=\int\mathrm{d}\@vec{\Theta}_{\rm model}p\left(\@vec{\Theta}_{\rm
model}\right)p\left(I_{\rm{T}}^{i},I_{\rm{R}}^{i},\@vec{\Theta}_{\rm
meas}^{i}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm model}\right)$

and

$D_{i}=\int\mathrm{d}\@vec{\Theta}_{\rm model}p\left(\@vec{\Theta}_{\rm
model}\right)p\left(I_{\rm{R}}^{i},\@vec{\Theta}_{\rm
meas}^{i}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm model}\right)$

From these equations, $D_{i}$ may be interpreted as the number of objects
expected in the reference catalog with measured properties $\@vec{\Theta}_{\rm
meas}^{i}$; $N_{i}$ may be interpreted as the number of objects expected in
both the reference and test catalogs with these same measured properties. In
this perspective, $P_{i}$ is a model for the fraction of reference objects
listed with some measured properties $\@vec{\Theta}_{\rm meas}^{i}$ that are
also present in the test sample.

These expressions may be expanded further by using selection variables
$\@vec{\Theta}_{\rm sel}$ for each of the surveys. We have:

$p\left(I_{\rm{T}}^{i},I_{\rm{R}}^{i},\@vec{\Theta}_{\rm
meas}^{i}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm
model}\right)=\int\mathrm{d}\@vec{\Theta}_{{\rm
sel},R}\mathrm{d}\@vec{\Theta}_{{\rm sel},T}\mathrm{d}\@vec{\alpha}\\\ \times
p\left(\@vec{\Theta}_{{\rm sel},T},\@vec{\Theta}_{{\rm
sel},R},\@vec{\alpha}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm
model}\right)\\\ \times
p\left(I_{\rm{R}}^{i}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{{\rm
sel},R}\right)p\left(I_{\rm{T}}^{i}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{{\rm
sel},T}\right)p\left(\@vec{\Theta}_{\rm
meas}^{i}\mathrel{}\middle|\mathrel{}\@vec{\alpha}\right)$

and similarly:

$p\left(I_{\rm{R}}^{i},\@vec{\Theta}_{\rm
meas}^{i}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm
model}\right)=\int\mathrm{d}\@vec{\Theta}_{{\rm
sel},R}\mathrm{d}\@vec{\alpha}p\left(\@vec{\Theta}_{{\rm
sel},R},\@vec{\alpha}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm
model}\right)\\\ \times
p\left(I_{\rm{R}}^{i}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{{\rm
sel},R}\right)p\left(\@vec{\Theta}_{\rm
meas}^{i}\mathrel{}\middle|\mathrel{}\@vec{\alpha}\right)$

The last two expressions are valid as long as selection in the test survey
only depends on variable $\@vec{\Theta}_{{\rm sel},T}$ and selection in the
reference survey only depends on variable $\@vec{\Theta}_{{\rm sel},R}$. We
have introduced a new, intermediate variable $\@vec{\alpha}$ to reflect the
dependence of the measurement vector. For instance, if the measurement vector
is X-ray flux, $\@vec{\alpha}$ may represent the cluster luminosity and
redshift.

The interest of this decomposition lies in it involving the selection
functions $p\left(I\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm
sel}\right)$, which are tabulated for each survey. It also accounts for
covariances in the selection and measurement processes, through the
multivariate distribution of $(\@vec{\Theta}_{{\rm sel},T},\@vec{\Theta}_{{\rm
sel},R},\@vec{\alpha})$ given the values of $\@vec{\Theta}_{\rm model}$.

The probability $P_{i}$ (Eq. 8) is the main outcome of this model. It is
calculated independently for each object listed in the reference catalog. By
comparing this probability to the actual presence of a match in the test
catalog, we obtain powerful diagnostics on the validity of the models, and the
presence of outliers in the population. In particular, we will denote by
‘surprising detections’ those objects listed in $R$ and detected in the test
catalog despite their low forecast probability ($P_{i}<0.025$ and
$I_{\rm{T}}^{i}=1$). We will denote ‘missed objects’ those listed in $R$ and
not detected in the test catalog despite their high forecast probability
($P_{i}>0.925$ and $I_{\rm{T}}^{i}=0$).

Our procedure to find matches between the reference and the test catalogs is
relies on simple positional match. For efficiency reason we trim the reference
catalog to the approximate sky footprint of the test sample and conversely. We
then search for symmetric 2-arcmin matches between the samples, keeping only
the best match, that is the closest distance. Such a procedure does not
require redshifts to be identical in both catalogs. Our choice of the maximal
separation corresponds to the most sensible cross-matching radius for surveys
and instruments considered in this study (Bulbul et al., 2022).

Beyond the interpretation of individual $P_{i}$ values, we also base our
evaluation on three aggregated diagnostics. First, from repeated Bernoulli
samples of $P_{i}$ values, we obtain a distribution for the expected number of
matching entries between both catalogs, and its associated uncertainty.
Comparing the actual number of matches to this distribution provides a bulk
validation test of the model. The second diagnostic consists in grouping
clusters by their $P_{i}$ values in bins of finite width, and in computing the
actual fraction of matches in each bin. This diagnostics informs on the
reliability of the $P_{i}$ values, hence on the underlying model. Our third
diagnostic empirically tests for leakage in our model, by introducing two
numbers representing the fraction of objects that should be in the test
catalog, but they are not ($\delta_{1}$) and the fraction of objects that
should not be in the test catalog, but they are present ($\delta_{2}$). We
compute posterior distributions for $\delta_{1}$ and $\delta_{2}$ by sampling
the following likelihood (product of independent Bernoulli likelihoods):

$\mathcal{L}(\delta_{1},\delta_{2})=\prod_{i\,\in\,R}\left[P_{i}(1-\delta_{1})+(1-P_{i})\delta_{2}\right]^{I_{\rm{T}}^{i}}\\\
\times\left[P_{i}\delta_{1}+(1-P_{i})(1-\delta_{2})\right]^{1-I_{\rm{T}}^{i}}$
(9)

The product runs over all clusters in the reference catalog, and
$I_{\rm{T}}^{i}$ equals 1 if the cluster is matched in the test catalog, it
equals 0 otherwise. If the posterior distribution for $\delta_{1}$
significantly excludes zero, it indicates that part of the clusters predicted
to be in the test catalog actually escape our model and are not detected for
some reason. Similarly, a distribution of $\delta_{2}$ significantly away from
zero hints towards our model not capturing properly the number of undetected
objects. In this computation we assume flat prior distribution for
$\delta_{1}$ and $\delta_{2}$ bounded in the [0, 1] interval. We sample the
posterior with a Monte-Carlo Markov Chain algorithm provided in the emcee
package (Foreman-Mackey et al., 2013).

Using the formalism depicted in this section, we have run several study cases:

  1. 1.

eFEDS as reference, eRASS1-cosmo as test (Sect. 6.2);

  2. 2.

eRASS1-cosmo as reference, eFEDS as test (Sect. 6.3);

  3. 3.

eFEDS as reference, eRASS1-primary as test (Sect. 6.4);

  4. 4.

SPTpol-100d as reference, eRASS1-cosmo as test (Sect. 6.5);

  5. 5.

eRASS1-cosmo as reference, SPTpol-100d as test (Sect. 6.6).

  6. 6.

SPTpol-ECS as reference, eRASS1-cosmo as test (Sect. 6.7).

  7. 7.

eRASS1-cosmo as reference, SPTpol-ECS as test (Sect. 6.8).

Table 2 provides a summary of the setup for each test, along with the number
of matched entries between the catalogs.

Table 2: Setup summary of the external validation tests. The last column
indicates the expected number of matching entries $\langle N\rangle$ predicted
by our model, together with its uncertainty.

Referenceaa$a$Number in parentheses indicate the number of entries in the reference catalog selected in the approximate common sky footprint. | Test |  $N_{\rm match}$bb$b$The number of matching entries $N_{\rm match}$ results from $2\arcmin$ coordinate positional matching. | $\@vec{\Theta}_{{\rm sel},R}$ | $\@vec{\Theta}_{{\rm sel},T}$ | $\@vec{\Theta}_{\rm meas}$ | $\alpha$ | $\langle N\rangle$
---|---|---|---|---|---|---|---
eFEDS (531) | Cosmo | 35 | $L_{X},z,EM_{0},T_{\rm exp}$ | $CR,z,\mathcal{H}$ |  F_300kpc, $z_{\rm eFEDS}$ | $L_{X},z$ | $19.5\pm 3.2$
Cosmo (85) | eFEDS | 36 | $CR,z,\mathcal{H}$ | $L_{X},z,EM_{0},T_{\rm exp}$ |  F500_0520, $z_{\lambda}$ | $L_{X},z$ | $35.7\pm 1.3$
eFEDS (531) | Main | 61 | $L_{X},z,EM_{0},T_{\rm exp}$ | $CR,z,\mathcal{H}$ |  F_300kpc, $z_{\rm eFEDS}$ | $L_{X},z$ | $30.5\pm 4.2$
SPTpolcc$c$SPTpol is the SPTpol-100d catalog presented in Huang et al. (2020), keeping only entries above $z=0.25$. (65) | Cosmo | 8 | $\zeta_{\rm SPTpol}$ | $CR,z,\mathcal{H}$ |  $\xi$, $z_{\rm SPT}$ | $L_{X},z,\zeta_{\rm SPTpol}$ | $11.4\pm 2.4$
Cosmo (19) | SPTpol | 8 | $CR,z,\mathcal{H}$ | $\zeta_{\rm SPTpol}$ |  F500_0520, $z_{\lambda}$ | $L_{X},z,\zeta_{\rm SPTpol}$ | $7.0\pm 1.3$
SPTecsdd$d$SPTecs is the SPTpol-ECS catalog presented in Bleem et al. (2020), keeping only entries with $\xi>5$ and $z>0.25$. (163) | Cosmo | 103 | $\zeta_{\rm ECS},\gamma_{\rm field}$ | $CR,z,\mathcal{H}$ |  $\xi$, $z_{\rm SPT}$ | $L_{X},z,\zeta_{\rm ECS}$ | $103\pm 4$
Cosmo (968) | SPTecs | 101 | $\zeta_{\rm ECS},\gamma_{\rm field}$ | $CR,z,\mathcal{H}$ |  F500_0520, $z_{\lambda}$ | $L_{X},z,\zeta_{\rm ECS}$ | $65.4\pm 5.8$

333

###  6.2 Testing the eRASS1 cosmological cluster sample against eFEDS

The deep eFEDS survey consists of a 140 deg2 area scanned by eROSITA during
the CalPV phase with exposure about 10 times that of the average eRASS1
survey. Specifically, the average net exposure of eFEDS is about 1200 s, while
in eRASS1 it amounts to a median 80 s in the same field (in both cases
accounting for vignetting).

We opt for a population model described by $\@vec{\Theta}_{\rm
model}=(M_{500},z)$ representing the (true) mass within $R_{500c}$ and the
(true) redshift of halos. The prior distribution $p\left(\@vec{\Theta}_{\rm
model}\right)$ is the halo mass function
$\mathrm{d}n/\mathrm{d}M_{500}\mathrm{d}z$. Its exact shape in this specific
experiment is not of critical relevance, given the wide sensitivity difference
between eFEDS and eRASS1. We assume it follows the fit of Tinker et al.
(2008).

The eFEDS cluster sample (Liu et al., 2022) comprises 542 candidate galaxy
groups and clusters, among them 531 are in the eRASS1 catalog footprint. Each
entry is associated to a measured redshift $z_{\rm eFEDS}$ (column z in the
catalog) and to an estimated $0.5-2$ keV flux within a 300 kpc aperture
(column F_300kpc, expressed in units erg s-1 cm-2). These two elements
constitute our vector of measured features $\@vec{\Theta}_{\rm meas}$. We
construct a generative model for flux and redshift that depends on
$\@vec{\alpha}=(L_{X},z,M_{500})$, where $L_{X}$ is a cluster (true) 0.5-2 keV
luminosity integrated in a cylinder of radius $R_{500}$. Our generative model
further assumes independence of flux and redshift measurements, i.e.:

$p\left(\@vec{\Theta}_{\rm
meas}\mathrel{}\middle|\mathrel{}\@vec{\alpha}\right)=p\left(z_{\rm
eFEDS}\mathrel{}\middle|\mathrel{}z\right)\times
p\left(\ln\texttt{F\\_300kpc}{}\mathrel{}\middle|\mathrel{}L_{X},z,M_{500}\right)$
(10)

We take redshift measurements to be unbiased, with uniform Gaussian scatter at
fixed $z$, that is supposed equal to $0.01(1+z)$:

$p\left(z_{\rm
eFEDS}\mathrel{}\middle|\mathrel{}z\right)\sim\mathcal{N}\left(z\mathrel{};\mathrel{}0.01(1+z)\right)$
(11)

We suppose a log-normally distributed flux F_300kpc given the values of
$\@vec{\alpha}$, that is:

$p\left(\ln\texttt{F\\_300kpc}{}\mathrel{}\middle|\mathrel{}L_{X},z,M_{500}\right)\sim\mathcal{N}\left(\ln\mu_{\texttt{F\\_300kpc}{}}\mathrel{};\mathrel{}\sigma_{\ln\texttt{F\\_300kpc}{}}\right)$
(12)

In this expression, both $\mu_{\texttt{F\\_300kpc}{}}$ and
$\sigma_{\ln\texttt{F\\_300kpc}{}}$ depend on $(L_{X},z,M_{500})$ as follows.
We first convert any set of values $(L_{X},z)$ into a $0.5-2$ keV flux within
$R_{500}$ by assuming an APEC plasma model with temperature set by a standard
luminosity-temperature relation (Bulbul et al., 2019). Foreground absorption
is taken into account using a position-dependent value of the hydrogen column
density (HI4PI Collaboration et al., 2016). Aperture correction (from
$R_{500}$ to 300 kpc) requires knowledge of an emissivity profile for the ICM.
We take an isothermal, isometallicity gas density profile obtained from local
X-COP clusters (Ghirardini et al., 2019). The aperture correction thus depends
solely on the value of $R_{500}$ (itself deriving from $M_{500}$ and $z$). We
assume unbiased flux measurements and set $\mu_{\texttt{F\\_300kpc}{}}$ to the
value of this aperture-corrected flux. We build a simple model for measurement
errors by setting a power-law model, with an extra term $\epsilon^{i}$:

$\sigma_{\ln\texttt{F\\_300kpc}{}}=0.14\left(\frac{\mu_{\texttt{F\\_300kpc}{}}}{10^{-13}\,\mathrm{erg\,s}^{-1}\,\mathrm{cm}^{-2}}\right)^{-0.47}+\epsilon^{i}$
(13)

Numerical coefficients in this expression were obtained by fitting a power-law
model to the flux uncertainties reported in Liu et al. (2022). In order to
account for cluster-to-cluster variability of the error, the extra term
$\epsilon^{i}$ is specific to each cluster (labelled by $i$). As an
assumption, $\epsilon^{i}$ is the deviation of the reported error to the value
predicted by the empirical power-law model. The eFEDS cluster catalog contains
some clusters for which only an upper limit is reported on the measured 300
kpc flux. For these objects, we still assume a log-normal model. Its central
value $\mu_{\texttt{F\\_300kpc}{}}$ is taken to be half of the upper limit
value; the dispersion is such that $\epsilon^{i}$ is the deviation of the
reported upper limit to the fixed value of $2.5\times 10^{-14}$ erg s-1 cm-2
(median of all upper limits in the catalog). Although imperfect, the
generative model depicted here should capture the main trends underlying the
flux measurement process. It is aimed at predicting what would be the reported
flux F_300kpc and redshift $z_{\rm eFEDS}$ of an eFEDS-detected cluster, for
any set of (true) values $\@vec{\alpha}=L_{X},z,M_{500}$.

The reference survey selection function
$p\left(I_{\rm{R}}^{i}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{{\rm
sel},R}\right)$ is the eFEDS selection function
$p\left(I_{\rm{eFEDS}}\mathrel{}\middle|\mathrel{}L_{X},z,EM_{0},T_{\rm
exp}\right)$. It is constructed as described in Liu et al. (2022), with
addition of the morphological parameter $EM_{0}$ characterizing the central
emissivity of the clusters (App. E). The exposure time $T_{\rm exp}$ varies
from one cluster to the other and its value is read from the eFEDS exposure
map (Brunner et al., 2022).

The eRASS1 selection
$p\left(I_{\rm{cosmo}}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{{\rm
sel},T}\right)$ is the main model ingredient we wish to validate in this work.
We present results for the model used in the eRASS1 cosmological analysis
(Ghirardini et al., submitted, see also Fig. 10), involving count-rate $CR$,
redshift $z$ as well as position-dependent, cluster-independent, selection
features $\mathcal{H}^{i}$. These features are galactic absorption column
density, local background brightness and exposure time and their values are
read from static sky maps at each cluster position. The eRASS1 cosmology
sample is restricted to clusters with eRASS1-measured redshift
$0.1<z_{\lambda}<0.8$. We assume Gaussian-distributed errors for the eRASS1
redshifts with standard deviation $\sigma_{z}/(1+z)=0.015$, see also Eq. 14.
The criterion IN_ZVLIM=True, recommended by Kluge et al. (submitted) is
present in our model through the implementation of a sky mask.

The last ingredient required in the model is the distribution
$p\left(\@vec{\Theta}_{{\rm sel},T},\@vec{\Theta}_{{\rm
sel},R},\@vec{\alpha}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm
model}\right)$. Given the assumptions enumerated in this section, it reduces
to $p\left(L_{X}\mathrel{}\middle|\mathrel{}M_{500},z\right)$. We take the
scaling relation of Bulbul et al. (2019) relating the core-included $0.5-2$
keV luminosity within $R_{500}$ to the $M_{500}$ mass and redshift $z$. It
follows a log-normal distribution with scatter 0.27 in the luminosity
direction. In practice, this operation implies computing a probabilistic
weight on a three-dimensional grid representing the values of mass, redshift
and luminosity.

Putting together the model ingredients, we are able to predict for each
cluster in eFEDS, associated to a given flux and redshift measurement, what
would be its probability to appear in the eRASS1 cosmological catalog of
clusters. By performing a 2 arcmin match between the eFEDS and eRASS1 catalog,
we find 35 actual clusters in common between both samples (see also Bulbul et
al., in press), and 90% of the matches are separated by less than $30\arcsec$.

Fig. 15 summarizes the outcome of this exercise. As expected, low-flux
clusters are associated to a low $P_{i}$ value. Reassuringly, none of the
‘bright’ (rigorously, the high-$P_{i}$) eFEDS clusters is missed in the eRASS1
cosmology sample. The model predicts $19.5\pm 3.2$ eFEDS clusters to be in the
eRASS1 cosmology sample, while the actual number is 35 (that is, a 5-$\sigma$
underestimate). We find the empirical parameters $\delta_{1}$ and $\delta_{2}$
have posterior distributions compatible with zero . However, the $P_{i}$
values appear systematically lower than the actual fraction of matches in each
bin of $P_{i}$ values (Fig. 16). Overall these results point to a slight
mismatch between the predicted detection rates and the actual number of
matches. We will discuss this issue later in the paper.

Figure 15: Testing the models with eRASS1-cosmo as a test catalog, eFEDS as a
reference catalog. Each dot corresponds to one of the 531 eFEDS clusters in
the plane of measured 300 kpc flux (F_300kpc in units erg s-1 cm-2) and
measured redshift ($z_{\rm eFEDS}$). eFEDS clusters with only flux upper
limits are placed at the bottom with star symbols. Black circles indicate
clusters also found in the eRASS1 cosmology sample. Colour reflects the model
probability $P_{i}$ (Eq. 8) computed for each eFEDS cluster. Figure 16:
Testing the models with eRASS1-cosmo as test catalog, eFEDS as reference
catalog. The 531 eFEDS entries are grouped in bins of the model output
probabilities $P_{i}$ (horizontal axis). In each bin, the actual fraction of
matches in the eRASS1 catalog is reported on the vertical axis. Error bars
represent the estimated 68% confidence range (App. C). The dotted line
indicates the one-to-one relation, that would follow a perfectly reliable
model. In general, the model predicts slightly lower probabilities than
actually observed.

###  6.3 Testing the eFEDS cluster sample against eRASS1 cosmology sample

By exchanging the role of eFEDS and eRASS1 in the test presented above, we can
use the eRASS1 cosmology catalog as a reference sample and test the eFEDS
selection model. The interest of this operation is limited to sanity checks,
because eRASS1 is 15 times shallower than eFEDS in their overlapping area. We
consider only the 85 eRASS1 clusters extracted from the cosmology sample
located in the fiducial eFEDS area bounded in Right Ascension by 8 h and 10 h
and in Declination by $-5$ and $8\deg$.

For describing each of the eRASS1 detections, we construct a model for the
flux and redshift measurements of each entry. We assume independence of the
estimated redshift $z_{\lambda}$ (column Z_LAMBDA in the eRASS1 cosmology
catalog) and the estimated $R_{500}$ flux in the band 0.5-2 keV (F500_0520,
units erg s-1 cm-2) and we form the model:

$p\left(z_{\lambda}\mathrel{}\middle|\mathrel{}z\right)\sim\mathcal{N}\left(z\mathrel{};\mathrel{}0.015(1+z)\right)$
(14)

This model expresses normally distributed redshift uncertainties with error
increasing with redshift. For flux measurements, we assume a log-normal
distribution:

$p\left(\ln\texttt{F500\\_0520}{}\mathrel{}\middle|\mathrel{}L_{X},z\right)\sim\mathcal{N}\left(\ln\mu_{\texttt{F500\\_0520}{}}\mathrel{};\mathrel{}\sigma_{\ln\texttt{F500\\_0520}{}}\right)$
(15)

In this equation, $\mu_{\texttt{F500\\_0520}{}}$ and
$\sigma_{\texttt{F500\\_0520}{}}$ depend on $L_{X}$ the true luminosity within
$R_{500}$ and on the true redshift $z$. Assuming unbiased flux measurements,
$\mu_{\texttt{F500\\_0520}{}}$ is obtained by converting the cluster rest-
frame $0.5-2$ keV luminosity into the observer-frame $0.5-2$ keV flux,
assuming an isothermal ICM and a position-dependent hydrogen column density
taken from the HI4PI survey (HI4PI Collaboration et al., 2016). As for the
uncertainty model, we proceed by fitting a simple model linking the catalog
flux and errors reported in Bulbul et al. (in press). We take:

$\sigma_{\ln\texttt{F500\\_0520}{}}=0.48\left(\frac{\mu_{\texttt{F500\\_0520}{}}}{10^{-13}\,\mathrm{erg\,s}^{-1}\,\mathrm{cm}^{-2}}\right)^{-0.51}+\epsilon^{i}.$
(16)

The quantity $\epsilon_{i}$ is specific to each cluster in the eRASS1 sample
and equals the deviation of the catalog uncertainty relative to the power-law
model.

Fig. 17 illustrates the outcome of this comparison exercise. The model
predicts $35.7\pm 1.3$ eRASS1 clusters to be found in eFEDS. The actual number
of matches is 36444Border effects due to more or less restrictive sky masking
lead to a different number of matches than mentioned in the previous section
(35 matches). This effect is explicitly taken into account in our model.. Note
that a number of eRASS1 clusters fall in the non-exposed eFEDS area, hence
their value of $P_{i}=0$. Most of the clusters covered by eFEDS have very high
probabilities of being detected ($P_{i}\simeq 1$), as a consequence of its
much deeper exposure. Reassuringly the eFEDS survey does not miss any eRASS1
cluster, highlighting in turn the low level of spurious contamination in the
eRASS1 cosmology sample. In particular, it is not considered problematic that
both clusters 1eRASS J094023.3+022824 and 1eRASS J084147.8-031154 have no
corresponding match in eFEDS despite their respective predicted values
$P_{i}=0.79$ and $0.55$. At their location on sky, the eFEDS survey is very
shallow (with 140 s and 10 s effective exposure time, respectively) and the
eFEDS selection model is not well calibrated in this regime.

Figure 17: Figure similar to Fig. 15, now taking eRASS1-cosmo as the reference
sample (85 entries) and eFEDS as the test sample. The vertical axis represents
F500_0520, the $0.5-2$ keV flux within $R_{500}$ reported in the eRASS1
catalog (units $10^{-14}$ erg s-1 cm-2). The horizontal axis is the catalog
redshift $z_{\lambda}$. Clusters falling in zones of zero eFEDS exposure have
vanishing probability $(P_{i}=0)$. Other eRASS1 clusters are very likely to be
detected in eFEDS ($P_{i}\simeq 1$); in fact they are, as shown by the black
circles. The isolated point at the bottom of the figure only has a flux upper
limit $1.5\times 10^{-13}$ erg s-1 cm-2.

The model explains correctly the outcome of matching eFEDS to eRASS1 cosmology
catalog. There is no evidence for a bias in the value of the probabilities
$P_{i}$; most of them are close to zero or one, in agreement with the cluster
being absent or present in the eFEDS catalog (respectively). The values of
$\delta_{1}$ and $\delta_{2}$ are themselves compatible with zero, indicating
again consistency between the model and the present catalogs being tested.

###  6.4 Testing the eRASS1-main sample against eFEDS

We now turn to a test using the extended version of the eRASS1 catalog,
comprising all X-ray clusters with $\mathcal{L}_{\rm ext}>3$ (Bulbul et al.,
in press). In contrast to the cosmology sample this catalog does not apply any
redshift subselection. The model we construct is very similar to that
presented in Sect. 6.2. In particular the distributions of eFEDS flux and
redshift measurements are identical to that depicted in Eq. 10 to 13. The main
modification consists in using a model selection function
$p\left(I_{\rm{T}}^{i}=I_{\rm{main}}^{i}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{{\rm
sel},T}\right)$. We take as selection variables the count-rate $CR$, true
redshift $z$ and sky position-dependent parameters $\mathcal{H}$.

Our model predicts $30.5\pm 4.2$ matches, while the actual number is 61 (hence
a larger than 7-$\sigma$ difference). We show one diagnostic plot on Fig. 18,
comparing the agreement between the $P_{i}$ values output of the model and the
actual fraction of matches in the eRASS1 catalog. This validation experiment
using the eRASS1-main sample highlights a higher discrepancy than that
obtained in the experiment described in Sect. 6.2, using the eRASS1 cosmology
catalog. The model predicts less clusters to be found in eRASS1 than actually
observed, due to underestimated $P_{i}$ values. This is confirmed with the
posterior distribution of $(\delta_{1},\delta_{2})$ shown in Fig. 19,
significantly away from the zero-point. Seven objects are qualified as
‘surprises’ ($P_{i}<0.025$ and matched). According to our model, they are too
faint to appear in the eRASS1 sample. For two cases (1eRASS J093024.6-020635
and 1eRASS J084558.5+031340), deeper eFEDS data reveals point-sources near the
cluster emission, that were correctly excised to provide a eFEDS flux
measurement, however in the eRASS1 data they are undistinguishable from the
cluster. We do not find any cluster missed by eRASS1 according to this model.
Overall, this result indicates our selection function for the primary cluster
sample is less well understood than for the cosmology sample. Alternatively,
our generative model for eFEDS flux F_300kpc may fail to exactly reproduce the
measurement process described in the eFEDS cluster catalog paper (Liu et al.,
2022). We will return to this issue in Sect. 7.

Figure 18: Figure similar as Fig. 16, using the primary sample of clusters
eRASS1-main as a test sample, instead of the eRASS1 cosmology sample.  Figure
19: Posterior distribution for the pair of parameters
$(\delta_{1},\delta_{2})$ introduced in Eq. 9, for the experiment using
eRASS1-main as a test and eFEDS as reference sample. The histograms represent
the marginalized posterior distributions, with mean and 68% range overprinted.
The contours represent the equivalent 0.5-, 1-, 1.5- and 2-$\sigma$
distribution of the joint posterior. Clearly the parameter $\delta_{2}$
departs from zero, indicating that a fraction of the systems we predict not be
detected leak into the actual eRASS1 sample.

###  6.5 Testing the eRASS1 cosmological cluster sample against SPTpol-100d

Huang et al. (2020) present a catalog of 89 galaxy clusters discovered via the
Sunyaev-Zeldovich (SZ) effect in the SPTpol-100d field, located at R.A. =
23h30m and Dec = $-55\deg$. In the area bounded by 23h ¡ R.A. ¡ 0h and
$-60\deg$ ¡ Dec ¡ $-55\deg$, there are 65 SPTpol clusters at $z>0.25$. This
lower redshift limit is chosen to account for growing incompleteness of the
SPT catalog in the low-$z$ regime. In all that follows we will force the
SPTpol-100d selection function to take value zero at redshifts below $z=0.25$.

The SPTpol-100d selection function depends on one single parameter $\zeta_{\rm
SPTpol}$. We write:

$p\left(I_{\rm{SPTpol}}\mathrel{}\middle|\mathrel{}\zeta_{\rm
SPTpol}\right)=\int^{u(\zeta)}_{-\infty}\left(2\pi\right)^{-1/2}e^{-x^{2}/2}\mathrm{d}x$
(17)

The upper integration bound is a function of $\zeta_{\rm SPTpol}$, namely
$u(\zeta_{\rm SPTpol})=\left[\left(3+\zeta_{\rm
SPTpol}^{2}\right)^{1/2}-4.5\right]$. This expression reflects the SPTpol-100d
selection being a mere thresholding on the signal-to-noise parameter
$\xi>4.5$. The parameter $\xi$ is biased and normally distributed with unit
standard deviation (e.g. Bleem et al., 2020; Huang et al., 2020):

$p\left(\xi\mathrel{}\middle|\mathrel{}\zeta_{\rm
SPTpol}\right)\sim\mathcal{N}\left(\sqrt{\zeta_{\rm
SPTpol}^{2}+3}\mathrel{};\mathrel{}1\right).$ (18)

The measurement vector $\@vec{\Theta}_{\rm meas}$ is taken from the
SPTpol-100d catalog, namely the xi and redshift columns, standing respectively
for the measured, uncorrected signal-to-noise $\xi$ and the measured redshift
$z_{\rm SPT}$. We again assume approximate independence of both measurements
and we take as model:

$p\left(z_{\rm
SPT}\mathrel{}\middle|\mathrel{}z\right)\sim\mathcal{N}\left(z\mathrel{};\mathrel{}0.02(1+z)\right)$

This expression broadly reflects the redshift measurement uncertainty. We use
again Eq. 18 to generate the distribution for $\xi$, now seen as a measurement
feature.

Involving the selection model for the eRASS1 cosmology sample (Sect. 6.2)
requires incorporation of the $0.5-2$ keV luminosity $L_{X}$ as a latent
variable in the model. Consequently we need to construct a joint model for the
distribution of luminosities and $\zeta_{\rm SPTpol}$. We combine the
$L_{X}-M_{500}-z$ scaling relation mentioned previously (Bulbul et al., 2019)
and the $\zeta-M_{500}-z$ relation555We multiply their value of $A_{SZ}$ by
$\gamma_{\rm field}=2.66$, as prescribed by the authors. We also account for a
misprinted $h$ in their Eq. 10 (S. Bocquet, priv. comm.) of Huang et al.
(2020) into a bivariate log-normal distribution with covariance matrix:

$\mathrm{Cov}\left(\ln L_{X},\ln\zeta_{\rm
SPTpol}\right)=\pmatrix{\sigma}_{\ln L_{X}}^{2}&\rho\sigma_{\ln
L_{X}}\sigma_{\ln\zeta}\\\ \rho\sigma_{\ln
L_{X}}\sigma_{\ln\zeta}\sigma_{\ln\zeta}^{2}$

Here, $\sigma_{\ln L_{X}}$ and $\sigma_{\ln\zeta}$ are the log-normal scatters
of both distributions at fixed mass $M_{500}$ and redshift $z$ (with values
respectively 0.27 and 0.18). The value of $\rho$ ranges between $-1$ and $1$,
it reflects correlated scatter between the SZ signal-to-noise and X-ray
luminosity at fixed mass and redshift. We take a fiducial value $\rho=0.2$, we
will investigate the impact of varying this parameter in Sect. 7.

The model predicts $11.4\pm 2.4$ matches, while the actual number of matches
is 8. Fig. 20 shows the distribution of values $P_{i}$ output of the model, in
the plane of the measured SPT mass and redshift. A small number of SPTpol-100d
clusters have their $P_{i}$ high enough to be likely detectable in eRASS1. In
fact, the matched instances are those with high $P_{i}$ values in general.

Figure 20: Figure similar as Fig. 15, using the cosmology sample of eRASS1
clusters as a test sample, and the SPTpol-100d catalog as a reference sample
(65 objects at $z>0.25$, coloured dots). The vertical axis represents the
column xi in the SPTpol catalog, standing for the measured cluster signal-to-
noise (unitless). The horizontal axis represents the measured redshift in the
SPT catalog. The model predicts only a handful of eRASS1 detections among
those objects, consistently with the observed number of 8 matches.

Fig. 21 demonstrates the fair agreement between the 65 $P_{i}$ values outcome
of the model and the actual fraction of matches (in bins of $P_{i}$ values).
From these diagnostics we conclude there is no evidence for inconsistencies
either in the catalogs nor in the model. In particular, the selection function
models account correctly for the observed matches between SPTpol-100d and
eRASS1-cosmo.

Figure 21: Figure similar as Fig. 16, using the cosmology sample of eRASS1
clusters as a test sample, and the SPTpol-100d cluster sample as reference.
Within the large uncertainties, there is good agreement between the predicted
eRASS1 detection probability and the actual fraction of matched systems.

###  6.6 Testing the SPTpol-100d cluster sample against eRASS1-cosmo

The previous test is reverted by exchanging the roles of the SPTpol-100d
catalog with that of the eRASS1 cosmology sample. We thus test the SPTpol
catalog, taking eRASS1 as reference. There are 19 eRASS1 clusters in the
common footprint between both surveys. Construction of the model involves both
selection functions and the same joint probability distribution
$p\left(L_{X},\zeta_{\rm SPTpol}\mathrel{}\middle|\mathrel{}M_{500},z\right)$
as described in Sect. 6.5. We also use a similar generative model as in Sect.
6.3 for the flux F500_0520 and redshift $z_{\lambda}$ listed in the eRASS1
catalog. The 19 output probabilities $P_{i}$ returned by the model are shown
in Fig. 22. The model predicts $7.0\pm 1.3$ matches, while the actual number
is 8, indicating consistency between the results and the predictions.

Figure 22: Figure similar as Fig. 17, using the cosmology sample of eRASS1
clusters as a reference sample and the SPTpol-100d cluster sample as a test
sample. Two eRASS1 clusters with high-$P_{i}$ values and absent from the
SPTpol-100d catalog (hollow magenta dots) are discussed in the text. The star
symbol at the bottom of the figure corresponds to a eRASS1 entry with upper
flux limit $7\times 10^{-14}$ erg s-1 cm-2.

Fig. 23 demonstrates adequacy between the 19 values of $P_{i}$ predicted by
the model and the actual fraction of matches in each bin. Low number
statistics prevents from gaining more insight into the model, calling for
further experiments with larger samples of clusters.

Figure 23: Figure similar as Fig. 16, using the cosmology sample of eRASS1
clusters as a reference sample and the SPTpol-100d cluster sample as a test
sample. The low number statistics (only 19 systems) prevents from finer
binning in $P_{i}$ values (horizontal axis). Nevertheless, there is
satisfactory agreement with the dotted (one-to-one) line.

As apparent from Fig. 22, one cluster is not found in the SPTpol-100d survey,
although our model predicts $P_{i}=0.90$. The eRASS1 cluster 1eRASS
J230029.2-510022 is measured at $z_{\lambda}=0.37$ with
$\texttt{F500\\_0520}=3.4\times 10^{-13}$ erg s-1 cm-2 and should be massive
enough to appear in the SPTpol survey. It is missed by SPTpol, very likely due
to its location close to the very edge of the surveyed field (at 23h Right
Ascension).

Similarly, the eRASS1 cluster 1eRASS J231306.7-550417 measured with
$z_{\lambda}=0.31$ with $\texttt{F500\\_0520}=3.4\times 10^{-13}$ erg s-1 cm-2
is not detected by SPTpol despite the model predicted $P_{i}=0.72$.

###  6.7 Testing the eRASS1 cosmological cluster sample against SPTpol-ECS

The model described for the SPTpol-100d catalog is straightforwardly adapted
to the SPTpol-ECS catalog presented in Bleem et al. (2020). The ECS sample is
collected over $2770\deg^{2}$ area in two regions of the Southern Hemisphere.
A total of 163 clusters with significance $\xi>5$ and located at redshift
$z>0.25$ fall within the footprint of the eRASS1 cosmology sample. It is
shallower than the SPTpol-100d survey, hence we expect a majority of the
clusters to be found in the eRASS1 sample, up to $z=0.8$.

Our model follows exactly the same steps as in Sect. 6.5, with following
modifications:

  * •

the measurement features $\xi$ and $z_{\rm SPT}$ are taken from columns XI and
REDSHIFT of the SPT-ECS catalog;

  * •

the unbiased significance parameter ruling SZ detections is denoted
$\zeta_{\rm ECS}$. Its logarithm scales as a power-law with mass, as described
in Eq. 4 of Bleem et al. (2020). In particular we multiply the normalisation
$A_{SZ}$ of this scaling relation by a factor $\gamma_{\rm field}$ specific to
each sub-field of the SPT-ECS survey. An additional multiplicative factor
$\gamma_{\rm ECS}=1.124$ rescales all values of $A_{SZ}$, thus reflecting
calibration updates (Table 1 and Sect. 5 in Bleem et al., 2020);

  * •

we only select entries with $\xi$ greater than 5.

In the following we keep a fixed value $\rho=0.2$ for the correlated scatter
between $\ln L_{X}$ and $\ln\zeta_{\rm ECS}$. There are 103 matching instances
between both catalogs, based on a $2\arcmin$ positional match. Our model
predicts $103\pm 4$ matches, in excellent agreement with the observed number.
Fig. 24 shows the distribution of the 163 entries in the (measured)
significance and redshift plane. Most of the SPT-ECS clusters at $z<0.8$ are
detected by eRASS1 (black circles). There is no obvious missing cluster among
the 163 SPT detections. However we identified one ‘surprising’ detection, SPT-
CLJ0333-3707 (black ‘X’ in the figure, $P_{i}\simeq 0$). It is listed in the
ECS catalog with significance 5.1 and photometric redshift $z_{\rm
SPT}=1.05\pm 0.04$. The matching instance in eRASS1 is 1eRASS
J033323.9-370744, listed at $z_{\lambda}=0.52\pm 0.01$. This cluster is one of
the 3% ambiguous cases discussed in Kluge et al. (submitted), where two
unrelated clusters overlap along the same line of sight.

Figure 24: Figure similar as Fig. 20, using 163 clusters in SPTpol-ECS catalog
as a reference sample and the eRASS1 cosmology sample as test.

Fig. 25 demonstrates the reliability of the $P_{i}$ values as predicted by the
model. We conclude to satisfactory agreement between the model and the
observed matches and non-matches in the two catalogs.

Figure 25: Figure similar as Fig. 21, obtained by using the SPT-ECS sample as
reference (163 entries) and using the cosmology sample of eRASS1 clusters as
test.

###  6.8 Testing the SPTpol-ECS sample against eRASS1 cosmology sample

We conclude this series of tests by reversing the experiment of Sect. 6.7.
There are 958 eRASS1 clusters in the cosmology catalog falling within the sky
footprint of SPTpol-ECS survey. The model takes as measurement features
F500_0520 and $z_{\lambda}$ and it assumes identical generative distributions
as in Sect. 6.3.

Fig. 26 shows the outcome of this calculation. Our model predicts $65.4\pm
5.8$ matches, significantly less than the observed 101 matches observed in the
common footprint. In particular we find three unexpected SPT-ECS detections
($P_{i}<0.025$). Secondary components along the line of sight of these three
clusters may be responsible for the discrepant mass estimates, either by
boosting the SZ significance or by decreasing the measured X-ray flux. They
are illustrated in Fig. 27:

  * •

1eRASS J021731.8-320002 is listed at $z_{\lambda}=0.352\pm 0.009$ with flux
$2.1\pm 0.4\times 10^{-14}$ erg s-1 cm-2 and extent likelihood 14.4. The
estimated X-ray mass is $M_{500,X}\simeq 3\times 10^{14}\,M_{\odot}$. The
matched SPT-ECS detection is at $z_{\rm SPT}=0.354\pm 0.009$ (from
Magellan/FourStar) and $\xi=5.1$. Its estimated SZ mass is $M_{500,{\rm
SPT}}\simeq 4\times 10^{14}\,M_{\odot}$.

  * •

1eRASS J034733.3-333351 is listed at $z_{\lambda}=0.465\pm 0.009$ with flux
$8_{-3}^{+10}\times 10^{-14}$ erg s-1 cm-2 and extent likelihood 7.5. The
estimated X-ray mass is $M_{500,X}\simeq 2\times 10^{14}$ M⊙. The matched SPT-
ECS detection is at $z_{\rm SPT}=0.45\pm 0.01$ (photometric, from
Magellan/FourStar) and $\xi=5.4$. Its estimated mass is $M_{500,{\rm
SPT}}\simeq 4\times 10^{14}$ M⊙.

  * •

1RASS J045420.6-373708 is listed at $z_{\lambda}=0.513\pm 0.008$ without any
flux measurement, nor an upper limit and extent likelihood 6.3. The
corresponding SPT detection is at $z_{\rm SPT}=0.509\pm 0.006$ (photometric)
and significance $\xi=5.9$, with estimated mass $M_{500,{\rm SPT}}\simeq
4.8\times 10^{14}$ M⊙.

Figure 26: Figure similar as Fig. 22, using 958 eRASS1 clusters from the
cosmology sample as reference, and SPTpol-ECS catalog as a test sample. Star
symbols at the bottom of the figure indicate eRASS1 objects with no flux
measurement, possibly with an upper flux limit.
---
Figure 27: Optical cutouts (DR10 Legacy survey) and eRASS1 X-ray contours
overlaid. These eRASS1 clusters are present in the cosmology sample and are
unexpectedly matched to SZ detected clusters in the SPT-ECS sample. Our model
indeed predicts a probability value $P_{i}=p\left(I_{\rm{\rm SPT-
ECS}}\mathrel{}\middle|\mathrel{}I_{\rm{\rm
cosmo}},\texttt{F500\\_0520},z_{\lambda}\right)<0.025$ for each of these
entries.

Fig. 28 assesses the reliability of the $P_{i}$ values as returned by the
model. The predicted probabilities are generally lower than the actual
fraction of matches. Several causes for this shift may include: our fiducial
$L_{X}-M_{500}-z$ relation (here taken from Bulbul et al., 2019) predicting
too high values of X-ray luminosity at fixed halo mass; or an underestimate of
eRASS1 flux measurements; or an overestimate of some SPT significance values
$\xi$ due, e.g. to Malmquist bias.

Figure 28: Figure similar as Fig. 23, obtained by using the eRASS1 cosmology
sample (958 entries) as reference and using the SPTpol-ECS as test sample.

##  7 Discussion

###  7.1 Linking our models with $\mathcal{L}_{\rm ext}$ and $\mathcal{L}_{\rm
det}$

In processing real and simulated data, eSASS assigns several measured
quantities to each detected source. Among the most relevant for understanding
selection effects are detection ($\mathcal{L}_{\rm det}$) and extent
($\mathcal{L}_{\rm ext}$) likelihoods. Applying thresholds on these quantities
leads to defining a new sample, characterized by a balance between
contamination and completeness (e.g. Seppi et al., 2022). In the specific case
of the eRASS1 cluster catalog, a source is detected if it shows
$\mathcal{L}_{\rm det}$ greater than 5. The threshold on $\mathcal{L}_{\rm
ext}$ for selecting this source as a cluster is either 3 or 6, the latter
being chosen for the cosmology sample.

Such parameters are absent from all models discussed in this paper, neither as
an input nor as an output of their predictions. It may seem at odds with
studies (e.g. Rix et al., 2021) suggesting to model selection function based
on catalog parameters cut only (here collectively designed under a generic
name, $\xi_{X}$). In fact, let us acknowledge that:

$p\left(I\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm
model}\right)=p\left(\xi_{X}>\xi_{X}^{\rm
min}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm
model}\right)=\int_{\xi_{X}^{\rm
min}}^{+\infty}\mathrm{d}\xi_{X}p\left(\xi_{X}\mathrel{}\middle|\mathrel{}\@vec{\Theta}_{\rm
model}\right)$

Building a selection model then translates into constructing a distribution of
$\xi_{X}$ given a set of model parameters. This may involve auxiliary
variables (such as observed flux, see e.g. Grandis et al., 2020) and a
combination of empirical scaling laws between $\xi_{X}$ and those variables,
involving some parameterized form of scatter. An appealing method consists in
letting all or a subset of these parameters free in the scientific analysis,
together with the cosmological parameters (e.g. Vanderlinde et al., 2010).
These scaling relation parameters would then be ‘learnt’ directly from data
themselves; they would self-adjust to catalog parameters and contribute in
increasing the global cosmology likelihood, via calibration through weak
gravitational lensing or other mass proxies. In this approach, the exact shape
of distributions (e.g. log-normal or Gaussian scatter) is a strong prior
assumption. In absence of other information (e.g. external or simulated
catalogs), we extrapolate our knowledge on brighter systems to fainter ones,
potentially leading to interpretation errors close to the detection limit
(Gallagher & Smeenk, 2023). We may even imagine extreme situations in which
these assumptions are not verifiable, e.g. if the unimodal scatter in scaling
relations does not account for a separate class of faint, flat surface-
brightness systems that escapes detection in our catalog.

A safer alternative may consist in constructing a model $\xi_{X}$ directly
from $\mathcal{L}_{\rm ext}$ and $\mathcal{L}_{\rm det}$ as recorded in
simulated datasets. There are also difficulties associated with this approach.
Most notably, only detected (resp. selected) objects have an associated value
of $\mathcal{L}_{\rm det}$ (resp. $\mathcal{L}_{\rm ext}$). A replacement
value needs to be assigned to those simulated systems missing a detection. One
may also attempt to run the source detection algorithm with lower thresholds
(detlikemin and extlikemin in eSASS) in order to assign a value to fainter
objects; however this has a number of non-trivial side-effects from the
algorithmic point of view, such as masking, source splitting, background
estimates, etc. As a consequence, there is no guarantee that the samples
selected with $\mathcal{L}_{\rm det}>5$ in both runs (the one with
detlikemin=5 and the other with detlikemin¡5) are truly identical.

Interestingly, the logistic regression formalism presented in Sect. 4 uses a
real-valued function $f$, whose expression is given in Eq. 5. In our analysis
$f$ is a latent variable, without predefined physical meaning. The selection
probability monotonically increases with $f$, through the sigmoid function of
Eq. 4. As a useful cross-check, we show the correlation between $f$ and the
value of $\mathcal{L}_{\rm det}$ and $\mathcal{L}_{\rm ext}$ for the logistic
regression involving 35 features (Fig. 29). Reminding that only binary flags
where exposed to the model in the training phase, it may be reassuring that
the probability output is driven by a quantity behaving similarly as detection
and extent likelihood.

|
---|---
Figure 29: Relation between the latent variable
$f(\eta)=w_{0}+\sum_{j}w_{j}\eta_{j}$ calculated for each simulated cluster in
the test sample after performing logistic regression with 35 surface
brightness features, and the values of $\mathcal{L}_{\rm det}$ and
$\mathcal{L}_{\rm ext}$ as provided by eSASS for those matched to a detection.
Each dot stands for a simulated cluster, density contours ease visualisation
in crowded regions. Undetected clusters have no associated $\mathcal{L}_{\rm
det}$ value (red histogram in the lower panel). Clusters detected but not
categorised as extended have zero value for $\mathcal{L}_{\rm ext}$ (orange
histogram). The apparent correlation between $f$ and both measurements
indicate that the model has ‘discovered’ the importance of $\mathcal{L}_{\rm
det}$ and $\mathcal{L}_{\rm ext}$ in the selection process.

We find a similar result with the mean of the latent variable $\@vec{f}$
governing the predictions output of the Gaussian Process classifier (see App.
D). Fig. 30 in particular shows the a posteriori correlation between $\langle
f\rangle$ and the detection and extent likelihood obtained in a test sample of
simulated clusters.

|
---|---
Figure 30: Similar as Fig. 29, replacing the x-axis by the mean of the latent
variable implicit in the Gaussian Process classifier $p\left(I_{\rm{\rm
cosmo}}\mathrel{}\middle|\mathrel{}CR,z,\mathcal{H}\right)$.

###  7.2 Underestimate of the eRASS1 selection model?

Our external validation results in Sect. 6 suggest an inconsistency between
the predicted number of eRASS1-detected objects among the eFEDS sample and the
actual number of matching entries. The difference is about 5-$\sigma$ when
considering the eRASS1 cosmology sample and above 7-$\sigma$ with the primary
eRASS1 cluster sample. Our validation against the much brighter SPT-selected
samples does not reveal such inconsistency. This suggests that our model
somehow fails at describing the fainter population of eFEDS clusters.

####  7.2.1 Uncertainties in the generative model

We may first incriminate our simple generative model for eFEDS-measured
fluxes. Our model for F_300kpc (Eq. 12 and 13) assumes unbiased flux
measurements regardless of the brightness of the source. If the reported eFEDS
fluxes (Liu et al., 2022) of faint systems were known to be biased low for
some reason, our models would predict a higher number of matched entries in
the eRASS1 sample and we would conclude to no discrepancy. On the one hand,
eFEDS luminosity measurements appear to agree with _XMM-Newton_ measurements
(e.g. Turner et al., 2022), as also shown by the luminosity functions in Liu
et al. (2022). On the other hand, we have found a 15% offset between eRASS-
measured fluxes and _Chandra_ measurements, the latter showing higher value,
likely attributed to calibration issues (Bulbul et al., in press).

We obtained very similar predictions, and thus similarly discrepant results,
when modelling the measured $R_{500}$ X-ray luminosity of eFEDS clusters
(column Lx500 in the catalog of Bahar et al., 2022) in lieu of the measured
flux. This operation somewhat simplifies our generative model, and it
suppresses any potential pitfall in deriving aperture corrections and
spectrally-dependent K-corrections. There we assumed unbiased, normally
distributed luminosity measurements, with an empirical uncertainty model
following:

$\sigma_{\texttt{Lx500}}=0.2\times 10^{43}\,{\rm
erg\,s}^{-1}\left(\frac{L_{X}}{10^{43}\,{\rm
erg\,s}^{-1}}\right)^{0.9}+\epsilon^{i}$ (19)

In this expression $\epsilon^{i}$ is a constant for each eFEDS cluster indexed
by $i$. We obtained $\epsilon^{i}$ by subtracting the value predicted by the
power-law model (without $\epsilon^{i}$) in the above equation to the
uncertainty reported in the catalog of Bahar et al. (2022).

####  7.2.2 Uncertainties in trained models

Our selection models inherently embed some uncertainty, due to the limited
size of the training sample and to the complex algorithm leading to optimize
many hundreds of Gaussian Process hyperparameters at once (App. D). We have
taken benefit of our multiple Poisson realizations of the eRASS1 twin
simulations to produce seven selection models, each relying on a distinct
training set as large as eight times the actual eRASS1 sky. We folded these
selection functions through our eFEDS/eRASS1 comparison algorithm. In all
cases we retrieved the under-prediction of eRASS1 clusters, the discrepancy
varying between 5 and 7-$\sigma$ among the seven models. The variability and
uncertainty attached to the training phase are thus too small to explain the
observed discrepancy.

####  7.2.3 Uncertainties in the training sample

We now discuss whether our eRASS1 training set may be biased in some way,
which would offset the models from the truth. Such offsets would not be
visible in our internal validation tests (Sect. 5.4) since the test sample is
constructed in exactly the same way as the training sample.

Our baseline procedure associating simulated clusters to detected sources
leads to fewer selected halos than with the positional matching procedure (see
Sect. 3.3). It is especially visible at cluster fluxes around $10^{-14}$ erg
s-1 cm-2, where the number of matched sources can differ up to a factor two
(Fig. 1). We have trained a selection function model with a new set of
detection flags constructed with the alternative matching technique. This
model indeed predicts slightly higher detection probabilities at low flux, in
accordance with the update in the training set. However, the $5-7$-$\sigma$
discrepancy remains in the comparison with the eRASS1 clusters detected in
eFEDS. This is not surprising, since the eRASS1 flux limit in the eFEDS field
is around $10^{-13}$ erg s-1 cm-2, well above the regime where the two
matching methods depart from each other (Fig. 1).

After completion of this work, we found minor differences in the eSASS
software configuration parameters used in processing the simulations (version
eSASSusers_201009, Seppi et al., 2022) and in processing the actual eRASS1
data (version 010, Merloni et al., 2024). These changes mainly impact the
background map creation and source splitting into multiple entries. In
particular, in real data the ermldet algorithm would attempt to split
detections into two sources only if the number of counts is greater than 25;
while in our simulations no specific threshold has been set on this value. We
have tried and processed a subset of the twin simulation with the updated
algorithm. We find a tentative higher detection probability than in the
original mock catalog, which may explain part of the deficit in the predicted
number of eRASS1 sources. However, this finding relies only on a subset of the
entire twin simulations, insufficient to firmly close the issue.

Our final test consists in using a selection model that explicitly involves
the cluster sizes $R_{500}$ in their features (Fig. 35). By using a selection
model in the form $p\left(I_{\rm{\rm
main}}\mathrel{}\middle|\mathrel{}f_{X},R_{500},\mathcal{H}\right)$, we do not
rely any longer on the distribution of sizes imprinted in the simulation.
Rather, we assume our $L_{X}-M_{500}$ is correct and we distribute our
modelled flux within the radius associated to $M_{500}$. Folding this new
model into our eFEDS/eRASS1 analysis, we find better agreement between
predictions and actual number of matches. Hence, using this alternative
selection model instead of our baseline model, the discrepancy reduces from 7-
to 5-$\sigma$ for the eRASS1-main catalog and from 5- to 3-$\sigma$ for the
cosmology catalog. This finding corroborates Fig. 6 in Comparat et al. (2020),
showing that the average $R_{500}$ radius of our simulated clusters at fixed
luminosity is larger than in several published samples666Different
cosmological models between samples and simulations entering this figure may
induce deviations of a few percent in the luminosity and angular distance
values, which is not enough to account for the shift..

####  7.2.4 Sensitivity of the external validation test

In this section we have attempted to address the apparent mismatch between the
number of eRASS1 systems predicted in the eFEDS field and the actual number of
matches. We could not identify one single cause for the discrepancy, rather we
have found several potential difficulties. Overall, our experiments reveal the
high sensitivity of our external validation tests to small variations in the
model assumptions. They also highlight the amount of effort needed to properly
understand selection effects at low fluxes and close to the detection
threshold.

###  7.3 X-ray and SZ correlated detection

In Sect. 6 we have introduced a joint distribution for the X-ray luminosity
and SZ significance $p\left(\ln
L_{X},\ln\zeta\mathrel{}\middle|\mathrel{}M_{500},z\right)$ as a bivariate
normal distribution with correlation coefficient $\rho=0.2$. We now
investigate its impact on the results. We varied $\rho$ between $-0.8$ and
$0.8$, while repeating the experiment with the SPT-ECS sample as reference $R$
and the eRASS1-cosmo sample as test $T$ (Sect. 6.7). We form the quantity:

$B=-2\ln\mathcal{L}=-2\sum_{i\in R}I_{\rm{T}}^{i}\ln
P_{i}+(1-I_{\rm{T}}^{i})\ln(1-P_{i})$ (20)

Its mimimum $B_{\rm min}$ is obtained by varying $\rho$, all other model
parameters held fixed at their (supposedly) true value. We removed the cluster
SPT-CLJ0333-3707 whose redshift measurements disagree with the matched entry
1eRASS J033323.9-370744 (Sect. 6.7). There are 162 entries entering the sum.
Following Cash (1979), $B-B_{\rm min}$ approximately follows a $\chi^{2}$
distribution with one degree of freedom. Fig. 31 shows the curve $B-B_{\rm
min}$ for various values of $\rho$, from which we derive an approximate 68%
confidence interval $-0.4<\rho<0.8$. Consistently with previous studies (e.g.
de Haan et al., 2016; Dietrich et al., 2019) our data does not favour strongly
correlated or anti-correlated scatter between the X-ray and SZ signals.

Figure 31: Value of the combined Bernoulli negative log-likelihood (Eq. 20) in
the comparison of the SPT-ECS sample (reference) and the eRASS1 cosmology
sample (test) for various values of the correlation coefficient $\rho$. This
coefficient expresses the level of covariance between $\ln L_{X}$ and
$\ln\zeta$ at fixed halo mass and redshift. Assuming all other model
parameters are correct, the constraint obtained from the catalog comparison
exercise depicted in this paper is relatively loose: $\rho\simeq 0.2\pm 0.6$
(68% confidence level) as deduced from $(B-B_{\rm min})\sim\chi^{2}_{1}$.

The fiducial value $\rho=0.2$ chosen in Sect. 6 is in fact the best-fit, with
$B_{\rm min}=99.9$. We have checked the goodness of fit by multiple random
Bernoulli draws from the $P_{i}$ values returned by the best-fit model; only
24% of the simulations provide a value of $B$ smaller than $B_{\rm min}$.
Given the many other assumptions made (e.g. fixed scaling relation parameters,
etc.), we cannot firmly exclude specific values of $\rho$ from our dataset. A
comparison involving larger samples may put stronger constraints on $\rho$
through this method.

It is interesting to note that $\rho$ is intimately related to the ICM
distribution for clusters in the reference survey. A physical link should
exist between $L_{X}$ and $\zeta$ that allows computing $\rho$ or at least a
prior distribution, without using the halo mass $M_{500c}$. However, our
implementation of the SPT-SZ selection function through $\zeta$ imposes to
rely on an empirical scaling relation with mass, which prevents from
investigating further in this direction.

##  8 Conclusions

This work presents an investigation of X-ray selection function models aimed
at describing the completeness of the primary and cosmological eRASS1 cluster
samples. Driven by our astrophysical and cosmological objectives, we have
developed a framework inspired by classification problems in statistics and
especially machine learning research. Such empirical models are well-suited to
our purpose, since selection of sources in astronomical catalogs is (at least
in part) resulting from human decisions, whose detailed physical modelling may
be cumbersome. The models we selected are intentionally simple, enabling their
interpretability and explainability.

Our main findings are summarized as follows:

  1. 1.

We have performed extensive twin simulations of the eRASS1_DE sky under
certain cosmological and physical assumptions. These simulations reproduce to
a large extent the chain of selection filters along the path from sources to
final catalogs.

  2. 2.

We have identified which simulated halos are detected and selected by means of
custom matching algorithms. One algorithm relies on tracking the origin of
simulated photons, the other algorithm on Bayesian positional matching. We
have highlighted relevant differences in the outcome of the two algorithms,
mostly below cluster fluxes a few times $10^{-14}$ erg s-1 cm-2.

  3. 3.

By means of a simple logistic (linear) model, we have asked which surface
brightness features govern the detectability and selection of clusters. In
particular, we find the central $20\arcsec$ counts help in detecting a
cluster, but not so in characterizing it as extended. Cluster counts beyond
approximately $1.5\arcmin$ from the cluster centre play a minor role in the
selection of clusters. This model may be seen as a simplified emulator of the
eSASS source detection algorithm.

  4. 4.

Acknowledging that selection function based on surface brightness features
require important modelling effort on the cosmological likelihood side, we
have developed models relying on intermediate variables. We identified a
hierarchy of models, some strongly depending on assumptions made in the twin
simulations, e.g. $p\left(I\mathrel{}\middle|\mathrel{}M_{500},z\right)$,
others enabling less dependence on such assumptions e.g.
$p\left(I\mathrel{}\middle|\mathrel{}CR,z,T_{\rm exp},N_{H},{\rm bkg}\right)$.
Using the latter class of models comes at a cost, namely a heavier modelling
effort in scientific analyses.

  5. 5.

We have tested all our models with left-over simulation samples. In particular
we established a ranking of models based on their precision, their sensitivity
(also known as recall), and the reliability of their predicted probabilities.
Models accounting for at least one morphological feature (e.g. the cluster
$R_{500}$ radius or its central emission measure $EM_{0}$) perform better.

  6. 6.

Using real data, we have tested our models by means of a custom population
model. Comparing to eFEDS or SZ-based surveys, we find no obvious missing
cluster from our eRASS1 samples. Conversely, there is no eRASS1 cluster in the
cosmology sample that is strongly missing in other tested surveys, in
agreement with the known low contamination of the sample (Seppi et al., 2022;
Ghirardini et al., submitted).

  7. 7.

While comparison with the (brighter) SPT-SZ population reveals no specific
issue, we find more eRASS1 clusters in the eFEDS sample than our models would
predict. We have highlighted a series of putative causes for the discrepancy.
We will pursue work to deepen our understanding of the selection of faint
clusters.

In future work we will refine the simulation set with updated cluster models
and halo simulations. An alternative assesment of selection systematics may be
gained through the use of ‘inject-and-retrieve’ methods (e.g. Pacaud et al.,
2006; Suchyta et al., 2016; Kong et al., 2020), where mock cluster images are
distributed in real data. Moreover, ultimate validation of selection models
should rely on a performance metric that relates directly to cosmological
parameters, especially since our purpose is to constrain those parameters from
data. Developing multiple eRASS simulations spanning a range of cosmological
models is a promising route in this direction, although computationally
intensive. Interestingly, forward models of the observed large scale structure
using Effective field theory (EFT) have progressed tremendously (e.g. Jasche &
Lavaux, 2019; Schmidt, 2021a). These models predict the cosmological matter
density field at percent (and even sub percent) accuracy. In the future, these
models will enable tight constraints on cosmological parameters by using not
only linear but also mildly non-linear information contained in the density
field (Schmidt, 2021b; Kostić et al., 2023) and by marginalizing fully over
the selection function. Finally, pursuing the comparison exercise depicted in
this paper with multiple other surveys across a range of wavelengths will
provide interesting constraints on model parameters and reveal any significant
population outliers worth a deep follow-up study.

###### Acknowledgements.

This work is based on data from eROSITA, the soft X-ray instrument aboard SRG,
a joint Russian-German science mission supported by the Russian Space Agency
(Roskosmos), in the interests of the Russian Academy of Sciences represented
by its Space Research Institute (IKI), and the Deutsches Zentrum für Luft- und
Raumfahrt (DLR). The SRG spacecraft was built by Lavochkin Association (NPOL)
and its subcontractors, and is operated by NPOL with support from the Max
Planck Institute for Extraterrestrial Physics (MPE). The development and
construction of the eROSITA X-ray instrument was led by MPE, with
contributions from the Dr. Karl Remeis Observatory Bamberg & ECAP (FAU
Erlangen-Nuernberg), the University of Hamburg Observatory, the Leibniz
Institute for Astrophysics Potsdam (AIP), and the Institute for Astronomy and
Astrophysics of the University of Tübingen, with the support of DLR and the
Max Planck Society. The Argelander Institute for Astronomy of the University
of Bonn and the Ludwig Maximilians Universität Munich also participated in the
science preparation for eROSITA. The eROSITA data shown here were processed
using the eSASS software system developed by the German eROSITA consortium.
N. Clerc acknowledges financial support from CNES. E. Bulbul, A. Liu, V.
Ghirardini, C. Garrel, S. Zelmer, and X. Zhang acknowledge financial support
from the European Research Council (ERC) Consolidator Grant under the European
Union’s Horizon 2020 research and innovation program (grant agreement CoG
DarkQuest No 101002585).
Some of the results in this paper have been derived using the following
packages: scikit-learn (Pedregosa et al., 2011),
MOCpy777https://github.com/cds-astro/mocpy/ (Fernique et al., 2014),
GPy888https://github.com/SheffieldML/GPy (GPy, since 2012), astropy (Astropy
Collaboration et al., 2022), matplotlib (Hunter, 2007), healpy and
HEALPix999http://healpix.sourceforge.net (Górski et al., 2005; Zonca et al.,
2019), HEASoftpy101010https://heasarc.gsfc.nasa.gov/lheasoft/heasoftpy/,
climin111111https://github.com/BRML/climin(Bayer et al., 2015).

## References

  * Agresti & Coull (1998) Agresti, A. & Coull, B. A. 1998, The American Statistician, 52, 119
  * Anderson et al. (2015) Anderson, M. E., Gaspari, M., White, S. D. M., Wang, W., & Dai, X. 2015, MNRAS, 449, 3806
  * Artis et al. (in prep.) Artis, E., Ghirardini, V., Bulbul, E., et al. in prep.
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167
  * Bahar et al. (2022) Bahar, Y. E., Bulbul, E., Clerc, N., et al. 2022, A&A, 661, A7
  * Bayer et al. (2015) Bayer, J., Osendorfer, C., Diot-Girard, S., Rueckstiess, T., & Urban, S. 2015, TUM, Tech. Rep.
  * Bleem et al. (2020) Bleem, L. E., Bocquet, S., Stalder, B., et al. 2020, ApJS, 247, 25
  * Blei et al. (2017) Blei, D. M., Kucukelbir, A., & McAuliffe, J. D. 2017, Journal of the American Statistical Association, 112, 859
  * Böhringer & Chon (2015) Böhringer, H. & Chon, G. 2015, A&A, 574, L8
  * Boubert & Everall (2020) Boubert, D. & Everall, A. 2020, Monthly Notices of the Royal Astronomical Society, 497, 4246
  * Boubert & Everall (2021) Boubert, D. & Everall, A. 2021, Monthly Notices of the Royal Astronomical Society, 510, 4626
  * Brunner et al. (2022) Brunner, H., Liu, T., Lamer, G., et al. 2022, A&A, 661, A1
  * Bulbul et al. (2019) Bulbul, E., Chiu, I. N., Mohr, J. J., et al. 2019, ApJ, 871, 50
  * Bulbul et al. (in press) Bulbul, E., Liu, A., Kluge, M., et al. in press, A&A
  * Bulbul et al. (2022) Bulbul, E., Liu, A., Pasini, T., et al. 2022, A&A, 661, A10
  * Cash (1979) Cash, W. 1979, ApJ, 228, 939
  * Chuang et al. (2019) Chuang, C.-H., Yepes, G., Kitaura, F.-S., et al. 2019, MNRAS, 487, 48
  * Clerc & Finoguenov (2023) Clerc, N. & Finoguenov, A. 2023, in Handbook of X-ray and Gamma-ray Astrophysics. Edited by Cosimo Bambi and Andrea Santangelo, 123
  * Clerc et al. (2018) Clerc, N., Ramos-Ceja, M. E., Ridl, J., et al. 2018, A&A, 617, A92
  * Comparat et al. (2020) Comparat, J., Eckert, D., Finoguenov, A., et al. 2020, The Open Journal of Astrophysics, 3, 13
  * Comparat et al. (2019) Comparat, J., Merloni, A., Salvato, M., et al. 2019, MNRAS, 487, 2005
  * Comparat et al. (2022) Comparat, J., Truong, N., Merloni, A., et al. 2022, A&A, 666, A156
  * Dauser et al. (2019) Dauser, T., Falkner, S., Lorenz, M., et al. 2019, A&A, 630, A66
  * de Haan et al. (2016) de Haan, T., Benson, B. A., Bleem, L. E., et al. 2016, ApJ, 832, 95
  * Debackere et al. (2022) Debackere, S. N. B., Hoekstra, H., Schaye, J., Heitmann, K., & Habib, S. 2022, MNRAS, 515, 3383