* Dietrich et al. (2019) Dietrich, J. P., Bocquet, S., Schrabback, T., et al. 2019, MNRAS, 483, 2871
  * Fernique et al. (2014) Fernique, P., Boch, T., Donaldson, T., et al. 2014, MOC - HEALPix Multi-Order Coverage map Version 1.0, IVOA Recommendation 02 June 2014
  * Finoguenov et al. (2020) Finoguenov, A., Rykoff, E., Clerc, N., et al. 2020, A&A, 638, A114
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Gallagher & Smeenk (2023) Gallagher, S. C. & Smeenk, C. 2023, What’s in a Survey? Simulation-Induced Selection Effects in Astronomy, ed. N. Mills Boyd, S. De Baerdemaeker, K. Heng, & V. Matarese (Cham: Springer International Publishing), 207–222
  * Garrel et al. (2022) Garrel, C., Pierre, M., Valageas, P., et al. 2022, A&A, 663, A3
  * Ghirardini et al. (submitted) Ghirardini, V., Bulbul, E., Artis, E., et al. submitted, A&A
  * Ghirardini et al. (2019) Ghirardini, V., Eckert, D., Ettori, S., et al. 2019, A&A, 621, A41
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759
  * GPy (since 2012) GPy. since 2012, GPy: A Gaussian process framework in python, http://github.com/SheffieldML/GPy
  * Grandis et al. (2020) Grandis, S., Klein, M., Mohr, J. J., et al. 2020, MNRAS, 498, 771
  * Hensman et al. (2015) Hensman, J., Matthews, A., & Ghahramani, Z. 2015, in Proceedings of Machine Learning Research, Vol. 38, Proceedings of the Eighteenth International Conference on Artificial Intelligence and Statistics, ed. G. Lebanon & S. V. N. Vishwanathan (San Diego, California, USA: PMLR), 351–360
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116
  * Hoyle et al. (2012) Hoyle, B., Jimenez, R., Verde, L., & Hotchkiss, S. 2012, J. Cosmology Astropart. Phys., 2012, 009
  * Hu & Kravtsov (2003) Hu, W. & Kravtsov, A. V. 2003, ApJ, 584, 702
  * Huang et al. (2020) Huang, N., Bleem, L. E., Stalder, B., et al. 2020, AJ, 159, 110
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Jasche & Lavaux (2019) Jasche, J. & Lavaux, G. 2019, A&A, 625, A64
  * Kluge et al. (submitted) Kluge, M., Comparat, J., Liu, A., et al. submitted, A&A
  * Kong et al. (2020) Kong, H., Burleigh, K. J., Ross, A., et al. 2020, MNRAS, 499, 3943
  * Kostić et al. (2023) Kostić, A., Nguyen, N.-M., Schmidt, F., & Reinecke, M. 2023, J. Cosmology Astropart. Phys., 2023, 063
  * Lindholm et al. (2021) Lindholm, V., Finoguenov, A., Comparat, J., et al. 2021, A&A, 646, A8
  * Liu et al. (2022) Liu, A., Bulbul, E., Ghirardini, V., et al. 2022, A&A, 661, A2
  * Liu et al. (2023) Liu, A., Bulbul, E., Ramos-Ceja, M. E., et al. 2023, A&A, 670, A96
  * Mantz (2019) Mantz, A. B. 2019, MNRAS, 485, 4863
  * Mantz et al. (2015) Mantz, A. B., von der Linden, A., Allen, S. W., et al. 2015, MNRAS, 446, 2205
  * Marulli et al. (2018) Marulli, F., Veropalumbo, A., Sereno, M., et al. 2018, A&A, 620, A1
  * Matthews (1975) Matthews, B. 1975, Biochimica et Biophysica Acta (BBA) - Protein Structure, 405, 442
  * Merloni et al. (2024) Merloni, A., Lamer, G., Liu, T., et al. 2024, A&A, 682, A34
  * Newcombe (1998) Newcombe, R. G. 1998, Statistics in Medicine, 17, 857
  * Pacaud et al. (2006) Pacaud, F., Pierre, M., Refregier, A., et al. 2006, MNRAS, 372, 578
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13
  * Predehl et al. (2021) Predehl, P., Andritschke, R., Arefiev, V., et al. 2021, A&A, 647, A1
  * Rasmussen & Williams (2006) Rasmussen, C. E. & Williams, C. K. I. 2006, Gaussian Processes for Machine Learning
  * Rix et al. (2021) Rix, H.-W., Hogg, D. W., Boubert, D., et al. 2021, AJ, 162, 142
  * Salvato et al. (2018) Salvato, M., Buchner, J., Budavári, T., et al. 2018, MNRAS, 473, 4937
  * Schmidt (2021a) Schmidt, F. 2021a, J. Cosmology Astropart. Phys., 2021, 033
  * Schmidt (2021b) Schmidt, F. 2021b, J. Cosmology Astropart. Phys., 2021, 032
  * Seppi et al. (2022) Seppi, R., Comparat, J., Bulbul, E., et al. 2022, A&A, 665, A78
  * Suchyta et al. (2016) Suchyta, E., Huff, E. M., Aleksić, J., et al. 2016, MNRAS, 457, 786
  * Tinker et al. (2008) Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709
  * Turner et al. (2022) Turner, D. J., Giles, P. A., Romer, A. K., et al. 2022, MNRAS, 517, 657
  * Vanderlinde et al. (2010) Vanderlinde, K., Crawford, T. M., de Haan, T., et al. 2010, ApJ, 722, 1180
  * Veronica et al. (2023) Veronica, A., Reiprich, T. H., Pacaud, F., et al. 2023, arXiv e-prints, arXiv:2311.07488
  * Vollset (1993) Vollset, S. E. 1993, Statistics in medicine, 12, 809
  * Wilson (1927) Wilson, E. B. 1927, Journal of the American Statistical Association, 22, 209
  * Zeiler (2012) Zeiler, M. D. 2012 [arXiv:1212.5701]
  * Zonca et al. (2019) Zonca, A., Singer, L., Lenz, D., et al. 2019, Journal of Open Source Software, 4, 1298

##  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 $\@vec{\eta}_{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:

$p\left(I_{i}\mathrel{}\middle|\mathrel{}\@vec{\eta}_{i}\right)=\varphi\left(f(\@vec{\eta}_{i})\right)=\varphi\left(\@vec{w}^{\prime\intercal}\@vec{\eta}_{i}\right)$
(21)

We used $\varphi(t)=\left(1+e^{-t}\right)^{-1}$ and
$\@vec{w}=(w_{0},w_{1},w_{2},\dotsc,w_{m})$ the vector of coefficients
($w_{0}$ is the intercept of the underlying linear model). We have defined
$\@vec{\eta}_{i}^{\prime}=(1,\eta_{i,1},\eta_{i,2},\dotsc,\eta_{i,m})$.

The fit is performed by minimizing the following cost function (known as log-
loss) over all possible values for $\@vec{w}$, where we have written
$q_{i}=p\left(I_{i}\mathrel{}\middle|\mathrel{}\@vec{\eta}_{i}\right)$:

$\mathcal{C}(\@vec{w})=\sum_{i=1}^{n}\left(-I_{i}\ln(q_{i})-(1-I_{i})\ln(1-q_{i})\right)$
(22)

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

$\partial q_{i}/\partial w_{j}=q_{i}(1-q_{i})\eta^{\prime}_{i,j}$ (23)

We obtain:

$H_{jk}=\frac{\partial\mathcal{C}}{\partial w_{j}\partial
w_{k}}=\sum_{i=1}^{n}q_{i}(1-q_{i})\eta^{\prime}_{i,j}\eta^{\prime}_{i,k}$
(24)

We obtain an estimate of the covariance matrix of best-fit coefficients
$\@vec{w}$ by inverting the $(m+1)\times(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 characterized as extended), given
$7\times 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 32 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. 32 differ notably through the coefficients associated to the
cluster emission (black crosses). The central $0-20\arcsec$ 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.

Figure 32: Similar figure as Fig. 4, but for a logistic model
$p\left(I_{\rm{eRASS1}}\mathrel{}\middle|\mathrel{}\rm counts\right)$, 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\leq n\leq N$. We estimated the 68%
confidence range $[L,U]$ by means of the score interval (Wilson 1927) with
continuity correction (Vollset 1993):

$\displaystyle L$ $\displaystyle=$
$\displaystyle\frac{n_{-}+\lambda^{2}/2}{N+\lambda^{2}}-\lambda\frac{\sqrt{n_{-}-n_{-}^{2}/N+\lambda^{2}/4}}{N+\lambda^{2}}$
(25) $\displaystyle U$ $\displaystyle=$
$\displaystyle\frac{n_{+}+\lambda^{2}/2}{N+\lambda^{2}}+\lambda\frac{\sqrt{n_{+}-n_{+}^{2}/N+\lambda^{2}/4}}{N+\lambda^{2}}$
(26)

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 $\lambda$ is the $1-\alpha/2$ quantile of
the normal distribution, in our case we set $\lambda=1$ for $\alpha=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 emphasize 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 since 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(\@vec{x})$, function of the $m$ features
$\@vec{x}=\\{x_{i}\\}_{1\dotsc 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 $\Phi(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
$\mu(\@vec{x})$ and the covariance function $k$ suffice to fully specify the
Gaussian Process. We write:

$f\sim\mathcal{GP}\left(\mu(\@vec{x})\mathrel{};\mathrel{}k(\@vec{x},\@vec{x^{\prime}})\right)$
(27)

Using $\mathbb{E}$ to denote the expectation value of a random variable, we
have defined:

$\mu(\@vec{x})=\mathbb{E}\left[f(\@vec{x})\right]$ (28)
$k(\@vec{x},\@vec{x^{\prime}})=\mathbb{E}\left[\left(f(\@vec{x})-\mu(\@vec{x})\right)\left(f(\@vec{x^{\prime}})-\mu(\@vec{x^{\prime}})\right)\right]$
(29)

In our setup the mean function $\mu$ 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:

$k(\@vec{x},\@vec{x^{\prime}})=\sigma^{2}\exp\left[-\frac{1}{2}\sum_{i=1}^{m}\left(\frac{x_{i}-x^{\prime}_{i}}{l_{i}}\right)^{2}\right]$
(30)

In this expression $\sigma$ 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 $\sigma^{2}$. 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^{*}=(\@vec{x}^{*}_{1},\@vec{x}^{*}_{2},\dotsc,\@vec{x}^{*}_{n^{*}})$ one
may generate a random vector $\@vec{f^{*}}$ by sampling the multivariate
normal distribution
$\@vec{f^{*}}\sim\mathcal{N}\left(\@vec{0}\mathrel{};\mathrel{}K(X^{*},X^{*})\right)$.
Here $K(X^{*},X^{*})$ is a shorthand notation for the $n^{*}\times n^{*}$
covariance matrix whose coefficients are the
$k(\@vec{x}^{*}_{p},\@vec{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=\\{\@vec{x}_{i}\\}_{i=1\dotsc n}$ is
associated to $n$ binary observations $\@vec{I}=\\{I_{i}\\}_{i=1\dotsc n}$,
with $I_{i}$ taking either value 0 or 1 (compliant with our notations in Table
1). The ensemble $(X,\@vec{I})$ constitutes the training set. We assume a
Bernoulli likelihood to describe the data given the model values, i.e.:

$p\left(I_{i}\mathrel{}\middle|\mathrel{}f_{i}\right)=\Phi(f_{i})^{I_{i}}\left(1-\Phi(f_{i})\right)^{1-I_{i}}$
(31)

Given a new test point $\@vec{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:

$p\left(I^{*}\mathrel{}\middle|\mathrel{}X,\@vec{I},\@vec{x}^{*}\right)=\int\Phi(f^{*})p\left(f^{*}\mathrel{}\middle|\mathrel{}X,\@vec{I},\@vec{x}^{*}\right)\mathrm{d}f^{*}$
(32)

It is therefore required to obtain a distribution of the latent function
$f^{*}$ at the point $\@vec{x}^{*}$. We have that:

$p\left(f^{*}\mathrel{}\middle|\mathrel{}X,\@vec{I},\@vec{x}^{*}\right)=\int
p\left(f^{*}\mathrel{}\middle|\mathrel{}X,\@vec{x}^{*},\@vec{f}\right)p\left(\@vec{f}\mathrel{}\middle|\mathrel{}X,\@vec{I}\right)\mathrm{d}\@vec{f}$
(33)

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\left(\@vec{f},f^{*}\mathrel{}\middle|\mathrel{}X,\@vec{x}^{*}\right)$ on
the latent function variables $\@vec{f}$:

$f^{*}\mathrel{}|\mathrel{}X,\@vec{x}^{*},\@vec{f}\sim\mathcal{N}\left(\@vec{k}^{*\intercal}K^{-1}\@vec{f}\mathrel{};\mathrel{}k(\@vec{x}^{*},\@vec{x}^{*})-\@vec{k}^{*\intercal}K^{-1}\@vec{k}^{*}\right)$
(34)

We have defined $K\equiv K(X,X)$ and $\@vec{k}^{*}\equiv K(\@vec{x}^{*},X)$,
both are obtained through the RBF kernel, see Eq. 30. If the number of
training points $n$ is large, the inversion of the $n\times n$ matrix $K$ may
be cumbersome and this represents a common numerical bottleneck in standard GP
analysis. The term
$p\left(\@vec{f}\mathrel{}\middle|\mathrel{}X,\@vec{I}\right)=p\left(\@vec{I}\mathrel{}\middle|\mathrel{}\@vec{f}\right)p\left(\@vec{f}\mathrel{}\middle|\mathrel{}X\right)/p\left(\@vec{I}\mathrel{}\middle|\mathrel{}X\right)$
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\left(\@vec{I}\mathrel{}\middle|\mathrel{}X\right)=\int
p\left(\@vec{I}\mathrel{}\middle|\mathrel{}\@vec{f},X\right)p\left(\@vec{f}\mathrel{}\middle|\mathrel{}X\right)\mathrm{d}\@vec{f}$
plays an important role. The hyperparameters of the model (including those
from the kernel) are optimally found by maximizing 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 $(\@vec{Z},\@vec{u})$.
The $k\times m$ coordinates of these points, called inducing inputs, will be
optimized together with other hyperparameters during the training phase.
Loosely speaking, these points will approximately ‘summarize’ 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:

$\ln
p\left(\@vec{I}\mathrel{}\middle|\mathrel{}X\right)\geq\sum_{i=1}^{n}\mathbb{E}_{q(f_{i})}\left[\ln
p\left(I_{i}\mathrel{}\middle|\mathrel{}f_{i}\right)\right]-{\rm
KL}\left[q(\@vec{u})\mathrel{}\middle|\middle|\mathrel{}p(\@vec{u})\right]$
(35)

In this expression ${\rm KL}$ represents the Kullback-Leibler divergence, it
measures the proximity between the two densities $q(\@vec{u})$ and
$p(\@vec{u})$. The distribution $q(\@vec{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
optimization problem. In our case the selected family is constituted by the
normal multivariate distributions, hence
$q(\@vec{u})\sim\mathcal{N}\left(\@vec{m}\mathrel{};\mathrel{}S\right)$. This
formulation requires the introduction of an additional set of hyperparameters
indexing the members in the family: the vector $\@vec{m}$, with $k$
coefficients, and the symmetric matrix $S$, with $k(k+1)/2$ independent
coefficients. The algorithm specifically uses lower Cholesky factorization
$S=LL^{\top}$ for numerical stability reasons. The expectation under the sum
sign is taken with respect to the marginals $q(f_{i})$ of the distribution
$q(\@vec{f})$, defined such that $q(\@vec{f})\equiv\int
p\left(\@vec{f}\mathrel{}\middle|\mathrel{}\@vec{u}\right)q(\@vec{u})\mathrm{d}\@vec{u}$.
The computation of
$p\left(\@vec{f}\mathrel{}\middle|\mathrel{}\@vec{u}\right)$ follows directly
from the GP prior assumption, in a similar way as in Eq. 34 (replacing
$\@vec{f}$ with $\@vec{u}$ and $f^{*}$ with $\@vec{f}$).

In the training phase, the SVGP algorithm consists in finding the set of
hyperparameters maximizing the bound on the marginal likelihood, that is the
right-hand side of Eq. 35. Gradient descent algorithms are well-suited to such
optimization tasks. Our algorithm uses stochastic optimization 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\ll n$) puts much of
the cost on the calculation of the expectation
$\mathbb{E}_{q(f_{i})}\left[\dotsc\right]$ and of the gradients with respect
to the hyperparameters. In any case, there is no requirement for inverting of
a $n\times n$ matrix any longer. Model predictions are also numerically
efficient, thanks to the approximation of the posterior
$p\left(\@vec{f},\@vec{u}\mathrel{}\middle|\mathrel{}X,\@vec{I}\right)\rightarrow
p\left(\@vec{f}\mathrel{}\middle|\mathrel{}\@vec{u}\right)q(\@vec{u})$. In
contrast to the standard GPC method (Eq. 33), computing a prediction does not
need knowledge of the values of the training set $(X,\@vec{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 set121212This number reduces to
$n=4.8\times 10^{5}$ when training models for the cosmology sample, as we
restrict simulated clusters to the range $0.05<z<0.85$., $n=5.7\times 10^{5}$.
The number of features may take values ranging from $m=2$ (e.g. in the case of
models based on $M_{500}$ and $z$) to $m=6$ (e.g. for models with
$L_{X},z,EM_{0},N_{H},T_{\rm 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 optimizing 651 hyperparameters. This number decomposes
into: 6 for the kernel ($\sigma$ and the five length-scales), $30\times 5$
coordinates of the inducing points $\@vec{Z}$, 30 for the values of the
$\@vec{m}$ parameter and $30\times 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
optimization does provide satisfactory results only if mini-batches are of
large enough size. In our case we choose random mini-batches containing
$2^{16}=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^{4}$)
before stopping the training phase. In the course of our model production, we
have noticed some models converge to undesired solutions, as apparent from
e.g. 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 finalized, 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\left(I^{*}\mathrel{}\middle|\mathrel{}\@vec{x}^{*}\right)$ is a very quick
operation, even for large numbers of test points $\@vec{x}^{*}$ (Eq. 32). This
is in fact a requirement of the cosmological modelling (Ghirardini et al.
submitted, 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 $\pm 1\sigma$ enveloppe of the latent
function posterior (given by the variance in Eq. 34) into the link function
$\Phi$.

##  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_{500}$ flux in the 0.5–2 keV band is an alternative to the
unabsorbed count-rate. It has units erg s-1 cm-2 and it is available for all
simulated clusters from the twin simulation. Fig. 33 shows slices through this
model $p\left(I_{\rm{main}}\mathrel{}\middle|\mathrel{}f_{X},z,N_{H},T_{\rm
exp},{\rm bkg}\right)$. It has very similar behaviour as a selection expressed
with $CR$ (Fig. 10).

Figure 33: Representation of a model
$p\left(I_{\rm{main}}\mathrel{}\middle|\mathrel{}f_{X},z,N_{H},T_{\rm
exp},{\rm bkg}\right)$ for fixed values of cluster redshift $z=0.2$ and
$z=0.5$, local galactic absorption $N_{H}$, a nominal background level, two
values of exposure times $T_{\rm exp}$ (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, characterized by $R_{500}$ and its central emission measure,
characterized by $EM_{0}$, play a secondary role. We have (see also Comparat
et al. 2020):

$EM(R)=-\log_{10}\int_{-\infty}^{+\infty}n_{e}n_{p}\mathrm{d}l$ (36)

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-3).
For numerical stability reasons, the central emission measure $EM_{0}$ is
defined as the average of $EM(R)$ in the radial range $0<R<0.025\,R_{500}$ and
the integration bounds are set to $\pm 20\,R_{500}$ instead of infinity.

Morphological features intervene non-trivially in the selection function. At
fixed $R_{500}$ luminosity, a larger cluster size $R_{500}$ 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. in press). Similarly, low values of $EM_{0}$ indicate peaked
surface brightness profiles, decreasing the probability of selecting a cluster
as extended. Large values of $EM_{0}$ 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. 34 shows four slices through this multi-dimensional model, as
a function of the parameter $EM_{0}$ and for two different cluster
luminosities, two different exposure times.

Figure 34: Similar figure as Fig. 33, representing the outcome of a model
$p\left(I_{\rm{main}}\mathrel{}\middle|\mathrel{}L_{X},z,EM_{0},N_{H},T_{\rm
exp},{\rm bkg}\right)$ for fixed values of cluster redshift $z$, two different
values for the luminosity $L_{X}$, and local galactic absorption $N_{H}$, a
nominal background level and two values of exposure time $T_{\rm exp}$. The
shaded regions indicate the approximate 68% confidence range output of the
model.

Also on Fig. 35 we show the dependence on $R_{500}$ at fixed flux. The larger
the apparent cluster extent (here the projected $R_{500}$ 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.

Figure 35: Similar as Fig. 33, representing the outcome of a model
$p\left(I_{\rm{main}}\mathrel{}\middle|\mathrel{}f_{X},R_{500},N_{H},T_{\rm
exp},{\rm bkg}\right)$ for two different values for the flux $f_{X}$, fixed
local galactic absorption $N_{H}$, a nominal background level and two values
of exposure time $T_{\rm exp}$. The shaded regions indicate the approximate
68% confidence range output of the model.