© The Authors 2024

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 virialised 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), and for the investigation of extreme objects (e.g. Hoyle et al. 2012).

The eROSITA instrument on board 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 thousand galaxy clusters in the western Galactic hemisphere (Bulbul et al. 2024). Those extended sources are identified in multi-band optical surveys and their redshift is measured, hence providing a distance to us observers. The first catalogues of clusters discovered in the first eROSITA All-Sky Survey (eRASS1) western Galactic hemisphere (Bulbul et al. 2024; Kluge et al. 2024) support detailed individual cluster studies (e.g. Liu et al. 2023; Veronica et al. 2024) and population studies revealing for the first time the large-scale evolution of clusters and groups up to z ≃ 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 catalogues. The first cosmological results based on cluster number counts are presented in Ghirardini et al. (2024) and Artis et al. (2024). In this article, we explain how the selection function model supporting these results was constructed. In fact, virtually any cluster population study based on the published eRASS1 catalogues 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 the development of accurate selection models with powerful statistical methods. The recent Gaia survey is an illustrative instance where selection functions for the parent catalogue (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 the 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 uniquely relying on X-ray count 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 catalogues (X-ray and millimetre-band detections) in Sect. 6. We discuss our results in Sect. 7 and summarise our findings in Sect. 8.

Unless stated otherwise, we assume a flat ΛCDM cosmological model with parameters from Planck Collaboration XIII (2016). We use ln for natural logarithm, and the base ten logarithm is written as log₁₀. The notation 𝒩(μ;σ) indicates a normal distribution centred on μ with standard deviation σ. We use the ″ symbol for arcsecond units.

2. Motivation

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

$$\begin{array}{c}\hfill L(\mathrm{model})\propto p(\mathrm{data}|\mathrm{model}).\end{array}$$

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 Θ_model fully describing it (see Table 1):

$$\begin{array}{c}\hfill L({\mathbf{\Theta }}_{\mathrm{model}})\propto p\left(\mathrm{data}|{\mathbf{\Theta }}_{\mathrm{model}}\right).\end{array}$$(1)

Posterior distributions p(Θ_model|data) 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 Θ_bin: measured mass, redshift, flux, pairwise distances, etc. The likelihood functions L combine observed numbers ${\widehat{N}}_j$ and model-predicted numbers n_j by means of a statistics (j is an index for the bins). A famous example is the Poisson log-likelihood (see also Cash 1979), adequate for shot-noise dominated bins. The size of the bins is part of the likelihood design and it is usual practice to consider infinitesimally small bins containing ${\widehat{N}}_j=0$ or 1 object. In such a case, n_j transforms into dn/dΘ_bin. In presence of additional variance terms, the Poisson 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 catalogues. 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 parameterised 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(I|Θ_sel) the probability of selecting an object in the sample given the value of parameters Θ_sel. These values must be a prediction of the model. The expression for the modelled counts n_j may write formally:

$$\begin{array}{cc}\hfill n_j& ={\displaystyle \int \mathrm{d}}{\mathbf{\Theta }}_{\mathrm{bin}}\mathrm{d}{\mathbf{\Theta }}_{\mathrm{sel}}\mathrm{d}{\mathbf{\Theta }}_{\mathrm{model}}\mathbf{1}({\mathbf{\Theta }}_{\mathrm{bin}}\in j)\hfill \\ \hfill & \times p\left(I|{\mathbf{\Theta }}_{\mathrm{sel}}\right)p\left({\mathbf{\Theta }}_{\mathrm{bin}},{\mathbf{\Theta }}_{\mathrm{sel}}|{\mathbf{\Theta }}_{\mathrm{model}}\right)\frac{\mathrm{d}n}{\mathrm{d}{\mathbf{\Theta }}_{\mathrm{model}}}·\hfill \end{array}$$(2)

In the latter expression we have introduced dn/dΘ_model, the predicted number distribution of galaxy clusters under a certain model assumption. The indicator function 1 ensures only those objects with Θ_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 Θ_bin and Θ_sel so the model is accurate and the modelling effort is well balanced between p(I|Θ_sel) and p(Θ_bin, Θ_sel|Θ_model). Depending on how close Θ_sel is from purely observable quantities, selection function models involve a varying amount of astrophysical modelling. It is also important to recognise 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.

3. The eROSITA galaxy cluster survey and its simulated twin

3.1. The eRASS1 galaxy cluster samples

The primary eRASS1 cluster catalogue in the western Galactic hemisphere (hereafter eRASS1-primary or eRASS1-main) is described in detail in Bulbul et al. (2024) and Kluge et al. (2024). This catalogue 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 characterised using a point-spread function (PSF) fitting algorithm. The detection likelihood (DET_LIKE_0 in the parent catalogue, hereafter shortened to ℒ_det) and the extent likelihood (EXT_LIKE, hereafter ℒ_ext) are among the most relevant parameters for selecting extended sources. In particular, the eRASS1-primary catalogue is compiled based on a low threshold of ℒ_ext > 3 to maximise completeness. A rigorous cleaning process is applied to the X-ray source catalogue 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. (2024). It is compiled with a more conservative cut on the measured parameter ℒ_ext > 6 compared to the eRASS1-primary sample to maximise 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. 2024). The cosmology sample comprises 5259 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 catalogues 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 catalogues (see Ghirardini et al. 2024). Most of the discussion in the paper neglects this effect, and we refer to Kluge et al. (2024) 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 summarise here the main features of the eRASS1 digital twin depicted in Seppi et al. (2022). The parent halo catalogue 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, possibly biased against low-luminosity and low-mass objects. Therefore, low-mass groups M_500c < 5 × 10¹³ M_⊙ require a flux rescaling correction (see appendix of Seppi et al. 2022). This mass-dependent correction is derived by displacing simulated objects along the observed stellar mass–X-ray luminosity relation of groups and clusters (Anderson et al. 2015; Comparat et al. 2022; Zhang et al. 2024). Stellar masses for distinct haloes are predicted via the stellar to halo mass relation (and its scatter) from Moster et al. (2013), see Eq. (1) in Comparat et al. (2019). The correction factor so obtained is negligible at high mass, while it can reach down to 0.1 at 10¹³ M_⊙. This correction procedure preserves scatter in the X-ray luminosity at fixed halo mass. This approach is validated a posteriori by the good agreement between the observed and simulated galaxy cluster flux distributions.

The AGN model (Comparat et al. 2019) reproduces by construction the observed hard X-ray AGN luminosity function from Aird et al. (2015) and the number density as a function of AGN flux at redshifts 0 < z < 4. AGN are associated to dark matter halos and subhalos by combining together a stellar mass to halo mass relation with a duty cycle, and a abundance matching between stellar mass and X-ray luminosity. The observed rate of AGN in galaxies as a function of stellar mass and redshift (i.e. the AGN duty cycle, Georgakakis et al. 2017) is thus naturally reproduced. As the algorithm matches halos to AGN sources, a parametric scatter in stellar mass at fixed AGN luminosity is tuned to a value ensuring the distribution of specific accretion rates follows observed data. Since AGN and clusters are painted onto halos drawn from the same N-body simulation, these two populations are spatially correlated. Except for the link to stellar mass, our simulations make no account for environmental AGN dependence nor dependence on the star formation rate of galaxies (e.g. Aird et al. 2018). Despite being simple, this model produces AGN Halo Occupation Distributions in good agreement with eFEDS data in group- and cluster-size halos, and this is true in particular for central galaxies in those halos (Comparat et al. 2023). Brightest Central Galaxies (BCG) are not treated separately from other halos in our simulation. Because they are hosted in higher mass halos, the stellar mass associated to BCG is higher than that in average field or satellite galaxies; so is their associated AGN X-ray luminosity. In practice, due to the scatter in the abundance matching procedure and to the low value of the duty cycle of high-L_X AGN, it is rare to find a luminous AGN at the centre of a simulated cluster. In fact, only 0.4% of the M > 10¹⁴ M_⊙ simulated halos host a central AGN more X-ray luminous than the cluster itself.

A simulated half-sky twin contains more than 10⁶ simulated halos up to z ∼ 1.5. Besides AGN and clusters, foreground X-ray emitting stars and galactic foregrounds also have their share of simulated photons. A total of 373 316 simulated stars have a flux assigned according to the X-ray flux cumulative distribution measured in eFEDS (Schneider et al. 2022); their positions are selected at random among true Gaia DR2 positions to follow the observed density of stars. Diffuse foregrounds include mainly emission from the local hot bubble, eROSITA bubble, and the galactic halo (Predehl et al. 2020; Ponti et al. 2023; Locatelli et al. 2024; Zheng et al. 2024a,b). Normal galaxies whose X-ray emission is driven by hot gas and X-ray binaries, also contribute to the unresolved background component. The sky density of X-ray emitting galaxies is well below the number of AGN in stars at the limiting flux of eRASS1 (Lehmer et al. 2012), hence we do not simulate these sources individually, rather their contribution is entirely considered as part of the eRASS1 unresolved background component. Our practical implementation relies on resampling real eRASS1 background events obtained after conservatively masking each eRASS1 detected source region. This approach enables the generation of spatially varying background (at the angular resolution of 0.9 deg) that correctly reproduces actual eRASS1 data.

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 × 3.6 deg² 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 ℒ_det and ℒ_ext.

All clusters simulated in a twin mock are associated to a set of properties, such as mass M₅₀₀, true redshift z, X-ray luminosity L_X, flux f_X or count-rate CR (all within R₅₀₀, 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 catalogue 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_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 catalogues. 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-σ error circle, and an angular extent (we take the core-radius of the best-fit β-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 catalogues in presence of positional uncertainties. We first take the input catalogue as reference and look for matches in the detection catalogue. 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_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_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 catalogues enclose complex configurations with multiple input sources projected along neighbouring 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_any. Our experiments showed that the result is rather insensitive to their exact values. For this first pass (input catalogue as reference) we kept only pairs with p_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 catalogues. For this second pass (detection catalogue as reference) we do not input any prior and we chose p_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_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 illustrates the comparison between both matching procedures. They agree well on a wide range of fluxes. Deviations are visible at a low flux, that is, below a few 10⁻¹⁴ erg s⁻¹ cm⁻². 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 catalogues 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.

Fig. 1. Comparing the outcome of two matching procedures relating simulated halos (plain black histogram) to detected sources (photon-based matching is a dashed blue line, and position-based matching is a thin orange line). A simulated halo is flagged as selected (Imain = 1) whenever it is (solidly) matched to an extended source in the detection catalogue. The details of the matching algorithms impact the flux distribution of detected sources, especially at flux below a few 10−14 erg s−1 cm−2, with deviations up to a factor 2 in certain bins.

Fig. 2. Example case of blending, where a bright AGN drives the detection and the presence of a faint extended source favours 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 ∼5 × 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 (Imain = 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 realisation of the twin eRASS1 sky comprises about 10⁶ halos, with less than 4% being detected (I_eRASS1 = 1) and less than 1% being selected (I_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 realisations, and randomly shuffling their content. Although the realisations 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 × 2/3 ≃ 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(I|Θ_sel), this is the main interest of this paper. The choice of Θ_sel usually results from a compromise between precision of the selection function and complexity of the model p(Θ_bin,Θ_sel|Θ_model).

A natural choice for the selection parameters Θ_sel is the collection of pixel values that form the images of clusters. Then p(I|image) is obtained by feeding eSASS with an image and applying thresholds on the values (ℒ_ext, EXT) returned by the algorithm; because eSASS is deterministic, p(I|image) 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′×10′ 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⋅hours, which is unrealistically long for a cosmological analysis. To this cost must be added the time to compute a model cluster image given Θ_model.

Realising 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 emphasise that the number of counts should be a prediction of the model p(counts|Θ_model). 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″ cluster 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″. We transform the counts (integer values N) with the following:

$$\begin{array}{c}\hfill \eta =ln(1+\frac{\mathrm{\Sigma }}{c}).\end{array}$$(3)

We set Σ = N/(πR²) and c = 0.0002 counts arcsec⁻². In our training set, N ranges from 0 to 3194 counts within 90″ aperture. Figure 3 displays the histogram of η associated to all halos in the training set, and for those flagged as selected (I_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 (ℒ_ext > 3) given the value η. 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₀ and the slope w_η of the linear model f(η) = w₀ + w_ηη. The model probabilities are such that

$$\begin{array}{c}\hfill p\left(I_{\mathrm{main}}|\eta \right)={(1+e^{-f(\eta )})}^{-1}.\end{array}$$(4)

Fig. 3. Distribution of clusters in the training set as a function of Σ (units counts arcsec−2), and the average surface brightness in the 90″ radius around their centre. The x-axis is rescaled with c = 2 × 10−4 arcsec−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 Appendix C. The vertical dashed line indicates the transition $\widehat{\eta }=1.94$ in the logistic model described by Eqs. (3) and (4) – corresponding to N ≃ 30 counts.

The cost function to minimise during the fitting procedure is the log-loss 𝒞(w₀, w_η), its expression is shown in Eq. (A.2).

Thanks to the one-to-one relation between η and N (Eq. (3)), we thus obtain a model p(I_main|N). We can interpret the values of the best-fit coefficients as follows: w_η governs the sharpness of the transition in the S-shaped curve, $\widehat{\eta }=(-{\mathit{w}}_0/{\mathit{w}}_{\eta })$ is the value at which this transition occurs. In our experiment, we find w_η = 4.59 ± 0.04 and w₀ = −8.92 ± 0.06, leading to $\widehat{\eta }=1.94$. This corresponds to a transition of the S-shaped curve taking place at N ≃ 30 cluster counts in the 90″ circular region. Positiveness of w_η reflects the fact that detection probability increases with number of counts N. The computation of coefficient uncertainties is detailed in Appendix 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 seven 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 × 7 = 35-element feature vector η = {η_j}_{j = 1…35} and perform logistic regression, involving a total of 36 coefficients {w_j}_{j = 0…35} (one for each feature, plus the intercept). Introducing Σ_j the average surface brightness in an annulus, we write the following:

$$\begin{array}{c}\hfill f(\eta )=w_0+{\displaystyle \sum _{j=1}^{35}}{w}_j{\eta }_j,with{\eta }_j=ln(1+\frac{{\mathrm{\Sigma }}_j}{c}).\end{array}$$(5)

One may interpret the value of the coefficients w_j as the sensitivity of the detection rate p(I_main|η) to a small variation of the surface brightness of one component in one annulus. Indeed we have

$$\begin{array}{c}\hfill {\frac{\partial p}{\partial {\eta }_j}|}_{{\eta }_k}=p(1-p)w_j,\end{array}$$

where the derivative is performed at constant η_k (k ≠ j). It yields the following:

$$\begin{array}{c}\hfill {\frac{\partial p}{\partial {\mathrm{\Sigma }}_j}|}_{{\mathrm{\Sigma }}_k}=\frac{p(1-p)}{{\mathrm{\Sigma }}_j+c}{w}_j.\end{array}$$(6)

In other terms, everything else maintained fixed, a small relative variation ϵ = ΔΣ_j/Σ_j of the surface brightness of one component in one given annulus leads to a variation of the detection probability proportional to w_jϵ. In particular, positive coefficients relate to those features that marginally contribute to selecting a cluster, and conversely for negative coefficients. Figure 4 conveniently represents the value of the 35 coefficients w_j (j ≠ 0) split into each component and annulus; the last coefficient not shown on this figure is w₀ = −8.82 ± 0.08.

Fig. 4. Representation of the 35 coefficients wj of a logistic regression model p(Imain|counts), trained to predict whether a cluster was 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″ 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″ surface brightness is not surprising. A detected source (I_eRASS1 = 1) must be classified as extended in order to be selected (I_main = 1); from this perspective, it is much more profitable to increase the counts beyond the eRASS point-spread function radius (∼30″). This argument may appear more clearly from Fig. B.1, 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″ 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′ 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. Figure 5 represents the so-called “precision-recall” curve. This curve is obtained by scanning values of a threshold set to convert probabilities p(I_main|counts) into a binary classification I_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 recall¹ 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 × 5) models. This agrees with common sense, in that adding more detailed information provides more precise models.

Fig. 5. Performance of three different logistic models p(Imain|counts) 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″. The dashed grey line stands for a model using as input only the 90″ 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).

A second test is shown on Fig. 6. This evaluation is performed by binning test clusters by their value of p(I_main|counts) 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 η ≃ 2.5 − 4 (corresponding to N ≃ 60 − 270 counts in the 90″ aperture).

Fig. 6. Reliablity of two different logistic models p(Imain|counts) 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 (Appendix 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 × 5 parameters as input, grey dots and dashed lines for the model using only cluster counts within 90″ 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 clarity 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 idealised 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(counts|Θ_model) represents a computational bottleneck.

5. Selection models with intermediate variables

We may instead choose Θ_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 Θ_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 2012). Details on the procedure and algorithm are discussed in Appendix D. The very nature of GPC allows to issue not only an estimate of p(I|Θ_sel), 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 Θ_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. 2024). In Appendix E we describe several other models, notably those including cluster morphology in their input parameters.

5.1. Mass and redshift

Setting Θ_sel = {M_500c, z} allows to gain insights on the selection process in the fundamental parameter space of cluster cosmology. Figure 7 is a visualisation 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 ≃ 5 × 10¹⁴ M_⊙.

Fig. 7. Representation of the model p(Imain|M500c, z) predicting a cluster to be detected and selected with an extent likelihood above three. The explanatory variables (features) are labels attached to simulated clusters in the twin simulation, standing for galaxy cluster M500c 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 M500c = 3 × 1014 M⊙ (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 XIII 2016) imprinted in the simulation.

5.2. Luminosity and redshift

Setting Θ_sel = {L_X, z} advantageously removes the strong dependence of the model upon the M₅₀₀ → L_X relation that is imprinted in the twin simulation. It is therefore up to the cosmologist to model the relation p(L_X|Θ_model) with (most likely) parameterised 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.

Fig. 8. Representation of the model p(Imain|LX, z) predicting a cluster to be detected and selected with an extent likelihood above three. The explanatory variables (features) are labels attached to simulated clusters in the twin simulation, standing for galaxy cluster luminosity LX measured in the 0.5–2 keV energy band at the cluster rest-frame (units erg s−1, rescaled in log10) 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. Figure 9 is a partial visualisation of the model output, all following quantities but exposure time T_exp being fixed at z = 0.3, N_H = 3 × 10²⁰ cm⁻², background brightness 5.2 × 10⁻¹⁵ erg s⁻¹ cm ⁻² arcmin⁻² in the 0.3–2.3 keV energy band, L_X = 10⁴³ or 10⁴⁴ erg s⁻¹. 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.

Fig. 9. Representation of the models p(Imain|LX, z, NH, Texp, bkg) and p(Icosmo|LX, z, NH, Texp, bkg) for fixed values of cluster redshift z, two different values for the luminosity LX, and local galactic absorption NH, in addition to a nominal background level and a range of exposure times Texp (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 Θ_sel as a five-component feature, namely: the cluster redshift z, sky position-dependent quantities ℋ≡ (N_H, T_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₅₀₀ aperture and unaffected by shot noise. It derives from the formula:

$$\begin{array}{c}\hfill CR=\frac{L_X}{4\pi d_L{(z)}^2\mathrm{ECF}(z,k_{B}T,N_H=0)}({\mathrm{units}\mathrm{s}}^{-1}).\end{array}$$(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_BT, 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. Figure 10 represents a few slices through the multi-dimensional model p(I_cosmo|CR, z, N_H, T_exp, bkg), 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.

Fig. 10. Representation of a model p(Icosmo|CR,z,NH,Texp,bkg) for two values of cluster redshift z = 0.2 and z = 0.5, fixed local galactic absorption NH, a nominal background level, two values of exposure times Texp (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(CR,z|Θ_model) fast and easy in likelihood (Eq. (2)). However, this simplicity comes at a cost, since all variations of the selection due to cluster morphology, for example, are marginalised 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 Θ_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.

Figure 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 Θ_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 ℋ = (N_H, T_exp, 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₀ parameter (as presented in Appendix E). The average precision increases to 71%, closer to that obtained with the 7 × 5 surface brightness model depicted in Fig. 5. Finally, a model involving (absorbed) flux, apparent R₅₀₀ radius and ℋ shows good performance with an average precision reaching 68%.

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₀, see Appendix E) provides a noticeable enhancement of the model precision (average precision 74%).

Fig. 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 the legend. A model involving luminosity, redshift, morphological features and local sky information performs best (yellow dot-dashed curve).

Fig. 12. Similar to Fig. 11, but 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. (2024).

We now assess the reliability of the probabilistic output of the models. In Fig. 13 we compare the same five models for the primary catalogue 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 catalogue. The result for the cosmology sample classifiers appear in Fig. 14 and they all are very well aligned along the one-to-one line.

Fig. 13. Evaluating the reliability of five models predicting the probability of a cluster to appear in the eRASS1 primary cluster catalogue (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.

Fig. 14. Similar to Fig. 13, but 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 catalogue R. In general both catalogues 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 catalogue an associated probability:

$$\begin{array}{c}\hfill P_i=p(I_T^i|I_R^i,{\mathbf{\Theta }}_{\mathrm{meas}}^i).\end{array}$$(8)

There, I_T and I_R are random boolean variables indicating whether an object is listed in catalogue T and R, respectively. Θ_meas is a collection of measured properties attached to the cluster: its measured flux, redshift, luminosity, extent, etc. are likely candidates. Using Bayes’ rule, P_i writes as the ratio of two joint probability distributions, the denominator being interpreted as the number of objects expected in the reference catalogue with measured properties ${\mathbf{\Theta }}_{\mathrm{meas}}^i$; the numerator as the number of objects expected in both the reference and test catalogues with these same measured properties. In this perspective, P_i is a model for the fraction of reference objects listed with some measured properties ${\mathbf{\Theta }}_{\mathrm{meas}}^i$ that are also present in the test sample. Introducing Θ_model and Θ_sel enables hierarchical Bayesian expansion of the joint probabilities into elementary factors. The numerator specifically involves the conditional probability distribution:

$$\begin{array}{cc}\hfill p(I_T^i,I_R^i,{\mathbf{\Theta }}_{\mathrm{meas}}^i|{\mathbf{\Theta }}_{\mathrm{model}})& ={\displaystyle \int \mathrm{d}}{\mathbf{\Theta }}_{\mathrm{sel},R}\mathrm{d}{\mathbf{\Theta }}_{\mathrm{sel},T}\mathrm{d}\mathit{\alpha }\hfill \\ \hfill & \times p({\mathbf{\Theta }}_{\mathrm{sel},T},{\mathbf{\Theta }}_{\mathrm{sel},R},\mathit{\alpha }|{\mathbf{\Theta }}_{\mathrm{model}})\hfill \\ \hfill & \times p\left(I_R^i|{\mathbf{\Theta }}_{\mathrm{sel},R}\right)p\left(I_T^i|{\mathbf{\Theta }}_{\mathrm{sel},T}\right)p\left({\mathbf{\Theta }}_{\mathrm{meas}}^i|\mathit{\alpha }\right),\hfill \end{array}$$

and similarly for the denominator:

$$\begin{array}{cc}\hfill p(I_R^i,{\mathbf{\Theta }}_{\mathrm{meas}}^i|{\mathbf{\Theta }}_{\mathrm{model}})& ={\displaystyle \int \mathrm{d}}{\mathbf{\Theta }}_{\mathrm{sel},R}\mathrm{d}\mathit{\alpha }p\left({\mathbf{\Theta }}_{\mathrm{sel},R},\mathit{\alpha }|{\mathbf{\Theta }}_{\mathrm{model}}\right)\hfill \\ \hfill & \times p\left(I_R^i|{\mathbf{\Theta }}_{\mathrm{sel},R}\right)p\left({\mathbf{\Theta }}_{\mathrm{meas}}^i|\mathit{\alpha }\right).\hfill \end{array}$$

The last two expressions are valid as long as selection in the test survey only depends on variable Θ_sel, T and selection in the reference survey only depends on variable Θ_sel, R. We have introduced a new, intermediate variable α to reflect the dependence of the measurement vector. For instance, if the measurement vector is X-ray flux, α may represent the cluster luminosity and redshift. The interest of this decomposition lies in it involving the selection functions p(I|Θ_sel), which are tabulated for each survey. It also accounts for covariances in the selection and measurement processes, through the multivariate distribution of (Θ_sel, T, Θ_sel, R, α) given the values of Θ_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 catalogue. By comparing this probability to the actual presence of a match in the test catalogue, 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 catalogue despite their low forecast probability (P_i < 0.025 and $I_T^i=1$). We will denote “missed objects” those listed in R and not detected in the test catalogue despite their high forecast probability (P_i > 0.925 and $I_T^i=0$).

Our procedure to find matches between the reference and the test catalogues relies on simple positional match. For efficiency reason we trim the reference catalogue 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 catalogues. 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 catalogues, 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 catalogue, but they are not (δ₁) and the fraction of objects that should not be in the test catalogue, but they are present (δ₂). We compute posterior distributions for δ₁ and δ₂ by sampling the following likelihood (product of independent Bernoulli likelihoods):

$$\begin{array}{cc}\hfill \mathcal{L}({\delta }_1,{\delta }_2)& ={\displaystyle \prod _{i\in R}}{[P_i(1-{\delta }_1)+(1-P_i){\delta }_2]}^{I_T^i}\hfill \\ \hfill & \times {[P_i{\delta }_1+(1-P_i)(1-{\delta }_2)]}^{1-I_T^i}.\hfill \end{array}$$(9)

The product runs over all clusters in the reference catalogue, and $I_T^i$ equals 1 if the cluster is matched in the test catalogue, it equals 0 otherwise. If the posterior distribution for δ₁ significantly excludes zero, it indicates that part of the clusters predicted to be in the test catalogue actually escape our model and are not detected for some reason. Similarly, a distribution of δ₂ 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 δ₁ and δ₂ 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. eFEDS as reference, eRASS1-cosmo as test (Sect. 6.2);

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

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

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

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

6. SPTpol-ECS as reference, eRASS1-cosmo as test (Sect. 6.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 catalogues.

6.2. Testing the eRASS1 cosmological cluster sample against eFEDS

The deep eFEDS survey consists of a 140 deg² 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 Θ_model = (M₅₀₀, z) representing the (true) mass within R_500c and the (true) redshift of halos. The prior distribution p(Θ_model) is the halo mass function dn/dM₅₀₀dz. 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 catalogue footprint. Each entry is associated to a measured redshift z_eFEDS (column z in the catalogue) and to an estimated 0.5–2 keV flux within a 300 kpc aperture (column F_300kpc, expressed in units erg s⁻¹ cm⁻²). These two elements constitute our vector of measured features Θ_meas. We construct a generative model for flux and redshift that depends on α = (L_X, z, M₅₀₀), where L_X is a cluster (true) 0.5–2 keV luminosity integrated in a cylinder of radius R₅₀₀. Our generative model further assumes independence of flux and redshift measurements, that is:

$$\begin{array}{c}\hfill p\left({\mathbf{\Theta }}_{\mathrm{meas}}|\mathit{\alpha }\right)=p\left(z_{\mathrm{eFEDS}}|z\right)\times p(ln\mathtt{F}\_\mathtt{300}\mathtt{kpc}|L_X,z,M_{500}).\end{array}$$(10)

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

$$\begin{array}{c}\hfill p\left(z_{\mathrm{eFEDS}}|z\right)\sim \mathcal{N}(z;0.01(1+z)).\end{array}$$(11)

We suppose a log-normally distributed flux F_300kpc given the values of α, that is:

$$\begin{array}{c}\hfill p(ln\mathtt{F}\_\mathtt{300}\mathtt{kpc}|L_X,z,M_{500})\sim \mathcal{N}(ln{\mu }_{\mathtt{F}\_\mathtt{300}\mathtt{kpc}};{\sigma }_{ln\mathtt{F}\_\mathtt{300}\mathtt{kpc}}).\end{array}$$(12)

In this expression, both μ_{F_300kpc} and σ_{lnF_300kpc} depend on (L_X, z, M₅₀₀) as follows. We first convert any set of values (L_X, z) into a 0.5–2 keV flux within R₅₀₀ 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 2016). Aperture correction (from R₅₀₀ 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₅₀₀ (itself deriving from M₅₀₀ and z). We assume unbiased flux measurements and set μ_{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 ϵⁱ:

$$\begin{array}{c}\hfill {\sigma }_{ln\mathtt{F}\_\mathtt{300}\mathtt{kpc}}=0.14{\left(\frac{{\mu }_{\mathtt{F}\_\mathtt{300}\mathtt{kpc}}}{{10}^{-13}{\mathrm{erg}\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}}\right)}^{-0.47}+{ϵ}^i.\end{array}$$(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 ϵⁱ is specific to each cluster (labelled by i). As an assumption, ϵⁱ is the deviation of the reported error to the value predicted by the empirical power-law model. The eFEDS cluster catalogue 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 μ_{F_300kpc} is taken to be half of the upper limit value; the dispersion is such that ϵⁱ is the deviation of the reported upper limit to the fixed value of 2.5 × 10⁻¹⁴ erg s⁻¹ cm⁻² (median of all upper limits in the catalogue). 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_eFEDS of an eFEDS-detected cluster, for any set of (true) values α = L_X, z, M₅₀₀.

The reference survey selection function $p\left(I_R^i|{\mathit{\Theta }}_{\text{sel},R}\right)$ is the eFEDS selection function p(I_eFEDS|L_X,z,EM₀,T_exp). It is constructed as described in Liu et al. (2022), with addition of the morphological parameter EM₀ characterising the central emissivity of the clusters (Appendix E). The exposure time T_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(I_cosmo|Θ_sel, T) 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. 2024, see also Fig. 10), involving count-rate CR, redshift z as well as position-dependent, cluster-independent, selection features ℋⁱ. 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_λ < 0.8. We assume Gaussian-distributed errors for the eRASS1 redshifts with standard deviation σ_z/(1 + z) = 0.015, see also Eq. (14). The criterion IN_ZVLIM = True, recommended by Kluge et al. (2024) is present in our model through the implementation of a sky mask.

The last ingredient required in the model is the distribution p(Θ_sel, T,Θ_sel, R,α|Θ_model). Given the assumptions enumerated in this section, it reduces to p(L_X|M₅₀₀,z). We take the scaling relation of Bulbul et al. (2019) relating the core-included 0.5–2 keV luminosity within R₅₀₀ to the M₅₀₀ 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 catalogue of clusters. By performing a 2 arcmin match between the eFEDS and eRASS1 catalogue, we find 35 actual clusters in common between both samples (see also Bulbul et al. 2024), and 90% of the matches are separated by less than 30″.

Figure 15 summarises 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 ± 3.2 eFEDS clusters to be in the eRASS1 cosmology sample, while the actual number is 35 (that is, a 5-σ underestimate). We find the empirical parameters δ₁ and δ₂ 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 (Sect. 7.2).

Fig. 15. Testing the models with eRASS1-cosmo as a test catalogue, and eFEDS as a reference catalogue. 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 (zeFEDS). 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 Pi (Eq. (8)) computed for each eFEDS cluster.

Fig. 16. Testing the models with eRASS1-cosmo as a test catalogue, and eFEDS as a reference catalogue. The 531 eFEDS entries are grouped in bins of the model output probabilities Pi (horizontal axis). In each bin, the actual fraction of matches in the eRASS1 catalogue is reported on the vertical axis. Error bars represent the estimated 68% confidence range (Appendix 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 catalogue 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 8deg.

To describe each of the eRASS1 detections, we constructed a model for the flux and redshift measurements of each entry. We assume independence of the estimated redshift z_λ (column Z_LAMBDA in the eRASS1 cosmology catalogue) and the estimated R₅₀₀ flux in the band 0.5–2 keV (F500_0520, units erg s⁻¹ cm⁻²) and we form the model:

$$\begin{array}{c}\hfill p\left(z_{\lambda }|z\right)\sim \mathcal{N}(z;0.015(1+z)).\end{array}$$(14)

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

$$\begin{array}{c}\hfill p(ln\mathtt{F}\mathtt{500}\_\mathtt{0520}|L_X,z)\sim \mathcal{N}(ln{\mu }_{\mathtt{F}\mathtt{500}\_\mathtt{0520}};{\sigma }_{ln\mathtt{F}\mathtt{500}\_\mathtt{0520}}).\end{array}$$(15)

In this equation, μ_{F500_0520} and σ_{F500_0520} depend on L_X the true luminosity within R₅₀₀ and on the true redshift z. Assuming unbiased flux measurements, μ_{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 2016). As for the uncertainty model, we proceed by fitting a simple model linking the catalogue flux and errors reported in Bulbul et al. (2024). We take:

$$\begin{array}{c}\hfill {\sigma }_{ln\mathtt{F}\mathtt{500}\_\mathtt{0520}}=0.48{\left(\frac{{\mu }_{\mathtt{F}\mathtt{500}\_\mathtt{0520}}}{{10}^{-13}{\mathrm{erg}\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}}\right)}^{-0.51}+{ϵ}^i.\end{array}$$(16)

The quantity ϵ_i is specific to each cluster in the eRASS1 sample and equals the deviation of the catalogue uncertainty relative to the power-law model.

Figure 17 illustrates the outcome of this comparison exercise. The model predicts 35.7 ± 1.3 eRASS1 clusters to be found in eFEDS. The actual number of matches is 36². 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 ≃ 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.

Fig. 17. Similar to Fig. 15, but 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 R500 reported in the eRASS1 catalogue (units 10−14 erg s−1 cm−2). The horizontal axis is the catalogue redshift zλ. Clusters falling in zones of zero eFEDS exposure have vanishing probability (Pi = 0). Other eRASS1 clusters are very likely to be detected in eFEDS (Pi ≃ 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 × 10−13 erg s−1 cm−2.

The model correctly explains the outcome of matching eFEDS to the eRASS1 cosmology catalogue. 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 catalogue (respectively). The values of δ₁ and δ₂ are themselves compatible with zero, indicating again consistency between the model and the present catalogues being tested.

6.4. Testing the eRASS1-main sample against eFEDS

We now turn to a test using the extended version of the eRASS1 catalogue, comprising all X-ray clusters with ℒ_ext > 3 (Bulbul et al. 2024). In contrast to the cosmology sample this catalogue 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 Eqs. (10)–(13). The main modification consists in using a model selection function $p\left(I_T^i=I_{\text{main}}^i|{\mathit{\Theta }}_{\text{sel},T}\right)$. We take as selection variables the count-rate CR, true redshift z and sky position-dependent parameters ℋ.

Our model predicts 30.5 ± 4.2 matches, while the actual number is 61 (hence a larger than 7-σ 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 catalogue. 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 catalogue. 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 (δ₁, δ₂) 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 catalogue paper (Liu et al. 2022). We will return to this issue in Sect. 7.

Fig. 18. Similar to Fig. 16, but using the primary sample of clusters eRASS1-main as a test sample, instead of the eRASS1 cosmology sample.

Fig. 19. Posterior distribution for the pair of parameters (δ1 and δ2) introduced in Eq. (9), for the experiment using eRASS1-main as a test and eFEDS as a reference sample. The histograms represent the marginalised posterior distributions, with mean and 68% range overprinted. The contours represent the equivalent 0.5-, 1-, 1.5- and 2-σ distribution of the joint posterior. Clearly the parameter δ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 catalogue of 89 galaxy clusters discovered via the Sunyaev–Zeldovich (SZ) effect in the SPTpol-100d field, located at RA = 23h30m and Dec = −55deg. In the area bounded by 23h < RA < 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 catalogue 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. At higher redshifts, the selection function model takes an analytic form involving the SZ signal-to-noise parameter ξ and the unbiased significance ζ_SPTpol (see Appendix F). The measurement vector Θ_meas is taken from the SPTpol-100d catalogue, namely the xi and redshift columns, standing respectively for ξ and the measured redshift z_SPT. We again assume approximate independence of both measurements. We assume unbiased redshift measurements with Gaussian uncertainties 0.02(1 + z). We use again Eq. (F.2) to generate the distribution for ξ, 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. Consequently we need to construct a joint model for the distribution of luminosities and ζ_SPTpol. We combine the L_X − M₅₀₀ − z scaling relation mentioned previously (Bulbul et al. 2019) and the ζ − M₅₀₀ − z relation³ of Huang et al. (2020) into a bivariate log-normal distribution with covariance matrix:

$$\begin{array}{c}\hfill \mathrm{Cov}(ln{L}_X,ln{\zeta }_{\mathrm{SPTpol}})=\left(\begin{array}{cc}{\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\end{array}\right).\end{array}$$

Here, σ_lnLX and σ_lnζ are the log-normal scatters of both distributions at fixed mass M₅₀₀ and redshift z (with values respectively 0.27 and 0.18). The value of ρ 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 ρ = 0.2, we will investigate the impact of varying this parameter in Sect. 7.

The model predicts 11.4 ± 2.4 matches, while the actual number of matches is 8. Figure 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.

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

Figure 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 catalogues nor in the model. In particular, the selection function models account correctly for the observed matches between SPTpol-100d and eRASS1-cosmo.

Fig. 21. Similar to Fig. 16, but using the cosmology sample of eRASS1 clusters as a test sample, and the SPTpol-100d cluster sample as a 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 catalogue with that of the eRASS1 cosmology sample. We thus test the SPTpol catalogue, 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(L_X,ζ_SPTpol|M₅₀₀,z) 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_λ listed in the eRASS1 catalogue. Combining together the 19 p_i values returned by our model, we predict 7.0 ± 1.3 matches, while the actual number is 8, indicating consistency between the results and the predictions. We have also checked the adequacy between the values of P_i predicted by the model and the actual fraction of matches in each bin. Despite its statistical limitations, this experiment demonstrates good behaviour of the model in the high-redshift, intermediate-mass range.

Two cases deserve further examination. The eRASS1 cluster 1eRASS J230029.2−510022 is measured at z_λ = 0.37 with F500_0520 = 3.4 × 10⁻¹³ erg s⁻¹ cm⁻² and should be massive enough to appear in the SPTpol survey. In fact our models predicts P_i = 0.90. It is actually missed by SPTpol, 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_λ = 0.31 with F500_0520 = 3.4 × 10⁻¹³ erg s⁻¹ cm⁻² is not detected by SPTpol despite the model predicted P_i = 0.72. This is due to a sky mask covering a point source brighter than 6 mJy and located at a position coinciding with the galaxy cluster (Bleem, L., priv. comm.).

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

The model described for the SPTpol-100d catalogue is straightforwardly adapted to the SPTpol-ECS catalogue presented in Bleem et al. (2020). The ECS sample is collected over 2770 deg² area in two regions of the Southern Hemisphere. A total of 163 clusters with significance ξ > 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, except that the unbiased significance parameter ruling SZ detections is denoted ζ_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 γ_field specific to each sub-field of the SPT-ECS survey. An additional multiplicative factor γ_ECS = 1.124 rescales all values of A_SZ, thus reflecting calibration updates (Table 1 and Sect. 5 in Bleem et al. 2020). We keep a fixed value ρ = 0.2 for the correlated scatter between lnL_X and lnζ_ECS. There are 103 matching instances between both catalogues, based on a 2′ positional match. Our model predicts 103 ± 4 matches, in excellent agreement with the observed number. Figure 22 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 ≃ 0). It is listed in the ECS catalogue with significance 5.1 and photometric redshift z_SPT = 1.05 ± 0.04. The matching instance in eRASS1 is 1eRASS J033323.9−370744, listed at z_λ = 0.52 ± 0.01. This cluster is one of the 3% ambiguous cases discussed in Kluge et al. (2024), where two unrelated clusters overlap along the same line of sight.

Fig. 22. Similar to Fig. 20, but using 163 clusters in SPTpol-ECS catalogue as a reference sample and the eRASS1 cosmology sample as a test.

Figure 23 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 catalogues.

Fig. 23. Similar to Fig. 21, but obtained by using the SPT-ECS sample as a reference (163 entries) and using the cosmology sample of eRASS1 clusters as a 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 catalogue falling within the sky footprint of SPTpol-ECS survey. The model takes as measurement features F500_0520 and z_λ and it assumes identical generative distributions as in Sect. 6.3.

Figure 24 shows the outcome of this calculation. Our model predicts 65.4 ± 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. These cases are discussed in Appendix G.

Fig. 24. Similar to Fig. 17, but using 958 eRASS1 clusters from the cosmology sample as a reference, and the SPTpol-ECS catalogue 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 25 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₅₀₀ − 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 ξ due for instance to Malmquist bias. As discussed in Sect. 7.3, progress in this direction requires implementation of a self-consistent ICM model propagated to the X-ray and SZ selection functions. We therefore defer detailed discussion of this apparent offset to future work, in particular through the use of deeper SPT-SZ cluster samples (e.g. SPT-3G).

Fig. 25. Similar to Fig. 16, but obtained by using the eRASS1 cosmology sample (958 entries) as a reference and using SPTpol-ECS as a test sample.

7. Discussion

7.1. Linking our models with ℒ_ext and ℒ_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 (ℒ_det) and extent (ℒ_ext) likelihoods. Applying thresholds on these quantities leads to defining a new sample, characterised by a balance between contamination and completeness (e.g. Seppi et al. 2022). In the specific case of the eRASS1 cluster catalogue, a source is detected if it shows ℒ_det greater than 5. The threshold on ℒ_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 catalogue parameters cut only (here collectively designed under a generic name, ξ_X). In fact, let us acknowledge that:

$$\begin{array}{c}\hfill p\left(I|{\mathbf{\Theta }}_{\mathrm{model}}\right)=p({\xi }_X>{\xi }_X^{\min }|{\mathbf{\Theta }}_{\mathrm{model}})={\displaystyle {\int }_{{\xi }_X^{\min }}^{+\infty }}\mathrm{d}{\xi }_{X}p\left({\xi }_X|{\mathbf{\Theta }}_{\mathrm{model}}\right).\end{array}$$

Building a selection model then translates into constructing a distribution of ξ_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 ξ_X and those variables, involving some parameterised 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 catalogue 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 catalogues), 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, for example 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 catalogue.

A safer alternative may consist in constructing a model ξ_X directly from ℒ_ext and ℒ_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 ℒ_det (resp. ℒ_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 ℒ_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 ℒ_det and ℒ_ext for the logistic regression involving 35 features (Fig. 26). 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.

Fig. 26. Relation between the latent variable f(η) = w0 + ∑jwjηj calculated for each simulated cluster in the test sample after performing logistic regression with 35 surface brightness features, and the values of ℒdet and ℒ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 ℒdet value (red histogram in the lower panel). Clusters detected but not categorised as extended have zero value for ℒext (orange histogram). The apparent correlation between f and both measurements indicate that the model has “discovered” the importance of ℒdet and ℒext in the selection process.

We find a similar result with the mean of the latent variable f governing the predictions output of the Gaussian process classifier (see Appendix D). Figure 27 in particular shows the a posteriori correlation between ⟨f⟩ and the detection and extent likelihood obtained in a test sample of simulated clusters.

Fig. 27. Similar to Fig. 26, but replacing the x-axis by the mean of the latent variable implicit in the Gaussian process classifier p(Icosmo|CR,z,ℋ).

7.2. Potential 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-σ when considering the eRASS1 cosmology sample and above 7-σ 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 (Eqs. (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. 2024).

We obtained very similar predictions, and thus similarly discrepant results, when modelling the measured R₅₀₀ X-ray luminosity of eFEDS clusters (column Lx500 in the catalogue 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:

$$\begin{array}{c}\hfill {\sigma }_{\mathtt{Lx}\mathtt{500}}=0.2\times {10}^{43}{\mathrm{erg}\mathrm{s}}^{-1}{\left(\frac{L_X}{{10}^{43}{\mathrm{erg}\mathrm{s}}^{-1}}\right)}^{0.9}+{ϵ}^i.\end{array}$$(17)

In this expression ϵⁱ is a constant for each eFEDS cluster indexed by i. We obtained ϵⁱ by subtracting the value predicted by the power-law model (without ϵⁱ) in the above equation to the uncertainty reported in the catalogue 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 optimise many hundreds of Gaussian process hyperparameters at once (Appendix D). We have taken benefit of our multiple Poisson realisations 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-σ 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⁻¹⁴ erg s⁻¹ cm⁻², 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-σ 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⁻¹³ erg s⁻¹ cm⁻², 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 catalogue, 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₅₀₀ in their features (Fig. E.3). By using a selection model in the form p(I_main|f_X,R₅₀₀,ℋ), we do not rely any longer on the distribution of sizes imprinted in the simulation. Rather, we assume our L_X − M₅₀₀ is correct and we distribute our modelled flux within the radius associated to M₅₀₀. 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-σ for the eRASS1-main catalogue and from 5- to 3-σ for the cosmology catalogue. This finding corroborates Fig. 6 in Comparat et al. (2020), showing that the average R₅₀₀ radius of our simulated clusters at fixed luminosity is larger than in several published samples⁴.

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(lnL_X,lnζ|M₅₀₀,z) as a bivariate normal distribution with correlation coefficient ρ = 0.2. We now investigate its impact on the results. We varied ρ 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:

$$\begin{array}{c}\hfill B=-2ln\mathcal{L}=-2{\displaystyle \sum _{i\in R}}{I}_T^{i}ln{P}_i+(1-I_T^i)ln(1-P_i).\end{array}$$(18)

Its mimimum B_min is obtained by varying ρ, 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_min approximately follows a χ² distribution with one degree of freedom. Figure 28 shows the curve B − B_min for various values of ρ, from which we derive an approximate 68% confidence interval −0.4 < ρ < 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.

Fig. 28. Value of the combined Bernoulli negative log-likelihood (Eq. (18)) in the comparison of the SPT-ECS sample (reference) and the eRASS1 cosmology sample (test) for various values of the correlation coefficient ρ. This coefficient expresses the level of covariance between lnLX and lnζ at fixed halo mass and redshift. Assuming all other model parameters are correct, the constraint obtained from the catalogue comparison exercise depicted in this paper is relatively loose: ρ ≃ 0.2 ± 0.6 (68% confidence level) as deduced from $(B-B_{\text{min}})~{\chi }_1^2$.

The fiducial value ρ = 0.2 chosen in Sect. 6 is in fact the best-fit, with B_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_min. Given the many other assumptions made (e.g. fixed scaling relation parameters, etc.), we cannot firmly exclude specific values of ρ from our dataset. A comparison involving larger samples may put stronger constraints on ρ through this method.

It is interesting to note that ρ is intimately related to the ICM distribution for clusters in the reference survey. A physical link should exist between L_X and ζ that allows computing ρ or at least a prior distribution, without using the halo mass M_500c. However, our implementation of the SPT-SZ selection function through ζ 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 for our purpose since the selection of sources in astronomical catalogues results (at least in part) from human decisions, whose detailed physical modelling may be cumbersome. The models we selected are intentionally simple, enabling them to be interpreted and explained.

Our main findings are summarised as follows:

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 catalogues.

2. We have identified which simulated halos were 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⁻¹⁴ erg s⁻¹ cm⁻².

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″ counts help in detecting a cluster, but not so in characterising it as extended. Cluster counts beyond approximately 1.5′ 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. Acknowledging that the selection function based on surface brightness features requires an important modelling effort on the cosmological likelihood side, we have developed models relying on intermediate variables. We identified a hierarchy of models, with some strongly depending on assumptions made in the twin simulations, for example p(I|M₅₀₀,z), and others enabling less of a dependence on such assumptions, for example p(I|CR,z,T_exp,N_H,bkg). Using the latter class of models comes at a cost, namely a heavier modelling effort in scientific analyses.

5. We have tested all our models with leftover 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₅₀₀ radius or its central emission measure EM₀) perform better.

6. Using real data, we have tested our models by means of a custom population model. In a comparison run involving 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, which is in agreement with the known low contamination of the sample (Seppi et al. 2022; Ghirardini et al. 2024).

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.

The outcome of the numerical deliverables for this work have been used directly in other eRASS1 cosmology papers. The eRASS1 group population analysis (Bahar et al. 2024) considers the model depicted in Sect. 5.2. The baseline model presented in Sect. 5.3 feeds the main cosmological analysis based on cluster abundance (Ghirardini et al. 2024; Artis et al. 2024) and on cluster clustering (Seppi et al. 2024). This model also serves in estimating sample completeness and in assigning total masses to catalogued clusters (Bulbul et al. 2024). Other models presented in this paper have indirectly fed these analyses by providing a critical, quantitative assessment of our understanding of X-ray cluster selection in the context of eRASS1.

In future work we will refine the simulation set with updated cluster models and halo simulations. An alternative assessment 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, the ultimate validation of selection models should rely on a performance metric that directly relates 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 a 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 marginalising 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.

Acknowledgments

The authors thank the referee for a careful reading of the manuscript. 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), MOCpy (https://github.com/cds-astro/mocpy/) (Fernique et al. 2014), GPy (https://github.com/SheffieldML/GPy) (GPy 2012), astropy (Astropy Collaboration 2022), matplotlib (Hunter 2007), healpy and HEALPix (http://healpix.sourceforge.net) (Górski et al. 2005; Zonca et al. 2019), HEASoftpy (https://heasarc.gsfc.nasa.gov/lheasoft/heasoftpy/), climin (https://github.com/BRML/climin) (Bayer et al. 2015).

References

Appendix A: Uncertainties on logistic regression coefficients

Considering a galaxy cluster indexed by i among n simulated systems in the training set, we denote by η_i the vector of its m associated features. Simulations indicate whether this cluster is detected as extended (I_i = 1) or not (I_i = 0). Our implementation of logistic regression writes:

$$\begin{array}{c}\hfill p\left(I_i|{\mathit{\eta }}_i\right)=\phi (f({\mathit{\eta }}_i))=\phi \left({{\mathit{w}}^{\mathbf{\prime }}}^{⊺}{\mathit{\eta }}_i\right).\end{array}$$(A.1)

We used φ(t) = (1+e^−t)⁻¹ and w = (w₀,w₁,w₂,⋯,w_m) the vector of coefficients (w₀ is the intercept of the underlying linear model). We have defined ${{\mathit{\eta }}^{\mathbf{\prime }}}_i=(1,{\eta }_{i,1},{\eta }_{i,2},\cdots ,{\eta }_{i,m})$.

The fit is performed by minimising the following cost function (known as log-loss) over all possible values for w, where we have written q_i = p(I_i|η_i):

$$\begin{array}{c}\hfill \mathcal{C}(\mathit{w})={\displaystyle \sum _{i=1}^n}(-I_{i}ln(q_i)-(1-I_i)ln(1-q_i)).\end{array}$$(A.2)

A formula for the Hessian matrix H derives by noticing that:

$$\begin{array}{c}\hfill \partial q_i/\partial w_j=q_i(1-q_i){{\eta }^{\prime }}_{i,j}.\end{array}$$(A.3)

We obtain:

$$\begin{array}{c}\hfill H_{\mathit{jk}}=\frac{\partial \mathcal{C}}{\partial w_j\partial w_k}={\displaystyle \sum _{i=1}^n}{q}_i(1-q_i){{\eta }^{\prime }}_{i,j}{{\eta }^{\prime }}_{i,k}.\end{array}$$(A.4)

We obtain an estimate of the covariance matrix of best-fit coefficients w by inverting the (m + 1)×(m + 1) matrix H. This approximation is roughly valid as long as the best-fit model corresponds to the minimum of the cost function and in the limit of large n.

Appendix B: Logistic regression predicting detected clusters

We have shown results in Sect. 4 for a model trained to predict whether a cluster is selected (that is, detected and characterised as extended), given 7 × 5 surface brightness features (7 radial ranges and 5 components). We reiterate this exercise with a model trained to predict whether a cluster is detected. Figure B.1 shows the 35 coefficients associated to each of these features. As discussed in Sect. 4, these coefficients can be interpreted as the sensitivity of the detection to a small variation in a given feature. Fig. 4 and Fig. B.1 differ notably through the coefficients associated to the cluster emission (black crosses). The central 0 − 20″ surface brightness plays a major role in the detection model. Counts deposited at radii larger than 2 arcmin are relatively more important for detecting than for selecting a cluster.

Fig. B.1. Similar to Fig. 4, but for a logistic model p(IeRASS1|counts), trained to predict whether a cluster is detected among all simulated clusters.

Appendix C: Confidence range for proportions n/N

Some figures in the paper incorporate uncertainties on a certain ratio n/N with n and N both integer numbers, 0 ≤ n ≤ N. We estimated the 68% confidence range [L, U] by means of the score interval (Wilson 1927) with continuity correction (Vollset 1993):

$$\begin{array}{cc}& L=\frac{n_{-}+{\lambda }^2/2}{N+{\lambda }^2}-\lambda \frac{\sqrt{n_{-}-n_{-}^2/N+{\lambda }^2/4}}{N+{\lambda }^2},\hfill \end{array}$$(C.1)

$$\begin{array}{cc}& U=\frac{n_{+}+{\lambda }^2/2}{N+{\lambda }^2}+\lambda \frac{\sqrt{n_{+}-n_{+}^2/N+{\lambda }^2/4}}{N+{\lambda }^2},\hfill \end{array}$$(C.2)

where n₋ = n − 1/2 and n₊ = n + 1/2. The bounds are set to L = 0 if n = 0 and to U = 1 if n = N. The value for λ is the 1 − α/2 quantile of the normal distribution, in our case we set λ = 1 for α = 0.32 (hence, we show the 68% confidence intervals). Without continuity correction, we would have used n in place of n₋ and n₊ in those formulas.

These intervals should be adequate to represent uncertainties in the proportion n/N (Newcombe 1998), although the confidence interval may be too conservative in general (Agresti & Coull 1998). We emphasise that we do not use such uncertainties while computing our models, nor in performing calculations. They only serve for visualisation purposes.

Appendix D: Gaussian process classification

The selection function models presented in Sect. 5 rely on Gaussian process classification (GPC). To overcome the well-known computational bottleneck encountered in standard GPC we have applied the Stochastic Variational Gaussian Process (SVGP) algorithm (Hensman et al. 2015) implemented in the GPy library (GPy 2012). This section presents briefly the main characteristics of the method and some relevant references. For an application in the context of cluster astrophysics, see e.g. Debackere et al. (2022), in particular their Appendix B provides concise descriptions of the method for binomial likelihood (we use instead Bernoulli likelihood).

D.1. Defining a Gaussian process prior

We refer the reader to Rasmussen & Williams (2006) for a pedagogical introduction of GPC and for the notations followed in this Appendix. The principle of Gaussian process binary classification builds upon a latent variable f(x), function of the m features x = {x_i}_1⋯m and on the assumption of a Gaussian process (GP) prior on f. Values of the latent function are then mapped to probability values ranging between 0 and 1 by means of a bijective link function. In our specific case we choose the inverse probit for the link function, that is the cumulative Gaussian distribution Φ(f). Other sigmoid-shaped link functions would convene (e.g. the inverse logit).

By definition, any finite subset of random variables constituting a Gaussian process follows a joint multivariate normal distribution. The mean function μ(x) and the covariance function k suffice to fully specify the Gaussian process. We write:

$$\begin{array}{c}\hfill f\sim \mathcal{GP}(\mu (\mathit{x});k(\mathit{x},{\mathit{x}}^{\mathbf{\prime }})).\end{array}$$(D.1)

Using 𝔼 to denote the expectation value of a random variable, we have defined:

$$\begin{array}{cc}& \mu (\mathit{x})=\mathbb{E}[f(\mathit{x})],\hfill \end{array}$$(D.2)

$$\begin{array}{cc}& k(\mathit{x},{\mathit{x}}^{\mathbf{\prime }})=\mathbb{E}\left[(f(\mathit{x})-\mu (\mathit{x}))(f({\mathit{x}}^{\mathbf{\prime }})-\mu ({\mathit{x}}^{\mathbf{\prime }}))\right].\hfill \end{array}$$(D.3)

In our setup the mean function μ is constantly zero and the covariance function is the Radial Basis Function kernel (RBF), also known as stationary squared exponential kernel. This latter choice ensures extremely smooth properties of the underlying Gaussian process. One introduces m hyperparameters l_i entering the definition of the RBF kernel, defining the length-scales of variation along each feature dimension:

$$\begin{array}{c}\hfill k(\mathit{x},{\mathit{x}}^{\mathbf{\prime }})={\sigma }^{2}exp[-\frac{1}{2}{\displaystyle \sum _{i=1}^m}{\left(\frac{x_i-{x\prime }_i}{l_i}\right)}^2].\end{array}$$(D.4)

In this expression σ is an additional hyperparameter governing the amplitude of the correlation. The RBF correlation function is intuitively easy to understand: close pairs of points in the m-dimensional feature space are associated to strongly covariant random variables, the actual value of covariance being close to σ². Points far from each other are uncorrelated, with zero covariance. The length-scales govern the typical distance correlation cut-off along each dimension.

Given a GP prior and a set of n^* test points $X^{\ast }=({\mathit{x}}_1^{\ast },{\mathit{x}}_2^{\ast },\cdots ,{\mathit{x}}_{n^{\ast }}^{\ast }),$ one may generate a random vector f^* by sampling the multivariate normal distribution f^* ∼ 𝒩(0;K(X^*, X^*)). Here K(X^*, X^*) is a shorthand notation for the n^* × n^* covariance matrix whose coefficients are the $k({\mathit{x}}_p^{*},{\mathit{x}}_q^{*})$.

D.2. Predicting and training in standard GPC

In order to finely model our predictions, we wish to condition the GP on a set of known measurements (observations), also named training set in this paper. A set of n known input data points X = {x_i}_{i = 1⋯n} is associated to n binary observations I = {I_i}_{i = 1⋯n}, with I_i taking either value 0 or 1 (compliant with our notations in Table 1). The ensemble (X, I) constitutes the training set. We assume a Bernoulli likelihood to describe the data given the model values, that is:

$$\begin{array}{c}\hfill p\left(I_i|f_i\right)=\mathrm{\Phi }{(f_i)}^{I_i}{(1-\mathrm{\Phi }(f_i))}^{1-I_i}.\end{array}$$(D.5)

Given a new test point x^*, the model should provide a probabilistic prediction, that is interpreted as the cluster detection probability given our knowledge of the training dataset. Our desired model output writes:

$$\begin{array}{c}\hfill p\left(I^{\ast }|X,\mathit{I},{\mathit{x}}^{\ast }\right)={\displaystyle \int \mathrm{\Phi }}(f^{\ast })p\left(f^{\ast }|X,\mathit{I},{\mathit{x}}^{\ast }\right)\mathrm{d}{f}^{\ast }.\end{array}$$(D.6)

It is therefore required to obtain a distribution of the latent function f^* at the point x^*. We have that:

$$\begin{array}{c}\hfill p\left(f^{\ast }|X,\mathit{I},{\mathit{x}}^{\ast }\right)={\displaystyle \int p}(f^{\ast }|X,{\mathit{x}}^{\ast },\mathit{f})p(\mathit{f}|X,\mathit{I})\mathrm{d}\mathit{f}.\end{array}$$(D.7)

The first term in the integral derives from the GP prior directly. It can be expressed in a standard way by conditioning the joint Gaussian distribution p(f,f^*|X,x^*) on the latent function variables f:

$$\begin{array}{c}\hfill f^{\ast }|X,{\mathit{x}}^{\ast },\mathit{f}\sim \mathcal{N}({\mathit{k}}^{\ast ⊺}{K}^{-1}\mathit{f};k({\mathit{x}}^{\ast },{\mathit{x}}^{\ast })-{\mathit{k}}^{\ast ⊺}{K}^{-1}{\mathit{k}}^{\ast }).\end{array}$$(D.8)

We have defined K ≡ K(X, X) and k^* ≡ K(x^*, X), both are obtained through the RBF kernel, see Eq. D.4. If the number of training points n is large, the inversion of the n × n matrix K may be cumbersome and this represents a common numerical bottleneck in standard GP analysis. The term p(f|X,I) = p(I|f)p(f|X)/p(I|X) is the posterior over latent variables. It is generally not analytically tractable and often it must be approximated.

In GP classification, the marginal likelihood (or evidence) p(I|X) = ∫p(I|f,X)p(f|X)df plays an important role. The hyperparameters of the model (including those from the kernel) are optimally found by maximising the value of the marginal likelihood. Throughout this paper, we identify this operation with the ‘training phase’ of a GP classifier. In our case, an approximation is again needed as the marginal likelihood does not have a simple analytic expression.

D.3. Principles of the selected SVGP algorithm

Several works propose schemes to approximate the distributions and make the associated computations numerically tractable. We only briefly depict here the principles associated to the SVGP algorithm used in our modelling. Details and related references are found in Hensman et al. (2015).

It is useful to introduce a new set of hyperparameters made of k pairs of input data points and their associated output values (Z, u). The k × m coordinates of these points, called inducing inputs, will be optimised together with other hyperparameters during the training phase. Loosely speaking, these points will approximately ‘summarise’ the GP at locations that are most relevant. The second stage of approximations uses searching variational bounds (e.g. Blei et al. 2017, for a review), leading to the following inequality:

$$\begin{array}{c}\hfill lnp(\mathit{I}|X)\ge {\displaystyle \sum _{i=1}^n}{\mathbb{E}}_{q(f_i)}[lnp\left(I_i|f_i\right)]-\mathrm{KL}[q(\mathit{u})||p(\mathit{u})].\end{array}$$(D.9)

In this expression KL represents the Kullback-Leibler divergence, it measures the proximity between the two densities q(u) and p(u). The distribution q(u) belongs to a predefined family of distributions. Selecting a suitable family of distributions is at the core of variational inference, since it amounts to finding the most suitable q among this family. In practice, this procedure is equivalent to resolving an optimisation problem. In our case the selected family is constituted by the normal multivariate distributions, hence q(u)∼𝒩(m;S). This formulation requires the introduction of an additional set of hyperparameters indexing the members in the family: the vector m, with k coefficients, and the symmetric matrix S, with k(k + 1)/2 independent coefficients. The algorithm specifically uses lower Cholesky factorisation S = LL^⊤ for numerical stability reasons. The expectation under the sum sign is taken with respect to the marginals q(f_i) of the distribution q(f), defined such that q(f)≡∫p(f|u)q(u)du. The computation of p(f|u) follows directly from the GP prior assumption, in a similar way as in Eq. D.8 (replacing f with u and f^* with f).

In the training phase, the SVGP algorithm consists in finding the set of hyperparameters maximising the bound on the marginal likelihood, that is the right-hand side of Eq. D.9. Gradient descent algorithms are well-suited to such optimisation tasks. Our algorithm uses stochastic optimisation techniques, by drawing random mini-batches of the input data in a large number of successive steps.

The interest of the method depicted here consists in a reduced numerical complexity, as compared to standard GPC. Ensuring that the number of inducing points k is much smaller than the training set (i.e. k ≪ n) puts much of the cost on the calculation of the expectation ${\mathbb{E}}_{q(f_i)}[\cdots ]$ and of the gradients with respect to the hyperparameters. In any case, there is no requirement for inverting of a n × n matrix any longer. Model predictions are also numerically efficient, thanks to the approximation of the posterior p(f,u|X,I) → p(f|u)q(u). In contrast to the standard GPC method (Eq. D.7), computing a prediction does not need knowledge of the values of the training set (X, I); only the inducing points (along with all other hyperparameters) are required. Saving a trained model on disk is thus very memory efficient.

It is important to remind here that several approximations lead to the method outlined above, which results from a compromise between efficiency and accuracy. In particular, Blei et al. (2017) discuss open problems related to variational inference and the use of KL divergence to set a lower bound to the evidence. For our own application we rely on a posteriori validation of the models (e.g. Sect. 5.4). In general, the expectation values of the GPC models reflect well the cluster detection rates, but we find the associated variance to be slightly underestimated.

D.4. Setup of our GPC model construction

Let us now relate the previous notations to the setup of our numerical experiments, as described in Sect. 5. The number of training points is the number of simulated clusters in the training set⁵, n = 5.7 × 10⁵. The number of features may take values ranging from m = 2 (e.g. in the case of models based on M₅₀₀ and z) to m = 6 (e.g. for models with L_X, z, EM₀, N_H, T_exp and bkg). We choose a small number of inducing points k = 30, thus reaching an empirical trade-off between representativeness of the input data space, computational efficiency and smoothness of the model. The baseline model for cosmological studies has m = 5, this leads to optimising 651 hyperparameters. This number decomposes into: 6 for the kernel (σ and the five length-scales), 30 × 5 coordinates of the inducing points Z, 30 for the values of the m parameter and 30 × 31/2 for the coefficients of the L matrix. Before training, each feature in the training set X is rescaled to the [0, 1] interval using simple functions (e.g. linear or logarithmic), selected according to the dynamic range and distribution of each feature.

Given the large number of free hyperparameters, we used the ADADELTA algorithm (Zeiler 2012) as implemented in the climin library. Our training set is extremely unbalanced, with only a few percent of the simulated clusters listed as detected or selected. In this context, we have found that stochastic optimisation does provide satisfactory results only if mini-batches are of large enough size. In our case we choose random mini-batches containing 2¹⁶ = 65536 simulated clusters, among them a few hundreds are detected systems. Therefore only about 7 different mini-batches are available. We secured a sufficiently high number of iterations (typically of order 10⁴) before stopping the training phase. In the course of our model production, we have noticed some models converge to undesired solutions, as apparent from unusual best-fit length-scales. This behaviour was in many cases due to unlucky draws of the first mini-batches, containing an uncommon fraction of detected systems. In such failure cases, launching a new training with a fresh random seed would generally provide a satisfactory solution.

The training phase of a classifier took typically a few hours on a standard multi-core machine. Once finalised, each model went through the list of diagnostic tests listed in Sect. 5.4. Each satisfactory model was shared in the eRASS1 collaboration after dumping onto disk the values of the hyperparameters and of the factors rescaling the features to the [0, 1] interval. Prediction of the probabilistic values p(I^*|x^*) is a very quick operation, even for large numbers of test points x^* (Eq. D.6). This is in fact a requirement of the cosmological modelling (Ghirardini et al. 2024, and Eq. 2) to obtain computationally efficient response of the selection function models. A crude estimate of the uncertainties on the model outputs may be obtained by folding the ±1σ enveloppe of the latent function posterior (given by the variance in Eq. D.8) into the link function Φ.

Appendix E: Additional selection models with intermediate variables

This section presents additional models discussed throughout this paper, constructed from intermediate variables as discussed in Sect. 5.

E.1. Absorbed flux and redshift

The absorbed R₅₀₀ flux in the 0.5–2 keV band is an alternative to the unabsorbed count-rate. It has units erg s⁻¹ cm⁻² and it is available for all simulated clusters from the twin simulation. Fig. E.1 shows slices through this model p(I_main|f_X,z,N_H,T_exp,bkg). It has very similar behaviour as a selection expressed with CR (Fig. 10).

Fig. E.1. Representation of a model p(Imain|fX,z,NH,Texp,bkg) for fixed values of cluster redshift z = 0.2 and z = 0.5, local galactic absorption NH, a nominal background level, two values of exposure times Texp (units s) and a range of flux values (x-axis). The shaded regions indicate the approximate 68% confidence range output of the model.

E.2. Integrating morphological features

Closer inspection of twin simulations reveals that cluster flux (or count rate, or luminosity) is not the only determinant of the detection probability. Also its extent, characterised by R₅₀₀ and its central emission measure, characterised by EM₀, play a secondary role. We have (see also Comparat et al. 2020):

$$\begin{array}{c}\hfill EM(R)=-{log}_{10}{\displaystyle {\int }_{-\infty }^{+\infty }}{n}_e{n}_p\mathrm{d}l.\end{array}$$(E.1)

The integration is along the line of sight (coordinate l, expressed in units Mpc) at distance R from the cluster centre, n_e and n_p are the electron and proton densities in the galaxy cluster (expressed in units cm⁻³). For numerical stability reasons, the central emission measure EM₀ is defined as the average of EM(R) in the radial range 0 < R < 0.025 R₅₀₀ and the integration bounds are set to ±20 R₅₀₀ instead of infinity.

Morphological features intervene non-trivially in the selection function. At fixed R₅₀₀ luminosity, a larger cluster size R₅₀₀ dilutes the photons over a larger area, thus making detection of faint objects increasingly difficult for eSASS; very bright extended objects are affected differently, since detection is split into a number of smaller sources, not necessarily classified as extended by eSASS and requiring thorough cleaning (Bulbul et al. 2024). Similarly, low values of EM₀ indicate peaked surface brightness profiles, decreasing the probability of selecting a cluster as extended. Large values of EM₀ indicate very flat cores, and thus difficulties in detecting such faint objects against the background. It is clear however that these trends depend also on the background level, cluster flux, redshift, exposure time. The GPC formalism comes particularly handy in constructing a model without imposing strong priors on the exact role of these parameters. Fig. E.2 shows four slices through this multi-dimensional model, as a function of the parameter EM₀ and for two different cluster luminosities, two different exposure times.

Fig. E.2. Similar to Fig. E.1, but representing the outcome of a model p(Imain|LX,z,EM0,NH,Texp,bkg) for fixed values of cluster redshift z, two different values for the luminosity LX, and local galactic absorption NH, in addition to a nominal background level and two values of exposure time Texp. The shaded regions indicate the approximate 68% confidence range output of the model.

Also in Fig. E.3 we show the dependence on R₅₀₀ at fixed flux. The larger the apparent cluster extent (here the projected R₅₀₀ in units arcmin), the lower the probability of selecting the cluster in the sample. This is a consequence of diluting the photons emitted by the cluster over a larger area, hence decreasing the average surface brightness and therefore the ability to distinguish extended emission from background fluctuations.

Fig. E.3. Similar to Fig. E.1, but representing the outcome of a model p(Imain|fX,R500,NH,Texp,bkg) for two different values for the flux fX, fixed local galactic absorption NH, a nominal background level, and two values of exposure time Texp. The shaded regions indicate the approximate 68% confidence range output of the model.

Appendix F: Model for the SPT-SZ selection

The SPTpol-100d selection function depends on one single parameter ζ_SPTpol. We write:

$$\begin{array}{c}\hfill p\left(I_{\mathrm{SPTpol}}|{\zeta }_{\mathrm{SPTpol}}\right)={\displaystyle {\int }_{-\infty }^{u(\zeta )}}{\left(2\pi \right)}^{-1/2}{e}^{-x^2/2}\mathrm{d}x.\end{array}$$(F.1)

The upper integration bound is a function of ζ_SPTpol, namely $u({\zeta }_{\text{SPTpol}})=\left[{\left(3+{\zeta }_{\text{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 ξ > 4.5. The parameter ξ is biased and normally distributed with unit standard deviation (Huang et al. 2020):

$$\begin{array}{c}\hfill p\left(\xi |{\zeta }_{\mathrm{SPTpol}}\right)\sim \mathcal{N}(\sqrt{{{\zeta }_{\mathrm{SPTpol}}}^2+3};1).\end{array}$$(F.2)

The expression for the SPT-ECS selection function writes in a similar manner (Bleem et al. 2020). There, ζ_SPTpol is replaced by ζ_ECS, both quantities scaling with mass as power-laws, albeit with a different normalisation.

Appendix G: Unexpected SPT-ECS detections

The analysis presented in Sect. 6.8 has revealed three eRASS1 clusters having too low an X-ray flux for them to appear in the SPT-ECS sample. Our model indeed predicts a probability value P_i = p(I_{SPT − ECS}|I_cosmo,F500_0520,z_λ) < 0.025 for each of these entries:

• 1eRASS J021731.8-320002 (Fig. G.1) is listed at z_λ = 0.352 ± 0.009 with flux 2.1 ± 0.4 × 10⁻¹⁴ erg s⁻¹ cm⁻² and extent likelihood 14.4. The estimated X-ray mass is M_500, X ≃ 3 × 10¹⁴ M_⊙. The matched SPT-ECS detection is at z_SPT = 0.354 ± 0.009 (from Magellan/FourStar) and ξ = 5.1. Its estimated SZ mass is M_500, SPT ≃ 4 × 10¹⁴ M_⊙.

• 1eRASS J034733.3-333351 (Fig. G.2) is listed at z_λ = 0.465 ± 0.009 with flux $8_{-3}^{+10}\times {10}^{-14}$ erg s⁻¹ cm⁻² and extent likelihood 7.5. The estimated X-ray mass is M_500, X ≃ 2 × 10¹⁴ M_⊙. The matched SPT-ECS detection is at z_SPT = 0.45 ± 0.01 (photometric, from Magellan/FourStar) and ξ = 5.4. Its estimated mass is M_500, SPT ≃ 4 × 10¹⁴ M_⊙.

• 1RASS J045420.6-373708 (Fig. G.3) is listed at z_λ = 0.513 ± 0.008 without any flux measurement, nor an upper limit and extent likelihood 6.3. The corresponding SPT detection is at z_SPT = 0.509 ± 0.006 (photometric) and significance ξ = 5.9, with estimated mass M_500, SPT ≃ 4.8 × 10¹⁴ M_⊙.

Fig. G.1. Cluster 1eRASS J021731.8-320002. The image shows an optical cutout (DR10 Legacy survey) and eRASS1 X-ray contours overlaid.

Fig. G.2. Cluster 1eRASS J034733.3-333351.

Fig. G.3. Cluster 1eRASS J045420.6-373708.

All Tables

All Figures

Fig. 1. Comparing the outcome of two matching procedures relating simulated halos (plain black histogram) to detected sources (photon-based matching is a dashed blue line, and position-based matching is a thin orange line). A simulated halo is flagged as selected (Imain = 1) whenever it is (solidly) matched to an extended source in the detection catalogue. The details of the matching algorithms impact the flux distribution of detected sources, especially at flux below a few 10−14 erg s−1 cm−2, with deviations up to a factor 2 in certain bins.

In the text

Fig. 2. Example case of blending, where a bright AGN drives the detection and the presence of a faint extended source favours 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 ∼5 × 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 (Imain = 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.

In the text

Fig. 3. Distribution of clusters in the training set as a function of Σ (units counts arcsec−2), and the average surface brightness in the 90″ radius around their centre. The x-axis is rescaled with c = 2 × 10−4 arcsec−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 Appendix C. The vertical dashed line indicates the transition $\widehat{\eta }=1.94$ in the logistic model described by Eqs. (3) and (4) – corresponding to N ≃ 30 counts.

In the text

Fig. 4. Representation of the 35 coefficients wj of a logistic regression model p(Imain|counts), trained to predict whether a cluster was 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)).

In the text

Fig. 5. Performance of three different logistic models p(Imain|counts) 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″. The dashed grey line stands for a model using as input only the 90″ 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).

In the text

Fig. 6. Reliablity of two different logistic models p(Imain|counts) 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 (Appendix 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 × 5 parameters as input, grey dots and dashed lines for the model using only cluster counts within 90″ to make its prediction.

In the text

Fig. 7. Representation of the model p(Imain|M500c, z) predicting a cluster to be detected and selected with an extent likelihood above three. The explanatory variables (features) are labels attached to simulated clusters in the twin simulation, standing for galaxy cluster M500c 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 M500c = 3 × 1014 M⊙ (indicated with dotted lines).

In the text

Fig. 8. Representation of the model p(Imain|LX, z) predicting a cluster to be detected and selected with an extent likelihood above three. The explanatory variables (features) are labels attached to simulated clusters in the twin simulation, standing for galaxy cluster luminosity LX measured in the 0.5–2 keV energy band at the cluster rest-frame (units erg s−1, rescaled in log10) and for cosmological redshift z. Other details are similar as in Fig. 7.

In the text

Fig. 9. Representation of the models p(Imain|LX, z, NH, Texp, bkg) and p(Icosmo|LX, z, NH, Texp, bkg) for fixed values of cluster redshift z, two different values for the luminosity LX, and local galactic absorption NH, in addition to a nominal background level and a range of exposure times Texp (along the x-axis). The shaded regions indicate the approximate 68% confidence range output of the model.

In the text

	Fig. 10. Representation of a model p(Icosmo|CR,z,NH,Texp,bkg) for two values of cluster redshift z = 0.2 and z = 0.5, fixed local galactic absorption NH, a nominal background level, two values of exposure times Texp (units s), and a range of count-rate values (x-axis). The shaded regions indicate the approximate 68% confidence range output of the model.
In the text

	Fig. 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 the legend. A model involving luminosity, redshift, morphological features and local sky information performs best (yellow dot-dashed curve).
In the text

	Fig. 12. Similar to Fig. 11, but 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. (2024).
In the text

Fig. 13. Evaluating the reliability of five models predicting the probability of a cluster to appear in the eRASS1 primary cluster catalogue (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.

In the text

	Fig. 14. Similar to Fig. 13, but comparing three models predicting the presence of a cluster in the cosmology sample.
In the text

Fig. 15. Testing the models with eRASS1-cosmo as a test catalogue, and eFEDS as a reference catalogue. 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 (zeFEDS). 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 Pi (Eq. (8)) computed for each eFEDS cluster.

In the text

Fig. 16. Testing the models with eRASS1-cosmo as a test catalogue, and eFEDS as a reference catalogue. The 531 eFEDS entries are grouped in bins of the model output probabilities Pi (horizontal axis). In each bin, the actual fraction of matches in the eRASS1 catalogue is reported on the vertical axis. Error bars represent the estimated 68% confidence range (Appendix 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.

In the text

Fig. 17. Similar to Fig. 15, but 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 R500 reported in the eRASS1 catalogue (units 10−14 erg s−1 cm−2). The horizontal axis is the catalogue redshift zλ. Clusters falling in zones of zero eFEDS exposure have vanishing probability (Pi = 0). Other eRASS1 clusters are very likely to be detected in eFEDS (Pi ≃ 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 × 10−13 erg s−1 cm−2.

In the text

	Fig. 18. Similar to Fig. 16, but using the primary sample of clusters eRASS1-main as a test sample, instead of the eRASS1 cosmology sample.
In the text

Fig. 19. Posterior distribution for the pair of parameters (δ1 and δ2) introduced in Eq. (9), for the experiment using eRASS1-main as a test and eFEDS as a reference sample. The histograms represent the marginalised posterior distributions, with mean and 68% range overprinted. The contours represent the equivalent 0.5-, 1-, 1.5- and 2-σ distribution of the joint posterior. Clearly the parameter δ2 departs from zero, indicating that a fraction of the systems we predict not be detected leak into the actual eRASS1 sample.

In the text

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

In the text

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

	Fig. 22. Similar to Fig. 20, but using 163 clusters in SPTpol-ECS catalogue as a reference sample and the eRASS1 cosmology sample as a test.
In the text

	Fig. 23. Similar to Fig. 21, but obtained by using the SPT-ECS sample as a reference (163 entries) and using the cosmology sample of eRASS1 clusters as a test.
In the text

	Fig. 24. Similar to Fig. 17, but using 958 eRASS1 clusters from the cosmology sample as a reference, and the SPTpol-ECS catalogue 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.
In the text

	Fig. 25. Similar to Fig. 16, but obtained by using the eRASS1 cosmology sample (958 entries) as a reference and using SPTpol-ECS as a test sample.
In the text

Fig. 26. Relation between the latent variable f(η) = w0 + ∑jwjηj calculated for each simulated cluster in the test sample after performing logistic regression with 35 surface brightness features, and the values of ℒdet and ℒ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 ℒdet value (red histogram in the lower panel). Clusters detected but not categorised as extended have zero value for ℒext (orange histogram). The apparent correlation between f and both measurements indicate that the model has “discovered” the importance of ℒdet and ℒext in the selection process.

In the text

	Fig. 27. Similar to Fig. 26, but replacing the x-axis by the mean of the latent variable implicit in the Gaussian process classifier p(Icosmo|CR,z,ℋ).
In the text

Fig. 28. Value of the combined Bernoulli negative log-likelihood (Eq. (18)) in the comparison of the SPT-ECS sample (reference) and the eRASS1 cosmology sample (test) for various values of the correlation coefficient ρ. This coefficient expresses the level of covariance between lnLX and lnζ at fixed halo mass and redshift. Assuming all other model parameters are correct, the constraint obtained from the catalogue comparison exercise depicted in this paper is relatively loose: ρ ≃ 0.2 ± 0.6 (68% confidence level) as deduced from $(B-B_{\text{min}})~{\chi }_1^2$.

In the text

	Fig. B.1. Similar to Fig. 4, but for a logistic model p(IeRASS1|counts), trained to predict whether a cluster is detected among all simulated clusters.
In the text

	Fig. E.1. Representation of a model p(Imain|fX,z,NH,Texp,bkg) for fixed values of cluster redshift z = 0.2 and z = 0.5, local galactic absorption NH, a nominal background level, two values of exposure times Texp (units s) and a range of flux values (x-axis). The shaded regions indicate the approximate 68% confidence range output of the model.
In the text

	Fig. E.2. Similar to Fig. E.1, but representing the outcome of a model p(Imain|LX,z,EM0,NH,Texp,bkg) for fixed values of cluster redshift z, two different values for the luminosity LX, and local galactic absorption NH, in addition to a nominal background level and two values of exposure time Texp. The shaded regions indicate the approximate 68% confidence range output of the model.
In the text

	Fig. E.3. Similar to Fig. E.1, but representing the outcome of a model p(Imain|fX,R500,NH,Texp,bkg) for two different values for the flux fX, fixed local galactic absorption NH, a nominal background level, and two values of exposure time Texp. The shaded regions indicate the approximate 68% confidence range output of the model.
In the text

	Fig. G.1. Cluster 1eRASS J021731.8-320002. The image shows an optical cutout (DR10 Legacy survey) and eRASS1 X-ray contours overlaid.
In the text

	Fig. G.2. Cluster 1eRASS J034733.3-333351.
In the text

	Fig. G.3. Cluster 1eRASS J045420.6-373708.
In the text