# Characterising galaxy clusters' completeness function in Planck with hydrodynamical simulations

S. Gallo, M. Douspis, E. Soubrié, L. Salvati

► To cite this version:

S. Gallo, M. Douspis, E. Soubrié, L. Salvati. Characterising galaxy clusters' completeness function in Planck with hydrodynamical simulations. Astronomy and Astrophysics - A&A, 2024, 686, pp.A15. 10.1051/0004-6361/202347678. [hal-04587839](https://hal.science/hal-04587839)

HAL Id: hal-04587839

https://hal.science/hal-04587839

Submitted on 24 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Characterising galaxy clusters' completeness function in Planck with hydrodynamical simulations

S. Gallo<fnref n="1" />, M. Douspis<fnref n="1" />, E. Soubrié<fnref n="1" />,<fnref n="2" />✉, and L. Salvati<fnref n="1" />

Received 8 August 2023 / Accepted 22 November 2023

## ABSTRACT

Galaxy cluster number counts are an important probe with which to constrain cosmological parameters. One of the main ingredients of the analysis, along with accurate estimates of cluster masses, is the selection function, and in particular the completeness associated with the cluster sample under consideration. Incorrectly characterising this function can lead to biases in cosmological constraints. In this work, we want to study the completeness of the *Planck* cluster catalogue, estimating the probability of cluster detection in a realistic setting using hydrodynamical simulations. In particular, we probe the case in which the cluster model assumed in the detection method differs from the shapes and profiles of true galaxy clusters. We created around 9000 images of the Sunyaev–Zel'dovich effect from galaxy clusters from the IllustrisTNG simulation, and used a Monte Carlo injection method to estimate the completeness function. We studied the impact of having different cluster pressure profiles and complex cluster morphologies on the detection process. We find that the cluster profile has a significant effect on completeness, with clusters with steeper profiles producing a higher completeness than ones with flatter profiles. We also show that cluster morphology has a small impact on completeness, finding that elliptical clusters have a slightly lower probability of detection with respect to spherically symmetric ones. Finally, we investigate the impact of a different completeness function on a cosmological analysis with cluster number counts, showing a shift in the constraints on $\Omega_m$ and $\sigma_8$ that lies in the same direction as the shift driven by the mass bias.

**Key words.** methods: numerical – methods: statistical – galaxies: clusters: general – large-scale structure of Universe

## 1. Introduction

Galaxy clusters are the largest gravitationally bound objects in the Universe. They form from the highest peaks in the initial density fluctuations and grow through mergers and the accretion of smaller groups and galaxies, driven by their strong gravitational pull (Kravtsov & Borgani 2012). As a result, galaxy clusters provide valuable information about the growth of cosmic structures and help constrain parameters in the cosmological model (Oukbir & Blanchard 1992; Henry 1997; Allen et al. 2011). The abundance of galaxy clusters as a function of mass and redshift (cluster number counts), for example, is particularly sensitive to the cosmic matter density and density fluctuations ($\Omega_m$ and $\sigma_8$), as well as to the dark energy equation of state. Therefore, galaxy clusters are considered fundamental tools for understanding the Universe, and have been used as probes in numerous studies (e.g. Rozo et al. 2010; Planck Collaboration XXIV 2016; Pacaud et al. 2018; Bocquet et al. 2019; Costanzi et al. 2021).

For these reasons, there have been considerable efforts in recent years to build large galaxy cluster catalogues for cosmological analyses, exploiting their multi-component nature to detect them at different wavelengths: for example, in optical, the Dark Energy Survey (Abbott et al. 2020); in X-rays, the XXL survey (Pierre et al. 2016; Pacaud et al. 2018) and eROSITA (Liu et al. 2022); and at millimetre wavelengths, the Atacama Cosmology Telescope (ACT; Hilton et al. 2021), the South Pole telescope (SPT; Bleem et al. 2015), and the *Planck* survey (Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016).

In order to be able to extract cosmological information from a cluster survey, we must know its selection function. The selection function is indeed a key element for any statistical study carried out with a survey, because it connects the detected objects with the underlying true population in the survey area, characterising the relation between the two sets; it is a function of the cluster properties and depends on the characteristics of the survey, including the detection strategy. The selection function can be divided into two separate functions: the purity, which is the probability that a given detection corresponds to a real object, and the completeness, the probability that an object in the real population will be detected in the survey.

In particular, the selection function is one of the main ingredients for cosmological analyses with galaxy cluster number counts, because it provides – as a function of the cluster observables – an estimate of the fraction of objects detected over the total number of those objects present in the sky. This is very important information when comparing the number of clusters observed with the number predicted from theory. It is then clear why an accurate characterisation of the selection function is important: an incorrect estimation could lead to biases in cosmological parameter constraints. Another source of uncertainty when probing cosmology with galaxy clusters is the mass calibration. Cluster masses are not directly observable, and so we need to use other cluster properties that correlate with mass as proxies. The relation between the chosen observable property and the cluster mass is modelled via a statistical scaling relation, which is calibrated using multi-wavelength observations. Nonetheless, uncertainties or mischaracterisation

<footnote n="1">¹ Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405 Orsay, France
e-mail: stefano.gallo@universite-paris-saclay.fr</footnote>

<footnote n="2">² Institute of Applied Computing & Community Code (IAC³), UIB, Spain</footnote>
# Characterising galaxy clusters' completeness function in Planck with hydrodynamical simulations

S. Gallo<fnref n="1" />, M. Douspis<fnref n="1" />, E. Soubrié<fnref n="1" />,<fnref n="2" />✉, and L. Salvati<fnref n="1" />

Received 8 August 2023 / Accepted 22 November 2023

## ABSTRACT

Galaxy cluster number counts are an important probe with which to constrain cosmological parameters. One of the main ingredients of the analysis, along with accurate estimates of cluster masses, is the selection function, and in particular the completeness associated with the cluster sample under consideration. Incorrectly characterising this function can lead to biases in cosmological constraints. In this work, we want to study the completeness of the *Planck* cluster catalogue, estimating the probability of cluster detection in a realistic setting using hydrodynamical simulations. In particular, we probe the case in which the cluster model assumed in the detection method differs from the shapes and profiles of true galaxy clusters. We created around 9000 images of the Sunyaev–Zel'dovich effect from galaxy clusters from the IllustrisTNG simulation, and used a Monte Carlo injection method to estimate the completeness function. We studied the impact of having different cluster pressure profiles and complex cluster morphologies on the detection process. We find that the cluster profile has a significant effect on completeness, with clusters with steeper profiles producing a higher completeness than ones with flatter profiles. We also show that cluster morphology has a small impact on completeness, finding that elliptical clusters have a slightly lower probability of detection with respect to spherically symmetric ones. Finally, we investigate the impact of a different completeness function on a cosmological analysis with cluster number counts, showing a shift in the constraints on $\Omega_m$ and $\sigma_8$ that lies in the same direction as the shift driven by the mass bias.

**Key words.** methods: numerical – methods: statistical – galaxies: clusters: general – large-scale structure of Universe

## 1. Introduction

Galaxy clusters are the largest gravitationally bound objects in the Universe. They form from the highest peaks in the initial density fluctuations and grow through mergers and the accretion of smaller groups and galaxies, driven by their strong gravitational pull (Kravtsov & Borgani 2012). As a result, galaxy clusters provide valuable information about the growth of cosmic structures and help constrain parameters in the cosmological model (Oukbir & Blanchard 1992; Henry 1997; Allen et al. 2011). The abundance of galaxy clusters as a function of mass and redshift (cluster number counts), for example, is particularly sensitive to the cosmic matter density and density fluctuations ($\Omega_m$ and $\sigma_8$), as well as to the dark energy equation of state. Therefore, galaxy clusters are considered fundamental tools for understanding the Universe, and have been used as probes in numerous studies (e.g. Rozo et al. 2010; Planck Collaboration XXIV 2016; Pacaud et al. 2018; Bocquet et al. 2019; Costanzi et al. 2021).

For these reasons, there have been considerable efforts in recent years to build large galaxy cluster catalogues for cosmological analyses, exploiting their multi-component nature to detect them at different wavelengths: for example, in optical, the Dark Energy Survey (Abbott et al. 2020); in X-rays, the XXL survey (Pierre et al. 2016; Pacaud et al. 2018) and eROSITA (Liu et al. 2022); and at millimetre wavelengths, the Atacama Cosmology Telescope (ACT; Hilton et al. 2021), the South Pole telescope (SPT; Bleem et al. 2015), and the *Planck* survey (Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016).

In order to be able to extract cosmological information from a cluster survey, we must know its selection function. The selection function is indeed a key element for any statistical study carried out with a survey, because it connects the detected objects with the underlying true population in the survey area, characterising the relation between the two sets; it is a function of the cluster properties and depends on the characteristics of the survey, including the detection strategy. The selection function can be divided into two separate functions: the purity, which is the probability that a given detection corresponds to a real object, and the completeness, the probability that an object in the real population will be detected in the survey.

In particular, the selection function is one of the main ingredients for cosmological analyses with galaxy cluster number counts, because it provides – as a function of the cluster observables – an estimate of the fraction of objects detected over the total number of those objects present in the sky. This is very important information when comparing the number of clusters observed with the number predicted from theory. It is then clear why an accurate characterisation of the selection function is important: an incorrect estimation could lead to biases in cosmological parameter constraints. Another source of uncertainty when probing cosmology with galaxy clusters is the mass calibration. Cluster masses are not directly observable, and so we need to use other cluster properties that correlate with mass as proxies. The relation between the chosen observable property and the cluster mass is modelled via a statistical scaling relation, which is calibrated using multi-wavelength observations. Nonetheless, uncertainties or mischaracterisation

<footnote n="1">¹ Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405 Orsay, France
e-mail: stefano.gallo@universite-paris-saclay.fr</footnote>

<footnote n="2">² Institute of Applied Computing & Community Code (IAC³), UIB, Spain</footnote>
of scaling relation parameters can have an important impact on the results of a cosmological analysis. Therefore, a better understanding of these sources of bias and their possible correlations is particularly valuable, especially in light of the reported tension in the value of $\sigma_8$ (the amplitude of matter fluctuations), when computed using clusters or cosmic microwave background data (CMB; Planck Collaboration XX 2014; Planck Collaboration XXIV 2016; Bocquet et al. 2019).

The selection function is necessarily survey-specific. In this study, we aim to characterise the selection function for cosmological analyses carried out with galaxy clusters detected via the thermal Sunyaev–Zel'dovich (SZ) effect (Sunyaev & Zeldovich 1970, 1972, 1980) by the Planck satellite. For this reason, we take as reference the Planck MMF3 cosmological sample (Planck Collaboration XXVII 2016), where MMF3 is the detection algorithm used, which is based on the matched multi-filter technique (MMF, Herranz et al. 2002; Melin et al. 2006). This sample contains 439 clusters, with masses $\in[0.8, 14.7] \times 10^{14} M_{\odot}$ and redshifts $\in[0.01, 0.97]$, covering 65% of the sky. Due to its very high purity (>99.8%), in the rest of the article we focus exclusively on the completeness function.

The completeness associated with this survey was already studied by the Planck Collaboration and detailed with the various catalogue releases (Planck Collaboration VIII 2011; Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016). In these works, the completeness is first estimated assuming Gaussian noise of the SZ signal, obtaining a rather simple analytical form for the function. This was then compared with a more direct approach, which relies on mock observations obtained by injecting a population of simulated clusters into the real sky maps. With full knowledge of the “true” (injected) cluster population, it becomes possible to compare it to the output of the detection algorithm run on these mock maps. The aim of this strategy is to reproduce the (unknown) conditions of the real detection task as faithfully as possible in order to get the best possible estimate of the completeness, including all the contaminant effects that might be difficult to model analytically.

This approach was taken by the Planck Collaboration XXIX (2014), Planck Collaboration XXVII (2016) using simulated SZ signals that assumed spherical symmetry for the clusters, finding substantial agreement with the analytical completeness. The latter was then incorporated as the baseline estimate of the completeness in cosmological analyses of cluster number counts with the Planck SZ catalogue (e.g. Planck Collaboration XX 2014; Planck Collaboration XXIV 2016).

In this paper, we use these works as a starting point, and analyse the completeness focusing on the case in which the cluster model assumed in the detection algorithm is different from the “true” injected cluster signals. This situation is to be expected, to an extent, given that a template is by necessity a simplification constructed to match the “average” features of a selected cluster sample as closely as possible. In particular, the Planck cluster model is based on the assumption of spherical symmetry, and a single pressure profile is assumed.

However, galaxy clusters are known from both simulations and observations to be generally not spherical (Limousin et al. 2013), which is due to various dynamical effects such as mergers and asymmetric accretion through cosmic filaments (Gouin et al. 2020, 2022; Vallés-Pérez et al. 2020). Departure from spherical symmetry is a first clear difference from the detection template, which can possibly bias cluster detection, and therefore the effect of realistic morphology needs to be tested while characterising the completeness. A second difference between real clusters and the detection template might come from the pressure profile.

First of all, not all clusters have exactly the same profile; there are variations due, for example, to the clusters' dynamical state, and these induce a scatter around the average profile of the population. Moreover, the average profile might also differ from the one assumed in the template. This might happen, for example, if the profile assumed in the detection was measured from a biased sample of clusters.

Some tests to probe these aspects were performed in Planck Collaboration XXVII (2016). To investigate the effect of cluster morphology, a modest sample of hydrodynamically simulated clusters was used, with a fixed angular scale larger than the Planck beam, where the effect of cluster morphology is supposed to be most relevant. No significant difference was found in the completeness using either realistic or spherical morphologies. Regarding the profile scatter, Planck Collaboration XXVII (2016) show that the completeness computed from cluster images with pressure profiles scattered around the one assumed in the detection is generally consistent with the analytical completeness estimate, but report a widening effect in the completeness drop-off.

In this work, we study these effects in a comprehensive way, analysing the impact of having an imperfectly matching cluster model as template for the matched filter detection technique. To this end, we used a sample of clusters from a large-volume, state-of-the-art hydrodynamical simulation in a somewhat agnostic way. We extract clusters from the simulation at different redshifts, and produce images of their SZ signal as it would be seen on the sky based on each cluster's redshift and gas distribution. In this way, when computing the completeness, we automatically include any possible redshift dependence of the cluster properties and dynamical states, as well as the effect due to cluster morphologies, in order to increase the realism of the completeness estimation.

In Sect. 2, we present the sets of cluster images and the Planck sky maps used for the completeness analysis. In Sect. 3, we describe the MMF detection method and discuss the completeness function and the ways it is estimated in detail. The results on completeness from our realistic cluster images are presented in Sect. 4, where we investigate the impact of the different profiles and asymmetric morphologies of the clusters. We discuss these results and their impact on cosmological analyses in Sect. 5, together with the limitations of our methods. Finally, in Sect. 6, we draw conclusions and describe some future perspectives.

## 2. Data

First, we describe the creation of the mock galaxy cluster images that we used to estimate the completeness. Starting alternatively from publicly available simulation data and analytical pressure profiles, we computed the SZ signal from galaxy clusters as it would be seen by the Planck satellite. We then present the Planck sky maps, which we cleaned from real detections to serve as background for the injection of the mock cluster signals.

### 2.1. Cluster SZ images

At millimetre wavelengths, such as the ones observed by the Planck satellite, clusters of galaxies can be observed through the SZ effect. The SZ effect is a spectral distortion of the CMB due to the inverse Compton scattering of CMB photons off energetic electrons in the hot ionised gas in galaxy clusters.

While travelling through a galaxy cluster, CMB photons can be scattered by hot electrons in the intracluster gas and gain
of scaling relation parameters can have an important impact on the results of a cosmological analysis. Therefore, a better understanding of these sources of bias and their possible correlations is particularly valuable, especially in light of the reported tension in the value of $\sigma_8$ (the amplitude of matter fluctuations), when computed using clusters or cosmic microwave background data (CMB; Planck Collaboration XX 2014; Planck Collaboration XXIV 2016; Bocquet et al. 2019).

The selection function is necessarily survey-specific. In this study, we aim to characterise the selection function for cosmological analyses carried out with galaxy clusters detected via the thermal Sunyaev–Zel'dovich (SZ) effect (Sunyaev & Zeldovich 1970, 1972, 1980) by the Planck satellite. For this reason, we take as reference the Planck MMF3 cosmological sample (Planck Collaboration XXVII 2016), where MMF3 is the detection algorithm used, which is based on the matched multi-filter technique (MMF, Herranz et al. 2002; Melin et al. 2006). This sample contains 439 clusters, with masses $\in[0.8, 14.7] \times 10^{14} M_{\odot}$ and redshifts $\in[0.01, 0.97]$, covering 65% of the sky. Due to its very high purity (>99.8%), in the rest of the article we focus exclusively on the completeness function.

The completeness associated with this survey was already studied by the Planck Collaboration and detailed with the various catalogue releases (Planck Collaboration VIII 2011; Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016). In these works, the completeness is first estimated assuming Gaussian noise of the SZ signal, obtaining a rather simple analytical form for the function. This was then compared with a more direct approach, which relies on mock observations obtained by injecting a population of simulated clusters into the real sky maps. With full knowledge of the “true” (injected) cluster population, it becomes possible to compare it to the output of the detection algorithm run on these mock maps. The aim of this strategy is to reproduce the (unknown) conditions of the real detection task as faithfully as possible in order to get the best possible estimate of the completeness, including all the contaminant effects that might be difficult to model analytically.

This approach was taken by the Planck Collaboration XXIX (2014), Planck Collaboration XXVII (2016) using simulated SZ signals that assumed spherical symmetry for the clusters, finding substantial agreement with the analytical completeness. The latter was then incorporated as the baseline estimate of the completeness in cosmological analyses of cluster number counts with the Planck SZ catalogue (e.g. Planck Collaboration XX 2014; Planck Collaboration XXIV 2016).

In this paper, we use these works as a starting point, and analyse the completeness focusing on the case in which the cluster model assumed in the detection algorithm is different from the “true” injected cluster signals. This situation is to be expected, to an extent, given that a template is by necessity a simplification constructed to match the “average” features of a selected cluster sample as closely as possible. In particular, the Planck cluster model is based on the assumption of spherical symmetry, and a single pressure profile is assumed.

However, galaxy clusters are known from both simulations and observations to be generally not spherical (Limousin et al. 2013), which is due to various dynamical effects such as mergers and asymmetric accretion through cosmic filaments (Gouin et al. 2020, 2022; Vallés-Pérez et al. 2020). Departure from spherical symmetry is a first clear difference from the detection template, which can possibly bias cluster detection, and therefore the effect of realistic morphology needs to be tested while characterising the completeness. A second difference between real clusters and the detection template might come from the pressure profile.

First of all, not all clusters have exactly the same profile; there are variations due, for example, to the clusters' dynamical state, and these induce a scatter around the average profile of the population. Moreover, the average profile might also differ from the one assumed in the template. This might happen, for example, if the profile assumed in the detection was measured from a biased sample of clusters.

Some tests to probe these aspects were performed in Planck Collaboration XXVII (2016). To investigate the effect of cluster morphology, a modest sample of hydrodynamically simulated clusters was used, with a fixed angular scale larger than the Planck beam, where the effect of cluster morphology is supposed to be most relevant. No significant difference was found in the completeness using either realistic or spherical morphologies. Regarding the profile scatter, Planck Collaboration XXVII (2016) show that the completeness computed from cluster images with pressure profiles scattered around the one assumed in the detection is generally consistent with the analytical completeness estimate, but report a widening effect in the completeness drop-off.

In this work, we study these effects in a comprehensive way, analysing the impact of having an imperfectly matching cluster model as template for the matched filter detection technique. To this end, we used a sample of clusters from a large-volume, state-of-the-art hydrodynamical simulation in a somewhat agnostic way. We extract clusters from the simulation at different redshifts, and produce images of their SZ signal as it would be seen on the sky based on each cluster's redshift and gas distribution. In this way, when computing the completeness, we automatically include any possible redshift dependence of the cluster properties and dynamical states, as well as the effect due to cluster morphologies, in order to increase the realism of the completeness estimation.

In Sect. 2, we present the sets of cluster images and the Planck sky maps used for the completeness analysis. In Sect. 3, we describe the MMF detection method and discuss the completeness function and the ways it is estimated in detail. The results on completeness from our realistic cluster images are presented in Sect. 4, where we investigate the impact of the different profiles and asymmetric morphologies of the clusters. We discuss these results and their impact on cosmological analyses in Sect. 5, together with the limitations of our methods. Finally, in Sect. 6, we draw conclusions and describe some future perspectives.

## 2. Data

First, we describe the creation of the mock galaxy cluster images that we used to estimate the completeness. Starting alternatively from publicly available simulation data and analytical pressure profiles, we computed the SZ signal from galaxy clusters as it would be seen by the Planck satellite. We then present the Planck sky maps, which we cleaned from real detections to serve as background for the injection of the mock cluster signals.

### 2.1. Cluster SZ images

At millimetre wavelengths, such as the ones observed by the Planck satellite, clusters of galaxies can be observed through the SZ effect. The SZ effect is a spectral distortion of the CMB due to the inverse Compton scattering of CMB photons off energetic electrons in the hot ionised gas in galaxy clusters.

While travelling through a galaxy cluster, CMB photons can be scattered by hot electrons in the intracluster gas and gain
some energy from them. This produces a peculiar distortion in the spectrum of CMB radiation that can be distinguished and detected, and encodes information about the distribution of the hot gas. In particular, neglecting relativistic corrections, the amplitude of the SZ effect is proportional to the Comptonisation parameter $y$, which in turn is proportional to the integrated electron pressure $P_e$ along the line of sight:

$$y = \frac{\sigma_T}{m_e c^2} \int P_e(l) dl, \qquad (1)$$

where $\sigma_T$ is the Thomson cross-section, $m_e$ the electron mass, and $c$ is the speed of light.

Therefore, with the SZ effect, it is possible to probe the pressure distribution of the gas in the intracluster medium, and its signal can be used to detect galaxy clusters. In the following, we explain the production of a set of images of SZ signal, starting from the gas content in simulated galaxy clusters from the IllustrisTNG simulations.

### 2.1.1. Simulation

IllustrisTNG is a suite of cosmological magnetohydrodynamical simulations that follow the coupled evolution of dark and baryonic matter through cosmic time, starting from redshift $z = 127$ to the present time (Nelson et al. 2018, 2019; Pillepich et al. 2018; Springel et al. 2018; Naiman et al. 2018; Marinacci et al. 2018). These simulations are run using the moving-mesh code AREPO (Springel 2010), and assume cosmological parameters consistent with the results of Planck Collaboration XIII (2016). In this work, among the simulations in the suite, we focus on the realisation called TNG300-1, with a comoving volume of $205 \text{ Mpc h}^{-1}$ and a mass resolution for baryons of around $7.6 \times 10^6 \text{ M}_\odot \text{ h}^{-1}$.

The simulation output provides snapshots at different redshifts, including a catalogue of halos, created running a friends-of-friends (FoF) algorithm (Davis et al. 1985) on the dark matter particles, with linking length $b = 0.2$. Gas particles are then sorted according to the closest DM particle. This catalogue lists the positions of the halos (identified by the particle with the minimum gravitational potential energy), their masses $M_{500}$, and radii $R_{500}^1$, among other quantities.

As our goal is to characterise the performance of the MMF detection algorithm, and in particular the completeness, we need a set of clusters that provides good coverage of the region in mass and redshift where the detection becomes more difficult, and the completeness goes from 1 to 0. To get an estimate for this region, we take the Planck PSZ2 cluster sample (Planck Collaboration XXVII 2016, Fig. 26) and observe that the lowest mass of the detected clusters increases with redshift, from $\sim 10^{14} M_\odot$ at $z \sim 0.05$ to $\sim 3 \times 10^{14} M_\odot$ at $z \sim 0.3$. We want our sample of simulated clusters to extend to lower masses, sufficiently below the limit of the Planck detections, to be able to study the completeness down to 0. At the same time, we do not want to include too many clusters with almost zero probability of being detected, because those would not add any information to the completeness. Therefore, we select halos from the simulation imposing a redshift-dependent lower mass limit of $M_{500} \ge M_{\min} = \frac{6}{5}(4z + \frac{7}{15}) \times 10^{14} M_\odot$ in the redshift range $0.05 \le z \le 0.3$ (which corresponds to 18 snapshots of the simulation).

<figure><img src="image_1.png" /><figcaption>Fig. 1. Distribution of the galaxy clusters selected from the TNG300 simulation. Top panel: mass and redshift distribution of the galaxy clusters colour coded according to their $Y_{5R500}$. The violet line shows the lower mass limit imposed in the selection and described in the text. Bottom panel: cluster distribution in integrated $y$ signal, $Y_{5R500}$, and angular scale, $\theta_{500}$, colour coded according to cluster mass $M_{500}$.</figcaption></figure>

With this selection, we obtain a total of 1487 clusters, whose distribution in mass and redshift is shown in Fig. 1.

This choice of the redshift range, which is smaller than that of the Planck PSZ2 catalogue ($0.01 < z < 0.97$, Planck Collaboration XXVII 2016), is imposed by the limits of the simulation. Outside of this range, the cluster distribution in radius and SZ flux does not allow us to sample the completeness function properly. In particular, at redshift higher than 0.3, there are not enough high-mass halos, while at low redshift, the spacing between the snapshots leaves empty regions in the domain of the completeness function. This limit is discussed in greater detail in Sect. 5. Nevertheless, we verified that about 65% of the Planck PSZ2 cosmological sample falls inside the mass and redshift range covered by our simulation sample.

### 2.1.2. Simulation images

Having selected the sample of galaxy clusters, we proceeded by computing the projected images of the SZ effect to be used for the completeness estimation. First, we computed the electron pressure for each gas cell associated with the cluster halos, $P_e$. Starting from the cell densities, electron abundances, and internal energies provided in the simulation output, we computed

<footnote>$^1$ $R_{500}$ is the radius of a sphere centred on the halo within which the average density is 500 times the critical density $\rho_{\text{crit}}(z)$; $M_{500}$ is the mass contained inside $R_{500}$.</footnote>
some energy from them. This produces a peculiar distortion in the spectrum of CMB radiation that can be distinguished and detected, and encodes information about the distribution of the hot gas. In particular, neglecting relativistic corrections, the amplitude of the SZ effect is proportional to the Comptonisation parameter $y$, which in turn is proportional to the integrated electron pressure $P_e$ along the line of sight:

$$y = \frac{\sigma_T}{m_e c^2} \int P_e(l) dl, \qquad (1)$$

where $\sigma_T$ is the Thomson cross-section, $m_e$ the electron mass, and $c$ is the speed of light.

Therefore, with the SZ effect, it is possible to probe the pressure distribution of the gas in the intracluster medium, and its signal can be used to detect galaxy clusters. In the following, we explain the production of a set of images of SZ signal, starting from the gas content in simulated galaxy clusters from the IllustrisTNG simulations.

### 2.1.1. Simulation

IllustrisTNG is a suite of cosmological magnetohydrodynamical simulations that follow the coupled evolution of dark and baryonic matter through cosmic time, starting from redshift $z = 127$ to the present time (Nelson et al. 2018, 2019; Pillepich et al. 2018; Springel et al. 2018; Naiman et al. 2018; Marinacci et al. 2018). These simulations are run using the moving-mesh code AREPO (Springel 2010), and assume cosmological parameters consistent with the results of Planck Collaboration XIII (2016). In this work, among the simulations in the suite, we focus on the realisation called TNG300-1, with a comoving volume of $205 \text{ Mpc h}^{-1}$ and a mass resolution for baryons of around $7.6 \times 10^6 \text{ M}_\odot \text{ h}^{-1}$.

The simulation output provides snapshots at different redshifts, including a catalogue of halos, created running a friends-of-friends (FoF) algorithm (Davis et al. 1985) on the dark matter particles, with linking length $b = 0.2$. Gas particles are then sorted according to the closest DM particle. This catalogue lists the positions of the halos (identified by the particle with the minimum gravitational potential energy), their masses $M_{500}$, and radii $R_{500}^1$, among other quantities.

As our goal is to characterise the performance of the MMF detection algorithm, and in particular the completeness, we need a set of clusters that provides good coverage of the region in mass and redshift where the detection becomes more difficult, and the completeness goes from 1 to 0. To get an estimate for this region, we take the Planck PSZ2 cluster sample (Planck Collaboration XXVII 2016, Fig. 26) and observe that the lowest mass of the detected clusters increases with redshift, from $\sim 10^{14} M_\odot$ at $z \sim 0.05$ to $\sim 3 \times 10^{14} M_\odot$ at $z \sim 0.3$. We want our sample of simulated clusters to extend to lower masses, sufficiently below the limit of the Planck detections, to be able to study the completeness down to 0. At the same time, we do not want to include too many clusters with almost zero probability of being detected, because those would not add any information to the completeness. Therefore, we select halos from the simulation imposing a redshift-dependent lower mass limit of $M_{500} \ge M_{\min} = \frac{6}{5}(4z + \frac{7}{15}) \times 10^{14} M_\odot$ in the redshift range $0.05 \le z \le 0.3$ (which corresponds to 18 snapshots of the simulation).

<figure><img src="image_1.png" /><figcaption>Fig. 1. Distribution of the galaxy clusters selected from the TNG300 simulation. Top panel: mass and redshift distribution of the galaxy clusters colour coded according to their $Y_{5R500}$. The violet line shows the lower mass limit imposed in the selection and described in the text. Bottom panel: cluster distribution in integrated $y$ signal, $Y_{5R500}$, and angular scale, $\theta_{500}$, colour coded according to cluster mass $M_{500}$.</figcaption></figure>

With this selection, we obtain a total of 1487 clusters, whose distribution in mass and redshift is shown in Fig. 1.

This choice of the redshift range, which is smaller than that of the Planck PSZ2 catalogue ($0.01 < z < 0.97$, Planck Collaboration XXVII 2016), is imposed by the limits of the simulation. Outside of this range, the cluster distribution in radius and SZ flux does not allow us to sample the completeness function properly. In particular, at redshift higher than 0.3, there are not enough high-mass halos, while at low redshift, the spacing between the snapshots leaves empty regions in the domain of the completeness function. This limit is discussed in greater detail in Sect. 5. Nevertheless, we verified that about 65% of the Planck PSZ2 cosmological sample falls inside the mass and redshift range covered by our simulation sample.

### 2.1.2. Simulation images

Having selected the sample of galaxy clusters, we proceeded by computing the projected images of the SZ effect to be used for the completeness estimation. First, we computed the electron pressure for each gas cell associated with the cluster halos, $P_e$. Starting from the cell densities, electron abundances, and internal energies provided in the simulation output, we computed

<footnote>$^1$ $R_{500}$ is the radius of a sphere centred on the halo within which the average density is 500 times the critical density $\rho_{\text{crit}}(z)$; $M_{500}$ is the mass contained inside $R_{500}$.</footnote>
the electron number density $n_e$ and temperature $T_e$, and from these, the electron pressure as $P_e = k_b n_e T_e$. Then, we computed images of the Compton-y parameter from the clusters. For each cluster, we took six projections: three along the axis of the simulation box and three along the axis rotated by Euler angles $(\alpha, \beta, \gamma) = (45^\circ, 45^\circ, 45^\circ)$. Each image is 4 Mpc wide and is centred at the cluster position with a resolution of 256 × 256 pixels. In each pixel, the Compton-y parameter is computed as

$$y = \frac{\sigma_T}{m_e c^2} \int P_e dl \approx \frac{\sigma_T}{m_e c^2} \frac{\sum_i P_{e,i} V_i}{A_{\text{pix}}}, \quad (2)$$

where $P_{e,i}$ and $V_i$ are the electron pressure and the volume of the ith gas cell, respectively, and the index $i$ in the sum runs over all gas cells whose centre (provided by the simulation) falls inside the pixel area $A_{\text{pix}}$. This rather “crude” approximation of the line of sight integral only works when the typical size of the gas cells is smaller than the size of the pixels. In principle, this is true only in the centre of the clusters, where density is higher (along with the temperature and therefore pressure). However, in the outer parts of clusters, the cell sizes are larger, but at the same time, the pressure is low, making the error we commit by using this approximation negligible. We tested different pixel sizes and found no appreciable difference in the average cluster y profile (and integrated y signal).

Before being able to inject these images into the Planck frequency maps, they first need to be converted into angular coordinates, as they would be observed on the sky. They are rescaled according to their redshift $z$, using the relations:

$$\theta_{\text{img}} = \frac{4 \text{ Mpc}}{d_{\text{A}}(z)}, \qquad \theta_{500} = \frac{R_{500}}{d_{\text{A}}(z)}, \qquad (3)$$

where $d_{\text{A}}(z)$ is the angular diameter distance, $\theta_{\text{img}}$ is the angular size of the image, and $\theta_{500}$ is the equivalent of $R_{500}$ in angular coordinates and represents the scale of the cluster. For the rescaled images, we use a pixel size of 0.5 arcmin. This pixel size is therefore smaller than that of the Planck maps, of about 1.7 arcmin; this is done to avoid having excessively coarse images before convolving them with the Planck beams. Later, when the images are injected into the maps, their resolution is adapted to that of the maps.

To make the images consistent with Planck observations at the six frequencies of Planck HFI (100, 143, 217, 353, 545, 857 GHz), we convolved them with the corresponding beam, assumed to be circular Gaussian with FWHM taken from Planck Collaboration XXVII (2016), Planck Collaboration VII (2014). Finally, the Compton-y images were transformed into frequency images, multiplying the y parameter by the value of the thermal SZ effect at the six frequencies, neglecting relativistic corrections:

$$\frac{\Delta T}{T_{\text{CMB}}} = y \cdot g(v), \qquad (4)$$

where $T_{\text{CMB}}$ is the CMB temperature, and $g(v)$ is the spectral signature of the tSZ effect integrated over Planck frequency bandwidths taken from Planck Collaboration XXII (2016).

In this way, we obtained a set of 8922 cluster images at the six Planck frequencies, each associated with the cluster angular scale, $\theta_{500}$, and the integrated SZ signal within a radius of $5 \times \theta_{500}$ from the cluster centre, $Y_{5R500}$.

### 2.1.3. Circular images

As a complement to the set of images extracted from the simulation, we also generated different sets of spherically symmetric

<table><caption>**Table 1.** Generalised NFW pressure profile parameters of the different sets of spherical images.</caption><thead><tr><th>Name</th><th>c<sub>500</sub></th><th>α</th><th>β</th><th>γ</th></tr></thead><tbody><tr><td>Standard</td><td>1.177</td><td>1.051</td><td>5.4905</td><td>0.3081</td></tr><tr><td>Planck</td><td>1.81</td><td>1.33</td><td>4.13</td><td>0.31</td></tr><tr><td>PACT</td><td>1.18</td><td>1.08</td><td>4.30</td><td>0.31</td></tr><tr><td>Tramonte+23</td><td>2.1</td><td>2.2</td><td>5.3</td><td>0.31</td></tr><tr><td>Peaked</td><td>1.5</td><td>1.051</td><td>5.4905</td><td>0.3081</td></tr><tr><td>SimFit</td><td>5.1 × 10<sup>-3</sup></td><td>0.71</td><td>1.33</td><td>500</td></tr></tbody></table>

**Notes.** In order: Arnaud et al. (2010) profile (Standard), Planck Collaboration Int. V (2013) profile (Planck), Pointecouteau et al. (2021) profile (PACT), Tramonte et al. (2023) profile (Tramonte+23), profile obtained changing the $c_{500}$ of Arnaud et al. (2010) (Peaked), and profile obtained fitting the average profile from the simulation images (SimFit).

cluster images to test the consistency of the results obtained with the simulation set and explore their implications.

To best compare the completeness obtained from these sets of images with that from the simulation, we want their distribution in ($Y_{5R500}$, $\theta_{500}$) to be similar. Therefore, we constructed a catalogue of ($Y_{5R500}$, $\theta_{500}$), taking all the pairs from the simulation set and applying a random offset to the two values sampled from a Gaussian with a standard deviation of 5% of each value in order to sample the same region but not exactly the same values. These new pairs of values were then used to generate the spherical images.

The images in the different sets were all constructed starting from a spherical pressure distribution, which models the gas in an ideal galaxy cluster. The form of the pressure profile is the generalised Navarro-Frenk-White (gNFW; Navarro et al. 1997; Nagai et al. 2007; Arnaud et al. 2010):

$$p(x) \propto \frac{1}{(c_{500}x)^{\gamma}[1 + (c_{500}x)^{\alpha}]^{[\beta-\gamma]/\alpha}}, \qquad (5)$$

where $x = r/R_{500}$ is the radius in units of $R_{500}$, and $[c_{500}, \alpha, \beta, \gamma]$ are the parameters that determine the shape of the profile. We used six different sets of parameters to build the images, which are detailed in Table 1: four come from observational studies (Arnaud et al. 2010; Planck Collaboration Int. V 2013; Pointecouteau et al. 2021; Tramonte et al. 2023), and two are artificially constructed to approximate the average profile of the simulation images.

To obtain the SZ images for each of the different sets, we first integrated the pressure profile along one direction and transformed into a y map using Eq. (1). We then rescaled this map to match the various ($Y_{5R500}$, $\theta_{500}$) of the catalogue described above. Hereafter, the y images obtained were treated with the same steps as the simulation images; namely convolution with the Planck beams and transformation into frequency images using Eq. (4). The Compton-y profiles obtained from the pressure profiles of Table 1 are shown in Fig. 2.

### 2.2. Planck frequency maps

To construct the sky maps into which we injected our cluster images, we started from the six Planck HFI frequency maps from the second data release Planck Collaboration VIII (2016). These maps are given in HEALPix (Górski et al. 2005) pixelisation scheme, with $N_{\text{side}} = 2048$.
the electron number density $n_e$ and temperature $T_e$, and from these, the electron pressure as $P_e = k_b n_e T_e$. Then, we computed images of the Compton-y parameter from the clusters. For each cluster, we took six projections: three along the axis of the simulation box and three along the axis rotated by Euler angles $(\alpha, \beta, \gamma) = (45^\circ, 45^\circ, 45^\circ)$. Each image is 4 Mpc wide and is centred at the cluster position with a resolution of 256 × 256 pixels. In each pixel, the Compton-y parameter is computed as

$$y = \frac{\sigma_T}{m_e c^2} \int P_e dl \approx \frac{\sigma_T}{m_e c^2} \frac{\sum_i P_{e,i} V_i}{A_{\text{pix}}}, \quad (2)$$

where $P_{e,i}$ and $V_i$ are the electron pressure and the volume of the ith gas cell, respectively, and the index i in the sum runs over all gas cells whose centre (provided by the simulation) falls inside the pixel area $A_{\text{pix}}$. This rather “crude” approximation of the line of sight integral only works when the typical size of the gas cells is smaller than the size of the pixels. In principle, this is true only in the centre of the clusters, where density is higher (along with the temperature and therefore pressure). However, in the outer parts of clusters, the cell sizes are larger, but at the same time, the pressure is low, making the error we commit by using this approximation negligible. We tested different pixel sizes and found no appreciable difference in the average cluster y profile (and integrated y signal).

Before being able to inject these images into the Planck frequency maps, they first need to be converted into angular coordinates, as they would be observed on the sky. They are rescaled according to their redshift $z$, using the relations:

$$\theta_{\text{img}} = \frac{4 \text{ Mpc}}{d_{\text{A}}(z)}, \qquad \theta_{500} = \frac{R_{500}}{d_{\text{A}}(z)}, \qquad (3)$$

where $d_{\text{A}}(z)$ is the angular diameter distance, $\theta_{\text{img}}$ is the angular size of the image, and $\theta_{500}$ is the equivalent of $R_{500}$ in angular coordinates and represents the scale of the cluster. For the rescaled images, we use a pixel size of 0.5 arcmin. This pixel size is therefore smaller than that of the Planck maps, of about 1.7 arcmin; this is done to avoid having excessively coarse images before convolving them with the Planck beams. Later, when the images are injected into the maps, their resolution is adapted to that of the maps.

To make the images consistent with Planck observations at the six frequencies of Planck HFI (100, 143, 217, 353, 545, 857 GHz), we convolved them with the corresponding beam, assumed to be circular Gaussian with FWHM taken from Planck Collaboration XXVII (2016), Planck Collaboration VII (2014). Finally, the Compton-y images were transformed into frequency images, multiplying the y parameter by the value of the thermal SZ effect at the six frequencies, neglecting relativistic corrections:

$$\frac{\Delta T}{T_{\text{CMB}}} = y \cdot g(v), \qquad (4)$$

where $T_{\text{CMB}}$ is the CMB temperature, and $g(v)$ is the spectral signature of the tSZ effect integrated over Planck frequency bandwidths taken from Planck Collaboration XXII (2016).

In this way, we obtained a set of 8922 cluster images at the six Planck frequencies, each associated with the cluster angular scale, $\theta_{500}$, and the integrated SZ signal within a radius of $5 \times \theta_{500}$ from the cluster centre, $Y_{5R500}$.

### 2.1.3. Circular images

As a complement to the set of images extracted from the simulation, we also generated different sets of spherically symmetric

<table><caption>**Table 1.** Generalised NFW pressure profile parameters of the different sets of spherical images.</caption><thead><tr><th>Name</th><th>c<sub>500</sub></th><th>α</th><th>β</th><th>γ</th></tr></thead><tbody><tr><td>Standard</td><td>1.177</td><td>1.051</td><td>5.4905</td><td>0.3081</td></tr><tr><td>Planck</td><td>1.81</td><td>1.33</td><td>4.13</td><td>0.31</td></tr><tr><td>PACT</td><td>1.18</td><td>1.08</td><td>4.30</td><td>0.31</td></tr><tr><td>Tramonte+23</td><td>2.1</td><td>2.2</td><td>5.3</td><td>0.31</td></tr><tr><td>Peaked</td><td>1.5</td><td>1.051</td><td>5.4905</td><td>0.3081</td></tr><tr><td>SimFit</td><td>5.1 × 10<sup>-3</sup></td><td>0.71</td><td>1.33</td><td>500</td></tr></tbody></table>

**Notes.** In order: Arnaud et al. (2010) profile (Standard), Planck Collaboration Int. V (2013) profile (Planck), Pointecouteau et al. (2021) profile (PACT), Tramonte et al. (2023) profile (Tramonte+23), profile obtained changing the $c_{500}$ of Arnaud et al. (2010) (Peaked), and profile obtained fitting the average profile from the simulation images (SimFit).

cluster images to test the consistency of the results obtained with the simulation set and explore their implications.

To best compare the completeness obtained from these sets of images with that from the simulation, we want their distribution in ($Y_{5R500}$, $\theta_{500}$) to be similar. Therefore, we constructed a catalogue of ($Y_{5R500}$, $\theta_{500}$), taking all the pairs from the simulation set and applying a random offset to the two values sampled from a Gaussian with a standard deviation of 5% of each value in order to sample the same region but not exactly the same values. These new pairs of values were then used to generate the spherical images.

The images in the different sets were all constructed starting from a spherical pressure distribution, which models the gas in an ideal galaxy cluster. The form of the pressure profile is the generalised Navarro-Frenk-White (gNFW; Navarro et al. 1997; Nagai et al. 2007; Arnaud et al. 2010):

$$p(x) \propto \frac{1}{(c_{500}x)^{\gamma}[1 + (c_{500}x)^{\alpha}]^{[\beta-\gamma]/\alpha}}, \qquad (5)$$

where $x = r/R_{500}$ is the radius in units of $R_{500}$, and [$c_{500}$, $\alpha$, $\beta$, $\gamma$] are the parameters that determine the shape of the profile. We used six different sets of parameters to build the images, which are detailed in Table 1: four come from observational studies (Arnaud et al. 2010; Planck Collaboration Int. V 2013; Pointecouteau et al. 2021; Tramonte et al. 2023), and two are artificially constructed to approximate the average profile of the simulation images.

To obtain the SZ images for each of the different sets, we first integrated the pressure profile along one direction and transformed into a y map using Eq. (1). We then rescaled this map to match the various ($Y_{5R500}$, $\theta_{500}$) of the catalogue described above. Hereafter, the y images obtained were treated with the same steps as the simulation images; namely convolution with the Planck beams and transformation into frequency images using Eq. (4). The Compton-y profiles obtained from the pressure profiles of Table 1 are shown in Fig. 2.

### 2.2. Planck frequency maps

To construct the sky maps into which we injected our cluster images, we started from the six Planck HFI frequency maps from the second data release Planck Collaboration VIII (2016). These maps are given in HEALPix (Górski et al. 2005) pixelisation scheme, with $N_{\text{side}} = 2048$.
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Compton-$y$ profiles obtained from the gNFW profiles of Table 1, in units of $Y_{500}/R_{500}^2$, as a function of normalized radius, $R/R_{500}$. In orange: profile from Arnaud et al. (2010). In green: profile from Planck Collaboration Int. V (2013); in pink: profile from Pointecouteau et al. (2021). In red: profile from Tramonte et al. (2023). In purple: “Peaked” profile. In light blue: “SimFit” profile.</figcaption></figure>

We chose to use real Planck maps – following Planck Collaboration XXIX (2014) – to ensure the most realistic setting possible for the completeness analysis, including all sources of noise and contaminations present during the real detection process. To this end, we chose the maps that were used for the detection of the second Planck SZ cluster catalogue (PSZ2, Planck Collaboration XXVII 2016), as these contain the original cluster sample whose completeness we aim to characterise.

Given the fact that the MMF algorithm estimates the noise directly from the input maps, the injection of simulated cluster signals in addition to the real ones already present could change the properties of the noise, and therefore impact the detection. For this reason, we first subjected the maps to a cleaning process in order to try to remove the SZ signal from the real clusters. The way a given cluster is removed from the maps is the following: starting from the cluster's integrated SZ flux and angular scale, ($Y_{5R500}$, $\theta_{500}$), a circular image of the cluster SZ emission at the frequencies of the different maps is computed in the same way as in Sect. 2.1.3, that is, using the profile from Arnaud et al. (2010). This cluster emission is then subtracted from the maps at the position of the original detected cluster. The differences between the circular template and the real cluster signal leave a residual contribution, but its impact on the noise estimation is certainly smaller than that of the original cluster and is therefore negligible.

The cleaning proceeds in two steps: the first step consists in simply removing all the clusters from the PSZ2 catalogue, obtaining a first “cleaned” version of the maps. In the second step, the MMF detection is run on the new maps to check for additional signals identified as clusters beyond those contained in the PSZ2 catalogue, with a lower limit of 4.25 in signal-to-noise ratio (S/N). All new detections obtained in this way are then also removed from the frequency maps, which are now, in principle, free from any relevant SZ source down to $S/N \sim 4.5$. These are the final cleaned maps we use for the completeness analysis.

Associated with the clean frequency maps, we create a mask that covers the emission from the Galaxy and the Magellanic Cloud, as well as an area of five times the beam size around point

sources (from PCCS2 catalogue Planck Collaboration XXVI 2016) and regions of CO emission Planck Collaboration XIII (2014). The final unmasked sky fraction is about 78%.

# 3. Methods

## 3.1. MMF

The matched multi-frequency filter (MMF) algorithm (Herranz et al. 2002; Melin et al. 2006) is a commonly used algorithm for the detection of galaxy clusters through their SZ signal. Notably, it has been used to produce cluster catalogues from Planck data, along with ACT and SPT data (Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016; Hilton et al. 2021; Bleem et al. 2015). This method is designed to combine the prior knowledge of the SZ signal from galaxy clusters, namely the spectral signature and spatial characteristics, to produce an optimal filter that returns the maximal S/N in the presence of a galaxy cluster. While the spectral shape of the SZ signal is well known and universal (that is in the non-relativistic regime), one has to make some approximations for the spatial filter; the usual choice is to assume spherical symmetry for clusters, and model the radial pressure profile according to the average profile of observed clusters. However, galaxy clusters are of course known to be aspherical in general, and rather triaxial in shape, and most have even more complex morphologies due to anisotropic accretion, mergers, and astrophysical processes such as shocks and feedback mechanisms. All this variety in morphological features prevents a perfect match between the cluster SZ signal and the spatial template used, which, in principle, impacts the detection performance. This effect is, in practice, reduced by the smoothing induced by the instrument's beam, which tends to symmetrise the signal, especially for clusters with scales comparable to or smaller than the beam size.

More specifically, our implementation of the MMF algorithm is very similar to that used in the construction of the Planck SZ cosmological catalogue (Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016), namely MMF3. This latter uses a projected gNFW profile as a cluster model, with parameters from Arnaud et al. (2010), [$c_{500}$, $\gamma$, $\alpha$, $\beta$] = [1.177, 0.308, 1.051, 5.491], and two additional free parameters: the signal amplitude $Y_{5R500}$ and the cluster scale $\theta_{500}$. The spatial filters are therefore constructed from this model, varying the cluster size on a grid of 40 logarithmically spaced points from $\theta_{500} = 1.059$ to 41.195 arcmin. When run, the MMF algorithm first divides the Planck full-sky maps into 546 overlapping square patches of 10° on a side. Each patch is then filtered with the templates, resulting in a set of S/N maps. The peaks in these maps over a S/N threshold of 3 are our candidate cluster detections. The integrated cluster signal $Y_{5R500}$ is given by the filtered map at the peak position, while the cluster size is assumed to be the scale of the filter that maximises the S/N at the cluster location. Finally, the clusters are merged in a single full-sky catalogue by merging all detections with lower S/N that fall inside the $\theta_{500}$ of a certain detected cluster.

The final catalogue of detected sources therefore contains the position of the clusters, their S/Ns, and the estimated $\theta_{500}$ and $Y_{5R500}$. In addition, the algorithm returns an estimate of the noise for each patch and filter size. It is computed directly from the filtered maps, assuming the SZ effect to be small compared to the other astrophysical signals.
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Compton-y profiles obtained from the gNFW profiles of Table 1, in units of $Y_{500}/R_{500}^2$, as a function of normalized radius, $R/R_{500}$. In orange: profile from Arnaud et al. (2010). In green: profile from Planck Collaboration Int. V (2013); in pink: profile from Pointecouteau et al. (2021). In red: profile from Tramonte et al. (2023). In purple: “Peaked” profile. In light blue: “SimFit” profile.</figcaption></figure>

We chose to use real Planck maps – following Planck Collaboration XXIX (2014) – to ensure the most realistic setting possible for the completeness analysis, including all sources of noise and contaminations present during the real detection process. To this end, we chose the maps that were used for the detection of the second Planck SZ cluster catalogue (PSZ2, Planck Collaboration XXVII 2016), as these contain the original cluster sample whose completeness we aim to characterise.

Given the fact that the MMF algorithm estimates the noise directly from the input maps, the injection of simulated cluster signals in addition to the real ones already present could change the properties of the noise, and therefore impact the detection. For this reason, we first subjected the maps to a cleaning process in order to try to remove the SZ signal from the real clusters. The way a given cluster is removed from the maps is the following: starting from the cluster's integrated SZ flux and angular scale, ($Y_{5R500}$, $\theta_{500}$), a circular image of the cluster SZ emission at the frequencies of the different maps is computed in the same way as in Sect. 2.1.3, that is, using the profile from Arnaud et al. (2010). This cluster emission is then subtracted from the maps at the position of the original detected cluster. The differences between the circular template and the real cluster signal leave a residual contribution, but its impact on the noise estimation is certainly smaller than that of the original cluster and is therefore negligible.

The cleaning proceeds in two steps: the first step consists in simply removing all the clusters from the PSZ2 catalogue, obtaining a first “cleaned” version of the maps. In the second step, the MMF detection is run on the new maps to check for additional signals identified as clusters beyond those contained in the PSZ2 catalogue, with a lower limit of 4.25 in signal-to-noise ratio (S/N). All new detections obtained in this way are then also removed from the frequency maps, which are now, in principle, free from any relevant SZ source down to $S/N \sim 4.5$. These are the final cleaned maps we use for the completeness analysis.

Associated with the clean frequency maps, we create a mask that covers the emission from the Galaxy and the Magellanic Cloud, as well as an area of five times the beam size around point

sources (from PCCS2 catalogue Planck Collaboration XXVI 2016) and regions of CO emission Planck Collaboration XIII (2014). The final unmasked sky fraction is about 78%.

# 3. Methods

## 3.1. MMF

The matched multi-frequency filter (MMF) algorithm (Herranz et al. 2002; Melin et al. 2006) is a commonly used algorithm for the detection of galaxy clusters through their SZ signal. Notably, it has been used to produce cluster catalogues from Planck data, along with ACT and SPT data (Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016; Hilton et al. 2021; Bleem et al. 2015). This method is designed to combine the prior knowledge of the SZ signal from galaxy clusters, namely the spectral signature and spatial characteristics, to produce an optimal filter that returns the maximal S/N in the presence of a galaxy cluster. While the spectral shape of the SZ signal is well known and universal (that is in the non-relativistic regime), one has to make some approximations for the spatial filter; the usual choice is to assume spherical symmetry for clusters, and model the radial pressure profile according to the average profile of observed clusters. However, galaxy clusters are of course known to be aspherical in general, and rather triaxial in shape, and most have even more complex morphologies due to anisotropic accretion, mergers, and astrophysical processes such as shocks and feedback mechanisms. All this variety in morphological features prevents a perfect match between the cluster SZ signal and the spatial template used, which, in principle, impacts the detection performance. This effect is, in practice, reduced by the smoothing induced by the instrument's beam, which tends to symmetrise the signal, especially for clusters with scales comparable to or smaller than the beam size.

More specifically, our implementation of the MMF algorithm is very similar to that used in the construction of the Planck SZ cosmological catalogue (Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016), namely MMF3. This latter uses a projected gNFW profile as a cluster model, with parameters from Arnaud et al. (2010), [$c_{500}$, $\gamma$, $\alpha$, $\beta$] = [1.177, 0.308, 1.051, 5.491], and two additional free parameters: the signal amplitude $Y_{5R500}$ and the cluster scale $\theta_{500}$. The spatial filters are therefore constructed from this model, varying the cluster size on a grid of 40 logarithmically spaced points from $\theta_{500} = 1.059$ to 41.195 arcmin. When run, the MMF algorithm first divides the Planck full-sky maps into 546 overlapping square patches of 10° on a side. Each patch is then filtered with the templates, resulting in a set of S/N maps. The peaks in these maps over a S/N threshold of 3 are our candidate cluster detections. The integrated cluster signal $Y_{5R500}$ is given by the filtered map at the peak position, while the cluster size is assumed to be the scale of the filter that maximises the S/N at the cluster location. Finally, the clusters are merged in a single full-sky catalogue by merging all detections with lower S/N that fall inside the $\theta_{500}$ of a certain detected cluster.

The final catalogue of detected sources therefore contains the position of the clusters, their S/Ns, and the estimated $\theta_{500}$ and $Y_{5R500}$. In addition, the algorithm returns an estimate of the noise for each patch and filter size. It is computed directly from the filtered maps, assuming the SZ effect to be small compared to the other astrophysical signals.
## 3.2. Completeness

The completeness is the probability that a cluster with a certain ($Y_{5R500}$, $\theta_{500}$) in the true population will be detected, $P(d | Y_{5R500}$, $\theta_{500}$), given the survey and detection method one is considering (in our case, Planck and MMF, respectively). A first approximation, as explained in Planck Collaboration XXIX (2014), is to assume the noise on the Compton-$y$ parameter as Gaussian. In this case, the completeness can be determined analytically, and takes the form of an error function:

$$P(d | Y_{5R500}, \sigma_{Yi}(\theta_{500}), q) = \frac{1}{2} \left[ 1 + \operatorname{erf} \left( \frac{Y_{5R500} - q \sigma_{Yi}(\theta_{500})}{\sqrt{2} \sigma_{Yi}(\theta_{500})} \right) \right], \quad (6)$$

which depends on the integrated SZ signal, $Y_{5R500}$, the noise in the filtered maps at the scale $\theta_{500}$ in a given patch $i$, $\sigma_{Yi}(\theta_{500})$, and the threshold in S/N, $q$, over which a detection is accepted; here we use $q = 4.5$. The noise per patch per filter used for the ERF completeness is estimated here from the cleaned Planck maps in order to avoid any spurious contribution from the injected cluster signals (in the spirit of Zubeldia et al. 2023). Equation (6) is valid for one patch, and so to obtain the completeness for the full sky we average it over all patches, weighted by the area of the sky in each patch that is not covered by the mask.

A more direct way to estimate the completeness function is via Monte Carlo injection of simulated clusters directly in the sky maps (Melin et al. 2005; Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016). By injecting the simulated cluster images described in Sect. 2.1 into the Planck maps, and then running the MMF on these new maps, we are able to directly compare the “true” cluster set with the catalogue of detections returned by the algorithm. This method allows us to estimate the completeness considering all the possible sources of bias present in the detection (e.g. algorithmic effects, non-Gaussian noise, cluster morphology, etc.), which can be difficult to model analytically.

In order to obtain a sufficient number of detections to properly sample the completeness in the ($Y_{5R500}$, $\theta_{500}$) plane without altering the noise properties of the Planck maps, we create several realisations of injected sky maps, each containing 2000 cluster images. In this way, the average noise of the injected maps differs from that of the cleaned maps by less than 5% (in the scales of interest, below 10 arcmin).

We created 50 mock sky maps for the simulation images set and for the circular images with the Standard profile (Arnaud et al. 2010), while we created only 10 maps for each of the other sets with different gNFW profiles. For each mock map, we randomly selected 2000 cluster images, with their associated $Y_{5R500}$ and $\theta_{500}$. The images were then injected into randomly chosen positions in the cleaned Planck maps, uniformly distributed outside the Galactic and point-source mask, avoiding overlaps with other injected clusters. This results in there being, on average, about four clusters per ($10^\circ \times 10^\circ$) patch. The MMF detection algorithm was then run on the maps, obtaining a catalogue of candidate detections for each, with candidate positions, S/Ns, $\theta_{500}$, and $Y_{5R500}$. A threshold of $S/N > 4.5$ is imposed on the catalogues, the same as for the full Planck PSZ2 catalogue. We chose this threshold for all the completeness tests because of the larger number of detected clusters, which allows for better statistics for the completeness. Choosing a higher S/N threshold (e.g. $S/N > 6$, as for the Planck cosmology sample) yields equivalent results for the completeness, which justifies our choice a posteriori.

Comparing these catalogues with the respective catalogues of injected sources, we determined which clusters among the injected ones are correctly detected by the algorithm. The matching between the input and output catalogues is done in the following way: a cluster is considered detected if there is an entry in the MMF output catalogue within a radius of 5 arcmin around the position of the injected cluster. If there are multiple detections in that area, the one closest to the true position is chosen. Otherwise, the cluster is marked as undetected. This procedure is consistent with that of the Planck Collaboration XXIX (2014), Planck Collaboration XXVII (2016). The resulting single-sky matched catalogues, together with the detection state of their entries, are then stacked into a single catalogue, which is used to compute the completeness. This step is necessary to obtain good statistics in all ($Y_{5R500}$, $\theta_{500}$) bins.

Finally, the completeness for each cluster set is computed as the ratio between the number of detected and injected clusters in logarithmically spaced bins of ($Y_{5R500}$, $\theta_{500}$). If there are less than ten clusters in a bin, we deem it unreliable and discard it. The error bars on the completeness are computed via bootstrap resampling of the matched detections table<fnref n="2" />.

# 4. Results

Here, we present the completeness function estimated with the Monte Carlo injection method using the set of cluster images from the IllustrisTNG simulation, and discuss its departure from the analytical ERF completeness with the help of the other sets of circular cluster images with varying profiles, described in Sect. 2.1.3.

## 4.1. Completeness with the simulation images

The completeness function obtained from the simulation images set is shown in Fig. 3, as a function of $Y_{5R500}$, in six $\theta_{500}$ bins. It is compared with the analytical ERF completeness, and with the one obtained with the “Standard gNFW” images set that uses the Arnaud et al. (2010) profile (the same profile used to build the detection template). This set of circular images is meant to act as a benchmark for the injection method, and also as a test for the ideal case in which the cluster images match the detection template almost perfectly. Therefore, a priori, the matched filter should yield the maximum response for these images; and, as a consequence, the maximum completeness. We therefore expect the simulation images – which in addition to having a variety of different profiles also show complex morphologies that depart from spherical symmetry – to be more difficult to detect, which would mean a lower completeness.

We see that this prediction is not reflected in the results of Fig. 3. While the completeness of the “standard” images tends to agree well with the analytical ERF (which is in agreement with the results of Planck Collaboration XXIX 2014), the simulation images produce a completeness that is almost always higher than both other estimates, with differences of as high as 0.4 in some bins. This means that the simulation cluster images, despite their imperfect match with the detection template, have a higher detection probability than those for which the match is near perfect.

To help us understand this result, we can imagine splitting the contributions of the simulation images to the completeness

<footnote n="2">² This approach gives a useful estimation of the statistical error on the Monte Carlo mean completeness, but probably underestimates other sources of uncertainty, like the sky-by-sky variations.</footnote>
## 3.2. Completeness

The completeness is the probability that a cluster with a certain ($Y_{5R500}$, $\theta_{500}$) in the true population will be detected, $P(d | Y_{5R500}$, $\theta_{500}$), given the survey and detection method one is considering (in our case, Planck and MMF, respectively). A first approximation, as explained in Planck Collaboration XXIX (2014), is to assume the noise on the Compton-$y$ parameter as Gaussian. In this case, the completeness can be determined analytically, and takes the form of an error function:

$$P(d | Y_{5R500}, \sigma_{Yi}(\theta_{500}), q) = \frac{1}{2} \left[ 1 + \operatorname{erf} \left( \frac{Y_{5R500} - q \sigma_{Yi}(\theta_{500})}{\sqrt{2} \sigma_{Yi}(\theta_{500})} \right) \right], \quad (6)$$

which depends on the integrated SZ signal, $Y_{5R500}$, the noise in the filtered maps at the scale $\theta_{500}$ in a given patch $i$, $\sigma_{Yi}(\theta_{500})$, and the threshold in S/N, $q$, over which a detection is accepted; here we use $q = 4.5$. The noise per patch per filter used for the ERF completeness is estimated here from the cleaned Planck maps in order to avoid any spurious contribution from the injected cluster signals (in the spirit of Zubeldia et al. 2023). Equation (6) is valid for one patch, and so to obtain the completeness for the full sky we average it over all patches, weighted by the area of the sky in each patch that is not covered by the mask.

A more direct way to estimate the completeness function is via Monte Carlo injection of simulated clusters directly in the sky maps (Melin et al. 2005; Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016). By injecting the simulated cluster images described in Sect. 2.1 into the Planck maps, and then running the MMF on these new maps, we are able to directly compare the “true” cluster set with the catalogue of detections returned by the algorithm. This method allows us to estimate the completeness considering all the possible sources of bias present in the detection (e.g. algorithmic effects, non-Gaussian noise, cluster morphology, etc.), which can be difficult to model analytically.

In order to obtain a sufficient number of detections to properly sample the completeness in the ($Y_{5R500}$, $\theta_{500}$) plane without altering the noise properties of the Planck maps, we create several realisations of injected sky maps, each containing 2000 cluster images. In this way, the average noise of the injected maps differs from that of the cleaned maps by less than 5% (in the scales of interest, below 10 arcmin).

We created 50 mock sky maps for the simulation images set and for the circular images with the Standard profile (Arnaud et al. 2010), while we created only 10 maps for each of the other sets with different gNFW profiles. For each mock map, we randomly selected 2000 cluster images, with their associated $Y_{5R500}$ and $\theta_{500}$. The images were then injected into randomly chosen positions in the cleaned Planck maps, uniformly distributed outside the Galactic and point-source mask, avoiding overlaps with other injected clusters. This results in there being, on average, about four clusters per (10° × 10°) patch. The MMF detection algorithm was then run on the maps, obtaining a catalogue of candidate detections for each, with candidate positions, S/Ns, $\theta_{500}$, and $Y_{5R500}$. A threshold of $S/N > 4.5$ is imposed on the catalogues, the same as for the full Planck PSZ2 catalogue. We chose this threshold for all the completeness tests because of the larger number of detected clusters, which allows for better statistics for the completeness. Choosing a higher S/N threshold (e.g. $S/N > 6$, as for the Planck cosmology sample) yields equivalent results for the completeness, which justifies our choice a posteriori.

Comparing these catalogues with the respective catalogues of injected sources, we determined which clusters among the injected ones are correctly detected by the algorithm. The matching between the input and output catalogues is done in the following way: a cluster is considered detected if there is an entry in the MMF output catalogue within a radius of 5 arcmin around the position of the injected cluster. If there are multiple detections in that area, the one closest to the true position is chosen. Otherwise, the cluster is marked as undetected. This procedure is consistent with that of the Planck Collaboration XXIX (2014), Planck Collaboration XXVII (2016). The resulting single-sky matched catalogues, together with the detection state of their entries, are then stacked into a single catalogue, which is used to compute the completeness. This step is necessary to obtain good statistics in all ($Y_{5R500}$, $\theta_{500}$) bins.

Finally, the completeness for each cluster set is computed as the ratio between the number of detected and injected clusters in logarithmically spaced bins of ($Y_{5R500}$, $\theta_{500}$). If there are less than ten clusters in a bin, we deem it unreliable and discard it. The error bars on the completeness are computed via bootstrap resampling of the matched detections table<fnref n="2" />.

# 4. Results

Here, we present the completeness function estimated with the Monte Carlo injection method using the set of cluster images from the IllustrisTNG simulation, and discuss its departure from the analytical ERF completeness with the help of the other sets of circular cluster images with varying profiles, described in Sect. 2.1.3.

## 4.1. Completeness with the simulation images

The completeness function obtained from the simulation images set is shown in Fig. 3, as a function of $Y_{5R500}$, in six $\theta_{500}$ bins. It is compared with the analytical ERF completeness, and with the one obtained with the “Standard gNFW” images set that uses the Arnaud et al. (2010) profile (the same profile used to build the detection template). This set of circular images is meant to act as a benchmark for the injection method, and also as a test for the ideal case in which the cluster images match the detection template almost perfectly. Therefore, a priori, the matched filter should yield the maximum response for these images; and, as a consequence, the maximum completeness. We therefore expect the simulation images – which in addition to having a variety of different profiles also show complex morphologies that depart from spherical symmetry – to be more difficult to detect, which would mean a lower completeness.

We see that this prediction is not reflected in the results of Fig. 3. While the completeness of the “standard” images tends to agree well with the analytical ERF (which is in agreement with the results of Planck Collaboration XXIX 2014), the simulation images produce a completeness that is almost always higher than both other estimates, with differences of as high as 0.4 in some bins. This means that the simulation cluster images, despite their imperfect match with the detection template, have a higher detection probability than those for which the match is near perfect.

To help us understand this result, we can imagine splitting the contributions of the simulation images to the completeness

<footnote n="2">² This approach gives a useful estimation of the statistical error on the Monte Carlo mean completeness, but probably underestimates other sources of uncertainty, like the sky-by-sky variations.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Completeness as a function of $Y_{5R500}$ in six $	heta_{500}$ bins, computed with Monte Carlo injection for the simulation and circular sets of clusters (in blue and orange, respectively), and estimated with the ERF approximation, in green. The shaded areas represent the uncertainty on the completeness in each bin, computed via bootstrap resampling.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Average Compton-$y$ profile of simulation images (in blue, shaded area corresponds to standard deviation of profile sample), compared with the $y$ profile obtained from the Arnaud et al. (2010) pressure profile (in orange) and the “peaked” profile described in the text (in purple), in units of $R_{500}$ and $Y_{500}/R_{500}^2$.</figcaption></figure>

into two parts: the average $y$ profile and the non-spherical morphology of simulated clusters. In the following, we present a study of the impact of these two aspects on the completeness.

## 4.2. Impact of cluster profile

In Fig. 4, we compare the average $y$ profile of the simulation images with the integrated Arnaud et al. (2010) profile. We can see how the average simulation profile tends to be overall steeper than that from Arnaud et al. (2010), especially in the outer part from ~0.7 $R_{500}$, and is generally more concentrated, or peaked,

than the assumed profile. This trend could be the cause of the higher completeness of the simulation images compared to the Standard gNFW ones, given that the MMF algorithm, as it tries to find the optimal parameter to fit its flatter profile template to the cluster signal, might favour a smaller radius than the real one. This underestimation of the radius leads to an increase in S/N, given that the MMF noise estimate increases with filter radius. This picture is confirmed in the bottom left panel of Fig. 5, where the bias in the detected cluster radius is evident, with a clear underestimation, and a median ratio between detected and injected $\theta_{500}$ of the order of ~0.8. This bias mostly disappears when analysing the detected $\theta_{500}$ of the Standard gNFW set (with only about 1% median difference), whose profile is the same as that of the MMF template. This indicates that the performance of the detection process depends on the assumed cluster profile in the template. Looking at the distribution of the detected SZ signal $Y_{5R500}$ of the two sets of cluster images (Fig. 5, top row), we find that both tend to be overestimated with respect to the injected quantities, but with less difference between the two sets. We find an overestimation of about 9% for the Standard gNFW set, which is consistent with the result of Planck Collaboration XXIX (2014), while the Simulation set has a median overestimation of ~25%, which leads to an increased S/N that in turn contributes to the increase in completeness.

To test the impact of a different mean profile on the completeness, we used the Peaked gNFW set described in Sect. 2.1.3, which is built from a pressure profile that maintains the same parameters as the Arnaud et al. (2010) profile, except for the concentration parameter, $c_{500}$ = 1.5. This higher parameter determines a $y$ profile that is higher in the centre and lower around $R_{500}$ than the Arnaud et al. (2010) profile, and roughly reproduces the average profile of the simulated clusters at large radii, as can be seen in Fig. 4. Using this set of images, we compute the completeness via Monte Carlo injections again, and
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Completeness as a function of $Y_{5R500}$ in six $	heta_{500}$ bins, computed with Monte Carlo injection for the simulation and circular sets of clusters (in blue and orange, respectively), and estimated with the ERF approximation, in green. The shaded areas represent the uncertainty on the completeness in each bin, computed via bootstrap resampling.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Average Compton-$y$ profile of simulation images (in blue, shaded area corresponds to standard deviation of profile sample), compared with the $y$ profile obtained from the Arnaud et al. (2010) pressure profile (in orange) and the “peaked” profile described in the text (in purple), in units of $R_{500}$ and $Y_{500}/R_{500}^2$.</figcaption></figure>

into two parts: the average $y$ profile and the non-spherical morphology of simulated clusters. In the following, we present a study of the impact of these two aspects on the completeness.

## 4.2. Impact of cluster profile

In Fig. 4, we compare the average $y$ profile of the simulation images with the integrated Arnaud et al. (2010) profile. We can see how the average simulation profile tends to be overall steeper than that from Arnaud et al. (2010), especially in the outer part from ~0.7 $R_{500}$, and is generally more concentrated, or peaked,

than the assumed profile. This trend could be the cause of the higher completeness of the simulation images compared to the Standard gNFW ones, given that the MMF algorithm, as it tries to find the optimal parameter to fit its flatter profile template to the cluster signal, might favour a smaller radius than the real one. This underestimation of the radius leads to an increase in S/N, given that the MMF noise estimate increases with filter radius. This picture is confirmed in the bottom left panel of Fig. 5, where the bias in the detected cluster radius is evident, with a clear underestimation, and a median ratio between detected and injected $\theta_{500}$ of the order of ~0.8. This bias mostly disappears when analysing the detected $\theta_{500}$ of the Standard gNFW set (with only about 1% median difference), whose profile is the same as that of the MMF template. This indicates that the performance of the detection process depends on the assumed cluster profile in the template. Looking at the distribution of the detected SZ signal $Y_{5R500}$ of the two sets of cluster images (Fig. 5, top row), we find that both tend to be overestimated with respect to the injected quantities, but with less difference between the two sets. We find an overestimation of about 9% for the Standard gNFW set, which is consistent with the result of Planck Collaboration XXIX (2014), while the Simulation set has a median overestimation of ~25%, which leads to an increased S/N that in turn contributes to the increase in completeness.

To test the impact of a different mean profile on the completeness, we used the Peaked gNFW set described in Sect. 2.1.3, which is built from a pressure profile that maintains the same parameters as the Arnaud et al. (2010) profile, except for the concentration parameter, $c_{500}$ = 1.5. This higher parameter determines a $y$ profile that is higher in the centre and lower around $R_{500}$ than the Arnaud et al. (2010) profile, and roughly reproduces the average profile of the simulated clusters at large radii, as can be seen in Fig. 4. Using this set of images, we compute the completeness via Monte Carlo injections again, and
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Comparison of real cluster properties vs detected ones. The top row shows the total integrated SZ flux, $Y_{5R500}$, while the bottom row shows the cluster radius, $\theta_{500}$. The columns show two different sets of cluster images: simulation images and circular images with the Arnaud et al. (2010) profile.</figcaption></figure>

show the results in Fig. 6 for two $\theta_{500}$ bins as an example. The completeness computed using the images with the Peaked profile appears to agree more closely with the one obtained with the simulation images, in particular at larger cluster scales, while at low $\theta_{500}$ it appears to be more in between the simulation and the ERF completeness. This is confirmation that a change in the shape of the profile can lead to an increased detection probability, despite the imperfect match with the filter template.

To probe this result further, we tested the impact of using spherically symmetric images with three different gNFW profiles obtained by different studies, fitting the gNFW profile to real observed cluster samples. The chosen works are: Planck Collaboration Int. V (2013), where the profile was fitted on 62 SZ-selected clusters using Compton-$y$ data from Planck and X-ray observations from XMM-Newton; Pointecouteau et al. (2021), where the SZ signal from the combined map of Planck and ACT (Aghanim et al. 2019) of 31 clusters was used to constrain the gNFW parameters; and finally, Tramonte et al. (2023), where the authors fitted the $y$ profile obtained by stacking the Planck SZ signal of a large number of clusters from different surveys, in different redshift bins. From this latter study, we take the parameters obtained from 4421 clusters in the $z < 0.35$ redshift bin, which overlaps with our simulation sample. Of these profiles, shown in Fig. 7, the first two are flatter than the Arnaud et al. (2010) profile (with the second showing the greater difference), while the third is flatter in the inner cluster region (up to ~0.1 $R_{500}$) and then gets steeper than all the other profiles. With this choice, we tested the impact of the scatter in the observed profiles in both the flatter and steeper directions.

The results of this test are shown in Fig. 8 for two cluster scales. The effect of using different profile parameters is consistent with our previous results: a steeper profile tends to increase the completeness at all scales, while a flatter profile

<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Same as in Fig. 3, comparing the completeness obtained from the simulation images (blue) with that from the Peaked profile images (purple) and the ERF completeness (green).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Compton-$y$ profiles in units of $Y_{500}/R_{500}^2$ as a function of normalized radius, $R/R_{500}$. Blue: mean $y$ profile of the simulation images. Orange: profile from Arnaud et al. (2010). Green: profile from Planck Collaboration Int. V (2013). Pink: profile from Pointecouteau et al. (2021). Red: profile from Tramonte et al. (2023).</figcaption></figure>

generally decreases it. More specifically, we note that the profile of Planck Collaboration Int. V (2013), which most closely resembles the detection template, gives a completeness that is comparable with the theoretical estimate, although it tends to be slightly lower for the larger cluster scale. On a similar note,
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Comparison of real cluster properties vs detected ones. The top row shows the total integrated SZ flux, $Y_{5R500}$, while the bottom row shows the cluster radius, $\theta_{500}$. The columns show two different sets of cluster images: simulation images and circular images with the Arnaud et al. (2010) profile.</figcaption></figure>

show the results in Fig. 6 for two $\theta_{500}$ bins as an example. The completeness computed using the images with the Peaked profile appears to agree more closely with the one obtained with the simulation images, in particular at larger cluster scales, while at low $\theta_{500}$ it appears to be more in between the simulation and the ERF completeness. This is confirmation that a change in the shape of the profile can lead to an increased detection probability, despite the imperfect match with the filter template.

To probe this result further, we tested the impact of using spherically symmetric images with three different gNFW profiles obtained by different studies, fitting the gNFW profile to real observed cluster samples. The chosen works are: Planck Collaboration Int. V (2013), where the profile was fitted on 62 SZ-selected clusters using Compton-$y$ data from Planck and X-ray observations from XMM-Newton; Pointecouteau et al. (2021), where the SZ signal from the combined map of Planck and ACT (Aghanim et al. 2019) of 31 clusters was used to constrain the gNFW parameters; and finally, Tramonte et al. (2023), where the authors fitted the $y$ profile obtained by stacking the Planck SZ signal of a large number of clusters from different surveys, in different redshift bins. From this latter study, we take the parameters obtained from 4421 clusters in the $z < 0.35$ redshift bin, which overlaps with our simulation sample. Of these profiles, shown in Fig. 7, the first two are flatter than the Arnaud et al. (2010) profile (with the second showing the greater difference), while the third is flatter in the inner cluster region (up to ~0.1 $R_{500}$) and then gets steeper than all the other profiles. With this choice, we tested the impact of the scatter in the observed profiles in both the flatter and steeper directions.

The results of this test are shown in Fig. 8 for two cluster scales. The effect of using different profile parameters is consistent with our previous results: a steeper profile tends to increase the completeness at all scales, while a flatter profile

<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Same as in Fig. 3, comparing the completeness obtained from the simulation images (blue) with that from the Peaked profile images (purple) and the ERF completeness (green).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Compton-$y$ profiles in units of $Y_{500}/R_{500}^2$ as a function of normalized radius, $R/R_{500}$. Blue: mean $y$ profile of the simulation images. Orange: profile from Arnaud et al. (2010). Green: profile from Planck Collaboration Int. V (2013). Pink: profile from Pointecouteau et al. (2021). Red: profile from Tramonte et al. (2023).</figcaption></figure>

generally decreases it. More specifically, we note that the profile of Planck Collaboration Int. V (2013), which most closely resembles the detection template, gives a completeness that is comparable with the theoretical estimate, although it tends to be slightly lower for the larger cluster scale. On a similar note,
<figure><img src="image_8.png" /><figcaption>Fig. 8. Same as in Fig. 3, comparing the completeness from images with three different observed profiles. Green: Planck Collaboration Int. V (2013) profile. Pink: Pointecouteau et al. (2021) profile. Red: Tramonte et al. (2023) profile.</figcaption></figure>

the Pointecouteau et al. (2021) profile, the flattest of the three we use, has remarkably lower completeness than the ERF curve. On the other hand, the profile from Tramonte et al. (2023) gives a completeness that is generally higher than the ERF, the same effect observed with the simulation set. This further confirms the hypothesis that the shape of the pressure profile and the estimated completeness are related, in particular pointing towards a greater importance of the slope outside the very centre of the cluster, although this possibility requires further evaluation.

### 4.3. Impact of cluster asymmetry

We also studied the impact of the departure from spherical symmetry of the cluster images on the completeness. Following Gouin et al. (2022), we quantified the degree of asphericity using a method based on 2D multipole moment decomposition (first introduced by Schneider & Bartelmann 1997). The method consists of a decomposition of the projected cluster gas distribution $\Sigma(R, \phi)$ in harmonic modes $m$, integrated in a radial aperture:

$$Q_m(\Delta R) = \int_{R_{\min}}^{R_{\max}} \int_{0}^{2\pi} R \, dR \, d\phi \, e^{im\phi} \Sigma(R, \phi). \quad (7)$$

Each mode highlights different patterns in the gas distribution, with higher $m$ corresponding to smaller angular scales. Then, to

quantify the relevance of each mode with respect to the spherical symmetry, the multipolar ratio $\beta_m$ is defined (Gouin et al. 2022) as

$$\beta_m = \frac{|Q_m|}{|Q_0|} \quad (8)$$

Here, we applied this method to the Compton-$y$ images from the simulation (before the convolution with the Planck beams) in the radial aperture inside $\theta_{500}$. In particular, we concentrated on the $\beta_2$ multipolar ratio, which has been shown to correlate with the ellipticity of the cluster gas and is expected to be the leading non-spherical mode inside clusters (Gouin et al. 2022). We therefore decided to use it as a proxy for cluster asphericity.

After obtaining the value of $\beta_2$ for each simulation image, we selected two groups of images, the “more elliptical” and “more spherical” subsets, formed by the images with the 25% highest and 25% lowest values of $\beta_2$, respectively, for a total of 2231 images per group. We then computed the completeness for the two subsets, which are shown in Fig. 9 with the analytical ERF completeness for reference.

Comparing the completeness of the more and less spherical images, we see no appreciable difference between the two curves for cluster sizes below ~6 arcmin. This is expected, given that the Planck beam size is indeed about 6 arcmin on average (and goes up to 10 arcmin for the 100 GHz channel), effectively erasing most smaller-scale morphological differences and symmetrising the images. Above 6 arcmin, instead, we start to see a small difference between the two subsets, with the more elliptical images showing, on average, lower completeness than the more spherical ones, and the difference appears to become more important with larger cluster radii. This trend can be more clearly visualised in Fig. 10, where the difference between “more elliptical” and “more spherical” image completeness is shown in bins of $\theta_{500}$ and $Y_{5R500}$. Indeed, the completeness difference tends to increase with increasing $\theta_{500}$, a sign that for resolved clusters the morphology has an impact on the detection probability (i.e. the completeness); although this impact remains moderate, at least in the range of scales probed by our cluster sample.

A further test was performed to reveal the impact of realistic cluster morphologies on the completeness estimation, comparing the simulation set with a set of circular images with a gNFW profile fitted to the average profile from the simulation (named SimFit profile in Table 1)<fnref n="3" />. The two completeness functions are shown in Fig. 11. We see that the two curves overlap almost perfectly, suggesting that the $y$ profile is the dominant driver of completeness, while the small effect of the cluster morphology remarked in the previous paragraph is no longer visible, and is probably “washed out” by considering the full cluster population. Therefore, we are led to believe that, in the context of Planck and in the range of scales probed by our simulation cluster set, the departure of clusters from spherical symmetry does not induce an appreciable bias in the completeness function. Nevertheless, we show that there is an effect related to cluster morphology when separating the “most elliptical” images from the “most spherical” ones. This effect could be larger in a context where, for example, the beam size is smaller, which will be interesting to check in future work.

<footnote n="3">³ We checked for mass and redshift dependence of the Simulation set $y$ profiles and found only minor deviations much smaller than the scatter of the full sample.</footnote>
<figure><img src="image_8.png" /><figcaption>Fig. 8. Same as in Fig. 3, comparing the completeness from images with three different observed profiles. Green: Planck Collaboration Int. V (2013) profile. Pink: Pointecouteau et al. (2021) profile. Red: Tramonte et al. (2023) profile.</figcaption></figure>

the Pointecouteau et al. (2021) profile, the flattest of the three we use, has remarkably lower completeness than the ERF curve. On the other hand, the profile from Tramonte et al. (2023) gives a completeness that is generally higher than the ERF, the same effect observed with the simulation set. This further confirms the hypothesis that the shape of the pressure profile and the estimated completeness are related, in particular pointing towards a greater importance of the slope outside the very centre of the cluster, although this possibility requires further evaluation.

### 4.3. Impact of cluster asymmetry

We also studied the impact of the departure from spherical symmetry of the cluster images on the completeness. Following Gouin et al. (2022), we quantified the degree of asphericity using a method based on 2D multipole moment decomposition (first introduced by Schneider & Bartelmann 1997). The method consists of a decomposition of the projected cluster gas distribution $\Sigma(R, \phi)$ in harmonic modes $m$, integrated in a radial aperture:

$$Q_m(\Delta R) = \int_{R_{\min}}^{R_{\max}} \int_{0}^{2\pi} R \, dR \, d\phi \, e^{im\phi} \Sigma(R, \phi). \quad (7)$$

Each mode highlights different patterns in the gas distribution, with higher $m$ corresponding to smaller angular scales. Then, to

quantify the relevance of each mode with respect to the spherical symmetry, the multipolar ratio $\beta_m$ is defined (Gouin et al. 2022) as

$$\beta_m = \frac{|Q_m|}{|Q_0|} \quad (8)$$

Here, we applied this method to the Compton-$y$ images from the simulation (before the convolution with the Planck beams) in the radial aperture inside $\theta_{500}$. In particular, we concentrated on the $\beta_2$ multipolar ratio, which has been shown to correlate with the ellipticity of the cluster gas and is expected to be the leading non-spherical mode inside clusters (Gouin et al. 2022). We therefore decided to use it as a proxy for cluster asphericity.

After obtaining the value of $\beta_2$ for each simulation image, we selected two groups of images, the “more elliptical” and “more spherical” subsets, formed by the images with the 25% highest and 25% lowest values of $\beta_2$, respectively, for a total of 2231 images per group. We then computed the completeness for the two subsets, which are shown in Fig. 9 with the analytical ERF completeness for reference.

Comparing the completeness of the more and less spherical images, we see no appreciable difference between the two curves for cluster sizes below ~6 arcmin. This is expected, given that the Planck beam size is indeed about 6 arcmin on average (and goes up to 10 arcmin for the 100 GHz channel), effectively erasing most smaller-scale morphological differences and symmetrising the images. Above 6 arcmin, instead, we start to see a small difference between the two subsets, with the more elliptical images showing, on average, lower completeness than the more spherical ones, and the difference appears to become more important with larger cluster radii. This trend can be more clearly visualised in Fig. 10, where the difference between “more elliptical” and “more spherical” image completeness is shown in bins of $\theta_{500}$ and $Y_{5R500}$. Indeed, the completeness difference tends to increase with increasing $\theta_{500}$, a sign that for resolved clusters the morphology has an impact on the detection probability (i.e. the completeness); although this impact remains moderate, at least in the range of scales probed by our cluster sample.

A further test was performed to reveal the impact of realistic cluster morphologies on the completeness estimation, comparing the simulation set with a set of circular images with a gNFW profile fitted to the average profile from the simulation (named SimFit profile in Table 1)$^3$. The two completeness functions are shown in Fig. 11. We see that the two curves overlap almost perfectly, suggesting that the $y$ profile is the dominant driver of completeness, while the small effect of the cluster morphology remarked in the previous paragraph is no longer visible, and is probably “washed out” by considering the full cluster population. Therefore, we are led to believe that, in the context of Planck and in the range of scales probed by our simulation cluster set, the departure of clusters from spherical symmetry does not induce an appreciable bias in the completeness function. Nevertheless, we show that there is an effect related to cluster morphology when separating the “most elliptical” images from the “most spherical” ones. This effect could be larger in a context where, for example, the beam size is smaller, which will be interesting to check in future work.

<footnote>$^3$ We checked for mass and redshift dependence of the Simulation set $y$ profiles and found only minor deviations much smaller than the scatter of the full sample.</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Same as in Fig. 3, comparing the completeness of two subsets of the simulation images. Images with the 25% highest $\beta_2$ (more elliptical) are shown in red, while images with the 25% lowest $\beta_2$ (more spherical) are shown in blue.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Difference in completeness between the “more spherical” and “more elliptical” subsets of the simulation images, in bins of ($\theta_{500}$, $Y_{5R500}$). The bins in blue are the ones in which the completeness from the “more spherical” subset is higher, while those in red are those in which the completeness of the “more elliptical” subset dominates.</figcaption></figure>

### 4.4. Impact on cluster count cosmology

Finally, we tested the impact of changing the completeness function in the cosmological analysis. We choose two completeness examples from the different sets we obtained, one above and one below the ERF estimate: the simulation completeness and the PACT gNFW completeness. By changing only the completeness function and not the other elements of the analysis (such as the scaling relations or mass bias), we are actually testing a rather extreme scenario, in which the true cluster profile is different not only from the one assumed in the MMF templates, but also in the scaling relation calibration and mass bias estimations. Therefore,

the results we obtain here are to be taken as indicative of a trend of impact rather than providing actual values.

To simplify the analysis, we chose to approximate the completeness obtained from Monte Carlo injection with an adapted version of the ERF completeness. We used a simple prescription to approximate the completeness: we modify the functional form of the ERF completeness in Eq. (6) by adding a free multiplicative parameter in front of the per-patch noise of the original Planck maps $\sigma_{Yi}^{\text{new}}(\theta_{500}) = a_{\text{fit}}\sigma_{Yi}(\theta_{500})$, and optimise this parameter to fit the new ERF completeness to the completeness obtained from Monte Carlo injection. The fitted value of $a_{\text{fit}}$ for the simulation completeness is $a_{\text{fit}}^{\text{sim}} = 1.27$, and that for the PACT completeness is $a_{\text{fit}}^{\text{PACT}} = 0.65$. We show an example of the fitted ERF completeness compared to those obtained with Monte Carlo in Fig. 12. Above, we use a threshold in S/N of 4.5 to compute the completeness, but the Planck cosmological sample has a threshold of $S/N > 6$. In any case, we verified that the values of the noise shift we computed are the same also with a higher S/N threshold, and so we can simply set $q = 6$ in Eq. (6) to get the ERF completeness for the cosmological sample.

We performed a cosmological analysis of cluster number counts, following Planck Collaboration XXIV (2016). We used the Planck PSZ2 MMF3 cosmological cluster sample (Planck Collaboration XXVII 2016), adding observational constraints from Big Bang nucleosynthesis (BBN; Steigman 2008), baryon acoustic oscillation (BAO) measurements from SDSS-III DR12 (Alam et al. 2017), and a prior on $n_s$ from Planck Collaboration XVI (2014). For the mass-observable scaling relations and mass-bias parameters, we used the same priors as the baseline analysis in Planck Collaboration XXIV (2016). The constraints on the cosmological parameters are obtained with a Markov-chain Monte Carlo (MCMC) implemented in the code CosmOMC<fnref n="4" /> (Lewis & Bridle 2002). We performed the

<footnote n="4">⁴ https://cosmologist.info/cosmomc/</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Same as in Fig. 3, comparing the completeness of two subsets of the simulation images. Images with the 25% highest $\beta_2$ (more elliptical) are shown in red, while images with the 25% lowest $\beta_2$ (more spherical) are shown in blue.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Difference in completeness between the “more spherical” and “more elliptical” subsets of the simulation images, in bins of ($\theta_{500}$, $Y_{5R500}$). The bins in blue are the ones in which the completeness from the “more spherical” subset is higher, while those in red are those in which the completeness of the “more elliptical” subset dominates.</figcaption></figure>

### 4.4. Impact on cluster count cosmology

Finally, we tested the impact of changing the completeness function in the cosmological analysis. We choose two completeness examples from the different sets we obtained, one above and one below the ERF estimate: the simulation completeness and the PACT gNFW completeness. By changing only the completeness function and not the other elements of the analysis (such as the scaling relations or mass bias), we are actually testing a rather extreme scenario, in which the true cluster profile is different not only from the one assumed in the MMF templates, but also in the scaling relation calibration and mass bias estimations. Therefore,

the results we obtain here are to be taken as indicative of a trend of impact rather than providing actual values.

To simplify the analysis, we chose to approximate the completeness obtained from Monte Carlo injection with an adapted version of the ERF completeness. We used a simple prescription to approximate the completeness: we modify the functional form of the ERF completeness in Eq. (6) by adding a free multiplicative parameter in front of the per-patch noise of the original Planck maps $\sigma_{Yi}^{\text{new}}(\theta_{500}) = a_{\text{fit}}\sigma_{Yi}(\theta_{500})$, and optimise this parameter to fit the new ERF completeness to the completeness obtained from Monte Carlo injection. The fitted value of $a_{\text{fit}}$ for the simulation completeness is $a_{\text{fit}}^{\text{sim}} = 1.27$, and that for the PACT completeness is $a_{\text{fit}}^{\text{PACT}} = 0.65$. We show an example of the fitted ERF completeness compared to those obtained with Monte Carlo in Fig. 12. Above, we use a threshold in S/N of 4.5 to compute the completeness, but the Planck cosmological sample has a threshold of $S/N > 6$. In any case, we verified that the values of the noise shift we computed are the same also with a higher S/N threshold, and so we can simply set $q = 6$ in Eq. (6) to get the ERF completeness for the cosmological sample.

We performed a cosmological analysis of cluster number counts, following Planck Collaboration XXIV (2016). We used the Planck PSZ2 MMF3 cosmological cluster sample (Planck Collaboration XXVII 2016), adding observational constraints from Big Bang nucleosynthesis (BBN; Steigman 2008), baryon acoustic oscillation (BAO) measurements from SDSS-III DR12 (Alam et al. 2017), and a prior on $n_s$ from Planck Collaboration XVI (2014). For the mass-observable scaling relations and mass-bias parameters, we used the same priors as the baseline analysis in Planck Collaboration XXIV (2016). The constraints on the cosmological parameters are obtained with a Markov-chain Monte Carlo (MCMC) implemented in the code CosmOMC<fnref n="4" /> (Lewis & Bridle 2002). We performed the

<footnote n="4">⁴ https://cosmologist.info/cosmomc/</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 11.** Same as in Fig. 3, comparing the completeness obtained from the simulation images (in blue) with that from spherical images with the SimFit profile (light blue).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 12.** Comparison of the three completeness functions used in the cosmological analyses with the functions computed via Monte Carlo injection of simulation images and spherical images with the Pointecouteau et al. (2021) profile.</figcaption></figure>

analysis using the original ERF completeness and the approximated Simulation and PACT completeness described above. The results of the three analyses are compared in Fig. 13, for $H_0$, $\Omega_m$, and $\sigma_8$, and we list the constraints in Table 2.

From these results, we note how the change in the completeness moves the constraints on $\Omega_m$ and $\sigma_8$. In particular, a higher

<figure><img src="image_13.png" /><figcaption>**Fig. 13.** Constraints on cosmological parameters ($H_0$, $\Omega_m$, $\sigma_8$) from the Planck PSZ2 cluster number counts and BAO data, with three different completeness functions: the “standard” ERF completeness (same as in Planck Collaboration XXIV 2016), and two “fitted” versions obtained by rescaling the noise per patch by a constant in order to reproduce the completeness derived from the simulation clusters and the PACT (Pointecouteau et al. 2021) profile clusters. The filled contours represent 68% and 95% confidence regions.</figcaption></figure>

<table><caption>**Table 2.** Constraints on cosmological parameters obtained from the analysis of Planck PSZ2 cluster number counts and BAO, using three different completeness functions, as described in the text.</caption><thead><tr><th></th><th>Simulation completeness</th><th>ERF completeness</th><th>PACT completeness</th></tr></thead><tbody><tr><td>$H_0$</td><td>68.4<sup>+2.0</sup><sub>-2.4</sub></td><td>69.1<sup>+2.1</sup><sub>-2.4</sub></td><td>70.3<sup>+2.2</sup><sub>-2.5</sub></td></tr><tr><td>$\Omega_m$</td><td>0.301<sup>+0.019</sup><sub>-0.023</sub></td><td>0.314<sup>+0.020</sup><sub>-0.024</sub></td><td>0.337 ± 0.024</td></tr><tr><td>$\sigma_8$</td><td>0.742<sup>+0.025</sup><sub>-0.033</sub></td><td>0.770<sup>+0.027</sup><sub>-0.036</sub></td><td>0.823<sup>+0.030</sup><sub>-0.039</sub></td></tr></tbody></table>

completeness with respect to the ERF (in this case, the Simulation completeness) favours lower values of the two cosmological parameters, while a lower completeness (that of PACT) favours higher values. For $\Omega_m$, the shift amounts to about $0.6\sigma$ for the Simulation completeness, and about $1.1\sigma$ for that of PACT. For $\sigma_8$, the two shifts are $0.9\sigma$ and $1.8\sigma$, respectively. The relation between the completeness function and these two cosmological parameters can be understood in the following terms: as the completeness function informs us about the fraction of detected clusters over the total, keeping the number of detected objects fixed with a higher (respectively, lower) completeness means that the total number of “true” clusters gets lower (higher). This, in turn, favours a lower (higher) value for both $\Omega_m$ and $\sigma_8$. It is interesting to note how the direction of the shift in the ($\Omega_m$, $\sigma_8$) plane is the same as moving the prior on the mass bias, as shown in Fig. 7 of Planck Collaboration XXIV (2016).
<figure><img src="image_11.png" /><figcaption>**Fig. 11.** Same as in Fig. 3, comparing the completeness obtained from the simulation images (in blue) with that from spherical images with the SimFit profile (light blue).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 12.** Comparison of the three completeness functions used in the cosmological analyses with the functions computed via Monte Carlo injection of simulation images and spherical images with the Pointecouteau et al. (2021) profile.</figcaption></figure>

analysis using the original ERF completeness and the approximated Simulation and PACT completeness described above. The results of the three analyses are compared in Fig. 13, for $H_0$, $\Omega_m$, and $\sigma_8$, and we list the constraints in Table 2.

From these results, we note how the change in the completeness moves the constraints on $\Omega_m$ and $\sigma_8$. In particular, a higher

<figure><img src="image_13.png" /><figcaption>**Fig. 13.** Constraints on cosmological parameters ($H_0$, $\Omega_m$, $\sigma_8$) from the Planck PSZ2 cluster number counts and BAO data, with three different completeness functions: the “standard” ERF completeness (same as in Planck Collaboration XXIV 2016), and two “fitted” versions obtained by rescaling the noise per patch by a constant in order to reproduce the completeness derived from the simulation clusters and the PACT (Pointecouteau et al. 2021) profile clusters. The filled contours represent 68% and 95% confidence regions.</figcaption></figure>

<table><caption>**Table 2.** Constraints on cosmological parameters obtained from the analysis of Planck PSZ2 cluster number counts and BAO, using three different completeness functions, as described in the text.</caption><thead><tr><th></th><th>Simulation<br>completeness</th><th>ERF<br>completeness</th><th>PACT<br>completeness</th></tr></thead><tbody><tr><td>$H_0$</td><td>68.4<sup>+2.0</sup><sub>-2.4</sub></td><td>69.1<sup>+2.1</sup><sub>-2.4</sub></td><td>70.3<sup>+2.2</sup><sub>-2.5</sub></td></tr><tr><td>$\Omega_m$</td><td>0.301<sup>+0.019</sup><sub>-0.023</sub></td><td>0.314<sup>+0.020</sup><sub>-0.024</sub></td><td>0.337 ± 0.024</td></tr><tr><td>$\sigma_8$</td><td>0.742<sup>+0.025</sup><sub>-0.033</sub></td><td>0.770<sup>+0.027</sup><sub>-0.036</sub></td><td>0.823<sup>+0.030</sup><sub>-0.039</sub></td></tr></tbody></table>

completeness with respect to the ERF (in this case, the Simulation completeness) favours lower values of the two cosmological parameters, while a lower completeness (that of PACT) favours higher values. For $\Omega_m$, the shift amounts to about $0.6\sigma$ for the Simulation completeness, and about $1.1\sigma$ for that of PACT. For $\sigma_8$, the two shifts are $0.9\sigma$ and $1.8\sigma$, respectively. The relation between the completeness function and these two cosmological parameters can be understood in the following terms: as the completeness function informs us about the fraction of detected clusters over the total, keeping the number of detected objects fixed with a higher (respectively, lower) completeness means that the total number of “true” clusters gets lower (higher). This, in turn, favours a lower (higher) value for both $\Omega_m$ and $\sigma_8$. It is interesting to note how the direction of the shift in the ($\Omega_m$, $\sigma_8$) plane is the same as moving the prior on the mass bias, as shown in Fig. 7 of Planck Collaboration XXIV (2016).
**5. Discussion**

A clear limitation to our set of simulation images is the restricted range of scales we are able to analyse, approximately [3–14] arcmin, compared to the extent of the *Planck* catalogue (~[1–30] arcmin). This is due to the limits of the simulation: first, the limited volume of the simulation does not allow many high-mass clusters to form. This becomes a limiting factor at high redshift, because there are not enough high-mass clusters to sample the high-$Y_{5R500}$ tail of the completeness. This translates to a lower limit in $\theta_{500}$. The second limit is related to the spacing of the simulation snapshots in redshift, which is about 0.01 at low z. When translating the cluster distribution from the mass–redshift plane to the ($\theta_{500}$, $Y_{5R500}$) plane, we see large empty regions appear diagonally at high $\theta_{500}$, which prevent a good coverage of the completeness in that region, and effectively impose a higher limit on the cluster scale.

This limitation prevents us from probing the completeness in the extremes of the *Planck* cluster sample with the simulation set. This is particularly true regarding the effect of the departure from spherical symmetry of the cluster shapes, which in our sample is seen having an increased impact for bigger clusters. However, we believe that this limit has little impact on our main result, that is, the impact of the mean profile of the clusters on the completeness. Indeed, we performed a test using spherical clusters with a broad range in $\theta_{500}$ (reproducing the extent of the *Planck* catalogue) with different pressure profiles and observed a consistent shift in the completeness at both larger and smaller scales. Given the result of this test, we decided to use a reduced $\theta_{500}$ range for our sets of spherical images to ensure a fair comparison with the Simulation set.

Another aspect that might be deemed too simplistic is the use of the same pressure profile in all the injected spherical images for each of the sets we use. In this way, the cluster signals we inject are all just scaled versions of each other. This is clearly unrealistic, and neglects the scatter around the mean profile present in the data. The main reason is that, given that we already incorporate a realistic treatment with the Simulation set, we treat the sets of spherically symmetric images as idealised test cases, and use them to understand the response of the detection algorithm to different profile shapes, without purporting to carry out a realistic test. Moreover, the impact of profile scatter between spherical images has already been detailed in Planck Collaboration XXIX (2014), Planck Collaboration XXVII (2016), who found a widening effect in the completeness drop-off region. This effect is qualitatively different from the shift we observe when using different mean profiles. We therefore believe that including a scatter in the profile would not significantly impact our conclusions, and for this reason we decided not to include it.

A third possible criticism to the work presented so far is the use of just one simulation. Indeed, the specific shape and values of the completeness function we find for the simulation set are certainly simulation-dependent to an extent. First, the chosen simulation is realistic enough in terms of pressure profile of the cluster gas, as the mean profile of the simulation clusters is not far from that observed in *Planck* and ACT by Tramonte et al. (2023). Also, we believe that the main result of this work is the dependence of the completeness on additional cluster characteristics to those of the ERF estimate (i.e. cluster size and total integrated signal), notably the steepness of the Compton-$y$ profile, which, as we demonstrate, has an important impact on the completeness.

**6. Conclusion**

In this work, we analysed the completeness function for galaxy clusters detected in SZ by *Planck* with the MMF algorithm, focusing on the case where the cluster model assumed in the detection method differs from the “true” one, both in terms of shape and pressure profile. To this end, we performed a Monte Carlo injection of SZ signal from simulated clusters into cleaned *Planck* sky maps. Comparing the injected sources with those detected by the MMF algorithm, we computed the completeness in bins of injected ($Y_{5R500}$, $\theta_{500}$) (i.e. total SZ signal and radius of clusters). We used a set of realistic cluster images from the IllustrisTNG simulation, as well as images of spherically symmetric clusters generated from different pressure profiles.

The main result of our analysis is an apparent relation between the completeness function and the shape of the pressure profile of the injected clusters with respect to the assumed cluster model. This effect is present when using both simulation images and spherical images. More specifically, we observe that cluster sets with steeper profiles (e.g. the average simulation one, or that from Tramonte et al. 2023) produce completeness functions that are higher than what we obtain using clusters with the same profile as in the detection template (our Standard gNFW set). On the other hand, cluster sets with flatter profiles (e.g. Pointecouteau et al. 2021) lead to lower completeness than the Standard gNFW set. The completeness of the Standard gNFW set is consistent with the theoretical ERF estimate (based on the assumption of Gaussian noise in the SZ signal), as expected. We explain this relation by observing a bias in the estimated properties of detected clusters, in particular the cluster radius $\theta_{500}$: when the cluster profile is steeper (respectively, flatter) than the template, the MMF algorithm tends to underestimate (overestimate) its $\theta_{500}$, which in turn leads to an increase (decrease) in its S/N, given the noise dependence on the scale.

We also tested the impact of the asymmetric cluster morphology of our simulated set on the completeness. Comparing the most and least spherical subsets, we found that the effect of the asphericity is generally small, especially below the *Planck* beam size, but there is a tendency of the “more elliptical” subset to be systematically less complete than the “more spherical” one. This effect tends to become more important as $\theta_{500}$ increases, in particular above the size of the *Planck* beam, at least in the range covered by our cluster set.

Finally, we studied the effect of modifying the completeness function in the cosmological analysis with galaxy cluster number counts. Taking two completeness examples that deviate substantially from the ERF estimate, we find that the impact on the cosmological parameters ($\Omega_m$, $\sigma_8$) is a shift in the constraints of roughly $1\sigma$, with higher completeness favouring lower values of the cosmological parameters, and vice versa. We recall that what we tested here is an extreme scenario in which the true clusters differ from those used to constrain the detection template, scaling relations, mass bias, and so on. Interestingly, we note that the shift due to changes in the completeness goes in the same direction as the shift due to the mass bias, and therefore the two effects should probably be taken into account together. This work highlights the fact that any bias and uncertainty on the completeness should be propagated into the galaxy cluster number count likelihood for robust and accurate analyses.

*Acknowledgements.* The authors acknowledge the fruitful discussions and comments from Nabila Aghanim and the other members of the Cosmology team at IAS. S.G. acknowledges financial support from the Ecole Doctorale d'Astronomie et d'Astrophysique d'Ile-de-France (ED AAIF). This
**5. Discussion**

A clear limitation to our set of simulation images is the restricted range of scales we are able to analyse, approximately [3–14] arcmin, compared to the extent of the *Planck* catalogue (~[1–30] arcmin). This is due to the limits of the simulation: first, the limited volume of the simulation does not allow many high-mass clusters to form. This becomes a limiting factor at high redshift, because there are not enough high-mass clusters to sample the high-$Y_{5R500}$ tail of the completeness. This translates to a lower limit in $\theta_{500}$. The second limit is related to the spacing of the simulation snapshots in redshift, which is about 0.01 at low z. When translating the cluster distribution from the mass–redshift plane to the ($\theta_{500}$, $Y_{5R500}$) plane, we see large empty regions appear diagonally at high $\theta_{500}$, which prevent a good coverage of the completeness in that region, and effectively impose a higher limit on the cluster scale.

This limitation prevents us from probing the completeness in the extremes of the *Planck* cluster sample with the simulation set. This is particularly true regarding the effect of the departure from spherical symmetry of the cluster shapes, which in our sample is seen having an increased impact for bigger clusters. However, we believe that this limit has little impact on our main result, that is, the impact of the mean profile of the clusters on the completeness. Indeed, we performed a test using spherical clusters with a broad range in $\theta_{500}$ (reproducing the extent of the *Planck* catalogue) with different pressure profiles and observed a consistent shift in the completeness at both larger and smaller scales. Given the result of this test, we decided to use a reduced $\theta_{500}$ range for our sets of spherical images to ensure a fair comparison with the Simulation set.

Another aspect that might be deemed too simplistic is the use of the same pressure profile in all the injected spherical images for each of the sets we use. In this way, the cluster signals we inject are all just scaled versions of each other. This is clearly unrealistic, and neglects the scatter around the mean profile present in the data. The main reason is that, given that we already incorporate a realistic treatment with the Simulation set, we treat the sets of spherically symmetric images as idealised test cases, and use them to understand the response of the detection algorithm to different profile shapes, without purporting to carry out a realistic test. Moreover, the impact of profile scatter between spherical images has already been detailed in Planck Collaboration XXIX (2014), Planck Collaboration XXVII (2016), who found a widening effect in the completeness drop-off region. This effect is qualitatively different from the shift we observe when using different mean profiles. We therefore believe that including a scatter in the profile would not significantly impact our conclusions, and for this reason we decided not to include it.

A third possible criticism to the work presented so far is the use of just one simulation. Indeed, the specific shape and values of the completeness function we find for the simulation set are certainly simulation-dependent to an extent. First, the chosen simulation is realistic enough in terms of pressure profile of the cluster gas, as the mean profile of the simulation clusters is not far from that observed in *Planck* and ACT by Tramonte et al. (2023). Also, we believe that the main result of this work is the dependence of the completeness on additional cluster characteristics to those of the ERF estimate (i.e. cluster size and total integrated signal), notably the steepness of the Compton-$y$ profile, which, as we demonstrate, has an important impact on the completeness.

**6. Conclusion**

In this work, we analysed the completeness function for galaxy clusters detected in SZ by *Planck* with the MMF algorithm, focusing on the case where the cluster model assumed in the detection method differs from the “true” one, both in terms of shape and pressure profile. To this end, we performed a Monte Carlo injection of SZ signal from simulated clusters into cleaned *Planck* sky maps. Comparing the injected sources with those detected by the MMF algorithm, we computed the completeness in bins of injected ($Y_{5R500}$, $\theta_{500}$) (i.e. total SZ signal and radius of clusters). We used a set of realistic cluster images from the IllustrisTNG simulation, as well as images of spherically symmetric clusters generated from different pressure profiles.

The main result of our analysis is an apparent relation between the completeness function and the shape of the pressure profile of the injected clusters with respect to the assumed cluster model. This effect is present when using both simulation images and spherical images. More specifically, we observe that cluster sets with steeper profiles (e.g. the average simulation one, or that from Tramonte et al. 2023) produce completeness functions that are higher than what we obtain using clusters with the same profile as in the detection template (our Standard gNFW set). On the other hand, cluster sets with flatter profiles (e.g. Pointecouteau et al. 2021) lead to lower completeness than the Standard gNFW set. The completeness of the Standard gNFW set is consistent with the theoretical ERF estimate (based on the assumption of Gaussian noise in the SZ signal), as expected. We explain this relation by observing a bias in the estimated properties of detected clusters, in particular the cluster radius $\theta_{500}$: when the cluster profile is steeper (respectively, flatter) than the template, the MMF algorithm tends to underestimate (overestimate) its $\theta_{500}$, which in turn leads to an increase (decrease) in its S/N, given the noise dependence on the scale.

We also tested the impact of the asymmetric cluster morphology of our simulated set on the completeness. Comparing the most and least spherical subsets, we found that the effect of the asphericity is generally small, especially below the *Planck* beam size, but there is a tendency of the “more elliptical” subset to be systematically less complete than the “more spherical” one. This effect tends to become more important as $\theta_{500}$ increases, in particular above the size of the *Planck* beam, at least in the range covered by our cluster set.

Finally, we studied the effect of modifying the completeness function in the cosmological analysis with galaxy cluster number counts. Taking two completeness examples that deviate substantially from the ERF estimate, we find that the impact on the cosmological parameters ($\Omega_m$, $\sigma_8$) is a shift in the constraints of roughly $1\sigma$, with higher completeness favouring lower values of the cosmological parameters, and vice versa. We recall that what we tested here is an extreme scenario in which the true clusters differ from those used to constrain the detection template, scaling relations, mass bias, and so on. Interestingly, we note that the shift due to changes in the completeness goes in the same direction as the shift due to the mass bias, and therefore the two effects should probably be taken into account together. This work highlights the fact that any bias and uncertainty on the completeness should be propagated into the galaxy cluster number count likelihood for robust and accurate analyses.

**Acknowledgements.** The authors acknowledge the fruitful discussions and comments from Nabila Aghanim and the other members of the Cosmology team at IAS. S.G. acknowledges financial support from the Ecole Doctorale d'Astronomie et d'Astrophysique d'Ile-de-France (ED AAIF). This
research has made use of the computation facility of the Integrated Data
and Operation Center (IDOC, https://idoc.ias.u-psud.fr) at the Institut
d'Astrophysique Spatiale (IAS), as well as the SZ-Cluster Database (https:
//szdb.osups.universite-paris-saclay.fr). This project was carried
out using the Python libraries matplotlib (Hunter 2007), numpy (Harris et al.
2020), Astropy (Astropy Collaboration 2013, 2018, 2022).

References

Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2020, Phys. Rev. D, 102, 023509

Aghanim, N., Douspis, M., Hurier, G., et al. 2019, A&A, 632, A47

Alam, S., Ata, M., Bailey, S., et al. 2017, MNRAS, 470, 2617

Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, ARA&A, 49, 409

Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92

Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167

Bleem, L. E., Stalder, B., de Haan, T., et al. 2015, ApJS, 216, 27

Bocquet, S., Dietrich, J. P., Schrabback, T., et al. 2019, ApJ, 878, 55

Costanzi, M., Saro, A., Bocquet, S., et al. 2021, Phys. Rev. D, 103, 043522

Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Gouin, C., Aghanim, N., Bonjean, V., & Douspis, M. 2020, A&A, 635, A195

Gouin, C., Gallo, S., & Aghanim, N. 2022, A&A, 664, A198

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Henry, J. P. 1997, ApJ, 489, L1

Herranz, D., Sanz, J. L., Barreiro, R. B., & Martínez-González, E. 2002, ApJ, 580, 610

Hilton, M., Sifón, C., Naess, S., et al. 2021, ApJS, 253, 3

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Kravtsov, A. V., & Borgani, S. 2012, ARA&A, 50, 353

Lewis, A., & Bridle, S. 2002, Phys. Rev. D, 66, 103511

Limousin, M., Morandi, A., Sereno, M., et al. 2013, Space Sci. Rev., 177, 155

Liu, A., Bulbul, E., Ghirardini, V., et al. 2022, A&A, 661, A2

Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113

Melin, J. B., Bartlett, J. G., & Delabrouille, J. 2005, A&A, 429, 417

Melin, J. B., Bartlett, J. G., & Delabrouille, J. 2006, A&A, 459, 341

Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1

Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624

Nelson, D., Springel, V., Pillepich, A., et al. 2019, Comput. Astrophys. Cosmol., 6, 2

Oukbir, J., & Blanchard, A. 1992, A&A, 262, L21

Pacaud, F., Pierre, M., Melin, J. B., et al. 2018, A&A, 620, A10

Pierre, M., Pacaud, F., Adami, C., et al. 2016, A&A, 592, A1

Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648

Planck Collaboration VIII. 2011, A&A, 536, A8

Planck Collaboration VII. 2014, A&A, 571, A7

Planck Collaboration XIII. 2014, A&A, 571, A13

Planck Collaboration XVI. 2014, A&A, 571, A16

Planck Collaboration XX. 2014, A&A, 571, A20

Planck Collaboration XXIX. 2014, A&A, 571, A29

Planck Collaboration VIII. 2016, A&A, 594, A8

Planck Collaboration XIII. 2016, A&A, 594, A13

Planck Collaboration XXII. 2016, A&A, 594, A22

Planck Collaboration XXIV. 2016, A&A, 594, A24

Planck Collaboration XXVI. 2016, A&A, 594, A26

Planck Collaboration XXVII. 2016, A&A, 594, A27

Planck Collaboration Int. V. 2013, A&A, 550, A131

Pointecouteau, E., Santiago-Bautista, I., Douspis, M., et al. 2021, A&A, 651, A73

Rozo, E., Wechsler, R. H., Rykoff, E. S., et al. 2010, ApJ, 708, 645

Schneider, P., & Bartelmann, M. 1997, MNRAS, 286, 696

Springel, V. 2010, MNRAS, 401, 791

Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676

Steigman, G. 2008, arXiv e-prints [arXiv:0807.3004]

Sunyaev, R. A., & Zeldovich, Y. B. 1970, Ap&SS, 7, 3

Sunyaev, R. A., & Zeldovich, Y. B. 1972, Comments Astrophys. Space Phys., 4, 173

Sunyaev, R. A., & Zeldovich, I. B. 1980, ARA&A, 18, 537

Tramonte, D., Ma, Y.-Z., Yan, Z., et al. 2023, ApJS, 265, 55

Vallés-Pérez, D., Planelles, S., & Quilis, V. 2020, MNRAS, 499, 2303

Zubeldia, Í., Rotti, A., Chluba, J., & Battye, R. 2023, MNRAS, 522, 4766
research has made use of the computation facility of the Integrated Data
and Operation Center (IDOC, https://idoc.ias.u-psud.fr) at the Institut
d'Astrophysique Spatiale (IAS), as well as the SZ-Cluster Database (https:
//szdb.osups.universite-paris-saclay.fr). This project was carried
out using the Python libraries matplotlib (Hunter 2007), numpy (Harris et al.
2020), Astropy (Astropy Collaboration 2013, 2018, 2022).

References

Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2020, Phys. Rev. D, 102, 023509

Aghanim, N., Douspis, M., Hurier, G., et al. 2019, A&A, 632, A47

Alam, S., Ata, M., Bailey, S., et al. 2017, MNRAS, 470, 2617

Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, ARA&A, 49, 409

Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92

Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167

Bleem, L. E., Stalder, B., de Haan, T., et al. 2015, ApJS, 216, 27

Bocquet, S., Dietrich, J. P., Schrabback, T., et al. 2019, ApJ, 878, 55

Costanzi, M., Saro, A., Bocquet, S., et al. 2021, Phys. Rev. D, 103, 043522

Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Gouin, C., Aghanim, N., Bonjean, V., & Douspis, M. 2020, A&A, 635, A195

Gouin, C., Gallo, S., & Aghanim, N. 2022, A&A, 664, A198

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Henry, J. P. 1997, ApJ, 489, L1

Herranz, D., Sanz, J. L., Barreiro, R. B., & Martínez-González, E. 2002, ApJ, 580, 610

Hilton, M., Sifón, C., Naess, S., et al. 2021, ApJS, 253, 3

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Kravtsov, A. V., & Borgani, S. 2012, ARA&A, 50, 353

Lewis, A., & Bridle, S. 2002, Phys. Rev. D, 66, 103511

Limousin, M., Morandi, A., Sereno, M., et al. 2013, Space Sci. Rev., 177, 155

Liu, A., Bulbul, E., Ghirardini, V., et al. 2022, A&A, 661, A2

Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113

Melin, J. B., Bartlett, J. G., & Delabrouille, J. 2005, A&A, 429, 417

Melin, J. B., Bartlett, J. G., & Delabrouille, J. 2006, A&A, 459, 341

Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1

Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624

Nelson, D., Springel, V., Pillepich, A., et al. 2019, Comput. Astrophys. Cosmol., 6, 2

Oukbir, J., & Blanchard, A. 1992, A&A, 262, L21

Pacaud, F., Pierre, M., Melin, J. B., et al. 2018, A&A, 620, A10

Pierre, M., Pacaud, F., Adami, C., et al. 2016, A&A, 592, A1

Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648

Planck Collaboration VIII. 2011, A&A, 536, A8

Planck Collaboration VII. 2014, A&A, 571, A7

Planck Collaboration XIII. 2014, A&A, 571, A13

Planck Collaboration XVI. 2014, A&A, 571, A16

Planck Collaboration XX. 2014, A&A, 571, A20

Planck Collaboration XXIX. 2014, A&A, 571, A29

Planck Collaboration VIII. 2016, A&A, 594, A8

Planck Collaboration XIII. 2016, A&A, 594, A13

Planck Collaboration XXII. 2016, A&A, 594, A22

Planck Collaboration XXIV. 2016, A&A, 594, A24

Planck Collaboration XXVI. 2016, A&A, 594, A26

Planck Collaboration XXVII. 2016, A&A, 594, A27

Planck Collaboration Int. V. 2013, A&A, 550, A131

Pointecouteau, E., Santiago-Bautista, I., Douspis, M., et al. 2021, A&A, 651, A73

Rozo, E., Wechsler, R. H., Rykoff, E. S., et al. 2010, ApJ, 708, 645

Schneider, P., & Bartelmann, M. 1997, MNRAS, 286, 696

Springel, V. 2010, MNRAS, 401, 791

Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676

Steigman, G. 2008, arXiv e-prints [arXiv:0807.3004]

Sunyaev, R. A., & Zeldovich, Y. B. 1970, Ap&SS, 7, 3

Sunyaev, R. A., & Zeldovich, Y. B. 1972, Comments Astrophys. Space Phys., 4, 173

Sunyaev, R. A., & Zeldovich, I. B. 1980, ARA&A, 18, 537

Tramonte, D., Ma, Y.-Z., Yan, Z., et al. 2023, ApJS, 265, 55

Vallés-Pérez, D., Planelles, S., & Quilis, V. 2020, MNRAS, 499, 2303

Zubeldia, Í., Rotti, A., Chluba, J., & Battye, R. 2023, MNRAS, 522, 4766