© The Authors 2025

1. Introduction

Galaxy clusters are excellent probes for both cosmological (Mantz et al. 2015; Planck Collaboration XXIV 2016; Costanzi et al. 2019; Marulli et al. 2021; Lesci et al. 2022a; Romanello et al. 2024; Seppi et al. 2024; Ghirardini et al. 2024) and astrophysical (Vazza et al. 2017; CHEX-MATE Collaboration 2021 Zhu et al. 2021; Sereno et al. 2021) studies. Counts and clustering represent the most powerful cosmological probes based on galaxy clusters (Sartoris et al. 2016; Fumagalli et al. 2024), followed by for example the two-halo term of cluster profiles (Giocoli et al. 2021; Ingoglia et al. 2022) and sparsity (Corasaniti et al. 2021). Accurate models describing such summary statistics can be attained through N-body dark matter simulations (Borgani & Kravtsov 2011; Angulo et al. 2012; Giocoli et al. 2012), that allow for the derivation of halo mass function (Sheth & Tormen 1999; Tinker et al. 2008; Despali et al. 2016; Euclid Collaboration: Castro et al. 2023) and halo bias (Sheth et al. 2001; Tinker et al. 2010; Euclid Collaboration: Castro et al. 2024a) models. As baryonic physics is expected to become one of the leading sources of systematic uncertainties in forthcoming cosmological analyses, several theoretical prescriptions describing its impact on cluster statistics have recently been proposed (Cui et al. 2012; Velliscig et al. 2014; Bocquet et al. 2016; Castro et al. 2021; Euclid Collaboration: Castro et al. 2024b).

Cluster mass measurements are the foundation of cosmological investigations based on such tracers, since mass is a fundamental cosmological quantity and several different cluster observational properties are good mass proxies (Teyssier et al. 2011; Pratt et al. 2019). In optical and near-infrared surveys, such mass proxies can be either the richness, which is the number of cluster member galaxies (Rykoff et al. 2014; Bellagamba et al. 2018; Maturi et al. 2019; Wen & Han 2024), the signal amplitude, in the case of matched-filter cluster detection algorithms (Bellagamba et al. 2018; Maturi et al. 2019), or the stellar mass of member galaxies (Palmese et al. 2020; Pereira et al. 2020; Doubrawa et al. 2023). Among the most reliable probes for measuring cluster masses we have weak gravitational lensing, which consists of the deflection of light rays from background sources due to the intervening cluster gravitational potential (see e.g. Bardeau et al. 2007; Okabe et al. 2010; Hoekstra et al. 2012; Melchior et al. 2015; Sereno et al. 2017; Schrabback et al. 2018; Stern et al. 2019; Giocoli et al. 2021; Ingoglia et al. 2022; Euclid Collaboration: Sereno et al. 2024; Grandis et al. 2024; Kleinebreil et al. 2025; Payerne et al. 2025). In fact, weak-lensing analyses do not rely on any assumptions about the dynamical state of the cluster, which is different from other methods used to estimate cluster masses based on the properties of the gas and member galaxies, such as cluster X-ray emission (Ettori et al. 2013; Eckert et al. 2020; Scheck et al. 2023), galaxy velocity dispersion (Kodi Ramanah et al. 2020; Ho et al. 2022; Biviano & Mamon 2023; Sereno et al. 2025), or the Sunyaev-Zeldovich (SZ) effect on the cosmic microwave background (CMB, Arnaud et al. 2010; Planck Collaboration XX 2014; Hilton et al. 2021). Given cluster weak-lensing measurements, mass estimates can be derived using robust models that describe the density profiles of the dark matter haloes that host galaxy clusters (Navarro et al. 1997; Baltz et al. 2009; Diemer & Kravtsov 2014), calibrated on N-body simulations.

In this work, we present the cosmological analysis based on the galaxy cluster catalogue by Maturi et al. (2025, referred to as M25 hereafter), which was built up through the use of the Adaptive Matched Identifier of Clustered Objects (AMICO, Bellagamba et al. 2018; Maturi et al. 2019) in the fourth data release of the Kilo-Degree Survey (referred to as KiDS-1000, Kuijken et al. 2019). To this end, we measured the stacked weak-lensing signal of 8730 clusters in bins of redshift and mass proxy, up to redshift z = 0.8 and based on the gold sample of galaxy shape measurements by Giblin et al. (2021). Then, we jointly modelled these measurements and cluster counts to derive constraints on fundamental cosmological parameters, such as the matter density parameter at z = 0, Ω_m, and the square root of the mass variance computed on a scale of 8 h⁻¹ Mpc at z = 0, σ₈. We note that this work presents an update of the mass calibration performed by Bellagamba et al. (2019, referred to as B19 hereafter), which was based on the AMICO cluster catalogue derived by Maturi et al. (2019) from the third KiDS data release (KiDS-DR3, de Jong et al. 2017). The need for a new mass calibration stems from the enhanced survey coverage and improved galaxy photometric redshift (photo-z) estimates. The latter, in fact, yield significant variations in the cluster mass proxy estimates between KiDS-DR3 and KiDS-1000 (M25).

The present work is part of a series of investigations that exploit the AMICO galaxy clusters in KiDS for both cosmological (Giocoli et al. 2021; Ingoglia et al. 2022; Lesci et al. 2022a,b; Busillo et al. 2023; Romanello et al. 2024) and astrophysical (Bellagamba et al. 2019; Sereno et al. 2020; Radovich et al. 2020; Puddu et al. 2021; Castignani et al. 2022, 2023; Tramonte et al. 2023; Giocoli et al. 2024) studies. The paper is organised as follows. In Sect. 2 we present the galaxy cluster and shear samples, while in Sect. 3 the selection function of the cluster sample is discussed. The measurement and stacking of the cluster weak-lensing profiles, along with the calibration of the background redshift distributions, are detailed in Sect. 4. Section 5 introduces the theoretical models adopted in the analysis, along with the likelihood function and parameter priors. In Sect. 6 we discuss our results and in Sect. 7 we draw our conclusions.

Throughout this paper, we assume a concordance flat Λ cold dark matter (ΛCDM) cosmological model. We refer to M₂₀₀ as the mass enclosed within the critical radius R₂₀₀, which is the distance from the cluster centre within which the mean density is 200 times the critical density of the Universe at the cluster redshift. The AB magnitude system is employed throughout this paper. The base 10 logarithm is referred to as log, while ln represents the natural logarithm. The statistical analyses presented in this paper are performed using the CosmoBolognaLib¹ (CBL, Marulli et al. 2016), a set of free software C++/Python numerical libraries for cosmological calculations. The linear matter power spectrum is computed with CAMB² (Lewis & Challinor 2011).

2. Dataset

This work is based on KiDS-1000 (Kuijken et al. 2019), which relies on observations carried out with the OmegaCAM wide-field imager (Kuijken 2011), installed on the European Southern Observatory (ESO) VLT Survey Telescope (VST, Capaccioli & Schipani 2011). The KiDS-1000 release covers 1006 tiles of about 1 deg² each. The 2 arcsec aperture photometry in the u, g, r, i bands is provided, with 5σ limiting magnitudes of 24.23, 25.12, 25.02 and 23.68 for the four aforementioned bands, respectively. In addition, aperture-matched Z, Y, J, H, K_s near-infrared photometry from the fully overlapping VISTA Kilo degree INfrared Galaxy survey (VIKING, Edge et al. 2013; Sutherland et al. 2015) is included. Based on these photometric bands, the KiDS-1000 galaxy sample includes photo-z estimates, derived with the Bayesian Photometric Redshift (BPZ, Benítez 2000) code, for more than 100 million objects, with a typical uncertainty of σ_z^gal/(1 + z) = 0.072 at full depth (Kuijken et al. 2019). As discussed in Sect. 2.2, the galaxy photo-zs adopted in this work differ from those delivered by Kuijken et al. (2019), though they maintain the same precision.

2.1. Galaxy cluster sample and mass proxy

The cluster catalogue on which this work is based, named AMICO KiDS-1000, was built applying the AMICO algorithm (Bellagamba et al. 2018; Maturi et al. 2019) to the KiDS-1000 photometric data (for a more detailed description, see M25). AMICO is based on an optimal matched filter and it is one of the two algorithms for cluster identification officially adopted by the Euclid mission (Euclid Collaboration: Adam et al. 2019), providing highly pure and complete samples of galaxy clusters based on photometric observations. To construct the AMICO KiDS-1000 cluster catalogue, all galaxies located in the regions affected by image artefacts, such as spikes and haloes around bright stars, have been removed. In addition, in this work we exclude the clusters that fall in the regions affected by photometric issues by imposing ARTIFACTS_FLAG < 2, and we neglect the tile with central coordinates RA = 196 and Dec = 1.5 due to its low data quality (see M25). This yields a final effective area of 839 deg², containing cluster detections up to z = 0.9. Furthermore, we consider clusters having S/N > 3.5 and in the redshift range z ∈ [0.1, 0.8], resulting in a sample of 22 396 objects. We conservatively exclude cluster detections at z > 0.8 due to low sample purity and completeness, which also results in poor matching with clusters in existing literature catalogues. Throughout this paper, we employ additional cuts based on the cluster mass proxy that yield 8730 and 7789 clusters in the weak-lensing and count analyses, respectively. The redshift selection is based on the cluster redshifts corrected for the bias derived by M25, relying on a comparison with a collection of spectroscopic redshifts available within KiDS (referred to as KiDZ, van den Busch et al. 2020, 2022), mainly based on the following surveys: Galaxy and Mass Assembly survey data release 4 (GAMA-DR4, Driver et al. 2011, 2022; Liske et al. 2015), Sloan Digital Sky Survey (SDSS, Alam et al. 2015), and 2-degree Field Lensing Survey (2dFLens, Blake et al. 2016). The origin of this bias lies in the systematic uncertainties associated with galaxy photo-zs. Throughout this paper, we assign these bias-corrected redshifts to each cluster. The inclusion of the statistical uncertainty on this correction leads to a final cluster redshift precision of σ_z/(1 + z) = 0.014. In the top panel of Fig. 1 we show the distribution of the corrected cluster redshifts, along with the redshift distribution of the AMICO KiDS-DR3 cluster sample adopted by B19 for weak-lensing mass calibration. In the latter sample, covering an effective area of 377 deg², detections with S/N > 3.5 and z ∈ [0.1, 0.6] were included, for a total of 6961 clusters. Thus, the KiDS-1000 sample covers a larger redshift range and contains more cluster detections thanks to the larger survey area. Furthermore, thanks to the improved photometry and the inclusion of VIKING bands, the significant incompleteness at z ∼ 0.35 reported by Maturi et al. (2019), caused by the 4000 Å break between g and r bands, is absent in the current sample. Indeed, the addition of such near-infrared bands improves the accuracy of galaxy photo-zs, reducing their dependency on the 4000 Å break. The cosmological sample used for the AMICO KiDS-DR3 cluster count analysis by Lesci et al. (2022a) contained fewer objects than the one used by B19, namely 3652. Thus, the number of clusters contributing to the counts in this work is more than twice that of KiDS-DR3.

Fig. 1. Redshift (top panel) and intrinsic richness (bottom panel) distributions of the AMICO galaxy clusters in KiDS-1000 (blue histograms), with S/N > 3.5 and z ∈ [0.1, 0.8], and in KiDS-DR3 (red hatched histograms), following the S/N and redshift selections employed in B19. In both panels, no λ* selections have been applied.

AMICO provides the estimate of two mass proxies for each cluster detection, namely the signal amplitude, A, and the intrinsic richness, λ^* (Bellagamba et al. 2018). The amplitude is obtained by convolving the observed galaxy 3D distribution and the cluster model adopted for the AMICO filter, for which the convolution of a Schechter luminosity function (Schechter 1976) and a Navarro et al. (1997, NFW) profile is assumed. Following Hennig et al. (2017), a concentration parameter of c = 3.59 is adopted in the NFW model (for more details, see M25). The contribution by the Brightest Cluster Galaxy (BCG) is not included in the AMICO filter, in order to prevent biases due to photo-z uncertainties. Indeed, the BCG signal often dominates over other cluster galaxies. The intrinsic richness is defined as follows:

$$\begin{array}{c}\hfill {\lambda }_j^{\ast }={\displaystyle \sum _i}{P}_i(j)\text{with}\{\begin{array}{c}{m}_i<m^{\ast }(z_j)+1.5\hfill \\ R_i(j)<R_{\max }(z_j)\hfill \end{array},\end{array}$$(1)

where P_i(j) is the probability assigned by AMICO to the ith galaxy of being a member of a given detection j, while z_j is the redshift of the jth detected cluster, m_i is the magnitude of the ith galaxy, R_i corresponds to the distance of the ith galaxy from the centre of the cluster, R_max(z_j) is the radius enclosing a mass of M₂₀₀ = 10¹⁴ h⁻¹ M_⊙, which is the typical mass of a galaxy cluster, and m^* corresponds to the r-band luminosity at the knee of the Schechter function used in the cluster model assumed for the AMICO filter. Thus, λ^* represents the sum of the galaxy membership probabilities, constrained by the conditions given in Eq. (1). As shown by Bellagamba et al. (2018), Maturi et al. (2023), Toni et al. (2024, 2025), the sum of the membership probabilities is an excellent estimator of the true number of member galaxies. We note that the enclosed mass M₂₀₀ is kept fixed in the AMICO detection process to reduce the noise in λ^* measurements.

In this work, we focused on the calibration of the log λ^* − log M₂₀₀ relation, using it to model cluster abundance and derive constraints on cosmological parameters. Compared to the amplitude A, λ^* is less dependent on the cluster model assumed for AMICO detection (Maturi et al. 2019). Furthermore, the weight assigned to each galaxy in Eq. (1) is by construction lesser than or equal to 1, given the definition of probability. This implies that while some galaxies with especially high weights might increase the value of A, their impact on λ^* would be milder. Nevertheless, given that A proved to be an effective proxy of the number of cluster galaxies in simulations (Bellagamba et al. 2018), we will thoroughly investigate the performance of this mass proxy in future studies based on KiDS data. In the bottom panel of Fig. 1 we show the λ^* distribution of the cluster sample. As detailed in M25, the λ^* distribution of KiDS-DR3 clusters is remarkably different from the one in KiDS-1000, due to the improvement of galaxy photo-zs in the latter survey thanks to the inclusion of VIKING photometry. Specifically, the better galaxy photo-zs led to the increase of λ^* for all KiDS-1000 galaxy clusters. We also note that, as a consequence of the S/N cut applied in this analysis, the clusters with the lowest λ^* values shown in Fig. 1 are predominantly located at low redshifts (M25). As discussed in the following, these objects are not included in the analysis due to the low sample purity at small λ^*.

2.2. Shear sample

To estimate the weak-lensing signal of the AMICO KiDS-1000 galaxy clusters, we based our analysis on the KiDS-1000 gold shear catalogue (Wright et al. 2020; Hildebrandt et al. 2021; Giblin et al. 2021). This sample comprises about 21 million galaxies, covering an effective area of 777.4 deg², with a weighted number density of $n_{\mathrm{eff}}^{\mathrm{gold}}=6.17$ arcmin⁻². It includes weak-lensing shear measurements from the deep KiDS r-band observations, as these are the ones with better seeing properties and yield the highest source density. The estimator used in this analysis is lensftit (Miller et al. 2007; Fenech Conti et al. 2017), a likelihood-based model fit method that has been used successfully in the analyses of other datasets, such as the Canada France Hawaii Telescope Lensing Survey (CFHTLenS, Miller et al. 2013) and the Red Cluster Survey (Hildebrandt et al. 2016). The photo-z calibration performed for cosmic shear studies (Asgari et al. 2021), based on deep spectroscopic reference catalogues that are weighted with the help of self-organising maps (SOMs, Kohonen 1982; Wright et al. 2020), has been performed by Hildebrandt et al. (2021) in the redshift range z ∈ [0.1, 1.2].

As in Kuijken et al. (2019), M25 used BPZ (Benítez 2000) to derive the galaxy photo-zs for cluster detection. However, the prior redshift probability distribution used in BPZ by Kuijken et al. (2019) appears to generate a remarkable redshift bias for bright, low-redshift galaxies, despite the fact that it reduces uncertainties and catastrophic failures for faint galaxies at higher redshifts. For this reason, M25 chose the same prior on redshift probabilities adopted in the KiDS-DR3 analysis (de Jong et al. 2017). This leads to a galaxy photo-z uncertainty of σ_z^gal/(1 + z) = 0.074 given the r magnitude limit r < 23, in agreement with Kuijken et al. (2019). In this analysis, we adopted the galaxy photo-zs derived in this way. As we shall discuss later in Sect. 4.4, we calibrated the galaxy redshift distribution using the SOM algorithm provided by Wright et al. (2020).

3. Cluster selection function and blinding

In this section, we delve into the description of the cluster selection function estimates and summarise the blinding of the statistical analysis presented in this paper. To derive the selection function of the galaxy cluster sample, namely its completeness and the uncertainties associated with cluster observables, along with its purity, we used the data-driven mock cluster catalogue produced by M25. This mock sample was constructed using the Selection Function extrActor (SinFoniA) code (Maturi et al. 2019, 2025). SinFoniA created mock galaxy catalogues based on the KiDS-1000 photometric data, randomly extracting field galaxies and cluster members according to their measured membership probabilities using a Monte Carlo approach. These mock catalogues maintain the autocorrelations of galaxies and clusters, the cross-correlation between clusters and galaxies, the survey footprint, and masks, along with the local variations in data properties caused by the varying depth of the survey. The mock catalogues are then processed using AMICO, allowing for a comparison between the detections and the clusters injected into the mock catalogues, ultimately deriving the final statistical properties of the sample.

We selected all AMICO detections in this mock catalogue that have a S/N > 3.5, which corresponds to the S/N threshold used to define the cluster sample described in Sect. 2.1. Then, we measured the cluster purity, 𝒫_clu, defined as the number of AMICO detections having a true cluster counterpart over the total number of AMICO detections, as a function of observed intrinisc richness and redshift, namely λ_ob^* and z_ob, respectively. Figure 2 displays the values of 𝒫_clu related to each bin considered for the cluster count and weak-lensing measurements detailed in the following. We remind that the impurities of the cluster sample yield abundance measurements that are biased high, while the stacked weak-lensing measurements are biased low. Indeed, weak-lensing profiles derived around random positions in the sky are expected to be statistically consistent with zero. We account for these effects in the models described in Sect. 5.

Fig. 2. Purity of the cluster sample as a function of λob*, in the redshift bins z ∈ [0.1, 0.3) (top panel), z ∈ [0.3, 0.45) (middle panel), and z ∈ [0.45, 0.8] (bottom panel). The pink histograms show the purity derived in the λob* bins used for cluster weak-lensing measurements, while the black hatched histograms display the purity in the bins used for cluster counts.

To estimate the uncertainties on the observed intrinisc richness, λ_ob^*, in the mock catalogue we measure the probability density function P(Δx | Δλ_tr^*, Δz_tr), where x = (λ_ob^* − λ_tr^*)/λ_tr^*. Here, λ_tr^* is the true intrinisc richness, while Δλ_tr^* and Δz_tr represent bins of true intrinsic richness and true redshift, respectively. We find that P(Δx | Δλ_tr^*, Δz_tr) can be accurately modelled as a Gaussian for any Δz_tr and Δλ_tr^*. Through a joint modelling of this distribution in different bins of z_tr and λ_tr^*, performed by means of a Bayesian analysis based on a Markov chain Monte Carlo (MCMC), assuming a Gaussian likelihood, Poisson uncertainties, and large uniform priors on the free parameters presented in the following, we find that the mean of P(Δx | Δλ_tr^*, Δz_tr) can be expressed as

$$\begin{array}{c}\hfill {\mu }_x=exp[-{\lambda }_{\mathrm{tr}}^{\ast }(A_{\mu }+B_{\mu }{z}_{\mathrm{tr}})],\end{array}$$(2)

where A_μ = 0.198 ± 0.005 and B_μ = −0.179 ± 0.007, while its standard deviation has the following form

$$\begin{array}{c}\hfill {\sigma }_x=A_{\sigma }exp(-B_{\sigma }{\lambda }_{\mathrm{tr}}^{\ast }),\end{array}$$(3)

where A_σ = 0.320 ± 0.008 and B_σ = 0.011 ± 0.001. Equations (2) and (3) are represented in Fig. 3. We notice that, differently from μ_x, σ_x depends only on λ_tr^*. Indeed, we do not find any significant trend with z_tr for this quantity. Due to masking, blending of cluster detections, and projection effects, σ_x is up to 35% larger than the Poisson relative uncertainty on λ_tr^*. In addition, we note that μ_x decreases with λ_tr^* and increases with z_tr mainly due to projection effects (see e.g. Myles et al. 2021; Cao et al. 2025). We also observe that, as optical cluster finders preferentially select haloes embedded within filaments aligned with the line of sight, the λ^* bias due to projection effects is positively correlated with a cluster shear bias, amounting to up to 25% in the two-halo regime (Sunayama et al. 2020; Wu et al. 2022; Sunayama 2023; Zhou et al. 2024). Nonetheless, as will be discussed in Sect. 6.3, this weak-lensing bias has a negligible impact on our results thanks to our choice of the cluster-centric radial range. Following an approach analogous to the one used to estimate P(x | λ_tr^*, z_tr), we find that P(z_ob | λ_tr^*, z_tr) is Gaussian and its mean and standard deviation do not significantly depend on λ_tr^*. We also remark that, as mentioned in Sect. 2.1, throughout this paper we use the cluster redshift uncertainties derived by M25 from the comparison of the AMICO sample with the KiDZ spectroscopic redshifts.

Fig. 3. Top panel: Mean of P(x | λtr*, ztr), namely μx, as a function of λtr* and in different bins of ztr (see legend). Bottom panel: Grey band represents the 68% confidence level of the σx model, while the dashed orange line shows the Poisson relative uncertainty on λtr*.

While P(x | λ_tr^*, z_tr) is a probability density that allows one to quantify how many clusters fall in a nearby λ_ob^* bin due to observational uncertainties, the amount of clusters recovered from the underlying true distribution is given by the completeness of the sample 𝒞_clu(λ_tr^*, z_tr). It is defined as the number of cluster detections with S/N > 3.5 having a true counterpart divided by the number of true clusters, as a function of λ_tr^* and z_tr. Figure 4 displays the measured completeness and its smoothed version based on Chebyshev polynomials, as discussed in M25, for some values of z_tr. Although we evaluated 𝒞_clu in seven redshift intervals ranging from 0.1 to 0.9, only four are displayed to avoid overcrowding the plot. These smoothed measurements are then interpolated as a function of z_tr and λ_tr^*, in order to include 𝒞_clu in the theoretical models described in Sect. 5. We note that 𝒞_clu is estimated up to z_tr = 0.9 to account for uncertainties in the observed redshift, z_ob, which is capped at z_ob = 0.8 (see the likelihood models in Sect. 5).

Fig. 4. Completeness of the cluster sample as a function of λtr*, for some values of ztr (see the legend). Dashed lines display the completeness measurements, while solid lines show the completeness smoothed with Chebyshev polynomials.

As detailed in M25, cluster completeness estimates remained blind until the complete definition and run of the pipeline presented in this paper. Specifically, GFL received three versions of 𝒞_clu, two of which were intentionally biased in relation to the true one by MM, and built up the pipeline using only one of them. Except for MM, none of the authors knew the true value of 𝒞_clu. More details on this blinding strategy and on its outcome are provided in Appendix B and M25.

We note that, due to the complexity of the cosmological analysis presented in this work, driven by both observational and theoretical systematic uncertainties, we do not select the cluster sample based on purity and completeness estimates. Instead, we adopt a data-driven approach, modelling the sample assuming alternative λ_ob^* and z_ob cuts to identify the combination that maximises the quality of the fit and the number of clusters included in the analysis. The impact of this sample selection is assessed in Sect. 6.

4. Stacked weak-lensing profiles measurement

In this section, we present the pipeline used to measure the weak-lensing profiles of the AMICO KiDS-1000 galaxy clusters. We detail the galaxy selections employed to define the background samples associated with the clusters in our sample. Then, we introduce a density profile estimator based on observed galaxy ellipticities and discuss the stacking of our measurements. Lastly, we assess the purity of the background galaxy samples and calibrate their redshift distributions. To perform the cluster weak-lensing measurements, we assume a flat Universe with H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_m = 0.3, where H₀ is the Hubble constant and Ω_m is the matter density parameter. We note that these cosmological assumptions do not affect the final results. Indeed, they induce geometric distortions in the weak-lensing measurements that are properly modelled in the cosmological analysis (see Sect. 5.3).

4.1. Selection of background sources

If galaxies that belong to the clusters or in the foreground are mistakenly considered as background, the measured lensing signal can be significantly diluted (see e.g. Broadhurst et al. 2005; Medezinski et al. 2007; Dietrich et al. 2019; Rau et al. 2024). To minimise such contamination, we employed stringent selections based on the galaxy photo-z probability density function, denoted as p(z_g), where z_g represents the galaxy redshift, and on galaxy colours. This approach is in line with previous works, such as Sereno et al. (2017) and B19.

Specifically, to exclude galaxies with a significant probability of being at a redshift equal to or lower than the cluster redshift, we adopted the following photo-z selection

$$\begin{array}{c}\hfill z_{\mathrm{g},\min }>z+0.05,\end{array}$$(4)

where z_g, min is the minimum of the interval containing 95% of the probability around the first mode of p(z_g), namely ${\overline{z}}_{\mathrm{g}}$, while z is the mean redshift of the cluster. In Eq. (4), the 0.05 buffer is larger than the uncertainty associated with cluster redshifts, that is σ_z/(1 + z) = 0.014. A larger buffer, accounting for the photo-z quality of the background sources, is not expected to significantly impact our analysis. Indeed, as we will discuss in Sect. 4.3, the selection in Eq. (4) plays a marginal role in defining the background galaxy samples.

In addition to Eq. (4), we considered a colour selection in order to enhance the background sample completeness, similarly to Oguri et al. (2012), Medezinski et al. (2018a), Dietrich et al. (2019), Schrabback et al. (2021). Specifically, we relied on the colour selection calibrated by Euclid Collaboration: Lesci et al. (2024, referred to as L24 hereafter), based on griz photometry. The gri photometry used to calibrate this selection, corresponding to that employed in SDSS (Aihara et al. 2011), closely matches that used in KiDS (Kuijken et al. 2019), while the z band is redder in KiDS. Any biases deriving from this difference are accounted for in the self-organising map calibration presented in Sect. 4.4. When applied to reference photometric samples, such as COSMOS (Laigle et al. 2016; Weaver et al. 2022), this selection provides a background completeness ranging from 30% to 84% and a purity larger than 97% in the lens redshift range z ∈ [0.2, 0.8]. To account for cluster redshift uncertainties, we conservatively applied the colour selection at $\stackrel{\sim }{z}=z+0.05$, where 0.05 corresponds to the buffer used also in Eq. (4). For lenses at $\stackrel{\sim }{z}<0.2$, namely below the minimum of the domain covered by the L24 selection, we considered the selection valid for $\stackrel{\sim }{z}=0.2$. We also imposed that candidate background galaxies can pass the colour selection if the following condition is satisfied:

$$\begin{array}{c}\hfill {\overline{z}}_{\mathrm{g}}>z+0.05,\end{array}$$(5)

where the 0.05 buffer is the same as in Eq. (4). The complete background selection criterion is defined as follows

$$\begin{array}{c}\hfill \text{(photo-< Emphasis Type="Italic"> z< /Emphasis>}\text{selection) or}\text{(colour selection)}.\end{array}$$(6)

To assess the performance of photo-z and colour selections, we counted the number of galaxies that meet such selection criteria. We considered the same cluster redshift bins adopted for the weak-lensing measurements presented in the following analysis, namely z ∈ [0.1, 0.3), z ∈ [0.3, 0.45), and z ∈ [0.45, 0.8]. Figure 5 shows that the colour selection provides the largest number of background sources, namely N_bkg(z_g), at low z_g in all cluster redshift bins. This is expected, as the photo-z selection in Eq. (4) excludes a greater number of galaxies close to the clusters. On average, the photo-z selection appears to be more conservative compared to the colour selection, enhancing the number of background galaxies by at most 2% compared to the case of the colour selection alone. The purity of the background samples and the calibration of their redshift distributions are discussed in Sect. 4.4.

Fig. 5. Number of galaxies that satisfy the photo-z selection (red histograms) and the colour selection (blue hatched histograms) as a function of zg, in different cluster redshift bins, namely z ∈ [0.1, 0.3) (left panel), z ∈ [0.3, 0.45) (middle panel), and z ∈ [0.45, 0.8] (right panel). The grey histograms show the galaxies selected through the total selection criterion, defined in Eq. (6).

4.2. Measurement of the cluster profiles

The tangential shear, γ₊, is linked to the excess surface density profile of the lens, ΔΣ₊, via (see e.g. Sheldon et al. 2004)

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Sigma }}_{+}(R)=\overline{\mathrm{\Sigma }}(<R)-\mathrm{\Sigma }(R)={\mathrm{\Sigma }}_{\mathrm{crit}}{\gamma }_{+}(R),\end{array}$$(7)

where Σ(R) is the mass surface density, $\overline{\mathrm{\Sigma }}(<R)$ is its mean within the radius R, and Σ_crit is the critical surface density, expressed as (Bartelmann & Schneider 2001)

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{crit}}\equiv \frac{c^2}{4\pi G}\frac{D_{\mathrm{s}}}{D_{\mathrm{l}}{D}_{\mathrm{ls}}},\end{array}$$(8)

where D_s, D_l, and D_ls are the observer-source, observer-lens, and lens-source angular diameter distances, respectively.

The weak-lensing quantity directly related to galaxy ellipticities is the reduced shear, namely g = γ/(1 − κ), where κ ≡ Σ/Σ_crit is the convergence (Schneider & Seitz 1995). The observed galaxy ellipticity can be conveniently separated into a tangential component, e₊, and a cross component, e_×, as follows (see e.g. Viola et al. 2015):

$$\begin{array}{c}\hfill \left(\begin{array}{c}{e}_{+}\\ e_{\times }\end{array}\right)=\left(\begin{array}{cc}-cos(2\phi )& -sin(2\phi )\\ sin(2\phi )& -cos(2\phi )\end{array}\right)\left(\begin{array}{c}{e}_1\\ e_2\end{array}\right),\end{array}$$(9)

where φ is the position angle of the source with respect to the lens centre, while e₁ and e₂ are the observed ellipticity components. Therefore, the reduced tangential shear profile, g₊, of a cluster labelled as k can be estimated as follows:

$$\begin{array}{c}\hfill g_{+,k}(R_j)=\left({\displaystyle \frac{{\sum }_{i\in j}{w}_i{e}_{+,i}}{{\sum }_{i\in j}{w}_i}}\right)\frac{1}{1+{\mathcal{M}}_j},\end{array}$$(10)

where j is the radial annulus index, with an associated average projected radius R_j, corresponding to the central value of the radial bin. The latter is a good approximation of the so-called effective radius (see the Appendix in Sereno et al. 2017) and we adopted it because its measurement does not require any assumption on the lens properties. In addition, w_i is the statistical weight assigned to the measure of the source ellipticity of the ith galaxy (Sheldon et al. 2004), satisfying the background selection (Eq. (6)), falling in the jth radial annulus of the kth cluster. The term ℳ_j is the average correction due to the multiplicative noise bias in the shear estimate, defined as

$$\begin{array}{c}\hfill {\mathcal{M}}_j={\displaystyle \frac{{\sum }_{i\in j}{w}_i{m}_i}{{\sum }_{i\in j}{w}_i}},\end{array}$$(11)

where m_i is the multiplicative shear bias of the ith galaxy. For KiDS-1000, Giblin et al. (2021) derived m estimates in five redshift bins, which are consistent within 1σ with zero and have standard deviation ranging from 10⁻² to 2 × 10⁻², decreasing with redshift. Since the galaxy selection and photo-z calibration considered in this work differ from those adopted by Giblin et al. (2021), we extracted for each galaxy a value of m from the uniform distribution [ − 0.05, 0.05], covering the 2σ intervals of the m distributions derived by Giblin et al. (2021). This prior leads to ℳ_j ≃ 0. Thus, in practice, ℳ_j does not depend on the radial annulus. To propagate the statistical uncertainty on m, namely σ_m, into the final results, we included σ_m into the covariance matrix (see Sect. 5.4). Indeed, σ_m is also a relative uncertainty on g₊. We conservatively assumed σ_m = 0.02, corresponding to the largest m uncertainty value derived for the tomographic bins considered by Giblin et al. (2021). In principle, a term accounting for the shear additive bias, originated by detector-level effects such as inefficiencies in the charge transfer (Fenech Conti et al. 2017), should be included in Eq. (10). Giblin et al. (2021) showed that this bias is consistent with zero and has an associated statistical uncertainty of the order of 10⁻⁴, which is negligible in our analysis because it is one order of magnitude lower than the average statistical uncertainty associated with our stacked g₊ measurements.

4.3. Stacked profiles

For most of the clusters in the sample, the weak-lensing signal is too low to precisely measure their density profiles. For this reason, we derived average g₊ estimates by stacking the signal from ensembles of clusters, selected according to their mass proxy and redshift. Specifically, the reduced shear profile in the pth bin of λ^* and the qth bin of z, namely Δλ_p^* and Δz_q, respectively, is expressed as

$$\begin{array}{c}\hfill g_{+}(R_j,\mathrm{\Delta }{\lambda }_p^{\ast },\mathrm{\Delta }{z}_q)={\displaystyle \frac{{\sum }_{k\in \mathrm{\Delta }{\lambda }_p^{\ast },\mathrm{\Delta }{z}_q}{W}_{k,j}{g}_{+,k}(R_j)}{{\sum }_{k\in \mathrm{\Delta }{\lambda }_p^{\ast },\mathrm{\Delta }{z}_q}{W}_{k,j}}},\end{array}$$(12)

where g_+, k is given by Eq. (10), k runs over all clusters falling in the bins of λ^* and z, while W_k, j is the total weight for the jth radial bin of the kth cluster, estimated as

$$\begin{array}{c}\hfill W_{k,j}={\displaystyle \sum _{i\in j}}{w}_i,\end{array}$$(13)

where i runs over the background galaxies in the jth radial bin. As an alternative to R_j in Eq. (12), namely the central value of the jth radial bin, we computed also an effective radius, $R_j^{\mathrm{eff}}$, as follows (Umetsu et al. 2014; Sereno et al. 2017):

$$\begin{array}{c}\hfill R_j^{\mathrm{eff}}={\displaystyle \frac{{\sum }_{k\in \mathrm{\Delta }{\lambda }_p^{\ast },\mathrm{\Delta }{z}_q}{W}_{k,j}{R}_{k,j}^{\mathrm{eff}}}{{\sum }_{k\in \mathrm{\Delta }{\lambda }_p^{\ast },\mathrm{\Delta }{z}_q}{W}_{k,j}}},\end{array}$$(14)

where $R_{k,j}^{\mathrm{eff}}$ is the jth effective radius of the kth cluster, which is expressed as

$$\begin{array}{c}\hfill R_{k,j}^{\mathrm{eff}}={\displaystyle \frac{{\sum }_{i\in j}{w}_i{\stackrel{\sim }{R}}_i}{{\sum }_{i\in j}{w}_i}},\end{array}$$(15)

where ${\stackrel{\sim }{R}}_i$ is the projected distance of the ith galaxy from the centre of the kth cluster. We verified that $R_j^{\mathrm{eff}}\simeq R_j$, with differences of approximately 0.1%, and thus the use of $R_j^{\mathrm{eff}}$ as a radial estimate does not have a significant impact on the final results.

We measured the stacked reduced shear in 6 logarithmically-equispaced radial bins in the cluster-centric radial range R ∈ [0.4, 3.5] h⁻¹ Mpc. By excluding the central 400 h⁻¹kpc, the contamination due to cluster members is significantly reduced (see e.g. Medezinski et al. 2018b; Bellagamba et al. 2019). The analysis of the shear signal in the immediate vicinity of the cluster centre might also be affected by inaccuracies in the weak-lensing approximation and by the influence of the BCG on the matter distribution. Furthermore, the exclusion of scales larger than 3.5 h⁻¹ Mpc reduces dependence of mass calibration results on cosmological parameters, mainly originating from matter correlated with clusters (see e.g. Oguri & Hamana 2011). In addition, the impact of possible anisotropic boosts, which affect the correlation functions on large scales due to projection effects, is mitigated (Sunayama et al. 2020; Wu et al. 2022; Sunayama 2023; Park et al. 2023; Zhou et al. 2024).

In Fig. 6 we show the stacked g₊(R) measurements in the bins of redshift and mass proxy listed in Table 1. Specifically, we considered the clusters with λ^* > 20 for z ∈ [0.1, 0.3), λ^* > 25 for z ∈ [0.3, 0.45), and λ^* > 30 for z ∈ [0.45, 0.8]. This selection yields a total sample of 9049 clusters. Nonetheless, the number of clusters contributing to the stacked weak-lensing measurements is slightly lower, amounting to 8730, due to the smaller area covered by the shear sample compared to the cluster sample. These λ^* cuts were validated to ensure adequate model fit quality while maximising the number of clusters included in the analysis. The statistical part of the covariance matrix for each stack of g₊(R) was estimated with a bootstrap procedure with replacement, by performing 10 000 resamplings of the single cluster profiles contributing to the given stack. Given the large sample sizes, we do not expect significant differences compared to jackknife estimates. We used the bootstrap method to ensure robust off-diagonal covariance estimates, particularly to account for large-scale structure (LSS) contributions. Since we bootstrap the cluster profiles, the covariance arising from source galaxies shared by multiple clusters is not accounted for. However, our method effectively captures the intrinsic scatter in the observable-mass and concentration-mass relations, as well as miscentring effects. In addition, we note that in Stage-III surveys, the contribution from sources shared by multiple clusters is small compared to those from shape noise and LSS (McClintock et al. 2019). The inverted covariance matrix is corrected following Hartlap et al. (2007). Lastly, we did not consider the covariance between radial bins in different redshift and mass proxy bins, as its impact on the final results is known to be negligible (McClintock et al. 2019; Ingoglia et al. 2022).

Fig. 6. Stacked g+(R) profiles of the AMICO KiDS-1000 galaxy clusters in bins of z (increasing from top to bottom) and λ* (increasing from left to right). The black dots show the measures, while the error bars are the sum of bootstrap errors and statistical uncertainties coming from systematic errors (see Sect. 5.4). The blue bands represent the 68% confidence levels derived from the multivariate posterior of all the free parameters considered in the joint analysis of counts and weak lensing.

In Appendix A, we show that null tests exclude the presence of additive systematic uncertainties affecting weak-lensing measurements. Furthermore, we note that for cosmological purposes, g₊ measurements as a function of the angular projected separation from cluster centres are ideal, as they do not require any assumptions on the cosmological model. Nonetheless, in the case of g₊ measurements averaged over large redshift bins, a background galaxy selection based on physical projected separations, as the one adopted here, more effectively ensures the exclusion of undesired cluster regions. This selection requires the assumption of a fiducial cosmological model, with this leading to geometric distortions in the measured g₊. We accounted for these distortions in the weak-lensing modelling outlined in Sect. 5.3.

To assess the impact of alternative background selection criteria on the stacked measurements, we computed the average weak-lensing S/N in the qth cluster redshift bin as follows (Umetsu et al. 2020; Umetsu 2020):

$$\begin{array}{c}\hfill {⟨{(\mathrm{S}/\mathrm{N})}_{\mathrm{WL}}⟩}_{\mathrm{\Delta }{\lambda }^{\ast }}(\mathrm{\Delta }{z}_q)=\frac{1}{N_p}{\displaystyle \sum _{p=1}^{N_p}}\frac{{\displaystyle \sum _{j=1}^{N_R}}{g}_{+}(R_j,\mathrm{\Delta }{\lambda }_p^{\ast },\mathrm{\Delta }{z}_q)w_{g_{+}}(R_j,\mathrm{\Delta }{\lambda }_p^{\ast },\mathrm{\Delta }{z}_q)}{\sqrt{{\displaystyle \sum _{j=1}^{N_R}}{w}_{g_{+}}(R_j,\mathrm{\Delta }{\lambda }_p^{\ast },\mathrm{\Delta }{z}_q)}}·\end{array}$$(16)

Here, w_g+ = σ_g+⁻², where σ_g+ is the uncertainty on g₊, derived through bootstrap resampling. In addition, p runs over the number of λ^* bins falling within Δz_q, namely N_p, j runs over the number of radial bins within a stack, that is N_R. Figure 7 displays how ⟨(S/N)_WL⟩_Δλ^* and the number of background galaxies, N_bkg, are affected by the background selection choices. For this test, we also included a selection based on Eq. (5) only, referred to as a photo-z peak selection. The baseline background selection, namely the combination of colour and photo-z selections (Eq. (6)), leads to N_bkg and ⟨(S/N)_WL⟩_Δλ^* values that are compatible with those derived from the colour selection alone. This agrees with what discussed in Sect. 4.1 (see also Fig. 5). Consequently, as shown in Fig. 7, the photo-z selection alone provides a ⟨(S/N)_WL⟩_Δλ^* which is remarkably lower compared to what is derived from the baseline selection. A considerable enhancement of N_bkg is achieved at high z by adopting the photo-z peak selection. Nevertheless, due to the low purity of this selection (see Sect. 4.4), the ⟨(S/N)_WL⟩_Δλ^* is not appreciably larger than that obtained with the baseline selection. Thus, ⟨(S/N)_WL⟩_Δλ^* is not a good estimate of the relative contamination between alternative background selection schemes, because contaminants lower both the signal and the noise. In addition, we remark that Eq. (16) is an approximation of the S/N, as it neglects the cross-correlation between radial bins.

Fig. 7. Top panel: ⟨(S/N)WL⟩Δλ* defined in the cluster redshift bins adopted for the stacking, namely z ∈ [0.1, 0.3), z ∈ [0.3, 0.45), and z ∈ [0.45, 0.8]. Bottom panel: Number of background sources as a function of z. In both panels, quantities obtained from the combination of colour and photo-z selections (solid black lines), colour selection only (dashed green lines), photo-z selection only (dashed blue lines), and photo-z peak selection are shown. The black and green curves are almost overlapping.

4.4. SOM-reconstructed background distributions

Foreground and cluster galaxies dilute the weak-lensing signal, as their shapes are uncorrelated with the matter distribution of galaxy clusters. This holds strictly true only in the absence of intrinsic alignments. To assess the fraction of such contaminants, we rely on the spectroscopic galaxy sample developed by van den Busch et al. (2022) and extended by Wright et al. (2024), which includes 93 224 objects. Due to different selection effects, the redshift distribution of this spectroscopic sample differs remarkably from those of the photometric samples defined via the background selection (Eq. (6)). Thus, to reconstruct the true redshift distributions of the background galaxy samples, the magnitude and colour distributions of the spectroscopic sample have been weighted through the use of the SOM algorithm developed by Wright et al. (2020) and adopted by Hildebrandt et al. (2021) for the calibration of KiDS-1000 galaxy redshift distributions, based on the kohonen package (Wehrens & Buydens 2007; Wehrens & Kruisselbrink 2018). We remark that the background selection presented in Sect. 4.1 is applied galaxy-by-galaxy, while the SOM analysis provides a redshift calibration for a population of objects. For this reason, SOM results cannot be used in the background selection adopted for our g₊ measurements.

Based on our baseline background selection (Eq. (6)), we derived uncalibrated background redshift distributions given a cluster redshift, namely ${\widehat{n}}_{\mathrm{bkg}}(z_{\mathrm{g}}|z)$, in z bins having a width of δz = 0.05 and spanning the whole redshift range of the cluster sample. We remark that ${\widehat{n}}_{\mathrm{bkg}}(z_{\mathrm{g}}|z)$ is derived from the background galaxies associated with the clusters in the sample and falling within the cluster-centric projected radius range defined in Sect. 4.3. Then, we reconstructed the true background redshift distributions, n_bkg(z_g | z), following the methods described in Wright et al. (2020). Specifically, we provided the SOM code with the information on the 9 KiDS and VIKING photometric bands, along with all their colour combinations. We used SOMs with toroidal topology and 50 × 50 hexagonal cells, verifying that increasing the number of cells does not significantly affect the final results. As an illustration, the left panels in Fig. 8 show n_bkg(z_g | z) and ${\widehat{n}}_{\mathrm{bkg}}(z_{\mathrm{g}}|z)$ for z = 0.125 and z = 0.775. The median values of n_bkg(z_g | z) and ${\widehat{n}}_{\mathrm{bkg}}(z_{\mathrm{g}}|z)$, namely M(z) and $\widehat{M}(z)$, respectively, along with their percentage differences, are shown in the right panels of Fig. 8. As we can see, M is 8% larger than $\widehat{M}$ at low z, with this difference progressively reducing up to z = 0.5. For z > 0.5, $\widehat{M}$ becomes larger than M, with differences reaching about 10% at large z.

Fig. 8. Differences between the uncalibrated, $\widehat{n}(z_{\mathrm{g}}|z)$, and the SOM-weighted, n(zg | z), background redshift distributions. Left panels: n(zg | z) (blue histograms) and $\widehat{n}(z_{\mathrm{g}}|z)$ (red hatched histograms) for z = 0.125 (top) and z = 0.775 (bottom). The vertical dashed black lines represent the cluster redshifts. Right panels: In the top subpanel, the median values of n(zg | z) (solid blue curve) and $\widehat{n}(z_{\mathrm{g}}|z)$ (dashed red curve), namely M and $\widehat{M}$, respectively, are shown. The percentage differences between them are displayed in the bottom subpanel.

From n_bkg(z_g | z), we derived the purity of the background samples, namely 𝒫_bkg(z), defined as the ratio of the number of galaxies with z_g > z over the total number of selected galaxies. Figure 9 shows that 𝒫_bkg > 96% for z < 0.5, in agreement with L24. For z > 0.5, for which $\widehat{M}>M$, 𝒫_bkg decreases down to 90% at z = 0.75. Nevertheless, L24 predicted 𝒫_bkg > 97% for these values of z. This is mainly due to the uncertainties on KiDS photometry, which are larger compared to those affecting the galaxy samples considered by L24. To evaluate the overall impact of such impurities on the stacked cluster weak-lensing measurements, in Fig. 9 we also show ⟨𝒫_bkg(Δz_ob)⟩, that is the effective purity weighted over the cluster redshift distribution in a given Δz. We find that ⟨𝒫_bkg(Δz_ob)⟩ = 99% for z ∈ [0.1, 0.3), ⟨𝒫_bkg(Δz_ob)⟩ = 97% for z ∈ [0.3, 0.45), and ⟨𝒫_bkg(Δz_ob)⟩ = 92% for z ∈ [0.45, 0.8].

Fig. 9. Top panel: Background selection purity as a function of the cluster redshift. Bottom panel: Effective selection purity averaged over the cluster redshift distribution, in the cluster redshift bins z ∈ [0.1, 0.3), z ∈ [0.3, 0.45), and z ∈ [0.45, 0.8]. In both panels, the symbols are the same as those in Fig. 7. The black and green curves are almost overlapping.

As displayed in Fig. 9, 𝒫_bkg can be improved by up to four percentage points at high z by using the photo-z selection alone. Despite the fact that this could make the photo-z selection preferable for the last Δz bin, we do not expect this to have a significant impact on the final results. Indeed, as we shall discuss in the following, the final results on cosmology and cluster masses are dominated by low and intermediate redshift measurements (see also ⟨(S/N)_WL⟩_Δλ^* in Fig. 7). Furthermore, we note that the colour selection does not yield any background galaxies for cluster redshifts larger than z > 0.75, due to the buffer of 0.05 in the lens redshifts described in Sect. 4.1. Indeed, we remark that the colour selection by L24 used in this work is valid up to z = 0.8. Consequently, the combination of colour and photo-z selections yields the same results as the photo-z selection alone for z > 0.75, as shown in Fig. 9. Figure 9 also shows that the photo-z peak selection yields the lowest background purity at any z, with 𝒫_bkg ∼ 80% at z > 0.7. Together with the fact that this selection does not improve the weak-lensing S/N, as discussed in Sect. 4.3, we conclude that the combination of colour and photo-z selection is more robust for our analysis.

Furthermore, by reconstructing the background redshift distribution at a given z and cluster-centric projected distance, that is n_bkg(z_g | z, R), we find no significant dependence of 𝒫_bkg on R. This is expected, as the inner regions of the clusters, which may degrade the purity due to the contamination by cluster members, are excluded from the analysis. Consequently, we do not expect any significant dependence of 𝒫_bkg on λ^*. While previous studies (e.g. McClintock et al. 2019) have demonstrated significant cluster member contamination when selecting source galaxies based on their average redshift estimates, our background selection method (Sect. 4.1) minimises this contamination through two key features: a robust colour selection, and stricter galaxy photometric redshift requirements. The reduced shear unaffected by impurities, namely g₊^true, can be expressed as (see e.g. Dietrich et al. 2019)

$$\begin{array}{c}\hfill g_{+}^{\mathrm{true}}(z)=\frac{g_{+}(z)}{{\mathcal{P}}_{\mathrm{bkg}}(z)},\end{array}$$(17)

where g₊(z) is the measured reduced shear. Thus, in the following analysis we multiplied the theoretical model by 𝒫_bkg(z).

Using simulations, Wright et al. (2020) showed that the bias on the mean of the SOM-reconstructed redshift distributions is below 1%, on average, in the different tomographic bins considered for the KiDS-1000 cosmic shear analysis (see also Hildebrandt et al. 2021; Giblin et al. 2021). The statistical uncertainty on the mean redshift, accounting for photometric noise, shot-noise due to the limited sample size, spectroscopic selection effects and incompleteness, sample variance due to LSS, and inherent approximations in the simulations, amounts to at most 2% (Hildebrandt et al. 2021). Despite the differences between the galaxy photo-z estimates used in this work and those adopted for the KiDS-1000 cosmic shear analysis, stemming from the different priors used in BPZ by M25 (see Sect. 2.2), we expect that the uncertainties discussed above hold also for our sample. Since this uncertainty impacts the estimation of the ratio ℛ = D_ls/D_s (see Eq. (8)), we derived the mean ℛ, weighted by n_bkg(z_g | z), at the central values of the cluster redshift bins used in the analysis. We repeated this process by shifting n_bkg(z_g | z) by ±2% of its median. This results in a relative uncertainty on g₊ of 1% for the first two redshift bins and 4% for the last one. As discussed in Sect. 5.4, this uncertainty is propagated into the final posteriors.

5. Modelling

In this section, we detail the modelling used to simultaneously constrain the log λ^* − log M₂₀₀ scaling relation and cosmological parameters. Section 5.1 outlines the cluster count model. After introducing, in Sect. 5.2, the model describing a single galaxy cluster mass profile, we delve into the description of the theoretical expected values of the g₊ stacked measurements, in Sect. 5.3. Then, in Sect. 5.4, we present the joint Bayesian modelling of all the weak-lensing stacks derived in Sect. 4 and cluster abundance measurements.

5.1. Cluster count model

The expected number of clusters in a bin of observed intrinsic richness, Δλ_ob^*, and redshift, Δz_ob, is expressed as follows:

$$\begin{array}{cc}\hfill ⟨N(\mathrm{\Delta }{\lambda }_{\mathrm{ob}}^{\ast },\mathrm{\Delta }& z_{\mathrm{ob}})⟩=\frac{\mathrm{\Omega }}{{\mathcal{P}}_{\mathrm{clu}}(\mathrm{\Delta }{\lambda }_{\mathrm{ob}}^{\ast },\mathrm{\Delta }{z}_{\mathrm{ob}})}{\displaystyle {\int }_0^{\infty }}\mathrm{d}{z}_{\mathrm{tr}}\frac{{\mathrm{d}}^{2}V}{\mathrm{d}{z}_{\mathrm{tr}}\mathrm{d}\mathrm{\Omega }}\hfill \\ \hfill & \times {\displaystyle {\int }_0^{\infty }}\mathrm{d}M\frac{\mathrm{d}n(M,z_{\mathrm{tr}})}{\mathrm{d}M}{\mathcal{B}}_{\mathrm{HMF}}(M)\hfill \\ \hfill & \times {\displaystyle {\int }_0^{\infty }}\mathrm{d}{\lambda }_{\mathrm{tr}}^{\ast }{\mathcal{C}}_{\mathrm{clu}}({\lambda }_{\mathrm{tr}}^{\ast },z_{\mathrm{tr}})P({\lambda }_{\mathrm{tr}}^{\ast }|M,z_{\mathrm{tr}})\hfill \\ \hfill & \times {\displaystyle {\int }_{\mathrm{\Delta }{\lambda }_{\mathrm{ob}}^{\ast }}}\mathrm{d}{\lambda }_{\mathrm{ob}}^{\ast }P({\lambda }_{\mathrm{ob}}^{\ast }|{\lambda }_{\mathrm{tr}}^{\ast },z_{\mathrm{tr}}){\displaystyle {\int }_{\mathrm{\Delta }{z}_{\mathrm{ob}}}}\mathrm{d}{z}_{\mathrm{ob}}P(z_{\mathrm{ob}}|z_{\mathrm{tr}}),\hfill \end{array}$$(18)

where z_tr is the true redshift, V is the co-moving volume, Ω is the effective area covered by the cluster sample, M is the true mass, and dn(M, z_tr)/dM is the halo mass function, for which we adopted the model by Tinker et al. (2008). In Eq. (18), ${\mathcal{B}}_{\mathrm{HMF}}(M)$ is the halo mass function bias. Following Costanzi et al. (2019), ${\mathcal{B}}_{\mathrm{HMF}}(M)$ is expressed as follows:

$$\begin{array}{c}\hfill {\mathcal{B}}_{\mathrm{HMF}}(M)=slog\frac{M_{200\mathrm{m}}(M)}{M^{\ast }}+q,\end{array}$$(19)

where M^* is expressed in h⁻¹ M_⊙ and log M^* = 13.8, q and s are free parameters of the model (see Sect. 5.4), and M_200m(M) is the mass within a radius enclosing a mean density 200 times larger than the mean density of the Universe, at the halo redshift, given a generically defined mass M. Indeed, Costanzi et al. (2019) calibrated Eq. (19) assuming M_200m, while we adopted the critical overdensity to define masses throughout this paper. Thus, we estimated M_200m(M) using the mass conversion model by Ragagnin et al. (2021, Eq. (13)), testing our implementation against the HYDRO_MC code³. As we will discuss in Sect. 5.4, we propagated the uncertainties on the parameters entering Eq. (19) into the final results. This ensures that our cosmological constraints are robust against the assumption of alternative halo mass function models.

In Eq. (18), P(λ_tr^*|M, z_tr) is a log-normal probability density function, whose mean is given by the log λ^* − log M₂₀₀ scaling relation and its standard deviation is given by the intrinsic scatter, σ_intr:

$$\begin{array}{cc}\hfill P({\lambda }_{\mathrm{tr}}^{\ast }|M,z_{\mathrm{tr}})=& \frac{1}{ln(10){\lambda }_{\mathrm{tr}}^{\ast }\sqrt{2\pi }{\sigma }_{\mathrm{intr}}}exp(-\frac{{[log{\lambda }_{\mathrm{tr}}^{\ast }-\mu (M,z_{\mathrm{tr}})]}^2}{2{\sigma }_{\mathrm{intr}}^2}),\hfill \end{array}$$(20)

where μ(M, z_tr) is the mean of the distribution, expressed as

$$\begin{array}{c}\hfill \mu (M,z_{\mathrm{tr}})=\alpha +\beta log\frac{M}{M_{\mathrm{piv}}}+\gamma log\frac{H(z_{\mathrm{tr}})}{H(z_{\mathrm{piv}})}+log{\lambda }_{\mathrm{piv}}^{\ast },\end{array}$$(21)

where M_piv = 10¹⁴ h⁻¹ M_⊙, z_piv = 0.4, and λ_piv^* = 50 are the redshift, mass, and intrinsic richness pivots, respectively. The second last term in Eq. (21), including the Hubble function H(z), accounts for deviations in the redshift evolution from what is predicted in the self-similar growth scenario (Sereno & Ettori 2015). As detailed in Sect. 5.4, α, β, γ, and σ_intr in Eqs. (20) and (21) are free parameters in the modelling. Furthermore, 𝒫_clu and 𝒞_clu in Eq. (18) are the purity and completeness of the cluster sample, respectively, described in Sect. 3, while P(z_ob|z_tr) is a Gaussian probability density function with mean corresponding to z_tr and a standard deviation of 0.014(1 + z_tr) (see Sect. 2.1). Lastly, P(λ_ob^*|λ_tr^*, z_tr) is a Gaussian derived from the P(x|λ_tr^*, z_tr) distribution described in Sect. 3. Specifically, following Eq. (2), the mean of P(λ_ob^*|λ_tr^*, z_tr) is expressed as

$$\begin{array}{c}\hfill {\mu }_{{\lambda }^{\ast }}={\lambda }_{\mathrm{tr}}^{\ast }+{\lambda }_{\mathrm{tr}}^{\ast }exp[-{\lambda }_{\mathrm{tr}}^{\ast }({\overline{A}}_{\mu }+{\overline{B}}_{\mu }{z}_{\mathrm{tr}})],\end{array}$$(22)

while its standard deviation is based on Eq. (3) and is expressed as follows:

$$\begin{array}{c}\hfill {\sigma }_{{\lambda }^{\ast }}={\overline{A}}_{\sigma }{\lambda }_{\mathrm{tr}}^{\ast }exp(-{\overline{B}}_{\sigma }{\lambda }_{\mathrm{tr}}^{\ast }).\end{array}$$(23)

In Eqs. (22) and (23), ${\overline{A}}_{\mu }=0.198$, ${\overline{B}}_{\mu }=-0.179$, ${\overline{A}}_{\sigma }=0.320$, and ${\overline{B}}_{\sigma }=0.011$.

5.2. Basic lens model and miscentring

Here we present the models used to describe the mass profiles of single galaxy clusters. The expected value of a stacked weak-lensing profile is detailed in Sect. 5.3. Specifically, the 3D mass density profiles of single galaxy clusters are described by a truncated NFW profile (BMO, Baltz et al. 2009),

$$\begin{array}{c}\hfill \rho (r)=\frac{{\rho }_{\mathrm{s}}}{(r/r_{\mathrm{s}}){(1+r/r_{\mathrm{s}})}^2}{\left(\frac{r_{\mathrm{t}}^2}{r^2+r_{\mathrm{t}}^2}\right)}^2,\end{array}$$(24)

where ρ_s is the characteristic density, r_s is the scale radius, and r_t is the truncation radius. The scale radius is expressed as r_s = r₂₀₀/c₂₀₀, where r₂₀₀ is the radius enclosing a mass such that the corresponding mean density is 200 times the critical density of the Universe at that redshift, and c₂₀₀ is the concentration. As discussed in Sect. 5.3, c₂₀₀ is given by a log-linear relation with M₂₀₀ whose amplitude is constrained by our data. The truncation radius is defined as r_t = F_tr₂₀₀, where F_t is the truncation factor. The latter is a free parameter in the analysis, as detailed in Sect. 5.4. The BMO profile (Eq. (24)) provides more accurate mass estimates than a simple NFW parametrisation. Oguri & Hamana (2011) demonstrated that NFW-based fits to shear signals can introduce systematic biases on masses of about 10–15%.

We included a two-halo term in the model, whose contribution to the surface mass density is expressed as (Oguri & Takada 2011)

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{2\mathrm{h}}(\theta ,M,z)=\frac{{\rho }_{\text{m}}(z)b_{\text{h}}(M,z)}{{(1+z)}^3{D}_{\text{l}}^2(z)}{\displaystyle \int \frac{l\mathrm{d}l}{2\pi }}{J}_0(l\theta )P(k_{\text{l}},z),\end{array}$$(25)

where θ is the angular radius, J₀ is the 0th order Bessel function, k_l = l/(1 + z)/D_l(z), D_l is the lens angular diameter distance, b_h is the halo bias, for which the model by Tinker et al. (2010) was assumed, and P(k_l, z) is the linear matter power spectrum. Then the centred surface mass density is expressed as

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{cen}}(R)={\mathrm{\Sigma }}_{1\mathrm{h}}(R)+{\mathrm{\Sigma }}_{2\mathrm{h}}(R),\end{array}$$(26)

where Σ_1h(R), for which Oguri & Hamana (2011) proposed an analytic form, is derived from Eq. (24). The centred excess surface mass density has the following functional form:

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Sigma }}_{+,\mathrm{cen}}(R)=\frac{2}{R^2}{\displaystyle {\int }_0^R}\mathrm{d}rr{\mathrm{\Sigma }}_{\mathrm{cen}}(r)-{\mathrm{\Sigma }}_{\mathrm{cen}}(R).\end{array}$$(27)

From Eqs. (26) and (27), we express the tangential reduced shear component of a centred halo, g_+,cen, as in Seitz & Schneider (1997):

$$\begin{array}{c}\hfill g_{+,\mathrm{cen}}(R,M,z)=\frac{\mathrm{\Delta }{\mathrm{\Sigma }}_{+,\mathrm{cen}}(R,M,z)⟨{\mathrm{\Sigma }}_{\mathrm{crit}}^{-1}(z)⟩}{1-{\mathrm{\Sigma }}_{\mathrm{cen}}(R,M,z){⟨{\mathrm{\Sigma }}_{\mathrm{crit}}^{-1}(z)⟩}^{-1}⟨{\mathrm{\Sigma }}_{\mathrm{crit}}^{-2}(z)⟩}·\end{array}$$(28)

In Eq. (28), $⟨{\mathrm{\Sigma }}_{\mathrm{crit}}^{-\eta }⟩$ has the following expression

$$\begin{array}{c}\hfill ⟨{\mathrm{\Sigma }}_{\mathrm{crit}}^{-\eta }(z)⟩={\displaystyle \frac{{\int }_{z_{\mathrm{g}}>z}\mathrm{d}{z}_{\mathrm{g}}{\mathrm{\Sigma }}_{\mathrm{crit}}^{-\eta }(z_{\mathrm{g}},z)n(z_{\mathrm{g}}|z)}{{\int }_{z_{\mathrm{g}}>z}\mathrm{d}{z}_{\mathrm{g}}n(z_{\mathrm{g}}|z)}},\end{array}$$(29)

where Σ_crit is given by Eq. (8), while n(z_g | z) is the SOM-reconstructed background redshift distribution described in Sect. 4.4.

We also modelled the contribution to g₊ due to the miscentred population of clusters, as is done, for example, in Johnston et al. (2007) and Viola et al. (2015). We assumed that the probability of a lens being at a projected distance R_s from the chosen centre, namely P(R_s), follows a Rayleigh distribution, that is

$$\begin{array}{c}\hfill P(R_{\text{s}})=\frac{R_{\text{s}}}{{\sigma }_{\text{off}}^2}exp[-\frac{1}{2}(\frac{R_{\text{s}}}{{\sigma }_{\text{off}}}{)}^2],\end{array}$$(30)

where σ_off is the standard deviation of the halo misplacement distribution on the plane of the sky, expressed in units of h⁻¹ Mpc. The corresponding azimuthally averaged profile is given by (Yang et al. 2006)

$$\begin{array}{c}\hfill \mathrm{\Sigma }(R|R_{\text{s}})=\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }}{\mathrm{\Sigma }}_{\text{cen}}(\sqrt{R^2+R_{\text{s}}^2-2R{R}_{\text{s}}cos\theta })\mathrm{d}\theta ,\end{array}$$(31)

and the surface mass density distribution of a miscentred halo is expressed as

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\text{off}}(R)={\displaystyle \int P}(R_{\text{s}})\mathrm{\Sigma }(R|R_{\text{s}})\mathrm{d}{R}_{\text{s}}.\end{array}$$(32)

Analogously to Eq. (28), the tangential reduced shear component of a miscentred halo, g_+,off, has the following expression:

$$\begin{array}{c}\hfill g_{+,\mathrm{off}}(R,M,z)=\frac{\mathrm{\Delta }{\mathrm{\Sigma }}_{+,\mathrm{off}}(R,M,z)⟨{\mathrm{\Sigma }}_{\mathrm{crit}}^{-1}(z)⟩}{1-{\mathrm{\Sigma }}_{\mathrm{off}}(R,M,z){⟨{\mathrm{\Sigma }}_{\mathrm{crit}}^{-1}(z)⟩}^{-1}⟨{\mathrm{\Sigma }}_{\mathrm{crit}}^{-2}(z)⟩},\end{array}$$(33)

where Σ_off is given by Eq. (32), while ΔΣ_+,off is derived by replacing Σ_cen with Σ_off in Eq. (27).

5.3. Stacked weak-lensing model

The expected value of the stacked g₊ is expressed as follows:

$$\begin{array}{cc}\hfill ⟨g_{+}(R,\mathrm{\Delta }{\lambda }_{\mathrm{ob}}^{\ast },\mathrm{\Delta }{z}_{\mathrm{ob}})⟩=& (1-f_{\mathrm{off}})⟨g_{+,\mathrm{cen}}(R,\mathrm{\Delta }{\lambda }_{\mathrm{ob}}^{\ast },\mathrm{\Delta }{z}_{\mathrm{ob}})⟩\hfill \\ \hfill & +f_{\mathrm{off}}⟨g_{+,\mathrm{off}}(R,\mathrm{\Delta }{\lambda }_{\mathrm{ob}}^{\ast },\mathrm{\Delta }{z}_{\mathrm{ob}})⟩,\hfill \end{array}$$(34)

where f_off is the fraction of haloes that belong to the miscentred population, while ⟨g_+,cen⟩ and ⟨g_+,off⟩ are the average g_+,cen and g_+,off of the sample, respectively. Specifically, ⟨g_+,cen⟩ has the form

$$\begin{array}{cc}\hfill ⟨g_{+,\mathrm{cen}}(& R,\mathrm{\Delta }{\lambda }_{\mathrm{ob}}^{\ast },\mathrm{\Delta }{z}_{\mathrm{ob}})⟩=\frac{{\mathcal{P}}_{\mathrm{clu}}(\mathrm{\Delta }{\lambda }_{\mathrm{ob}}^{\ast },\mathrm{\Delta }{z}_{\mathrm{ob}})⟨{\mathcal{P}}_{\mathrm{bkg}}(\mathrm{\Delta }{z}_{\mathrm{ob}})⟩}{⟨n(\mathrm{\Delta }{\lambda }_{\mathrm{ob}}^{\ast },\mathrm{\Delta }{z}_{\mathrm{ob}})⟩}\hfill \\ \hfill & \times {\displaystyle {\int }_0^{\infty }\mathrm{d}{z}_{\mathrm{tr}}}\frac{{\mathrm{d}}^{2}V}{\mathrm{d}{z}_{\mathrm{tr}}\mathrm{d}\mathrm{\Omega }}{\displaystyle {\int }_0^{\infty }}\mathrm{d}M\frac{\mathrm{d}n(M,z_{\mathrm{tr}})}{\mathrm{d}M}{\mathcal{B}}_{\mathrm{HMF}}(M)\hfill \\ \hfill & \times g_{+,\mathrm{cen}}(R^{\mathrm{test}},M,z_{\mathrm{tr}}){\displaystyle {\int }_0^{\infty }}\mathrm{d}{\lambda }_{\mathrm{tr}}^{\ast }{\mathcal{C}}_{\mathrm{clu}}({\lambda }_{\mathrm{tr}}^{\ast },z_{\mathrm{tr}})P({\lambda }_{\mathrm{tr}}^{\ast }|M,z_{\mathrm{tr}})\hfill \\ \hfill & \times {\displaystyle {\int }_{\mathrm{\Delta }{\lambda }_{\mathrm{ob}}^{\ast }}}\mathrm{d}{\lambda }_{\mathrm{ob}}^{\ast }P({\lambda }_{\mathrm{ob}}^{\ast }|{\lambda }_{\mathrm{tr}}^{\ast },z_{\mathrm{tr}}){\displaystyle {\int }_{\mathrm{\Delta }{z}_{\mathrm{ob}}}}\mathrm{d}{z}_{\mathrm{ob}}P(z_{\mathrm{ob}}|z_{\mathrm{tr}}),\hfill \end{array}$$(35)

while ⟨g_+,off⟩ is obtained by replacing g_+,cen with g_+,off in Eq. (35). In Eq. (35), n(Δλ_ob^*, Δz_ob) is the expected density of observed haloes not corrected for sample impurities, having the following expression:

$$\begin{array}{cc}\hfill ⟨n(\mathrm{\Delta }{\lambda }_{\mathrm{ob}}^{\ast },\mathrm{\Delta }& z_{\mathrm{ob}})⟩={\displaystyle {\int }_0^{\infty }}\mathrm{d}{z}_{\mathrm{tr}}\frac{{\mathrm{d}}^{2}V}{\mathrm{d}{z}_{\mathrm{tr}}\mathrm{d}\mathrm{\Omega }}{\displaystyle {\int }_0^{\infty }}\mathrm{d}M\frac{\mathrm{d}n(M,z_{\mathrm{tr}})}{\mathrm{d}M}{\mathcal{B}}_{\mathrm{HMF}}(M)\hfill \\ \hfill & \times {\displaystyle {\int }_0^{\infty }}\mathrm{d}{\lambda }_{\mathrm{tr}}^{\ast }{\mathcal{C}}_{\mathrm{clu}}({\lambda }_{\mathrm{tr}}^{\ast },z_{\mathrm{tr}})P({\lambda }_{\mathrm{tr}}^{\ast }|M,z_{\mathrm{tr}})\hfill \\ \hfill & \times {\displaystyle {\int }_{\mathrm{\Delta }{\lambda }_{\mathrm{ob}}^{\ast }}}\mathrm{d}{\lambda }_{\mathrm{ob}}^{\ast }P({\lambda }_{\mathrm{ob}}^{\ast }|{\lambda }_{\mathrm{tr}}^{\ast },z_{\mathrm{tr}}){\displaystyle {\int }_{\mathrm{\Delta }{z}_{\mathrm{ob}}}}\mathrm{d}{z}_{\mathrm{ob}}P(z_{\mathrm{ob}}|z_{\mathrm{tr}}).\hfill \end{array}$$(36)

The cluster sample purity, 𝒫_clu, in Eq. (35) quantifies the suppression of the stacked weak-lensing signal due to false cluster detections, while ⟨𝒫_bkg⟩ is the effective background sample purity, described in Sect. 4.4. We remark that the contributions by halo triaxiality and projection effects were included in the covariance matrix, as discussed in Sect. 5.4.

The measurements presented in Sect. 4 were performed under the assumption of a fiducial cosmology. To address geometric distortions arising from this assumption, g_+,cen in Eq. (35) is evaluated at the test projected radius R^test, expressed as

$$\begin{array}{c}\hfill R^{\mathrm{test}}=\theta D_{\mathrm{l}}^{\mathrm{test}}=R^{\mathrm{fid}}\frac{D_{\mathrm{l}}^{\mathrm{test}}}{D_{\mathrm{l}}^{\mathrm{fid}}},\end{array}$$(37)

where θ is the angular separation from the cluster centre, R^fid is the projected radius in the fiducial cosmology, while $D_{\mathrm{l}}^{\mathrm{fid}}$ and $D_{\mathrm{l}}^{\mathrm{test}}$ represent the angular diameter distance of the cluster in the fiducial and test cosmologies, respectively.

Lastly, we assumed a log-linear model for the c₂₀₀ − M₂₀₀ relation, which is defined as

$$\begin{array}{c}\hfill log{c}_{200}=log{c}_0+c_{M}log\frac{M}{M_{\mathrm{piv}}}+c_{z}log\frac{1+z_{\mathrm{tr}}}{1+z_{\mathrm{piv}}},\end{array}$$(38)

where M_piv and z_piv are the same as those assumed in Eq. (21). As discussed in the following, log c₀ is a free parameter in the analysis, while c_M and c_z are fixed to fiducial values. We point out that Eq. (38) neglects the intrinsic scatter in the log c₂₀₀ − log M₂₀₀ relation. In Sect. 6 we address the impact of this scatter on our results.

5.4. Likelihood function and parameter priors

In this paper, we performed a joint Bayesian analysis of all the stacked weak-lensing and cluster count measurements, by means of an MCMC algorithm. The likelihood function is expressed as

$$\begin{array}{c}\hfill \mathcal{L}={\mathcal{L}}^{\mathrm{WL}}{\mathcal{L}}^{\mathrm{C}},\end{array}$$(39)

where ${\mathcal{L}}^{\mathrm{WL}}$ and ${\mathcal{L}}^{\mathrm{C}}$ are the cluster weak-lensing and count likelihoods, respectively. Notably, ${\mathcal{L}}^{\mathrm{WL}}$ takes the form

$$\begin{array}{c}\hfill {\mathcal{L}}^{\mathrm{WL}}={\displaystyle \prod _{i=1}^{N_{{\lambda }^{\ast }}^{\mathrm{WL}}}}{\displaystyle \prod _{j=1}^{N_z^{\mathrm{WL}}}}{\mathcal{L}}_{\mathit{ij}}^{\mathrm{WL}},\end{array}$$(40)

where $N_{{\lambda }^{\ast }}^{\mathrm{WL}}$ and $N_z^{\mathrm{WL}}$ are the numbers of cluster intrinsic richness and redshift bins used for the stacked measurements, respectively, while ${\mathcal{L}}_{\mathit{ij}}^{\mathrm{WL}}$ is a Gaussian likelihood function defined as

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathit{ij}}^{\mathrm{WL}}\propto exp(-{\chi }_{\mathit{ij}}^2/2),\end{array}$$(41)

with

$$\begin{array}{cc}& {\chi }_{\mathit{ij}}^2={\displaystyle \sum _{k=1}^{N_R}}{\displaystyle \sum _{l=1}^{N_R}}(g_{+,ijk}^{\mathrm{ob}}-g_{+,ijk}^{\mathrm{mod}})C_{\mathit{kl}}^{-1}(g_{+,ijl}^{\mathrm{ob}}-g_{+,ijl}^{\mathrm{mod}}).\hfill \end{array}$$(42)

Here, $g_{+}^{\mathrm{ob}}$ represents the observed stacked profile, given by Eq. (12), $g_{+}^{\mathrm{mod}}$ is the model, namely Eq. (34), while the indices k and l run over the number of radial bins, N_R, and $C_{\mathit{kl}}^{-1}$ is the inverse of the covariance matrix. In particular, C_kl is defined as

$$\begin{array}{c}\hfill C_{\mathit{kl}}=C_{\mathit{kl}}^{\mathrm{BT}}+C_{\mathit{kl}}^{\mathrm{sys}},\end{array}$$(43)

where $C_{\mathit{kl}}^{\mathrm{BT}}$ is estimated through a bootstrap resampling, as discussed in Sect. 4.3, while $C_{\mathit{kl}}^{\mathrm{sys}}$ accounts for residual uncertainties on systematic errors that are not included in the model (Eq. (35)). Specifically, $C_{\mathit{kl}}^{\mathrm{sys}}$ is written as

$$\begin{array}{c}\hfill C_{\mathit{kl}}^{\mathrm{sys}}=({\sigma }_m^2+{\sigma }_{\mathrm{SOM}}^2+{\sigma }_{\mathrm{OP}}^2)g_{+,k}^{\mathrm{ob}}{g}_{+,l}^{\mathrm{ob}},\end{array}$$(44)

where σ_m = 0.02 is the uncertainty on the multiplicative shear bias (Sect. 4.2), ${\sigma }_{\mathrm{SOM}}$ accounts for the uncertainty on the SOM-reconstructed background redshift distributions, amounting to ${\sigma }_{\mathrm{SOM}}=0.01$ for the first two cluster redshift bins and to ${\sigma }_{\mathrm{SOM}}=0.04$ for the last one (see Sect. 4.4), while ${\sigma }_{\mathrm{OP}}=0.03$ is the residual uncertainty due to orientation and projection effects. Indeed, optical cluster finders preferentially select haloes with the major axis aligned with the line of sight, as these objects produce a larger density contrast with respect to the distribution of the field galaxies. In stacked weak-lensing analyses, this implies a boost in the weak-lensing signal (Corless & King 2008; Dietrich et al. 2014; Osato et al. 2018; Zhang et al. 2023; Euclid Collaboration: Giocoli et al. 2024). An opposite mass bias comes from projections of secondary haloes aligned with the detected clusters, which are in fact blended in a single detection. This may cause an overestimation of the cluster optical mass proxies, and a consequent bias in the mass-observable relation (Myles et al. 2021; Wu et al. 2022). For the AMICO clusters in KiDS-DR3, Bellagamba et al. (2019) found that orientation effects counterbalance those caused by projections, leading to a negligible bias on the derived masses with a residual uncertainty of 3%. Also Simet et al. (2017) and Melchior et al. (2017), for galaxy clusters detected applying the red-sequence Matched-filter Probabilistic Percolation (redMaPPer, Rykoff et al. 2014, 2016) algorithm in the Dark Energy Survey (DES, The Dark Energy Survey Collaboration 2005), found that the combination of these bias sources is consistent with zero, deriving a mass corrective factor of 0.98 ± 0.03 (see also McClintock et al. 2019). Thus, we conservatively assumed a residual uncertainty on mass due to orientation and projections of 4%. To relate this uncertainty on mass to that on weak-lensing profiles, we express the logarithmic dependence of ΔΣ₊ on the mass M as (Melchior et al. 2017)

$$\begin{array}{c}\hfill \mathrm{\Gamma }=\frac{\mathrm{d}ln\mathrm{\Delta }{\mathrm{\Sigma }}_{+}(M)}{\mathrm{d}lnM},\end{array}$$(45)

with Γ ≃ 0.7 being a good approximation for a broad range of cluster masses, redshifts, and cluster-centric distances (Melchior et al. 2017; Sereno et al. 2017). From Eq. (45), we obtain

$$\begin{array}{c}\hfill \frac{\delta \mathrm{\Delta }{\mathrm{\Sigma }}_{+}}{\mathrm{\Delta }{\mathrm{\Sigma }}_{+}}=\mathrm{\Gamma }\frac{\delta M}{M}·\end{array}$$(46)

Based on Eq. (46), assuming a 4% uncertainty on mass due to halo orientation and projection effects results in a relative uncertainty of about 3%, that is ${\sigma }_{\mathrm{OP}}=0.03$, on ΔΣ₊ and, in turn, on g₊.

Following Lima & Hu (2004) and Lacasa & Grain (2019), the likelihood ${\mathcal{L}}^{\mathrm{C}}$ in Eq. (39) is expressed as the convolution of Poisson and Gaussian distributions:

$$\begin{array}{cc}\hfill {\mathcal{L}}^{\mathrm{C}}={\displaystyle \int \mathrm{d}}{\mathit{\delta }}_{\mathrm{b}}^{N_z^{\mathrm{C}}}& \left[{\displaystyle \prod _{i=1}^{N_{{\lambda }^{\ast }}^{\mathrm{C}}}}{\displaystyle \prod _{j=1}^{N_z^{\mathrm{C}}}}\text{Poiss}(N_{\mathit{ij}}^{\mathrm{ob}}|N_{\mathit{ij}}^{\mathrm{mod}}+\frac{\mathrm{d}{N}_{\mathit{ij}}^{\mathrm{mod}}}{\mathrm{d}{\delta }_{\mathrm{b}}}{\delta }_{\mathrm{b},j})\right]\hfill \\ \hfill & \times \mathcal{N}({\mathit{\delta }}_{\mathrm{b}}|0,S),\hfill \end{array}$$(47)

where N^ob is the observed number of clusters, N^mod corresponds to the cluster count model (Eq. (18)), while dN^mod/dδ_b is the response of the counts to the variation δ_b of the background matter density. In particular, the cluster count response is analogous to Eq. (18), where the integrand function is multiplied by the linear halo bias (see e.g. Lacasa & Grain 2019; Lesci et al. 2022a), for which we adopted the model by Tinker et al. (2010). In Eq. (47), 𝒩(δ_b|0, S) is a multivariate Gaussian probability density function describing the super-sample covariance (SSC) effects on cluster count measurements, which is a function of the background density fluctuation, δ_b, is centred on zero, and has covariance matrix S. Notably, ${\mathit{\delta }}_{\mathrm{b}}=\{{\delta }_{\text{b},1},\dots ,{\delta }_{\text{b},N_z^{\mathrm{C}}}\}$, where $N_z^{\mathrm{C}}$ is the number of redshift bins used for cluster count measurements. Consequently, the S matrix has dimension $N_z^{\mathrm{C}}\times N_z^{\mathrm{C}}$, and the matrix element S_pq has the following functional form (Lacasa & Grain 2019):

$$\begin{array}{c}\hfill S_{\mathit{pq}}=\frac{1}{8{\pi }^3{f}_{\mathrm{sky}}^2}{\displaystyle \sum _{\ell }}(2\ell +1){\mathcal{C}}_{\mathrm{m}}(\ell ){\displaystyle \int \mathrm{d}}k{k}^{2}P(k)\frac{U_p(k,\ell )}{I_p}\frac{U_q(k,\ell )}{I_q},\end{array}$$(48)

where f_sky = Ω/(4π) is the sky fraction covered by the survey, Ω is the survey area in steradians, 𝒞_m(ℓ) is the angular power spectrum of the footprint mask computed at the multipole ℓ, accounting for the masked regions and the survey geometry, and P(k) is the linear matter power spectrum computed at z = 0. In addition, U_p(k) in Eq. (48) is expressed as

$$\begin{array}{c}\hfill U_p(k,\ell )={\displaystyle \int \frac{\mathrm{d}V}{\mathrm{d}z}}\mathrm{d}z{W}_p(z)G(z)j_{\ell }[kr(z)],\end{array}$$(49)

where W_p(z) is a top-hat window function, G(z) is the linear growth factor, j_ℓ is the spherical Bessel function, and r(z) is the comoving distance, while I_p in Eq. (48) has the expression

$$\begin{array}{c}\hfill I_p={\displaystyle \int \frac{\mathrm{d}V}{\mathrm{d}z}}\mathrm{d}z{W}_p(z).\end{array}$$(50)

We note that S and the cluster count response in Eq. (47) did not vary in our analysis. Indeed, these quantities were derived at the parameter posterior medians obtained from the MCMC analysis performed without including the SSC contribution. As the SSC was included in the likelihood, the MCMC analysis was repeated using the updated parameter median values until the posteriors converged.

Fig. 10. Top panels: Counts of the AMICO KiDS-1000 galaxy clusters in the redshift bins adopted in the analysis (increasing from left to right). The black dots show the measures, with the error bars corresponding to Poisson uncertainties. The blue bands display the 68% confidence levels of the model, derived from the posterior of all the free parameters considered in the joint analysis of counts and weak lensing. Bottom panels: Pearson residuals. The horizontal dashed grey lines show the interval between −1 and 1.

In Table 2 we list the parameter priors adopted in the analysis. Wide uniform priors were assumed for the cold dark matter density parameter at z = 0, denoted as ${\mathrm{\Omega }}_{\mathrm{CDM}}$, and for the amplitude of the primordial matter power spectrum, A_s, such that Ω_m and σ₈ are derived parameters in our analysis. We adopted Gaussian priors on the baryon density parameter, Ω_b, and on the primordial spectral index, n_s, based on the results by Planck Collaboration VI (2020, Table 2, TT, TE + EE + lowE + lensing, referred to as Planck18 hereafter), assuming the same median values and imposing a standard deviation equal to 5σ. The prior on the normalised Hubble constant, h, is a Gaussian distribution centred on 0.7 with a standard deviation of 0.03. In addition, we assumed uniform priors on the log λ^* − log M₂₀₀ scaling relation parameters in Eq. (21), namely α, β, and γ, and on its intrinsic scatter, σ_intr. A uniform prior was assumed also for the amplitude of the log c₂₀₀ − log M₂₀₀ relation (Eq. (38)), namely log c₀. Following Duffy et al. (2008), we fixed c_M = −0.084 and c_z = −0.47. Through these priors on σ_intr and c₂₀₀, we expect to account for the theoretical uncertainties on the contribution due to baryonic matter to the galaxy cluster profiles (Schaller et al. 2015; Henson et al. 2017; Shirasaki et al. 2018; Lee et al. 2018; Beltz-Mohrmann & Berlind 2021; Shao & Anbajagane 2024). Furthermore, we assumed a Gaussian prior with mean 0.3 and standard deviation 0.1 on the fraction of miscentred clusters, f_off, imposing also f_off > 0. For the typical miscentring scale, we adopted the uniform prior σ_off ∈ [0, 0.5] h⁻¹ Mpc. These priors on the offset of galaxy-based cluster centres agree with the results from simulated (Yan et al. 2020; Sommer et al. 2024) and observed (Saro et al. 2015; Zhang et al. 2019; Seppi et al. 2023; Ding et al. 2025) galaxy, intra-cluster medium, and dark matter distributions. We considered a Gaussian prior on the truncation factor of the BMO profile, F_t, with mean equal to 3 and standard deviation of 0.5, in agreement with Oguri & Hamana (2011). In addition, we imposed F_t > 0.

In Eq. (19) we introduced the ${\mathcal{B}}_{\mathrm{HMF}}$ correction factor for the Tinker et al. (2008) halo mass function. Following Costanzi et al. (2019), we adopted a 2D Gaussian prior on the s and q parameters in Eq. (19)⁴, having mean ${\mu }_{\mathrm{HMF}}=(0.037,1.008)$ and covariance matrix $C_{\mathrm{HMF}}$ expressed as

$$\begin{array}{c}\hfill C_{\mathrm{HMF}}=\left(\begin{array}{cc}0.00019& 0.00024\\ 0.00024& 0.00038\end{array}\right).\end{array}$$(51)

Despite the fact that this prescription does not account for the uncertainties due to baryonic physics, the latter are subdominant compared to the precision on ${\mathcal{B}}_{\mathrm{HMF}}$ (see Costanzi et al. 2019, and references therein).

6. Results

In this section, we detail the constraints on cosmological parameters and on the log λ^* − log M₂₀₀ relation from the joint modelling of the stacked weak-lensing reduced shear profiles and counts of the AMICO KiDS-1000 clusters, carried out through a Bayesian analysis based on the likelihood function in Eq. (39) and assuming the priors listed in Table 2. We conducted a blind analysis by introducing biases in the completeness estimates of the cluster sample (see Sect. 3), with the results detailed in Appendix B. While the cluster redshift cut at z = 0.8 was driven by the low completeness and purity of the sample at higher redshifts (see Sect. 2.1), the λ^* selection criteria were optimised to maintain robust model fits while maximising the cluster sample size. Compared to the weak-lensing analysis (Sect. 4.3), we adopted more stringent cuts in λ^* for the counts, namely λ^* > 25 for z ∈ [0.1, 0.3), λ^* > 28 for z ∈ [0.3, 0.45), and λ^* > 30 for z ∈ [0.45, 0.8]. We verified that this sample selection has no significant impact on the final constraints and improves the quality of the cluster count fit. Differences between the lensing and count sample selections are expected, since low-λ^* count measurements are very precise and more prone to unmodelled systematics.

As mentioned in Sect. 5.4, we first run the pipeline without including the SSC. The median parameter values obtained from this initial run were then used to compute the S matrix and the cluster count response entering Eq. (47). This process was repeated iteratively until the posteriors converge. The resulting model values are overplotted to the stacked weak-lensing data in Fig. 6 and to the cluster count measurements in Fig. 10. By considering the weak-lensing measurements only, we obtained a reduced χ² of χ_red² = 1.1. Here, in the computation of the degrees of freedom, we neglected the parameters entering the mass function correction factor, along with Ω_b, n_s, and h, as their priors have no impact on the weak-lensing modelling. We found a good quality of the fit also for the counts, as shown by the Pearson residuals in the bottom panels of Fig. 10.

In Sect. 6.1 we discuss the constraints on cosmological parameters and their comparison with the results from literature analyses, while Sect. 6.2 focuses on the log λ^* − log M₂₀₀ relation. Finally, in Section 6.3, we assess the robustness of our findings against alternative modelling choices.

6.1. Cosmological constraints

From the initial MCMC run that does not include the SSC contribution, we obtained ${\mathrm{\Omega }}_{\mathrm{m}}=0.{21}_{-0.01}^{+0.02}$, ${\sigma }_8=0.{86}_{-0.03}^{+0.03}$, and $S_8\equiv {\sigma }_8{({\mathrm{\Omega }}_{\mathrm{m}}/0.3)}^{0.5}=0.{72}_{-0.02}^{+0.02}$. The inclusion of the SSC led to the following constraints:

$$\begin{array}{cc}\hfill {\mathrm{\Omega }}_{\mathrm{m}}& =0.{218}_{-0.021}^{+0.024},\hfill \\ \hfill {\sigma }_8& =0.{86}_{-0.03}^{+0.03},\hfill \\ \hfill S_8& =0.{74}_{-0.03}^{+0.03}.\hfill \end{array}$$(52)

We verified that the SSC iterative process converges after one iteration. These constraints, along with the ones on the other free model parameters considered in the analysis, are reported in Table 2. We find a 1σ agreement with the results by Lesci et al. (2022a), who obtained ${\mathrm{\Omega }}_{\mathrm{m}}=0.{24}_{-0.04}^{+0.03}$, ${\sigma }_8=0.{86}_{-0.07}^{+0.07}$, and $S_8=0.{78}_{-0.04}^{+0.04}$ by jointly modelling the cluster weak-lensing and count measurements of the AMICO clusters in KiDS-DR3. As we can see, thanks to the larger survey area and the extension in the covered cluster redshift range in our analysis (see Sect. 2.1), the statistical uncertainties on Ω_m and σ₈ are halved, while the one on S₈ is reduced by 25%. Defining the figure of merit (FoM) as

$$\begin{array}{c}\hfill \mathrm{FoM}=\frac{1}{\sqrt{\det (C({\mathrm{\Omega }}_{\mathrm{m}},{\sigma }_8))}},\end{array}$$(53)

where C(Ω_m, σ₈) is the Ω_m − σ₈ covariance matrix, we obtain FoM = 1680. Compared to FoM = 701 from the constraints by Lesci et al. (2022a), the FoM more than doubles in KiDS-1000.

Figure 11 displays the comparison with Lesci et al. (2022a) results in the Ω_m − σ₈ parameter space, along with the constraints by Planck18 based on observations of the CMB power spectrum. The difference with Planck18 derived by Lesci et al. (2022a) was of 2σ on Ω_m and of about 1σ on σ₈ and S₈⁵. Due to the lower median value and uncertainty, in this work we find a 4.5σ tension with Planck18 on Ω_m, while σ₈ agrees within 1.6σ with Planck18. The lower Ω_m compared to that obtained by Lesci et al. (2022a) leads to a 2.8σ tension on S₈ with Planck18. A tension on Ω_m, of about 3σ, is found also with the results from the latest baryon acoustic oscillation (BAO) measurements based on the Dark Energy Spectroscopic Instrument (DESI, Adame et al. 2025). Indeed, assuming a flat ΛCDM model, they found Ω_m = 0.295 ± 0.015.

Fig. 11. Constraints on Ωm and σ8, obtained from the joint modelling of cluster weak-lensing and count measurements described in this work (blue), by Lesci et al. (2022a) (grey), and by Planck18 (orange). The confidence contours correspond to 68% and 95%, while the bands over the 1D marginalised posteriors represent the interval between 16th and 84th percentiles.

Our results are compared with other literature analyses in Fig. 12. The S₈ constraint agrees within 1σ with the cluster count analysis by Costanzi et al. (2019), based on the photometric data from the SDSS Data Release 8 (SDSS-DR8, Aihara et al. 2011), as well as with the counts of South Pole Telescope SZ survey (SPT-SZ, de Haan et al. 2016) clusters by Bocquet et al. (2019), based on mass estimates from Chandra X-ray observations and weak-lensing measurements from the Hubble Space Telescope (HST) and the Magellan Telescope. We find agreement on S₈ also with the analysis by Bocquet et al. (2024a), based on SZ-selected clusters in SPT and weak-lensing data from DES and HST. The result by Abbott et al. (2020), derived by modelling the abundance of DES-Y1 (Drlica-Wagner et al. 2018) optically selected clusters, agrees within 2σ. This DES-Y1 result was confirmed by Aguena et al. (2023), who employed a different modelling of projection effects and a different MCMC sampler. Our S₈ result shows a 2.7σ difference against the latest DES-Y3 result by DES Collaboration (2025), which combines cluster and galaxy statistics. Here, a key methodological distinction is our inclusion of small-scale cluster lensing, which DES Collaboration (2025) exclude from their analysis. A large disagreement on S₈, of about 3.8σ, is found with the counts of X-ray-selected clusters in the first eROSITA All-Sky Survey (eRASS1, Bulbul et al. 2024) by Ghirardini et al. (2024), due to the 2.8σ difference in the Ω_m constraint. Indeed, our σ₈ constraint agrees within 1σ with the result from Ghirardini et al. (2024), as well as with the findings from the other cluster count analyses mentioned earlier. The latter also show an agreement of about 1σ with our result on Ω_m.

Fig. 12. Constraints on S8 (left panel), σ8 (middle panel), and Ωm (right panel) obtained, from top to bottom, in this work (blue), by Planck18 (orange), Lesci et al. (2022a) (grey), Costanzi et al. (2019) (green), Abbott et al. (2020) (red), DES Collaboration (2025) (dark blue), Bocquet et al. (2019) (brown), Bocquet et al. (2024a) (violet), Ghirardini et al. (2024) (cyan), Secco et al. (2022) (pink), Dalal et al. (2023) (dark green), Asgari et al. (2021) (black), and Wright et al. (2025a) (dark grey). The median as well as the 16th and 84th percentiles are shown.

Figure 12 shows that analyses based on optically selected clusters tend to yield lower values of Ω_m compared to Planck18. Following the DES-Y1 cosmological results by Abbott et al. (2020), extensive work was conducted to evaluate the impact of optical projection effects on cluster mass estimates (Sunayama et al. 2020; Wu et al. 2022; Sunayama 2023; Park et al. 2023; Zhou et al. 2024) and, in turn, on cosmological constraints. As discussed in Sect. 6.3, we are confident that our results are not subject to such biases.

Furthermore, as shown in Fig. 12, the S₈ result derived in this work is in excellent agreement with the cosmic shear analyses by Asgari et al. (2021), Secco et al. (2022), and Dalal et al. (2023), based on KiDS-1000 (Giblin et al. 2021), DES-Y3 (Gatti et al. 2021), and Hyper Suprime-Cam Year 3 (HSC-Y3, Li et al. 2022) galaxy shape catalogues, respectively. We find a 2σ difference with the S₈ constraint by Wright et al. (2025a), who analysed the cosmic shear measurements from the latest KiDS-Legacy data (Li et al. 2023; Wright et al. 2024). The same level of discrepancy on S₈ was found by Wright et al. (2025a) against the KiDS-1000 result by Asgari et al. (2021).

6.2. Mass calibration

Regarding the parameters of the log λ^* − log M₂₀₀ relation, we derive $\alpha ={-0.170}_{-0.020}^{+0.017}$, $\beta ={0.553}_{-0.021}^{+0.021}$, $\gamma ={0.32}_{-0.26}^{+0.28}$, and ${\sigma }_{\mathrm{intr}}={0.052}_{-0.015}^{+0.023}$ (see also Table 2 and Fig. 13). The result on σ_intr confirms the goodness of λ^* as a mass proxy, as already found for AMICO KiDS-DR3 clusters by Sereno et al. (2020) and Lesci et al. (2022a). In addition, this constraint agrees with literature results on the relation between cluster richness and weak-lensing mass based on real data (see e.g. Bleem et al. 2020; Ghirardini et al. 2024) and simulations (Giocoli et al. 2025). Given the constraints on the log λ^* − log M₂₀₀ relation and on cosmological parameters, the expected mass value for a single galaxy cluster can be obtained from the formula

$$\begin{array}{cc}\hfill ⟨M(& {\lambda }_{\mathrm{ob}}^{\ast },z_{\mathrm{ob}})⟩=\frac{1}{⟨n({\lambda }_{\mathrm{ob}}^{\ast },z_{\mathrm{ob}})⟩}\hfill \\ \hfill & \times {\displaystyle {\int }_0^{\infty }}\mathrm{d}{z}_{\mathrm{tr}}\frac{{\mathrm{d}}^{2}V}{\mathrm{d}{z}_{\mathrm{tr}}\mathrm{d}\mathrm{\Omega }}{\displaystyle {\int }_0^{\infty }}\mathrm{d}M\frac{\mathrm{d}n(M,z_{\mathrm{tr}})}{\mathrm{d}M}{\mathcal{B}}_{\mathrm{HMF}}(M)M\hfill \\ \hfill & \times {\displaystyle {\int }_0^{\infty }}\mathrm{d}{\lambda }_{\mathrm{tr}}^{\ast }{\mathcal{C}}_{\mathrm{clu}}({\lambda }_{\mathrm{tr}}^{\ast },z_{\mathrm{tr}})P({\lambda }_{\mathrm{tr}}^{\ast }|M,z_{\mathrm{tr}})P({\lambda }_{\mathrm{ob}}^{\ast }|{\lambda }_{\mathrm{tr}}^{\ast },z_{\mathrm{tr}})P(z_{\mathrm{ob}}|z_{\mathrm{tr}}),\hfill \end{array}$$(54)

Fig. 13. Constraints on the log λ* − log M200 relation parameters, obtained from the joint modelling of cluster weak-lensing and count measurements. The confidence contours correspond to 68% and 95%, while the bands over the 1D marginalised posteriors represent the interval between the 16th and 84th percentiles.

where

$$\begin{array}{cc}\hfill ⟨n& ({\lambda }_{\mathrm{ob}}^{\ast },z_{\mathrm{ob}})⟩={\displaystyle {\int }_0^{\infty }}\mathrm{d}{z}_{\mathrm{tr}}\frac{{\mathrm{d}}^{2}V}{\mathrm{d}{z}_{\mathrm{tr}}\mathrm{d}\mathrm{\Omega }}{\displaystyle {\int }_0^{\infty }}\mathrm{d}M\frac{\mathrm{d}n(M,z_{\mathrm{tr}})}{\mathrm{d}M}{\mathcal{B}}_{\mathrm{HMF}}(M)\hfill \\ \hfill & \times {\displaystyle {\int }_0^{\infty }}\mathrm{d}{\lambda }_{\mathrm{tr}}^{\ast }{\mathcal{C}}_{\mathrm{clu}}({\lambda }_{\mathrm{tr}}^{\ast },z_{\mathrm{tr}})P({\lambda }_{\mathrm{tr}}^{\ast }|M,z_{\mathrm{tr}})P({\lambda }_{\mathrm{ob}}^{\ast }|{\lambda }_{\mathrm{tr}}^{\ast },z_{\mathrm{tr}})P(z_{\mathrm{ob}}|z_{\mathrm{tr}}).\hfill \end{array}$$(55)

By computing Eq. (54) at each MCMC step for each cluster in the sample, considering the whole parameter space defined in Table 2, we obtained mean mass estimates and the corresponding standard deviations. The minimum cluster mass derived in this way is $M_{200}^{\min }=3\times {10}^{13}$h⁻¹ M_⊙ for z ∈ [0.1, 0.3), $M_{200}^{\min }=4\times {10}^{13}$h⁻¹ M_⊙ for z ∈ [0.3, 0.45), and $M_{200}^{\min }=5\times {10}^{13}$h⁻¹ M_⊙ for z ∈ [0.45, 0.8]. Furthermore, the median mass precision corresponds to about 8%. This precision shows a three percentage points enhancement compared to what derived by B19 from the AMICO KiDS-DR3 sample (Maturi et al. 2019). Thus, despite the fact that the galaxy weighted density in KiDS-1000, which amounts to 6.17 arcmin⁻² (Giblin et al. 2021), is lower compared to the one associated with the KiDS-DR3 sample, namely 8.53 arcmin⁻² (Hildebrandt et al. 2017), due to more stringent redshift and SOM-gold selections (Wright et al. 2020; Giblin et al. 2021), the larger survey area in KiDS-1000 implies a better precision on the mass scaling relation. We also remark that the analysis carried out in this work is not fully comparable with that performed by B19. The modelling procedure is different, as we derive the log λ^* − log M₂₀₀ relation directly from the stacked profiles, while B19 first derived the average masses of the stacks and then modelled their relation with cluster observables. The weak-lensing modelling carried out in this work has the advantage of allowing for a proper marginalisation of the cosmological posteriors over the local cluster properties. Moreover, the set of modelling parameters differs from that assumed by B19, since, for example, we marginalise the log λ^* − log M₂₀₀ relation over the halo truncation factor, F_t, and σ_intr. A thorough comparison between AMICO KiDS-1000 and KiDS-DR3 mass estimates is detailed in Appendix C. In addition, in Appendix D we discuss the mass calibration results in the case of alternative spherical overdensity definitions.

We find that F_t and the fraction of miscentred clusters, f_off, are not constrained, while we obtain ${\sigma }_{\mathrm{off}}={0.24}_{-0.11}^{+0.11}$h⁻¹ Mpc. This result on the miscentring scale agrees with literature results from both simulations and observations (Saro et al. 2015; Zhang et al. 2019; Yan et al. 2020; Sommer et al. 2024; Seppi et al. 2023). Furthermore, we verified that the constraint $log{c}_0={0.42}_{-0.13}^{+0.13}$ yields a log c₂₀₀ − log M₂₀₀ relation which is consistent within 1σ with the one by Duffy et al. (2008). Compared to the results by Dutton & Macciò (2014), Meneghetti et al. (2014), Child et al. (2018), Diemer & Joyce (2019), and Ishiyama et al. (2021), our log c₂₀₀ values are lower but remain consistent within 2σ.

6.3. Robustness of the results

Throughout this paper, we based our assumptions on simulations and prior observational constraints for parameters that cannot be fully determined by our data, as detailed in Sect. 5.4. In addition, the choice of sample selections is motivated solely by the goodness of the fit. For example, we excluded clusters at redshift z > 0.8 because, if included, they would degrade the quality of our fit. In this section we test how much these assumptions affect the final results. In the tests presented in this section, we did not include the SSC contribution, as it is not expected to affect our conclusions. Therefore, the comparisons below are made against the posteriors obtained from the baseline analysis where SSC is also excluded.

Firstly, to assess the impact of our assumptions on c₂₀₀, we set the mass evolution parameter of the concentration, namely c_M in Eq. (38), as free to vary. By assuming the large flat prior c_M ∈ [ − 1.5, 1.5], we obtain $c_M=0.{07}_{-0.19}^{+0.23}$, that is well consistent with the value assumed for the baseline analysis, namely c_M = −0.084. The statistical uncertainty on log c₀ is increased by 20% this case, and the constraint on log c₀ is consistent within 1σ with the one reported in Table 2. The log λ^* − log M₂₀₀ relation and the cosmological posteriors remain unchanged. Furthermore, as discussed in Sect. 5.3, we neglected the effect of the intrinsic scatter in the log c₂₀₀ − log M₂₀₀ relation. If assumed to be log-normal, the intrinsic scatter on the natural logarithm of the concentration given a mass amounts to about 35% (see e.g. Duffy et al. 2008; Bhattacharya et al. 2013; Diemer & Kravtsov 2015; Umetsu 2020). Applying this relative uncertainty to the median posterior value of log c₀ reported in Table 2, we obtain an intrinsic scatter of σ_intr^{log c₀} = 0.06, that is less than half the statistical uncertainty on log c₀, namely σ_stat^{log c₀} = 0.13. Treating σ_intr^{log c₀} and σ_stat^{log c₀} as uncorrelated, the total uncertainty on log c₀ amounts to σ_tot^{log c₀} = 0.14, which is very close to σ_stat^{log c₀}. Thus, we do not expect σ_intr^{log c₀} to have a significant role in our analysis.