# Weak lensing mass-richness relation of redMaPPer clusters in LSST DESC DC2 simulations

Constantin Payerne<fnref n="1" />,²,*, Zhuowen Zhang³, Michel Aguena⁴,⁵, Céline Combet², Thibault Guillemin⁶, Marina Ricci⁴, Nathan Amouroux⁶, Camille Avestruz⁷,⁸, Eduardo J. Barroso⁶, Arya Farahi⁹,¹⁰, Eve Kovacs¹¹, Calum Murray⁴,¹², Markus M. Rau¹³,¹⁴, Eli S. Rykoff¹⁵,¹⁶, Samuel J. Schmidt¹⁷, and the LSST Dark Energy Science Collaboration

¹ Université Paris-Saclay, CEA, IRFU, 91191, Gif-sur-Yvette, France

² Université Grenoble Alpes, CNRS, IN2P3, LPSC, 38000 Grenoble, France

³ Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA

⁴ Université Paris Cité, CNRS, IN2P3, APC, 75013 Paris, France

⁵ Italian National Institute of AstroPhysics, Osservatorio Astronomico di Trieste, Italy

⁶ Université de Savoie, CNRS, IN2P3, LAPP, Annecy-le-Vieux, France

⁷ Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA

⁸ Leinweber Center of Theoretical Physics, University of Michigan, Ann Arbor, MI 48109, USA

⁹ Department of Statistics and Data Sciences, The University of Texas at Austin, TX 78712, USA

¹⁰ The NSF-Simons AI Institute for Cosmic Origins, University of Texas at Austin, Austin, TX 78712, USA

¹<fnref n="1" /> HEP Division, Argonne National Laboratory, 9700 S. Cass Ave., Lemont, IL 60439, USA

¹² Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-sur-Yvette, France

¹³ School of Mathematics, Statistics and Physics, Newcastle University, Newcastle upon Tyne, NE17RU, UK

¹⁴ High Energy Physics Division, Argonne National Laboratory, Lemont, IL 60439, USA

¹⁵ Kavli Institute for Particle Astrophysics & Cosmology, PO Box 2450, Stanford University, Stanford, CA 94305, USA

¹⁶ SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA

¹⁷ Department of Physics and Astronomy, University of California, One Shields Avenue, Davis, CA 95616, USA

Received 11 February 2025 / Accepted 8 June 2025

## ABSTRACT

**Context.** Cluster scaling relations are key ingredients in cluster abundance-based cosmological studies. In optical cluster cosmology, where clusters are detected through their richness, cluster-weak gravitational lensing has proven to be a powerful tool to constrain the cluster mass-richness relation. This work is conducted as part of the Dark Energy Science Collaboration (DESC), which aims to analyze the Legacy Survey of Space and Time (LSST) of the Vera C. Rubin Observatory, starting in 2026.

**Aims.** Cluster properties inferred from weak lensing, such as mass, suffer from several sources of bias. In this paper, we aim to test the impact of modeling choices and observational systematics in cluster lensing on the inference of the mass-richness relation.

**Methods.** We constrained the mass-richness relation of 3600 clusters detected by the redMaPPer algorithm in the cosmoDC2 extra-galactic mock catalog of the LSST DESC DC2 simulation, covering 440 deg², using number count measurements and either stacked weak lensing profiles or mean cluster masses in several intervals of richness ($20 \le \lambda \le 200$) and redshift ($0.2 \le z \le 1$).

**Results.** We provide the first constraints on the redMaPPer cluster mass-richness relation detected in cosmoDC2. We find that for an LSST-like source galaxy density, our constraints are robust to changes in the concentration-mass relation, as well as the dark matter density profile modeling choices, when source redshifts and shapes are perfectly known. We find that photometric redshift uncertainties can introduce bias at the $1\sigma$ level, which could be mitigated by an overall correction factor fitted jointly with the scaling parameters. We find that including positive shear-richness covariance in the fit shifts the results by up to $0.5\sigma$. Our constraints also offer a fair comparison to a fiducial mass-richness relation, obtained from matching cosmoDC2 halo masses to redMaPPer-detected cluster richness results.

**Key words.** gravitational lensing: weak – methods: statistical – galaxies: clusters: general

## 1. Introduction

Galaxy clusters have been essential in the construction of the standard model of cosmology, providing some of the first evidence of dark matter (Zwicky 1937) through the motions of galaxies within clusters. Their spatial distribution has also provided evidence for the primordial origin of density fluctuations (Kaiser 1984). They originate from the gravitational collapse of large matter overdensities that are decoupled from the cosmic expansion and form galaxy clusters. They represent the most

massive structures in the Universe held together by gravity. Their formation history, as well as their spatial and mass distributions, are strongly dependent on the nature of gravity, the growth rate of large-scale structures, and the Universe's expansion history (Bartlett 1997; Allen et al. 2011; Kravtsov & Borgani 2012).

Over the years, cluster abundance has proven to provide competitive constraints<fnref n="1" />, complementary to other large-scale structure and geometrical probes, as well as with

<footnote n="1">¹ Especially on the total amount of matter in the Universe $\Omega_m$ and the amplitude of matter fluctuations $\sigma_8$.</footnote>

<footnote>* Corresponding author: constantin.payerne@gmail.com</footnote>
the analysis of the Cosmic Microwave Background, for instance, using X-ray clusters detected by the ROSAT<fnref n="2" /> All-Sky Survey (Mantz et al. 2015), XMM-Newton (Pacaud et al. 2018), or eROSITA<fnref n="3" /> Ghirardini et al. 2024); clusters detected through their Sunyaev-Zeldovich effect by the Planck satellite (Planck Collaboration XX 2014; Planck Collaboration XXIV 2016), South Pole Telescope (SPT; Bocquet et al. 2024; Vogt et al. 2024), or Atacama Cosmology Telescope (ACT; Hasselfield et al. 2013). Other studies have also used clusters detected through their galaxy member populations via, for instance, Dark Energy Survey (DES Collaboration 2020, 2025), Kilo-Degree Survey (KiDS; Lesci et al. 2022), Sloan Digital Sky Survey (SDSS; Fumagalli et al. 2024), or with shear-selected clusters from the Subaru Hyper Suprime-Cam (HSC; Chiu et al. 2024). Table 2 in Payerne et al. (2024) recaps the cluster abundance-based cosmological analyses prior to June 2024.

The constraining power of cluster number counts is currently limited by our understanding of the cluster scaling relations (e.g., Pratt et al. 2019; Costanzi et al. 2019; DES Collaboration 2020), which denote the statistical relationship between cluster observables and their total masses. At optical wavelengths, clusters are often detected through the color and brightness of their member galaxies (Rykoff et al. 2014). Still, clusters of galaxies can also be identified as high signal-to-noise ratio (S/N) peaks on weak-lensing aperture mass maps (Hetterscheidt et al. 2005; Chen et al. 2025).

Weak gravitational lensing has become a robust tool for constraining cluster masses (e.g., McClintock et al. 2019; Umetsu 2020; Murray et al. 2022; Mistele & Durakovic 2024; Grandis et al. 2024) through the coherent distortion of the shapes of background galaxies, caused by the bending of the light path due to the cluster's gravitational field. Cluster abundance and cluster weak lensing information (using either masses or profiles directly) are usually combined, as they display different degeneracies on scaling relation parameters and cosmology, which allows for tighter constraints to be placed on both cosmological parameters and cluster scaling relations (Mantz et al. 2015; Murata et al. 2019; Mulroy et al. 2019; DES Collaboration 2020; Lesci et al. 2022; Sunayama et al. 2024). However, the mapping between the measurements of weak gravitational lensing and cluster masses is not well understood, as it is affected by several sources of statistical noise and systematic uncertainties (e.g., Becker & Kravtsov 2011; Köhlinger et al. 2015; Grandis et al. 2021).

First, weak lensing cluster mass reconstruction is statistically limited by the finite galaxy number density. Ongoing surveys such as DES, KiDS, HSC, and the Euclid mission (Laureijs et al. 2011), as well as future wide optical surveys such as Legacy Survey of Space and Time<fnref n="4" /> (LSST; Abell et al. 2009) and Nancy Grace Roman Space Telescope (Spergel et al. 2015), are equipped to provide a large influx of data to reduce statistical uncertainties in lensing measurements to an unprecedented level. Other sources of statistical noise affect the measurement of the cluster lensing signal, including intrinsic variability of cluster morphology, line-of-sight correlated structures, and galaxy shape intrinsic variations.

Second, several systematic biases may affect cluster mass inference. Some are related to “theoretical uncertainties”, such as unknowns within the modeling of the dark matter density in the cluster field (Becker & Kravtsov 2011; Lee et al. 2018), or the modeling of the presence of baryons (Cromer et al. 2022). Other biases are more directly related to the observations themselves and may arise from shear calibration (Hernández-Martín et al. 2020), calibration of the photometric redshift distribution of the background galaxy sample (Wright et al. 2020), contamination of the source galaxy sample by foreground and cluster member galaxies (Varga et al. 2019), and miscentering, which is the possible offset between the cluster’s detected center and the center of its gravitational potential (Becker & Kravtsov 2011; Lee et al. 2018; Zhang et al. 2019; Sommer et al. 2022), as well as any selection and projection effects, when galaxies physically unassociated with clusters are counted as members by cluster finders (Myles et al. 2021; Lee et al. 2024; Myles et al. 2025). All of these sources of scatter and bias must be carefully controlled, given the unprecedented volume and quality of lensing data provided by the upcoming large-footprint lensing surveys.

The Vera Rubin Observatory in Chile will conduct LSST, a 10-year wide-field imaging survey that will cover 18 000 deg<fnref n="2" /> of the southern sky, starting at the end of 2025 (Abell et al. 2009). The LSST aims to measure the shapes and redshifts of a few billion galaxies up to z ~ 3, enabling the measurement of number counts and weak gravitational lensing for ~100 000 galaxy clusters up to z ~ 1 (Abell et al. 2009). In this context, the Dark Energy Science Collaboration (DESC) from LSST Dark Energy Science Collaboration (2012)<fnref n="5" /> is poised to deliver a cosmological analysis of upcoming LSST data.

In this work, we make use of the LSST DESC Data Challenge 2 (DC2; Korytov et al. 2019; Abolfathi et al. 2021) simulated dataset. We aim to infer the mass-richness relation of galaxy clusters detected by the redMaPPer cluster finder algorithm (Rykoff et al. 2014) in DC2, using a combination of cluster weak lensing and abundance.

This paper is organized as follows. In Sect. 2, we introduce the theoretical background of this work, describing how cluster number counts, stacked lensing profile, and mean mass of a cluster ensemble are intrinsically linked to the cluster mass-richness relation that our study is aimed at constraining. The DESC cosmoDC2 simulated dataset used in this work is presented in Sect. 3, describing in detail the galaxy catalog (used to estimate the weak lensing signal) and the DC2 redMaPPer cluster catalog; the fiducial mass-richness relation to which the results will be compared to is also given in that section. In Sect. 4, we present the methodology for building the data vectors (counts, stacked lensing profiles, or mean lensing masses) and the likelihoods that combine them with the models of Sect. 2 to infer the cluster mass-richness relation. In Sect. 5, we present the results of the inference procedure and discuss the stability of the derived redMaPPer-detected cluster mass-richness relation in cosmoDC2 when considering a variety of systematic effects. Finally, we summarize our findings and present our conclusions in Sect. 6.

## 2. Theoretical background

Since our aim is to constrain the mass-richness relation of redMaPPer-selected clusters in the LSST DESC DC2 simulations, we first present the general theoretical framework for mod-

<footnote n="2">² Short for Röntgensatellit (in German X-rays are called Röntgenstrahlen).</footnote>

<footnote n="3">³ Short for extended ROentgen Survey with an Imaging Telescope Array.</footnote>

<footnote n="4">⁴ http://www.lsst.org</footnote>

<footnote n="5">⁵ http://lsstdesc.org</footnote>
eling cluster number counts (Sect. 2.1) and cluster weak gravitational lensing (Sect. 2.2). We then introduce the parametrization of the mass-richness relation considered in this work (see Sect. 2.3).

## 2.1. Cluster number count

For clusters detected by their observed richness, the redshift-richness cluster number density is given by

$$ \frac{d^2 N(\lambda, z)}{d\Omega dz d\lambda} = \int_{m_{\min}}^{m_{\max}} dm \frac{dn(m, z)}{dm} \frac{d^2 V(z)}{dz d\Omega} \Phi(\lambda, m, z) P(\lambda|m, z), \quad (1) $$

where $\Omega$ is the survey sky area$^6$, $P(\lambda|m, z)$ corresponds to the mass-richness relation (see below), $dn(m, z)/dm$ is the halo mass function (which predicts the comoving number density of dark matter halos per mass interval), $d^2V(z)/dzd\Omega$ is the differential comoving volume, and $\Phi(\lambda, m, z)$ is the cluster selection function. The differential comoving volume is given by

$$ \frac{d^2 V(z)}{dz d\Omega} = d_H \frac{D_C^2(z)}{H(z)/H_0}, \quad (2) $$

where $H(z)$ is the Hubble parameter at redshift $z$ (with the Hubble constant given by $H_0 = H(z=0)$), $d_H = c/H_0$ is the Hubble distance, and $D_C$ is the radial comoving distance (Hogg 1999). The masses $m_{\min}$ and $m_{\max}$ correspond to the minimum and maximum mass accessible in the survey by the cluster finder, which are generally set by the survey strategy.

The selection function $\Phi$ denotes our ability to detect clusters. It is associated with the performance of the cluster finder algorithm and the survey strategy, and its calibration plays a crucial role in cluster-based cosmological analyses. The selection function of a cluster finder can be assessed by (i) using end-to-end simulated data, geometrically matching the observed cluster catalog to the underlying dark matter halo population (Euclid Collaboration 2019; Lesci et al. 2022; Bulbul et al. 2024); (ii) injecting mock clusters into the real dataset (Rykoff et al. 2014, 2016; Planck Collaboration XXVII 2016); or (iii) cross-validating the cluster finder's performance by comparing datasets at different wavelengths (Sadibekova et al. 2014; Saro et al. 2015). We follow the prescription presented in Aguena & Lima (2018), where the selection function is given by

$$ \Phi(\lambda, m, z) = \frac{c(m, z)}{p(\lambda, z)}. \quad (3) $$

First, the completeness, $c(m, z)$, gives the fraction of true underlying dark matter halos detected by the cluster finder algorithm. If the completeness is less than 1, then a portion of the underlying dark matter halo population is systematically missing from our dataset (Mantz 2019). Second, the purity $p(\lambda, z)$ denotes the fraction of “spurious” detections in the cluster catalog (false positives or misidentified structures along the line of sight), where $\lambda$ and $z$ are the cluster richness and redshift$^7$.

With all these elements, the predicted cluster number count, $N_{ij}$, in the $i$-th redshift and $j$-th richness bin is finally given by

$$ N_{ij} = \Omega \int_{z_i}^{z_{i+1}} dz \int_{\lambda_j}^{\lambda_{j+1}} d\lambda \frac{d^2 N(\lambda, z)}{d\Omega dz d\lambda}. \quad (4) $$

## 2.2. Cluster weak gravitational lensing

The observed ellipticity $\epsilon^{\text{obs}}$ of a lensed source galaxy is related to its intrinsic (unlensed) shape $\epsilon^{\text{int}}$ (Schneider et al. 1992) given by

$$ \epsilon^{\text{obs}} = \frac{\epsilon^{\text{int}} + g}{1 + g^* \epsilon^{\text{int}}}, \quad (5) $$

where $g = \gamma/(1 - \kappa)$ is the reduced shear, with $\gamma$ being the shear and $\kappa$ the convergence. In the weak lensing regime, we assume that $\kappa \ll 1$, so the first-order Taylor expansion in $\gamma$ of the observed ellipticity gives $\epsilon^{\text{obs}} \approx \gamma + \epsilon^{\text{int}}$. Taking the average of this expression and assuming that $\langle \epsilon^{\text{int}} \rangle = 0$ in the absence of large-scale galaxy intrinsic alignments, we find that the cluster's local shear can be estimated statistically through $\langle \epsilon^{\text{obs}} \rangle \approx \gamma$. In practice, we decompose galaxy shapes into tangential and cross components, defined respectively by

$$ \epsilon_+ + i\epsilon_- = -\epsilon^{\text{obs}} \exp\{-i2\phi\}, \quad (6) $$

where $\epsilon^{\text{obs}}$ is the observed galaxy's ellipticity and $\phi$ is the polar angle of the source galaxy relative to the center of the galaxy cluster$^8$. For any mass distribution, we find that $\langle \epsilon_+ \rangle_C = \gamma_+$, where $\gamma_+$ is the tangential shear, and $\langle \epsilon_- \rangle_C = \gamma_-$, where $\langle \cdot \rangle_C$ denotes the average along a closed (circular) loop (Bernstein & Nakajima 2009; Umetsu 2020).

The excess surface density, $\Delta\Sigma(R)$, is commonly used as the cluster weak-lensing shear estimator (Mandelbaum et al. 2005; Murata et al. 2019; McClintock et al. 2019). At a given projected radius, $R$, from the cluster center, the excess surface density is related to the local tangential shear $\gamma_+(R, z_l, z_s)$; where $z_l$ is the redshift of the cluster (denoted as “lens”) and $z_s$ is the redshift of the source galaxy via

$$ \Delta\Sigma(R) = \Sigma_{\text{crit}}(z_s, z_l) \gamma_+(R, z_l, z_s), \quad (7) $$

where the critical surface mass density $\Sigma_{\text{crit}}(z_s, z_l)$ is a geometrical factor given by

$$ \Sigma_{\text{crit}}(z_s, z_l) = \frac{c^2}{4\pi G} \frac{D_A(z_s)}{D_A(z_l)D_A(z_s, z_l)}, \quad (8) $$

where $D_A(z_l)$, $D_A(z_s)$, and $D_A(z_l, z_s)$ are, respectively, the physical angular diameter distances to the lens, to the source, and between the lens and the source (Hogg 1999).

To model the stacked lensing signal, $\Delta\Sigma(R)$ can be related to the surface mass density of the cluster $\Sigma(R)$. For a perfectly centered cluster, it is given by

$$ \Sigma(R) = \int_{-\infty}^{+\infty} dy \rho \left( \sqrt{y^2 + R^2} \right), \quad (9) $$

where $\rho$ is the three-dimensional dark matter density. Then, $\Delta\Sigma$ is given by

$$ \Delta\Sigma(R) = \Sigma(<R) - \Sigma(R), \quad (10) $$

where

$$ \Sigma(<R) = \frac{2}{R^2} \int_{0}^{R} dR' R' \Sigma(R'). \quad (11) $$

<footnote>$^6$ Here, we assume the survey depth is uniform across the sky.</footnote>

<footnote>$^7$ Purity is generally defined as $p(\lambda, z_{\text{obs}})$ instead of $p(\lambda, z)$, where $z_{\text{obs}}$ is the cluster finder-inferred redshift. For simplicity, we consider here that observed and true redshifts (i.e., the redshift of the underlying dark matter halo) are equal.</footnote>

<footnote>$^8$ The cluster center is typically identified from observational counterparts, such as the position of the brightest central galaxy. However, there may be an offset, called the miscentering radius, between the identified cluster center and the center of its gravitational potential.</footnote>
The total matter content around clusters is parametrized by the three-dimensional density $\rho(r)$ and, by extension, $\Sigma(R)$ through Eq. (9).

The matter density around clusters has two contributions. The first one is the one-halo term, denoting the intrinsic cluster matter density. We generally assume that the one-halo term is well described by the Navarro-Frank-White (NFW; Navarro et al. 1997) profile<fnref n="9" />, set by the spherical overdensity mass,

$$M_{200c} = \frac{4\pi r_{200c}^3}{3} 200\rho_c(z), \qquad (12)$$

where $r_{200c}$ is the radius of the sphere that contains an average matter density 200 times higher than the critical density $\rho_c(z)$. Such profiles also depend on a concentration parameter $c_{200c}$, which quantifies the level of concentration of mass in the innermost regions of the cluster (we discuss in detail in Sect. 5.2.1 the different concentration-mass relations used in the literature and their impact on the inference of the cluster mass). The one-halo term dominates the lensing signal for $R < 3.5$ Mpc, while the two-halo term – associated with the contribution to the matter density field from surrounding, correlated halos – becomes increasingly more important at larger scales (Oguri & Takada 2011).

The collapsed dark matter halos are not expected to be spherical, due to the non-spherical initial density peaks from which they form and also due to their complex individual accretion history in the cosmic web (Sheth et al. 2001). The underlying dark matter halos have been shown to have complex triaxial structures in simulations (Jing & Suto 2002; Schneider et al. 2012; Despali et al. 2014), as well as from the observation of their member galaxies (Binggeli 1982) or their weak lensing shear field (Oguri et al. 2003, 2010). In this analysis, we consider stacks of individual cluster lensing profiles. Not only does stacking profiles increase the S/N, which is particularly beneficial for low-mass clusters (where the gravitational lensing signal is weaker), but the stacking method also averages out the intrinsic triaxiality of individual halos and substructures (Corless & King 2009), effectively recovering a spherical profile.

As for the cluster count, accounting for the selection function of the cluster finder is essential to predict the average cluster lensing profile for an ensemble of clusters. In this work, we neglect the contributions from “spurious”<fnref n="10" /> redMaPPer detections (encoded in the purity) to the measured cluster lensing signal, as we consider it to be negligible compared to the cluster lensing induced by properly identified halos – encoded in the completeness. However, we recall that it is a strong assumption, since any redMaPPer-identified structure would inevitably induce positive lensing signal; however, this has not been explored extensively in the literature and this is not the purpose of this work, so we keep using the aforementioned assumption<fnref n="11" />. Then, the predicted stacked excess surface density profile

is given by

$$\begin{gathered} \Delta\Sigma_{ij}(R) = \frac{1}{N_{ij}} \int_{z_i}^{z_{i+1}} dz \int_{\lambda_j}^{\lambda_{j+1}} d\lambda \int_{m_{\min}}^{+\infty} dm \times \\ \frac{d^2 N(m, z)}{dz dm} c(m, z) P(\lambda|m, z) \Delta\Sigma(R|m, z), \end{gathered} \qquad (13)$$

where we use the completeness $c(m, z)$, namely, the fraction of true detected halos relative to the underlying halo population.

Alternatively, instead of using the stacked profile directly, we can use the mean mass of the clusters within the ij-th redshift-richness bin – derived from the corresponding excess surface density profile – to investigate the cluster scaling relation (we provide more details in Sect. 4). For this purpose, using the same arguments as above, the mean mass can be modeled by

$$\begin{gathered} M_{ij} = \frac{1}{N_{ij}} \int_{z_i}^{z_{i+1}} dz \int_{\lambda_j}^{\lambda_{j+1}} d\lambda \int_{m_{\min}}^{+\infty} dm \times \\ \frac{d^2 N(m, z)}{dz dm} c(m, z) P(\lambda|m, z) m. \end{gathered} \qquad (14)$$

In this analysis, either equation, Equations (13) or (14), will be used to model the cluster lensing information, in combination with cluster counts, to constrain the mass-richness relation (see the next section).

## 2.3. Cluster mass-richness relation

In this paper, we aim to infer the mass-richness relation of galaxy clusters detected by the redMaPPer (Rykoff et al. 2014) cluster finder in the LSST DESC DC2 simulations using a combination of cluster lensing and cluster counts. The observed cluster richness $\lambda$ traces the number of member galaxies inside a galaxy cluster. Due to (i) the complex cluster formation history and baryon physics and (ii) observational effects in richness measurement (see, e.g., Costanzi et al. 2018), the observed cluster richness, $\lambda$, is a statistical variable for a halo with a fixed mass, $m$. In this paper, we consider the log-normal scaling relation $P(\lambda|m, z)$ (see, e.g., Mantz et al. 2008; Evrard et al. 2014; Saro et al. 2015; Farahi et al. 2018; Murata et al. 2019; Anbajagane et al. 2020) given by

$$P(\ln \lambda | m, z) = \frac{1}{\sqrt{2\pi}\sigma_{\ln \lambda | m, z}} \exp \left\{ -\frac{[\ln \lambda - \langle \ln \lambda | m, z \rangle]^2}{2\sigma_{\ln \lambda | m, z}^2} \right\}, \qquad (15)$$

where

$$\langle \ln \lambda | m, z \rangle = \ln \lambda_0 + \mu_z \ln \left( \frac{1+z}{1+z_0} \right) + \mu_m \log_{10} \left( \frac{m}{m_0} \right), \qquad (16)$$

and

$$\sigma_{\ln \lambda | m, z} = \sigma_{\ln \lambda_0} + \sigma_z \ln \left( \frac{1+z}{1+z_0} \right) + \sigma_m \log_{10} \left( \frac{m}{m_0} \right). \qquad (17)$$

In the above equation, $\sigma_{\ln \lambda | m, z}$ is the total error, namely, including contributions from the richness measurement errors and the intrinsic sources of scatter (see, e.g., Murata et al. 2019). It is possible to split the two contributions by adding a Poisson variance term (Zhang et al. 2023) in Eq. (17) (see Appendix A). Let us note that the forward modeling $P(\lambda|m, z)$ is different from the backward formalism $P(m|\lambda, z)$, which has been used in several other works (see, e.g., Baxter et al. 2016; Melchior et al. 2017; Simet et al. 2017; Jimeno et al. 2018; McClintock et al. 2019).

<footnote n="9">⁹ The NFW profile, as well as other profiles used in the literature, are fitting functions calibrated on N-body simulations.</footnote>

<footnote n="10">¹⁰ Such as low-mass galaxy halos, misidentified structures along the line of sight.</footnote>

<footnote n="11">¹¹ In the latter, we measure the cluster lensing profile for $\lambda > 20$, such as the redMaPPer purity is ~99%.</footnote>
From the above, the cluster scaling relation has 6 free param-
eters: ln $\Lambda_0$, $\mu_z$, and $\mu_m$ denote respectively the offset richness,
the redshift dependence, and the mass dependence of the mean
cluster scaling relation, where the offset variance of the relation
is given by $\sigma_{ln\Lambda_0}$, and $\sigma_z$ and $\sigma_m$ denote the redshift and mass
dependencies of the variance.

3. Datasets

In this section we detail the datasets used in this paper. We first present the DESC DC2 (Sect. 3.1), from which the cosmoDC2 galaxy catalog and associated the photometric redshift catalogs (Sect. 3.2), and the catalog of clusters detected by redMaPPer (Rykoff et al. 2014) in cosmoDC2 (Sect. 3.3), are derived. Matching these clusters to dark matter halos in cosmoDC2 allows us to estimate a fiducial mass-richness relation (that will serve as reference for the remaining of the paper) and the selection function of the cosmoDC2 redMaPPer clusters, as detailed in Sect. 3.4.

3.1. *The simulated catalogs of the DESC Data Challenge 2*
(DC2)

DC2 is a vast simulated astronomical dataset covering ΩDC2 = 440 square degrees, designed to help develop and test the pipeline and analysis tools of the DESC for interpreting LSST data (see full details in Abolfathi et al. 2021).

The cosmoDC2 galaxy catalog detailed in the next Sect. 3.2 is built upon the OUTERRIM N-body (gravity-only) simulation (Heitmann et al. 2019). Each dark matter halo has been identified by a friend-of-friend (FoF) halo finder and assigned both a mass $M_{\text{FoF}}$ (the sum of the individual dark matter particles associated through the FoF algorithm) and a spherical overdensity mass $M_{200c}$ (obtained by fitting an NFW profile to the distribution of dark matter particles). Galaxies within halos were generated using the Galacticus semi-analytic model of galaxy formation (Benson 2012), and were placed into halos using GALSAMPLER (Hearin et al. 2020). The resulting galaxy catalog includes properties such as stellar mass, morphology, spectral energy distributions (SEDs), broadband filter magnitudes, and host halo information. Weak lensing shears and convergences at each galaxy position were computed using a ray-tracing algorithm applied to past light-cone particles from the simulation. We note that the performance of the ray-tracing procedure is degraded in the innermost regions of cluster fields (typically at $R \le 1$ Mpc); this resolution effect was observed in Korytov et al. (2019) and discussed in Kovacs et al. (2022) for galaxy-galaxy and cluster-galaxy lensing. In Sect. 4.1.2 we describe in more detail in how this issue is handled in our data analysis pipeline.

3.2. The cosmoDC2 extragalactic catalog and photometric redshift add-on catalogs

The cosmoDC2 extragalactic catalog (Korytov et al. 2019; Kovacs et al. 2022) contains ~2.26 billion galaxies, providing an inventory of ~550 properties per galaxy, including “true” galaxy attributes (e.g., true magnitudes in the six LSST bands, true redshift, true shapes, etc.) as well as ray-tracing quantities (shear and convergence). The catalog reaches a magnitude depth of 28 in the r-band and extends out to redshift $z \sim 3$. In this sense, it represents an idealized LSST dataset, as it includes only galaxies and does not account for observational contaminants such as dust extinction, stars, or instrumental systematics.

This work also uses two cosmoDC2 photometric redshift
add-on catalogs, produced by DESC using two representa-
tive, well-established approaches in the literature: FlexZBoost
(Izbicki & Lee 2017)<fnref n="12" /> and Bayesian photometric redshifts
(BPZ; Benítez 2011)<fnref n="13" />. FlexZBoost is an empirical, machine-
learning-based technique that has been shown to yield highly
accurate conditional photometric redshift density estimates in
prior studies (Schmidt et al. 2020). It uses spatially overlapping
data from LSST and accurate spectroscopic redshifts to estimate
the conditional redshift distribution p(z|m), where m corresponds
to LSST photometric magnitudes. BPZ is a template-fitting
method that formulates a likelihood of a galaxy’s observed col-
ors based on a set of redshifted SED models. This approach
incorporates informative priors, such as the expected redshift
distribution and galaxy type.

The DC2 simulation suite has two key limitations that affect the accuracy of photometric redshift estimates. First, the simulation does not model spectroscopic surveys or their galaxy selection functions. In cosmoDC2, FlexZBoost was trained on a complete subsample of galaxies down to $i < 25$. While this depth is representative of LSST, it extends beyond the range where real spectroscopic surveys are complete. As a result, the photometric redshift estimates obtained with FlexZBoost should be considered optimistic in terms of accuracy. Second, the color-redshift space of the simulated galaxies is not continuous, but discrete, due to the mock photometry construction process (Korytov et al. 2019)<fnref n="14" />. While template-fitting methods are generally less sensitive to incomplete calibration data, they depend more heavily on accurate SED models<fnref n="15" />. This artificial discreteness in the data leads to poor fits for template-based methods, which inherently assume that galaxy SEDs evolve smoothly with redshift<fnref n="16" />. Therefore, the accuracy of BPZ photometric redshifts in DC2 should be regarded as pessimistic.

We would expect realistic photometric redshift estimates to lie somewhere between the optimistic FlexZBoost and the pessimistic BPZ results from the DC2 runs. However, they are likely to be closer to FlexZBoost, given that the pessimism in BPZ primarily arises from nonphysical artifacts in the simulation.

3.3. The redMaPPer cluster catalog

As mentioned in the introduction, galaxy clusters can be identi-
fied through their member galaxies using galaxy catalogs. The
redMaPPer cluster finder (Rykoff et al. 2014) identifies clus-
ters via the presence of red-sequence galaxies, leveraging the
multiple optical bands available in the survey. redMaPPer has
already been widely used on SDSS (Abdullah et al. 2020) and
DES (DES Collaboration 2020), and has been run on cosmoDC2
using the six LSST-like bands *u*, *g*, *r*, *i*, *z*, and *y* (see Ricci
et al., in prep. for details). We note that since redMaPPer was

<footnote n="12">¹² https://github.com/rizbicki/FlexCoDE</footnote>

<footnote n="13">¹³ https://github.com/LSSTDESC/rail_bpz</footnote>

<footnote n="14">¹⁴ Specifically, a very low-resolution internal dust model was used, which introduced artificial “breaks” in some galaxy spectra, especially at high redshift.</footnote>

<footnote n="15">¹⁵ Template-fitting techniques are typically less affected by incomplete calibration data, but are more agnostic than empirical techniques, relying more on accurate SED modeling rather than the specificity of the calibration dataset.</footnote>

<footnote n="16">¹⁶ This is not a fundamental limitation of template-fitting methods, but stems from the strong prior assumptions imposed by empirical techniques when generalizing $p(z|m)$ from the calibration set to the LSST data.</footnote>
<figure><img src="image_1.png" /><figcaption>Fig. 1. Left: cosmoDC2 halo $M_{200c}$ masses (from the simulation) versus the redMaPPer cluster richnesses. The points are colored with the redMaPPer cluster redshift. The solid line is the best-fitted mean richness-mass relation in Eq. (16) at $z = 0.4$. The dashed lines represent the low mass and low richness cut used on the cosmoDC2-redMaPPer matched catalog for the fit of the fiducial scaling relation. Right: Histogram of richnesses in bins of mass, and best-fit probability distributions $P(\ln \lambda|M, z)$ from Eq. (15) (solid lines) for the different mass bins.</figcaption></figure>

applied directly to the cosmoDC2 galaxy catalog<fnref n="17" />, the resulting cluster catalog is free from atmospheric and instrumental systematics. The redMaPPer catalog provides galaxy cluster positions, redshifts, and richness values (computed as the sum of membership probabilities of galaxies around the cluster, see Rykoff et al. 2014), along with per-cluster membership galaxy catalogs, for 880 000 clusters within the ranges $5 < \lambda < 270$ and $0.1 < z < 1.15$.

## 3.4. The matched cosmoDC2 halo-redMaPPer cluster catalog

To infer the “true” cluster scaling relation, we can associate the redMaPPer-detected clusters with halos from the cosmoDC2 dark matter halo catalog using the DESC software C1EvaR<fnref n="18" />. The association is based on membership matching (see, e.g., Farahi et al. 2016), where we use the member galaxies<fnref n="19" /> of the detected redMaPPer clusters and those of the cosmoDC2 halos, respectively. Each halo is associated with all detected clusters with which it shares galaxies, and vice versa.

The one-to-one catalog is then obtained by selecting the associated system (redMaPPer cluster-cosmoDC2 halo pair) with the highest “membership fraction” (i.e., the overlap between the member galaxies of the cosmoDC2 halo and the redMaPPer cluster). We then retain only the reciprocal matches (i.e., two-way matches).

This matching procedure was performed using a cut of $M_{\text{FoF}} > 10^{13} M_{\odot}$ for the dark matter halo catalog and $\lambda > 5$ for the redMaPPer cluster catalog. We verified that applying a cut of $M_{200c} > 10^{13} M_{\odot}$ instead yields identical catalogs. The resulting redMaPPer-cosmoDC2 halo masses and richnesses are shown in Figure 1 (left panel). Then, each redMaPPer cluster

with richness, $\lambda_k$, was assigned a “true” spherical overdensity mass $M_k = M_{200c,k}$. The matched redMaPPer cluster-cosmoDC2 halo catalog objects in the mass-richness plane are shown in Figure 1 (left panel), color-coded according to their redshifts.

### 3.4.1. The fiducial mass-richness relation

From the cosmoDC2-redMaPPer matched catalog, we can infer a “fiducial” scaling relation. We consider the fiducial set to be the set of clusters with $\lambda > 10$ (horizontal dashed line in Figure 1, left panel), as well as $M > 4 \times 10^{13} M_{\odot}$ (vertical dashed line) and with $0.2 < z < 1$. The fiducial likelihood is given by

$$ \mathcal{L}_{\text{fid}} = \prod_{k=1}^{N_{\text{tot}}} P_{\lambda_k > 10}(\ln \lambda_k | M_k, z_k), \qquad (18) $$

where $P(\ln \lambda_k | M_k, z_k)$ is the mass-richness relation given in Eq. (15) for the k-th cluster-halo match (we consider the pivot redshift $z_0 = 0.5$ and pivot mass $m_0 = 10^{14.3} M_{\odot}$), and the subscript denotes truncated Gaussians with $\lambda_k > \lambda_{\min} = 10$. We use the emcee package (Foreman-Mackey et al. 2013), with flat priors inspired by Murata et al. (2019), as listed in Table 1 (second column). The best-fitted fiducial relation is shown in the third column of Table 1. These values can serve as a reference against which we compare the results obtained using cluster weak lensing and abundance.

The low-mass richness cuts applied to the matched catalog are deliberately conservative. These cuts are chosen to ensure that the completeness of the redMaPPer cluster catalog (see Sect. 3.4.2) does not influence the inference of the fiducial cluster scaling relation. In Appendix A we examine the impact of alternative cuts, motivated by the modeling of the mass-richness relation used in this work. While this topic warrants further investigation, for this study, we adopt the fiducial values listed in Table 1.

We show, over-plotted in black in the left panel of Figure 1, the mean scaling relation for the best-fit parameters at $z = 0.4$. The right panel of Figure 1 presents the normalized histograms of redMaPPer richness in three different mass bins (within the redshift range $z \in [0.2, 1]$). Over-plotted on these histograms is the best-fit probability density function, $P(\lambda|m_{\text{center}}, z_{\text{center}})$,

<footnote n="17">¹⁷ In addition to redMaPPer, other cluster detection methods have been applied to the cosmoDC2 data, including WaZP (Aguena et al. 2021a), AMICO (Bellagamba et al. 2018), and YOLO-CL (Grishin et al. 2023). The corresponding catalogs are currently available within DESC or are expected to be released soon.</footnote>

<footnote n="18">¹⁸ https://github.com/LSSTDESC/C1EvaR</footnote>

<footnote n="19">¹⁹ Specifically, we use the redMaPPer member galaxies for each cluster (i.e., the output of redMaPPer), and the cosmoDC2 galaxies flagged as “members” for each cosmoDC2 dark matter halo in the simulation.</footnote>
<table><caption>**Table 1.** Fiducial values of the parameters of the redMaPPer cluster mass-richness relation.</caption><thead><tr><th>Parameters</th><th>Priors</th><th>p ± Δp</th></tr></thead><tbody><tr><td>ln λ<sub>0</sub></td><td>U(0, +∞)</td><td>3.35 ± 0.01</td></tr><tr><td>μ<sub>z</sub></td><td>U(-2, 2)</td><td>0.06 ± 0.08</td></tr><tr><td>μ<sub>m</sub></td><td>U(0, +∞)</td><td>2.23 ± 0.05</td></tr><tr><td>σ<sub>ln λ<sub>0</sub></sub></td><td>U(0, +∞)</td><td>0.56 ± 0.01</td></tr><tr><td>σ<sub>z</sub></td><td>U(-2, 2)</td><td>-0.05 ± 0.04</td></tr><tr><td>σ<sub>m</sub></td><td>U(-2, 2)</td><td>0.10 ± 0.02</td></tr></tbody></table>

**Notes.** Second column: Prior we use in this analysis, in addition to $\sigma_{\ln \lambda_{m,z}} > 0$ (the dispersion is set to be positive). Third column: Posterior means for the six parameters of the scaling relation, by using individual dark matter halo masses and redMaPPer richnesses.

where $z_{\text{center}}$ and $m_{\text{center}}$ represent the center values of each redshift-mass bin. We observe that the predicted distribution matches the observed one quite well.

### 3.4.2. The redMaPPer selection function in DC2

To properly model cluster number counts, we must account for the selection function $\Phi$ of the redMaPPer algorithm (see, e.g., Euclid Collaboration 2019; Lesci et al. 2022). Defined in Eq. (3), the selection function reflects the fact that the cluster-finding algorithm may miss a fraction of true galaxy clusters (completeness) as well as detect “false” clusters that are not associated with underlying collapsed dark matter structures (purity).

From the matched catalog, we can disentangle the effects of purity and completeness within narrow mass-redshift-richness bins. The measured completeness $\tilde{c}_{ij}$ in the mass and redshift bins $[m_i, m_{i+1}]$ and $[z_j, z_{j+1}]$, and the measured purity $\tilde{p}_{kl}$ in the richness and observed redshift bins $[\lambda_k, \lambda_{k+1}]$ and $[z_{\text{obs},l}, z_{\text{obs},l+1}]$, are given by

$$ \tilde{c}_{ij} = \frac{N_{\text{match},ij}^{\text{halos}}}{N_{\text{tot},ij}^{\text{halos}}} \quad \text{and} \quad \tilde{p}_{kl} = \frac{N_{\text{match},kl}^{\text{halos}}}{N_{\text{obs},kl}^{\text{clusters}}}, \qquad (19) $$

where $N_{\text{tot},ij}^{\text{halos}}$ is the total number of halos in the ij mass-redshift bin, and $N_{\text{obs},kl}^{\text{clusters}}$ is the number of redMaPPer clusters in the kl richness-redshift bin. The quantity $N_{\text{match},ij}^{\text{halos}}$ (respectively $N_{\text{match},kl}^{\text{halos}}$) is the number of halos that have been matched to redMaPPer clusters in the ij mass-redshift bin (respectively in the kl richness-redshift bin). We model the observed completeness $\tilde{c}_{ij}$ and purity $\tilde{p}_{kl}$, allowing them to be included in the count/lensing formalism as detailed in Sect. 2. The performance of the DC2 redMaPPer cluster catalog is analyzed in detail in Ricci et al. (in prep.). We note that the redMaPPer cluster catalog is complete at the 80% level for $M_{200c} > 10^{14} M_{\odot}$, with a small redshift dependence, and is pure at the > 90% level for $\lambda > 12$.

## 4. Methodology to constrain the mass-richness relation

As shown in Sect. 2, the modeling of cluster number counts and/or cluster lensing (with stacked excess surface density profiles or mean cluster masses) in bins of redshift and richness depends on the mass-richness relation. To constrain the parameters of the latter using this observables requires building the

corresponding data vectors (Sect. 4.1) and defining the likelihood functions that link the data vectors, the covariances, and the models for each of these observables, as described in (Sect. 4.2). These likelihoods may then be used individually or in combination to constrain the parameters of the mass-richness relation. To implement the methodology described in this section, we rely on a variety of software tools that are briefly listed in Sect. 4.3.

### 4.1. Data vectors

We turn to the construction of the data vectors that will be used for the inference of the redMaPPer cluster mass-richness relation. For the stacking strategy we employ, we considered the redshift bin edges $z_i = 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1$ and the richness bin edges $\lambda_i = 20, 35, 70, 100, 200$. We use $\lambda > 20$ and $z > 0.2$, inspired by other redMaPPer-based cluster analyses (see, e.g., McClintock et al. 2019; DES Collaboration 2020). The $\lambda > 20$ cut ensures a high-purity cluster sample (Costanzi et al. 2019), close to 100% pure in our case, and restricting to $z > 0.2$ mimics the conservative cut used in DES cluster-based analyses (DES Collaboration 2020, 2025). This cut was implemented to prevent the degradation of redMaPPer performance at low redshifts$^{20}$, where the red-sequence galaxy population becomes harder to isolate due to the lack of u-band data in DES (which used g, r, i, and z bands). We note that robust detection of low-redshift clusters below $z = 0.2$ will be feasible with LSST, since u-band imaging will be available.

#### 4.1.1. Cluster number counts

We show in Figure 2 the measured count of redMaPPer clusters as a function of richness for different redshift bins. As expected, there are fewer clusters at higher richness than at lower richness.

#### 4.1.2. Stacked excess surface density profiles

Considering an ensemble of $N_l$ clusters (or “lenses”), where each cluster has $N_{ls}$ background source galaxies in the radial bin $[R - \Delta R/2, R + \Delta R/2]$, the maximum likelihood estimator of the stacked $\Delta\Sigma(R)$ profile is given by (Shirasaki & Takada 2018; Sheldon et al. 2004)

$$ \widehat{\Delta\Sigma}_+(R) = \frac{1}{\frac{N_l}{\sum_{l=1}^{N_l} \sum_{s=1}^{N_{ls}} w_{ls}}} \sum_{l=1}^{N_l} \sum_{s=1}^{N_{ls}} w_{ls} \widehat{\Sigma}_{\text{crit}}(z_s, z_l) \epsilon_+^{l,s}, \qquad (20) $$

where the sum runs over all lens-source pairs and only for sources located within the physical projected radius interval, $[R - \Delta R/2, R + \Delta R/2]$, from the lens $l$. Here, $\epsilon_+^{l,s}$ is the tangential ellipticity of the source galaxy, $s$, relative to the position of lens $l$, as given by Eq. (6). The quantity $\widehat{\Sigma}_{\text{crit}}(z_s, z_l)$ is the effective critical surface mass density of the lens-source system, averaged over the photometric redshift probability density function $p(z_s)$ of the galaxy with index $s$, such that $v_{\text{fill}}$

<footnote n="20">20 For $z < 0.2$ galaxies, the distinctive Balmer Break at 4000 Å falls in the u band.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 2. Measured count of redMaPPer cluster as a function of richness, for different richness bins. For each richness-redshift bin, the width of the shaded area represent the width of the richness bin, and the height correspond to the Poisson noise $\sqrt{N}$. The steeper slope of cluster counts at intermediate richness (between the second and third richness bin) for all redshift bins is primarily due to the choice of richness binning, which is not perfectly log-spaced. This steeper inner slope would be mitigated if a log-spaced richness binning were used.</figcaption></figure>

$$ \widehat{\Sigma}_{\text{crit}}(z_s, z_l)^{-1} = \int_{z_l}^{+\infty} dz_s p(z_s) \Sigma_{\text{crit}}(z_s, z_l)^{-1}. \quad (21) $$

The weights $w_{ls}$ maximize the S/N for this estimator (Sheldon et al. 2004) and can be written as the product $w_{ls} = w_{ls}^{\text{geo}} w_{ls}^{\text{shape}}$, such that

$$ w_{ls}^{\text{geo}} = \langle \widehat{\Sigma}_{\text{crit}}(z_s, z_l) \rangle^{-2}, \quad (22) $$

$$ w_{ls}^{\text{shape}} = \frac{1}{\sigma_{\text{rms}}^2(\epsilon_s^+) + \sigma_{\text{meas}}^2(\epsilon_s^+)}. \quad (23) $$

The quantity $\sigma_{\text{rms}}(\epsilon_s^+) = \sigma_{\text{SN}}$ is the intrinsic galaxy shape noise of the ellipticity's tangential component, arising from the fact that unlensed galaxies have an intrinsic distribution of ellipticities and position angles (see, e.g., Pranjal et al. 2023; Chang et al. 2013). In cosmoDC2, the intrinsic shape noise for each galaxy's ellipticity component is $\sigma_{\text{SN}} = 0.25^{21}$, which is comparable to what has been found in the literature (Chang et al. 2013). In contrast, $\sigma_{\text{meas}}(\epsilon_s^+)$ represents the uncertainty originating from galaxy ellipticity estimation on images (see, e.g., Hirata & Seljak 2003; Mandelbaum et al. 2005; Sheldon & Huff 2017). In the ideal case of cosmoDC2, where galaxy redshifts and shapes are perfectly measured, these weights reduce to $w_{ls}^{\text{geo}} = \Sigma_{\text{crit}}^{-2}(z_s, z_l)$ and $w_{ls}^{\text{shape}} = \sigma_{\text{rms}}^{-2}(\epsilon_s^+)$, where only the intrinsic variation of galaxy shapes remains.

The evaluation of the $\widehat{\Sigma}_{\text{+}}(R)$ data vector relies on a selection of background sources. The source selection is initially based on an r-band magnitude cut of $r < 28$ and an i-band magnitude cut of $i < 24.25$, aiming to reach a number density of galaxies $n_{\text{gal}} \approx 25 \text{ arcmin}^{-2}$, comparable to the number density expected in the context of LSST after 10 years of data (Chang et al. 2013). For each cluster in the redMaPPer catalog, we then extract the galaxy source catalog in a circular aperture of $R = 10$ Mpc, applying the

source selection $z_s > z_l + 0.2^{22}$, where $z_s$ is the true cosmoDC2 galaxy redshift and $z_l$ is the redMaPPer cluster redshift$^{23}$.

We constructed each source galaxy lensed ellipticity $\epsilon_{\text{obs}}$ from Eq. (5) using (i) its intrinsic shape $\epsilon_{\text{int}}$ and (ii) the shear and convergence values at the galaxy's location that are provided in the cosmoDC2 catalog. For each stack of $N_l$ clusters, we consider 10 log-spaced radial bins from 0.5 Mpc to 10 Mpc and estimate the stacked lensing profile using Eq. (20).

The corresponding excess surface density profiles are displayed in Figure 3 as a function of the distance to the cluster center, for different richness bins and different redshift bins. At fixed richness, we do not observe a particular trend with redshift at scales larger than 1 Mpc, as denoted by the vertical dashed lines.

However, for $R < 1$ Mpc, we observe that the stacked lensing profiles are progressively attenuated in the innermost regions with redshift. This attenuation has already been observed in DC2 galaxy-galaxy lensing (Korytov et al. 2019) and cluster-galaxy lensing (Kovacs et al. 2022) and is associated with the limited resolution of ray tracing used to compute the lensing shear and convergence at each galaxy position in the simulation. As a consequence, we chose to use only the $R > 1$ Mpc region for each stack in the analysis. This choice is valid across the full richness and redshift range. However, it also means that we cannot use the innermost regions with the largest S/N values. This intrinsic limitation will constrain the cluster-related forecasts of what could be achieved with the LSST data. We show in Figure 4 the S/N of some excess surface density profiles for a low (left panel) and a high (right panel) redshift bins. In the considered redshift

<footnote n="21">21 To infer the intrinsic cosmoDC2 galaxy shape noise, we compute $\sigma_{\text{SN}} = [\text{stddev}(\epsilon_1^{\text{int}}) + \text{stddev}(\epsilon_2^{\text{int}})]/2$, where $\epsilon_{[1,2]}^{\text{int}}$ are the two components of the intrinsic galaxy ellipticities, identified as ellipticity_{[1,2]}_true in the cosmoDC2 catalog.</footnote>

<footnote n="22">22 We adopted a redshift separation of $\Delta z = 0.2$ to mimic the common mitigation techniques of projection effects due to, e.g., contamination from foregrounds or cluster member galaxies. Additionally, since the lensing strength increases with the separation between the lens and the source, this threshold ensures that selected sources are more strongly lensed, thereby enhancing the S/N of the measurement. We do not explore less conservative cuts (e.g., $\Delta z = 0.1$ as in McClintock et al. 2019) in this work, as we expect them to have minimal impact on the results aside from increasing the per-cluster source galaxy number density, thereby enhancing the precision of the lensing measurement.</footnote>

<footnote n="23">23 In Sect. 5, we give other methodologies considered for source selection when using photometric redshifts.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Stacked excess surface density profile as a function of the distance to the cluster center, for different richness bins (colors) and different redshift bins (from top left to bottom right). The error bars are the diagonal elements of the bootstrap covariance matrices. The dashed black vertical line (resp. black filled vertical line) represents the $R > 1$ Mpc cut (resp. $R < 3.5$ Mpc cut).</figcaption></figure>

range, the S/N oscillates between 2 and 13 for the low-redshift bin, and between 0.5 and 8 for the high redshift bin. The profiles in the last richness bin have the lowest S/N, due to the low statistics in these bins (with ~10 times fewer clusters than in the first richness bin). Let us note that the stacking strategy allows us to significantly increase the S/N, as the S/N of individual cluster excess surface density profiles (not shown here) ranges from 0.1 to $3^{24}$.

### 4.1.3. Mean cluster masses

As mentioned earlier in this paper, instead of directly using the stacked profiles for weak-lensing information to constrain the mass-richness relation, we can use the corresponding mean cluster lensing masses as an alternative data vector. To do so, for each *ij*-th redshift-richness bin, we first need to infer the mean cluster mass $\log_{10} \hat{M}_{ij}$ and its corresponding error $\sigma(\log_{10} \hat{M}_{ij})$ from the corresponding stacked lensing profile $\Delta\hat{\Sigma}_{ij}$.

This is done by fitting, for each bin, a dark matter density profile$^{25}$ to the corresponding stacked lensing profile obtained in Sect. 4.1.2 using the Minuit$^{26}$ minimizer (James & Roos 1975), which is much faster than an MCMC. We ensured that the recovered mass and error from Minuit coincide at the <0.1σ level with the mean mass and errors inferred from an MCMC procedure.

We used a NFW profile (Navarro et al. 1997) with the concentration-mass relation of Duffy et al. (2008)$^{27}$. We show in the left panel of Figure 5 the inferred mean cluster masses as a function of richness for the different redshift bins after fitting the stacked profiles within the radial range$^{28}$ $R \in [1, 3.5]$ Mpc. The error bars on the *y*-axis correspond to the statistical error obtained from our fitting procedure of the cluster mass, while the *x*-axis error represents the dispersion of richness values in each bin.

For validation, we compared the recovered mean lensing mass to the true mean mass in each redshift-richness bin in the right panel of Figure 5, where the true masses are obtained by (i) matching the redMaPPer cluster catalog to dark matter halos (ii) computing the corresponding cosmoDC2 mass in each redshift-richness bin. We see a clear correlation between weak lensing mean masses and cosmoDC2 masses (we added the one-to-one line $x = y$ to allow for a visual comparison).

### 4.2. Likelihoods and priors

After building the data vectors, we turned our attention to inferring the parameters of the mass-richness relation defined in Sect. 2.3, namely $\theta = (\ln \lambda_0, \mu_z, \mu_m, \sigma_{\ln\lambda_0}, \sigma_z, \sigma_m)$. We adopted the pivot redshift $z_0 = 0.5$ and pivot mass $m_0 = 10^{14.3} M_{\odot}$ in the mass-richness relation in Eq. (16).

<footnote>$^{24}$ To go from stacked profile SNR to individual profile S/N, a good approximation is to rescale them by $1/\sqrt{N_l}$.</footnote>

<footnote>$^{25}$ Many parameterizations exist in the literature, such as the NFW profile (Navarro et al. 1997), as well as those by Einasto (1965) and Hernquist (1990).</footnote>

<footnote>$^{26}$ We use `minuit.migrad()` and `minuit.minos()` sequentially to (i) determine the best-fit mean cluster mass $\log_{10} M$ and (ii) compute the asymmetric error bars $\log_{10} M_{-\sigma_{-}}^{+\sigma_{+}}$. In Eq. (27), we adopt a Gaussian likelihood with symmetric errors, for which we use the mean error $(\sigma_{+} + \sigma_{-})/2$. However, asymmetric errors can be more accurately modeled using, for example, a split normal distribution.</footnote>

<footnote>$^{27}$ These modeling choices are referred to as the “baseline” analysis in the latter Section 5.1.</footnote>

<footnote>$^{28}$ We discuss the upper radial cut $R = 3.5$ Mpc in Sect. 4.2.2, which is set to use only the 1-halo term lensing signal.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 4. S/Ns of the stacked excess surface density profiles, for two different richness bins (left and right panel) and different richness bins (colors).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 5. Left: Mean cluster lensing mass as a function of mean richness within the stack, the different colors correspond to the different redshift bins. Right: Mean cluster lensing masses as a function of the DC2 matched masses from the cosmoDC2 dark matter halo catalog. The different colors correspond to the different redshift bins. The dashed line corresponds to y = x.</figcaption></figure>

### 4.2.1. Cluster count likelihood

We considered the cluster count likelihood (Hu & Kravtsov 2003) given by

$$ \mathcal{L}_N \propto |\Sigma_N|^{-1/2} \exp \left\{ -\frac{1}{2} \sum_{ijkl} (N - \tilde{N})_{ij} [\Sigma_N]^{-1} [ij, kl] (N - \tilde{N})_{kl} \right\} \quad (24) $$

In the above equation, the index *ij* (resp. *kl*) refers to the *i*-th redshift bin and the *j*-th richness bin (respectively the *k*-th redshift bin and the *l*-th richness bin). $\tilde{N}_{ij}$ (resp. $N_{ij}$) denotes the measured (resp. predicted) cluster count in the *ij*-th redshift-richness bin. For the count prediction, we adopted $m_{min} = 10^{12} M_{\odot}$ and $m_{max} = 10^{15.5} M_{\odot}$ as the minimum and maximum halo masses for the integration of the halo mass function of Despali et al. (2015) in Eq. (1). The best-fit selection function, $\Phi$, as described in Sect. 3.3, is used in Eq. (1). The quantity $\Sigma_N$ is the cluster count covariance matrix, which accounts for two contributions: the Poisson shot noise and the super-sample covariance (SSC). First, the Poisson noise affects the counting of clusters in uncorrelated bins (Poisson 1837). Second, the SSC denotes the contribution from the intrinsic fluctuations of the matter overdensity field, within and beyond the survey volume

(Hu & Kravtsov 2003). As a result, SSC introduces a covariance between the counts in different bins and increases their variance. SSC is particularly important in cluster abundance-based likelihoods, as it impacts the precision of the recovered cosmological parameters by approximately 20% for surveys such as the Vera C. Rubin Observatory's LSST or the Euclid mission (Fumagalli et al. 2021; Payerne et al. 2023, 2024), which are each expected to detect hundreds of thousands of clusters. The shot-noise and SSC contributions to the total covariance can, in principle, be estimated through resampling techniques directly on the data, such as jackknife resampling (Escoffier et al. 2016). However, due to the relatively small number of clusters detected in the DC2 galaxy catalog, jackknife estimates of the cluster count covariance matrix suffer from high noise. We therefore opt to use a theoretical prediction for the covariance instead, given by

$$ \Sigma_N[ij, kl] = N_{ij} \delta_{ik}^{K} \delta_{jl}^{K} + N_{ij} N_{kl} \langle b \rangle_{ij} \langle b \rangle_{kl} S_{ik}. \quad (25) $$

The first term in Eq. (25) represents the Poisson shot noise, with $\delta_{ik}^K$ denoting the Kronecker delta function. The second term corresponds to SSC, where $S_{ik} = \langle \delta_{m,i} \delta_{m,k} \rangle$ is the covariance of the smoothed matter overdensities $\delta_m$ in the *i*-th and
k-th redshift bins, respectively (Lacasa et al. 2018). The coupling between the matter overdensity field and the halo number density field is encoded by the halo bias $b(m, z)$ (Tinker et al. 2010), which depends on mass and redshift and is averaged over the ij-th (resp. kl-th) redshift-richness bin to yield $\langle b \rangle_{ij}$ (resp. $\langle b \rangle_{kl}$). Finally, $N_{ij}$ is the model prediction for the number of clusters in the ij-th redshift-richness bin, as defined in Eq. (4). This likelihood model, which incorporates both Poisson noise and SSC, has been adopted in previous cluster analyses, for instance, by DES Collaboration (2020), Costanzi et al. (2019), Sunayama et al. (2024).

### 4.2.2. Stacked lensing profile likelihood

We can constrain the mass-richness relation by using the stacked excess surface density profiles directly (Murata et al. 2019; Park et al. 2023; Sunayama et al. 2024). The likelihood for stacked excess surface density profiles, denoted $\mathcal{L}_{WLp}$, is given by

$$ \mathcal{L}_{\Delta\Sigma} \propto \exp\left(-\frac{1}{2} \sum_{ij} (\Delta\Sigma - \widehat{\Delta\Sigma})_{ij} \Sigma_{\Delta\Sigma}^{-1}[ij] (\Delta\Sigma - \widehat{\Delta\Sigma})_{ij}\right). \quad (26) $$

Here, $\widehat{\Delta\Sigma}_{ij}$ corresponds to the stacked excess surface density profile estimated in Eq. (20) for the ij-th redshift-richness bin, $\Delta\Sigma_{ij}$ is the corresponding theoretical prediction from Eq. (13), and $\Sigma_{\Delta\Sigma}[ij]$ denotes the associated covariance matrix. This covariance originates from a variety of phenomena, as discussed in Hoekstra (2003), Wu et al. (2019), Gruen et al. (2015), McClintock et al. (2019). A primary contribution is the intrinsic scatter in the shapes of background galaxies, combined with the limited number of clusters and source galaxies within each bin. This term, often referred to as shape noise, dominates the small-scale regime but can be reduced by increasing the sample size. Another significant contribution comes from uncorrelated large-scale structures along the line of sight, which introduce stochastic fluctuations in the lensing signal unrelated to the lensing cluster. In addition, correlated structures around the cluster – such as nearby halos and filaments – affect the lensing profile through variations in the two-halo term, thereby contributing an additive component to the covariance. At smaller scales, intrinsic variations in the properties of halos within each stack, including the scatter in concentration at fixed mass, halo ellipticity, and orientation, also contribute to the overall uncertainty. Finally, the spread in true halo masses among richness-selected clusters within a given bin introduces further variance in the measured lensing profiles.

As for cluster counts, the covariance of the cluster lensing profiles can be predicted analytically (see, e.g., Wu et al. 2019)<fnref n="29" /> or using semi-analytical approaches (see, e.g., McClintock et al. 2019). These methods require modeling a wide range of data properties, including the redshift distributions of the source and lens samples, the dispersion of intrinsic galaxy shapes, the dark matter halo density profile, and the level of scatter in concentration, miscentering, and morphology of clusters in the redMaPPer catalog. In our case, we choose a simpler, data-driven approach that estimates all these contributions simultaneously without relying on model assumptions, while also keeping computational costs relatively low. Specifically, we use bootstrap resampling (see, e.g., Simet et al. 2017; Parroni et al. 2017) with $N_{boot} = 400$

resampled cluster ensembles. We retain only the diagonal terms of the covariance matrices, as the off-diagonal elements are too noisy due to the limited cluster count statistics in each richness-redshift bin within the DC2 footprint. Wu et al. (2019) point out that neglecting off-diagonal terms at large scales ($R \gtrsim 5-6$ Mpc) can lead to an underestimation of parameter uncertainties inferred from weak-lensing profiles – such as cluster mass and concentration – especially when shape noise (i.e., uncorrelated sources of noise) is subdominant. However, since our analysis focuses exclusively on the $R < 3.5$ Mpc regime, this impact is expected to be minimal. Therefore, we assume diagonal covariance matrices throughout, as the off-diagonal terms remain extremely noisy (Phriksee et al. 2020).

### 4.2.3. Mean cluster mass likelihood

The cluster mass-richness relation can also be inferred using the information from mean mass estimates (McClintock et al. 2019; DES Collaboration 2020; Lesci et al. 2022). Then, the mean mass likelihood is given by

$$ \mathcal{L}_{M_{WL}} \propto \exp\left(-\frac{1}{2} \sum_{ij} \left( \frac{\log_{10} M_{ij} - \log_{10} \widehat{M}_{ij}}{\sigma(\log_{10} \widehat{M}_{ij})} \right)^2\right), \quad (27) $$

where $\log_{10} M_{ij}$ is the mass prediction that is given in Eq. (14) in the ij-th redshift-richness bin, and $\log_{10} \widehat{M}_{ij} \pm \sigma(\log_{10} \widehat{M}_{ij})$ is the measured mean mass.

It is important to highlight the differences between the one-step approach (described in Sect. 4.2.2) and the two-step approach (discussed in this section). The one-step approach offers greater flexibility. Specifically, it allows for more comprehensive modeling of various systematic effects that influence the lensing signal, including cluster miscentering, projection effects from both correlated and uncorrelated large-scale structures, contamination by cluster member galaxies, and deviations from spherical symmetry in the cluster mass distribution compared to the two-step approach. However, the one-step approach presents a significant computational challenge, as it necessitates marginalizing over numerous theoretical and observational effects (Aguena et al. 2023), compared to the two-step approach.

### 4.2.4. Total likelihood and priors

To infer the parameters of the mass-richness relation, the total likelihood $\mathcal{L}_{tot}$ is given by considering either counts, cluster lensing profiles, or cluster lensing masses, or the combination between counts and cluster lensing, namely,

$$ \mathcal{L}_{tot} = \begin{cases} \mathcal{L}_{N} & \text{if: abundance alone,} \\ \mathcal{L}_{\Delta\Sigma} & \text{if: lensing profiles alone,} \\ \mathcal{L}_{M_{WL}} & \text{if: lensing masses alone,} \\ \mathcal{L}_{\Delta\Sigma} \times \mathcal{L}_{N} & \text{if: lensing profiles + abundance,} \\ \mathcal{L}_{M_{WL}} \times \mathcal{L}_{N} & \text{if: lensing masses + abundance.} \end{cases} \quad (28) $$

It is important to note that the way the lensing and count likelihoods are combined assumes no covariance between the two. As shown by Costanzi et al. (2019), the cross-correlation between the abundance and weak lensing inferred quantities is consistent with zero.

We drew samples from the parameter posterior distribution given by

$$ \mathcal{P}(\theta|\text{data}) = \frac{\mathcal{L}_{\text{tot}}(\text{data}|\theta) \times \pi(\theta)}{\mathcal{L}_{\text{tot}}(\text{data})}, \quad (29) $$

<footnote n="29">²⁹ Wu et al. (2019) developed the cluster-lensing-cov package (https://github.com/hywu/cluster-lensing-cov), which provides analytical predictions for the excess surface density profile covariance.</footnote>
where $\theta = (\ln \lambda_0, \mu_z, \mu_m, \sigma_{\ln\lambda_0}, \sigma_z, \sigma_m)$, $\mathcal{L}_{tot}(data|\theta)$ is the total likelihood in Eq. (28) ($\mathcal{L}_{tot}(data)$ is the likelihood evidence) where “data” refers to either the stacked lensing profiles, the inferred mean lensing masses, the cluster counts, or combinations of abundance and lensing information.

To compute the prior $\pi(\theta)$, we consider the choices made for the fiducial mass-richness relation and summarized in the second column of Table 1. All fits are performed at a fixed cosmology, chosen to be the fiducial cosmology of the DC2 simulation, which is close to the seven-year Wilkinson Microwave Anisotropy Probe (WMAP) best-fit $\Lambda$ cold dark matter values (Komatsu et al. 2011) given by $h = H_0/100 = 0.71$ (the reduced Hubble constant $H_0$), $\Omega_{cdm}h^2 = 0.1109$ ($\Omega_{cdm}$ is the current fractional energy density associated with cold dark matter), $\Omega_b h^2 = 0.02258$ ($\Omega_b$ is the current fractional energy density associated with baryons), $n_s = 0.963$ (spectral index of the primordial power spectrum), $\sigma_8 = 0.8$ (dispersion of matter density fluctuations smoothed over $8 h^{-1}$ Mpc), and $w = -1$ (dark energy equation of state).

### 4.3. Software

This work allows us to demonstrate the use of some of the software tools developed in DESC (in combination with others) for concrete cluster analysis. For this work, we use several DESC software tools that are currently in development as part of the DESC pipeline for analyzing the upcoming LSST data. To extract the cluster, halo, and background galaxy catalogs from the DC2 dataset, we used the DESC package GCRCatalogs<fnref n="30" />, along with Qserv<fnref n="31" />, an open-source SQL database system (Massively Parallel Processing) originally designed to host LSST data. For the estimation of the stacked lensing profiles, we used the DESC Cluster Lensing Mass Modeling (CLMM; Aguena et al. 2021b) package<fnref n="32" />, which provides various tools for estimating cluster lensing profiles as well as for halo modeling. The Core Cosmology Library<fnref n="33" /> (CCL; Chisari et al. 2019) was used to predict the halo mass function, halo bias, and concentration-mass relations. A combination of these codes was produced within the LSSTDESC/CLCosmo Sim repository<fnref n="34" />. To compute the binned $S_{ij}$ terms for the cluster count covariance in Eq. (25), we use the PySSC<fnref n="35" /> package (Lacasa et al. 2018; Gouyou Beauchamps et al. 2022). Cluster scaling relation parameters were inferred using the Markov chain Monte Carlo (MCMC) implementation in the emcee package (Foreman-Mackey et al. 2013), with visualization performed via getdist (Lewis 2019). Finally, we also use the Minuit<fnref n="36" /> minimizer (James & Roos 1975) to fit the mean cluster mass per redshift-richness bin in the two-step procedure.

## 5. Results

All the ingredients are now in place to perform the analysis of the scaling relation of redMaPPer clusters. In Sect. 5.1, we begin by presenting the setup that serves as our baseline analysis. Subsequent sections present variations of this baseline, changing one ingredient at a time. Sect. 5.2 explores how the mod-

eling choices for the halo density and mass-concentration relation impact the results, while Sect. 5.3 focuses on observational systematic effects, specifically the impact of source photometric redshifts and the shear-richness covariance. Table 2 summarizes the tables and figures associated with the different study cases.

### 5.1. Defining the count and weak lensing baseline analyses: $N + M_{WL}$ versus $N + \Delta\Sigma$

We model the per-cluster excess surface density, $\Delta\Sigma(R|m, z)$, using a three-dimensional dark matter density profile given by

$$ \rho(r) = \rho_s \left(\frac{r}{r_s}\right)^{-1} \left(1 + \frac{r}{r_s}\right)^{-n} \quad (30) $$

For $n = 2$, the Eq. (30) corresponds to the NFW profile (Navarro et al. 1997). Then, $r_s = r_{200c}/c_{200c}$ is the scale radius, and $\rho_s$ is the scale density<fnref n="37" />. We use the concentration mass-relation of Duffy et al. (2008) (see 5.2.1 hereafter).

We restricted the fitting to $1 < R < 3.5$ Mpc. The typical value of $R_{max} = 3.5$ Mpc is inspired by Lee et al. (2018), Giocoli et al. (2021), Cromer et al. (2022), Murray et al. (2022), and Bocquet et al. (2024) and is typically the distance beyond which the two-halo term increasingly dominates the cluster lensing signal and cannot be neglected (this choice is discussed further in Appendix B and later in the text).

We show in the left panel of Figure 6 the 68% and 95% credible intervals on the cluster scaling relation when using (separately) the cluster abundance likelihood, stacked lensing profiles, or mean cluster masses. The vertical and horizontal lines represent the fiducial values of the redMaPPer cluster scaling relation inferred in Sect. 3.4.1. We observe that the abundance provides tighter constraints on the mean parameters, but has different degeneracies compared to lensing (either stacked profiles or mean masses). The contours exhibit characteristic degeneracies and significant overlap, which suggests they can be combined. Additionally, we note that the posteriors are fairly compatible with the fiducial constraints at the $1-2\sigma$ level.

The different degeneracy directions between the abundance and lensing contours provide tight final constraints on the scaling parameters when combined<fnref n="38" />. By combining abundance with lensing, we obtain the posteriors shown in the right panel of Figure 6, either with stacked profiles or mean masses<fnref n="39" />. We observe that using either mean masses or stacked profiles provides roughly the same constraints with comparable precision. The errors on the mean scaling parameters ($\ln \lambda_0, \mu_z$, and $\mu_m$) are reduced by a factor of ~7, while the errors on the scatter terms ($\sigma_{\ln\lambda_0}, \sigma_z$, and $\sigma_m$) are reduced by a factor of ~3–4. The one-dimensional posteriors obtained from our joint lensing and abundance likelihood fairly recover the fiducial values at the $<2\sigma$ level; however, we observe a $>2\sigma$ tension with the projected posterior on $\mu_z - \sigma_z$.

$$ \rho_s = \frac{200}{3} \frac{c_{200c}^3}{f(c_{200c})} \rho_c(z) \quad \text{where} \quad f(x) = \ln(1 + x) - x/(1 + x). \quad (31) $$

<footnote n="30">³⁰ https://github.com/LSSTDESC/gcr-catalogs</footnote>

<footnote n="31">³¹ https://qserv.lsst.io/index.html</footnote>

<footnote n="32">³² https://github.com/LSSTDESC/CLMM</footnote>

<footnote n="33">³³ https://github.com/LSSTDESC/CCL</footnote>

<footnote n="34">³⁴ https://github.com/LSSTDESC/CLCosmo Sim</footnote>

<footnote n="35">³⁵ https://github.com/fabienlacasa/PySSC</footnote>

<footnote n="36">³⁶ https://scikit-hep.org/iminuit/about.html</footnote>

<footnote n="37">³⁷ The scale density is given by</footnote>

<footnote n="38">³⁸ Again, we recall that the cosmological parameters are fixed to the DC2 values in the MCMC.</footnote>

<footnote n="39">³⁹ To compare joint constraints with abundance-alone or lensing-alone contours, we show in black in the left panel of Figure 6 the joint abundance and lensing profiles posterior.</footnote>
<table><caption>**Table 2.** Different study cases explored in Sect. 5.</caption><thead><tr><th>Analysis</th><th>Section</th><th>Numerical values</th><th>Posterior distribution</th></tr></thead><tbody><tr><td>Fiducial relation</td><td>Sect. 3.4.1</td><td>Table 1</td><td>—</td></tr><tr><td>Baseline</td><td>Sect. 5.1</td><td>Table E.1 (first block)</td><td>Figure 6</td></tr><tr><td>Modeling – c(M) relation</td><td>Sect. 5.2.1</td><td>Table E.1 (second block)</td><td>Figure 7 (left panel)</td></tr><tr><td>Modeling – dark matter density profile</td><td>Sect. 5.2.2</td><td>Table E.1 (third block)</td><td>Figure 7 (right panel)</td></tr><tr><td>Observational – photometric redshifts</td><td>Sect. 5.3.1</td><td>Table E.1 (fourth block)</td><td>Figure 9</td></tr><tr><td>Observational – shear-richness covariance</td><td>Sect. 5.3.2</td><td>Table E.1 (fifth block)</td><td>Figure 12</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Left: Posterior distribution of the scaling relation parameters using abundance alone (blue), lensing profiles alone (dashed red), or masses alone (solid green lines). The joint abundance and lensing profiles posterior is displayed in black. Right: Posterior distribution of scaling relation parameters when combining abundance and lensing profiles (solid blue lines) and combining abundance and cluster masses (orange-filled contours). The dashed contours are obtained by removing the last redshift bin for both abundance and lensing.</figcaption></figure>

We checked whether the stacked lensing profiles in the last redshift bin are still impacted by the ray-tracing shear resolution beyond $R = 1$ Mpc, despite our conservative radial cut, since the attenuation increases with redshift and is therefore maximal in the last bin. To assess this, we removed the last redshift bin from the estimation of the posterior (from both abundance and lensing profiles). This yields results compatible with our baseline analysis, albeit with wider error bars (the last redshift bin contains 30% of the full cluster catalog). This provides a stability check for our lensing and count pipeline, which offers a reasonable comparison to the “fiducial” constraints of the redMaP-Per mass-richness relation. Residual systematics inherent to our analysis – such as the calibration of the halo mass function and uncertainties in the selection function – may still induce some bias in the recovered parameters.

We also examine in Appendix B the impact of alternative $R_{max}$ cuts on the scaling relation parameters. We find that using the fully available radial range (up to $R_{max} = 10$ Mpc) without modeling the two-halo term yields results comparable to those obtained with the baseline $R_{max} = 3.5$ Mpc cut, showing only small biases (<1σ) and similar error bars. This is attributed to the noisy lensing measurements at larger scales, which diminish the impact of the two-halo term and do not provide

additional constraining power. We adopt the $R_{max} = 3.5$ Mpc cut throughout this work. However, the two-halo term cannot be neglected for larger datasets, which will provide a larger S/N in the two-halo regime. In the following, we consistently use the combination of lensing (either mean masses or stacked profiles) and cluster counts to probe the underlying cluster scaling relation, leveraging their complementarity in parameter space.

## 5.2. Impact of modeling choices (in the $N + M_{WL}$ framework)

The inference of the cluster scaling relation is impacted by several systematic effects, stemming from the many degrees of freedom in the modeling of cluster observables. The choice of a particular mass distribution model can significantly affect the inferred cluster lensing mass and, consequently, have a non-negligible impact on the mass-richness relation. We restricted our analysis to the one-halo regime, namely at scales below $R = 3.5$ Mpc, which is more sensitive to the internal halo properties (such as cluster mass and concentration) than the two-halo regime, which is primarily influenced by cosmological quantities through the power spectrum and halo bias (Oguri & Takada 2011).
<figure><img src="image_7.png" /><figcaption>Fig. 7. Left: Posterior distribution of the scaling relation parameters using count/mass likelihood. The concentration-free case is shown with the dashed lines, the other contours are obtained using various concentration-mass relations from the literature. Right: Posterior distribution of the scaling relation parameters using count/mass likelihood, without using a concentration-mass relation, and using different modeling of the cluster density profile, namely NFW (red filled contours), Einasto (solid lines), and Hernquist (dashed lines).</figcaption></figure>

### 5.2.1. Concentration-mass relation

As mentioned in Sect. 2, the cluster density profile depends on the concentration parameter $c_{200c}$. Cosmological simulations indicate that cluster concentration and mass are statistically correlated. This correlation is typically quantified by a fitting formula known as the concentration-mass relation, labeled $c(M)$. The concentration $c$ of an individual cluster with mass $M$ is then scattered around $c(M)$, with a dispersion denoted $\sigma_{c|M}$. Several $c(M)$ relations exist in the literature (see, e.g., Diemer & Kravtsov 2014; Duffy et al. 2008; Prada et al. 2012; Bhattacharya et al. 2013; Diemer & Joyce 2019). For example, our baseline choice in Sect. 5.1 adopts the $c(M)$ relation of Duffy et al. (2008), calibrated using N-body simulations based on the five-year WMAP cosmology (Komatsu et al. 2009), which differs by at most 5% from the seven-year WMAP cosmology (Komatsu et al. 2011) used in cosmoDC2. The alternative relations examined in this section (Prada et al. 2012; Bhattacharya et al. 2013; Diemer & Kravtsov 2014) have been shown to accurately describe the relationship between halo mass and concentration across a variety of Λ cold dark matter cosmologies. These relations, despite some intrinsic scatter, capture the general trend observed in simulations. Employing such a relation can reduce the number of free parameters in cluster-based analyses and thus improves constraints on other relevant parameters.

To test the impact of the $c(M)$ relation on the scaling relation parameters, we (i) fit the mean cluster masses using various $c(M)$ relations from the literature (where $M$ denotes the mean mass, considered as a free parameter), (ii) use these mean mass estimates in Eq. (27), and (iii) combine them with the cluster count likelihood from Eq. (24) to infer the scaling relation parameters. Alternatively, the concentration can be treated as a free parameter in each redshift-richness bin. In this case, we jointly fit the concentrations and mean masses

in each bin and use the resulting mass estimates in the mass likelihood.

In the left panel of Figure 7 the concentration-free case shows reasonable agreement with the fiducial constraints. The posteriors obtained using different concentration-mass relations are fully consistent with the free-concentration case. When adopting a fixed $c(M)$ relation, the inferred posteriors exhibit the same level of parameter correlation, with slightly reduced error bars due to the smaller number of free parameters in the fitting procedure. We note that the choice of concentration model is expected to have a larger impact in the inner regions ($R < 1$ Mpc), which are not probed here due to the limitations imposed by ray-tracing resolution.

### 5.2.2. Modeling of the matter density profile

Besides the NFW parameterization, Einasto (1965) suggested an empirical profile that has been shown to provide a slightly better fit to N-body simulations of galaxy clusters compared to the NFW profile – particularly at small scales and over a broad range of halo masses and redshifts (Klypin et al. 2016; Sereno et al. 2016; Wang et al. 2020). The corresponding three-dimensional density profile is given by

$$ \rho^{\text{ein}}(r) = \rho_{-2} \exp\left(-\frac{2}{\alpha}\left[\left(\frac{r}{r_{-2}}\right)^{\alpha} - 1\right]\right), \quad (32) $$

where $r_{-2} = r_{200c}/c_{200c}$ and the $\rho_{-2}$ is given in the footnote$^{40}$. In the above equation, $\alpha$ is the shape parameter. We adopt a typical

$$ \rho_{-2} = 200\rho_c(z) \frac{2}{3} \gamma^{-1} \left[ \frac{3}{\alpha}, \frac{2}{\alpha} c_{200c}^{\alpha} \right] \frac{c_{200c}^3}{e^{2/\alpha}} \left( \frac{2}{\alpha} \right)^{\frac{3-\alpha}{\alpha}}, \quad (33) $$

where $\gamma$ is the lower incomplete gamma function.

<footnote>$^{40}$ The density $\rho_{-2}$ is given by</footnote>
constant value of $\alpha = 0.25$, though it is worth noting that more sophisticated models incorporate mass, redshift, and cosmology dependence in $\alpha$, as observed in simulations (Gao et al. 2008). In such models, $\alpha$ typically ranges from 0.15 to 0.3 as a function of mass. Less commonly, the Hernquist density model (corresponding to $n = 1$ in Eq. (30)) has been used, for example, in Buote & Lewis (2004), Sanderson & Ponman (2009) to model the mass density of X-ray-detected galaxy clusters, though it has not been extensively used in the lensing literature.

Concentration-mass relations for the Einasto and Hernquist models have not been extensively studied in the literature. Therefore, we consider the fit without assuming a specific $c(M)$ relation, prioritizing the count/mass cluster likelihood instead. In the right panel of Figure 7, we show the posterior distribution of the scaling relation parameters when considering the count/mass likelihood, using the NFW, Einasto, and Hernquist profiles. The results are in good agreement with each other, likely due to the radial range used in the fit, where the profiles are similar and only start to diverge at smaller scales. By fitting the mean cluster mass from 1 Mpc to 3.5 Mpc, where the one-halo regime is dominant, and using either the NFW, Einasto, or Hernquist profile, we found that the inferred scaling relation remains fairly stable.

### 5.3. Observational systematic effects (in the $N + \Delta\Sigma$ framework)

We addressed several potential sources of bias arising from uncertainties in the modeling of the cluster lensing observable. Another major source of bias in cluster lensing analyses comes from the data itself, including the calibration of galaxy shapes (Hernández-Martín et al. 2020), the redshift distribution of the background galaxy sample (Wright et al. 2020), contamination of the source galaxy sample by foreground galaxies (Varga et al. 2019), and miscentering in the cluster catalog (Sommer et al. 2022).

Using the cosmoDC2 catalog, we rely on the true shapes of galaxies (i.e., accounting only for the intrinsic shape and local shear), and thus cannot investigate the effects of shape measurement errors with this dataset. In the following, we focus instead on the impact of the calibration of the background galaxy photometric redshift distribution on the cluster scaling relation (Sect. 5.3.1), as well as on the impact of shear-richness covariance (Sect. 5.3.2). With this simulated dataset and the analysis choices we made, the effects of miscentering and contamination of the source sample by cluster members, which generally need to be accounted for, are not significant; they are discussed in Appendices C.1 and C.2, respectively.

#### 5.3.1. Source galaxy photometric redshifts

In Sects. 5.1 and 5.2, we use the true galaxy redshifts to compute the geometrical lensing weights $w_{ls}^{\text{geo}}$. In contrast, with imaging data, each source galaxy behind clusters is assigned a photometric redshift distribution, derived from galaxy magnitudes measured in several optical bands. The uncertainty in the measured redshift can be incorporated into the lens-source geometrical weights through Eq. (21).

In this section, we test the impact of using the BPZ (Benítez 2011) and FlexZBoost (Izbicki & Lee 2017) output photometric redshifts on the inference of the cluster scaling relation. To do this, we modify the source selection, which is now based on photometric-related quantities. Figure 8 shows the mean photometric redshift and its dispersion as a function of the true cos-

moDC2 redshift (no quality cut is applied to the photometric-derived quantities). As FlexZBoost provides robust estimates up to $z = 3$, BPZ yields more dispersed and biased results, making the use of BPZ redshifts without quality cuts problematic above $z \sim 2$, and potentially biasing the constraint of the mass-richness relation. The increased dispersion and bias of BPZ estimates at higher redshifts is not unexpected, as the model photometry in the simulation included a nonphysical dust feature in some SEDs at high redshift that was not present in the SED templates used to compute theoretical fluxes by BPZ, leading to a bias and increased dispersion.

As discussed in Sect. 3, none of these estimations are realistic due to (i) the deep selection of the training sample up to $i < 25$ for FlexZBoost, which makes the results optimistic, and (ii) the “discreteness” in the color-redshift space of the modeled galaxies (Korytov et al. 2019), which negatively impacts the performance of BPZ. Updated redshifts are currently being produced by DESC. These updated redshifts were not available at the time of this work, so we proceeded with the first version of the photometric redshift catalogs. The actual LSST data will likely fall between the pessimistic BPZ and optimistic FlexZBoost runs we use. By considering both approaches, we can bracket the impact of photometric redshifts on the results.

For the source sample selection, we apply two cuts simultaneously: (i) based on the mean photometric redshift, namely, $\langle z_s \rangle > z_l + 0.2$, and (ii) based on the probability density function, with $P(z_s > z_l) > 0.8$, where

$$P(z_s > z_l) = \int_{z_l}^{+\infty} dz_s p(z_s) \qquad (34)$$

is the probability for a galaxy to be located behind the cluster. In the above equation, $p(z_s)$ is the galaxy photometric redshift distribution (obtained either from BPZ or FlexZBoost). Moreover, we do not apply any quality cut on the photometric redshifts, allowing us to test the robustness of the cluster scaling relation with respect to photometric redshifts in the “worst-case” scenario; namely, with no mitigation of the impact of outliers on the final results. We show in Figure 8 (right, upper panel) the distribution of true galaxy redshifts in cosmoDC2 after our magnitude cuts. We also show the background source sample after applying the $P(z_s > z_l)$ cut alone and after combining the $P(z_s > z_l)$ and mean redshift cuts on FlexZBoost redshifts, for a redMaPPer cluster at redshift $z = 0.38$. We see that our photometric source selection efficiently removes foreground and member galaxies. The same distributions are shown in Figure 8 (right, lower panel) for BPZ redshifts.

We show in Figure 9 the constraints on the scaling relation when considering photometric redshifts from either BPZ (left panel) or FlexZBoost (right panel). Using BPZ introduces a small negative bias of $\sim 1\sigma$ in the normalization of the scaling relation, as well as in the mass dependence, compared to the true redshift case.

We incorporate a simple model to account for a possible bias in photometric redshift calibration at the level of the stacked lensing profile model, as proposed by Simet et al. (2017). We correct the modeling of the stacked excess surface density profile in Eq. (13) in each redshift-richness bin by a unique multiplicative factor $(1 + b)$, common to all redshift-richness bins. This factor $(1 + b)$ is fitted jointly with the scaling relation parameters, using a flat prior of $[-0.5, 0.5]$ for $b$. The results are shown in Figure 9. For BPZ, we see that adding a free parameter to the fit increases the size of the posteriors as expected and helps decrease the tension between the photometric redshift case and the true redshift case. We recover a negative bias
<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Left: Mean photometric redshift against true cosmoDC2 redshift for BPZ and FlexZBoost methods. Right: Distribution of true cosmoDC2 redshifts after our magnitude cuts (black). Upper panel: The blue distribution represents the same distribution after a cut on the probability $P(z_s > z_i)$ using FlexZBoost photometric redshift, for a cluster at $z = 0.38$. The orange distribution is obtained after combining the $P(z_s > z_i)$ cut and the mean redshift cut. Lower panel: Same as upper panel but for BPZ redshifts. On both panels, the vertical red line indicates the position of the cluster redshift.</figcaption></figure>

$b = -0.02 \pm 0.03$, which is compatible with 0 at the $1\sigma$ level. For FlexZBoost (right panel), we first observe that without correction, the constraints are fully compatible with the true redshift case. Adding $(1 + b)$, we recover $b = 0.02 \pm 0.03$, and although the error bars increase, the posteriors do not shift significantly.

### 5.3.2. Impact of shear-richness covariance

In this section, we discuss a systematic bias arising from selection biases in cluster finder algorithms. For optically selected clusters, one relatively unexplored category of cluster systematics is the covariance between the two cluster observables, first, the shear estimated from background galaxies and second the richness, as an outcome of cluster finders. This covariance was left unconstrained in previous cluster cosmology analyses (McClintock et al. 2019; DES Collaboration 2020). This property covariance can induce additive biases that cannot be mitigated with increased sample size and reduced shape noise (Nord et al. 2008; Evrard et al. 2014; Farahi et al. 2018; Wu et al. 2019). To achieve percent-level mass calibration of clusters, this new category of systematics must be accurately and precisely quantified (Rozo et al. 2014; Zhou et al. 2024). In this work, we follow the prescription of Zhang et al. (2024) in modeling and quantifying the systematic uncertainty induced by the property covariance of optically selected clusters on the cluster scaling relation.

We make a note of the difference between intrinsic and extrinsic covariance. In Zhang et al. (2024) the authors measured the intrinsic covariance of cluster weak lensing observables by encircling galaxies within a three-dimensional physical radius around the halo to measure the “true” richness and second, measured the integrated dark matter density around the halo. This intrinsic contribution to the covariance originates from the galaxy assembly bias, the degree of which can be quantified

by secondary halo properties (e.g., concentration, mass accretion rate, kinetic-to-potential energy ratio) that relate to the formation history of the cluster-sized halo. On the other side, the extrinsic covariance arises from systematic biases introduced by the cluster-finding algorithm, such as projection effects (see, e.g., Myles et al. 2021, 2025). In this work, we measure the total covariance (intrinsic + extrinsic) using the redMaPPer cluster finder to determine the observed richness and the galaxy shear as a proxy for the surface density. This realistic mock catalog introduces observational systematics such as photometric redshift uncertainty (Graham et al. 2017) and shape noise (Wu et al. 2019) for the galaxy shear signal and projection and percolation effects (Costanzi et al. 2018), triaxiality bias (Zhang et al. 2023) and a set of other observational and modeling systematics (McClintock et al. 2019) related to a realistic cluster finder as redMaPPer.

Our goal is to model and correct for the systematic uncertainty in the estimated stacked lensing profiles in Eq. (20). From Evrard et al. (2014), Aihara et al. (2018), Wu et al. (2019), Zhang et al. (2024), the corrected stacked excess surface density profile is given by

$$ \Delta\Sigma_{ij}^{\text{corr}} = \Delta\Sigma_{ij} + \ln(10) \frac{[\beta_1]_{ij}}{\mu_m} \times \langle \text{Cov}(\Delta\Sigma, \ln \lambda|m, z) \rangle_{ij}, \quad (35) $$

where $\Delta\Sigma_{ij}$ is given in Eq. (13), and $[\beta_1]_{ij}$ is the average logarithmic slope of the Despali et al. (2015) halo mass function in the $ij$-th redshift-richness bin. For the latter, we compute $[\beta_1]_{ij}$ once before running MCMC, using the fiducial cluster scaling relation, and we consider the denominator $\mu_m$ to be free. $\langle \text{Cov}(\gamma, \ln \lambda|m, z) \rangle_{ij}$ is the averaged shear-richness covariance (see calculation details in Appendix D and also in Zhou et al. 2024), associated with potential selection bias in the cluster finder algorithms. The factor $\ln(10)$ is used to match the definition in Evrard et al. (2014). In the latter, we present the
<figure><img src="image_9.png" /><figcaption>Fig. 9. Impact of photometric redshifts of source galaxies on scaling relation parameter estimation. Left: The filled contours show the constraints when using true redshifts of source galaxies. Dashed contours are obtained when considering the BPZ photometric redshifts. Empty contours are obtained using BPZ photometric redshifts, after marginalizing over the (1 + b) correcting factor, considered as a free parameter, used to correct from photometric redshift calibration issues. Right: same as left panel, but for FlexZBoost.</figcaption></figure>

estimation of $\langle \text{Cov}(\Delta\Sigma, \ln \lambda|m, z) \rangle_{ij}$, to be used in the estimation of the cluster scaling relation parameters.

While Zhang et al. (2024) observed a radially dependent covariance at small scales $R < R_{200c}$ due to cluster formation physics, because of the attenuation of the cluster lensing signal at small scales due to the ray tracing resolution, we focus, instead, on the $R > 1$ Mpc scales. From mock-cluster studies analyzing selection effects such as Sunayama et al. (2020), Wu et al. (2022), Zhou et al. (2024), the total bias (intrinsic+extrinsic) is close to null at small scales and increases at larger scales for $R > R_{200c}$. For this reason, the covariance after applying the scale cut due to lensing attenuation may be sufficient for describing to impact of the total covariance.

We apply a kernel local linear regression (Farahi et al. 2022) to find the best fit local mean with mass as the independent variable and $\Delta\Sigma(R)$ and $\ln \lambda$ as the dependent variables. Specifically, the residuals of $\Delta\Sigma$ and $\ln \lambda$ are taken around their local mean quantities $\langle \Delta\Sigma | m, z, R \rangle$ and $\langle \ln \lambda | m, z \rangle$ such as

$$ \text{res}_{\Delta\Sigma}(m, z, R) = \Delta\Sigma - \langle \Delta\Sigma | m, z, R \rangle, \qquad (36) $$

$$ \text{res}_{\ln \lambda}(m, z) = \ln \lambda - \langle \ln \lambda | m, z \rangle, \qquad (37) $$

and the covariance is measured around the local residuals, namely,

$$ \langle \text{Cov}(\Delta\Sigma, \ln \lambda|m, z) \rangle_{ij} = \text{Cov}(\text{res}_{\Delta\Sigma}(m, z, R), \text{res}_{\ln \lambda}(m, z)). \qquad (38) $$

It is important to note that as we remove the residual dependence of the mean mass on the covariance, the size of the covariance itself can still nonetheless be a function of $(m, z)$, namely, the scatter and correlation are modeled as $\sigma_{\ln \lambda|M, z}(M, z)$ and $\sigma_{\Delta\Sigma|m, z}(m, z)$ and $r_{\ln \lambda-\Delta\Sigma|m, z}(m, z)$. For this reason, we test the covariance on small enough $(m, z)$ bins that the scatter can be considered homoskedastic, or in other words, approximated as a constant. In our benchmark tests, we choose to bin the covariance either in mass or redshift but not by both to retain enough

statistical constraint, especially for high mass, low redshift bins with a low cluster count.

We plot the covariance merged by mass and binned in redshift in Figure 10. Visually, we observe a slight positive covariance in most bins. A similar trend (not shown) is found in a benchmark test for covariances binned by richness and merged over $z \in [0.2, 0.7]$ (it can be shown that the covariances in $z \in [0.7, 1.0]$ are negligible). To explicitly show the mass and redshift dependence, we model the covariance as a constant bias across radius, which we denote as $b_{\text{Cov}}(M, z)$. The covariance trends with respect to richness (left panel) and redshift (right panel) are plotted in Figure 11. From the figure, we see that the amplitude of the effect is smaller by a factor of 10 compared to the amplitude of the stacked lensing profile (i.e., $10^{12}$ versus $10^{13}$).

To constrain the mass-richness relation while accounting for this effect, we modeled the stacked excess surface density profile using Eq. (35). We propagated the error of the shear-richness covariance (represented by the error bars in Figure 10) by adding the associated variance to the diagonal element of the total excess surface density profile covariance matrices, scaled by $[\log(10)\beta_1/\mu_m]^2$. We used the stacked profiles estimated with true background source redshifts and employed the inference setup defined in our baseline analysis, which involves a joint likelihood between counts and the stacked lensing profile. In Figure 12 (left panel), the filled (resp. dashed) contours show the posterior of the scaling relation parameters when accounting for (resp. not accounting for) the shear-richness covariance correction as calculated in Eq. (35). We observe that the “corrected” posterior is shifted by approximately $0.5\sigma$ compared to the baseline analysis. From this test, we find that the impact of property covariance significantly influences the scaling parameter inference. The size of this bias is larger than the one obtained using only FlexZBoost redshifts (see Sect. 5.3.1), but smaller than that from using only BPZ redshifts. Posterior shifts of roughly the same amplitude as the FlexZBoost-only (resp. BPZ-
<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Binned excess surface density-richness covariance in Eq. (38), as a function of the radius from the cluster center. We stack all clusters between $\lambda \in [20, 70)$ and subdivide by different redshift bins from the top left panel to the bottom right panel. The true redshift case is displayed in pink, the BPZ and FlexZBoost cases are represented in orange and blue, respectively.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 11.** Covariance in Eq. (38) modeled as a constant bias term $b_{\text{Cov}}(M, z)$ across radius when binned by richness (left panel) or redshift (right panel). The positive covariance at large scales at lower redshift and richness bins is consistent with expectations from projection effects.</figcaption></figure>

only) are obtained when using FlexZBoost (resp. BPZ) photo-
metric source redshifts in addition to shear-richness covariance,
as shown in the left panel of Figure 12. Finally, we show in the
right panel of Figure 12 all sources of observational systematic
effects that we tested in this work, namely the combination of
(i) source photometric redshifts (either FlexZBoost or BPZ) (ii)
the inclusion of the shear-richness covariance, and (iii) the addi-

tion of a free correcting photometric factor (1 + b) in the fit, to
be compared with the baseline case (blue filled contours: true
source redshift, and no shear-richness covariance).

As mentioned above, the amplitude of the shear-richness
covariance term is primarily influenced by the galaxy formation
history within clusters and, secondarily, by the method used to
 populate halos in the simulation. As mentioned in Korytov et al.
<figure><img src="image_12.png" /><figcaption>**Fig. 12.** Impact of shear-richness covariance on scaling relation parameter estimation. Left: Posterior distribution of the scaling relation parameters with (dashed blue) and without (filled blue) the shear-richness covariance correction in Eq. (35), considering true source galaxy redshifts. The red (resp. black) unfilled contours are obtained considering the shear-richness covariance and using FlexZBoost (resp. BPZ) photometric redshifts for source galaxies. Right: The blue filled contours are the same as in the left panel. The purple (resp. dashed green) unfilled contours are obtained considering the shear-richness covariance and using FlexZBoost (resp. BPZ) photometric redshifts, and letting the bias (1 + b) as a free parameter. These two effects explain somewhat the low amplitude (about 1 to 10%) of the shear-richness covariance (let us note that the ln(10)β₁/μm factor in Eq. (35) is the order of unity; it does not change our conclusions).</figcaption></figure>

(2019), after defining the number of galaxies per halo, galaxy positions are drawn randomly in the halo environment (according to a generic NFW profile), thus not following the halo complex shape. As a result, the simulated shear-richness correlation is likely suppressed relative to what would be expected in real data. SkySim5000 (Abolfathi et al. 2021), the extended version of cosmoDC2 covering 5000 deg², features improved ray-tracing resolution and a more realistic model for galaxy population in halos. In particular, galaxies are distributed based on the triaxial configuration of halo dark matter particles, better capturing the underlying mass distribution. Leveraging SkySim5000 in future work will allow for a more realistic characterization of shear-richness covariance, especially at smaller scales and over a larger cluster sample.

Second, the shear-richness covariance also depends on the performance of the cluster finder. In our analysis, we used the redMaPPer cluster catalog derived from mock datasets, where the detection process is rather idealized compared to what would be expected from real observational data.

# 6. Summary and conclusions

We conducted several analyses based on the lensing profiles and abundance of redMaPPer-detected clusters in the DC2 simulation to constrain their mass-richness relation. We defined a baseline analysis, with state-of-the-art modeling choices and estimators. On this basis, we were able to test the impact of other modeling choices and observational systematics on the inferred mass-richness relation, also comparing it to a “fiducial” scaling relation. To do that, we used either stacked cluster lensing profiles or mean cluster lensing masses in combination with cluster count to explore several analysis setups: fixing or varying the

cluster concentration, changing the dark matter density profile, adding shear-richness covariance.

## 6.1. Summary

For the various study cases explored in this paper, the best-fit values of the six-parameter redMaPPer mass-richness relation are summarized in Table E.1, with the corresponding study cases and their associated paper sections listed in Table 2.

In Figure 13, we illustrate the constraints on the three key parameters of the mean mass-richness relation defined in Eq. (16): ln λ₀, μz, and μm. The vertical bands correspond to the fiducial constraints reported in Table 1, indicating the 1σ and 2σ confidence intervals. We first performed the analysis under ideal conditions, using the true shapes of background galaxies (i.e., their intrinsic ellipticities sheared by the local shear field) and their true redshifts. While the limited ray-tracing resolution of the DC2 simulation affects the innermost cluster regions (R < 1 Mpc), we found that combining cluster counts and lensing observables (whether via mean mass measurements or stacked profiles) yields tighter constraints on the scaling relation than using each probe independently. This improvement arises from the complementary degeneracy directions of the two observables. Using a NFW (Navarro et al. 1997) profile with the concentration-mass relation of Duffy et al. (2008), and restricting the radial fit to the one-halo regime (1 < R < 3.5 Mpc), we obtain constraints that are consistent at the 1σ level with the fiducial values derived from the match between redMaPPer clusters and dark matter halos. This consistency validates our joint count-and-lensing approach as a reliable method for constraining the mass-richness relation. We further tested the robustness of our results by varying the modeling assumptions of the stacked weak lensing profiles. In particular, we showed that the inferred
<figure><img src="image_13.png" /><figcaption>Fig. 13. Constraints on the scaling relation parameters (only ln λ₀, μₓ, and μₘ) obtained in Sect. 5. The vertical shaded region in each subplot represents the fiducial constraints presented in Sect. 3.4.1. For clarity, each color used for the plot corresponds to a subsection in Sect. 5.</figcaption></figure>

posterkors are stable when adopting different concentration-mass relations, and remain consistent within 1σ with the results of the free-concentration case. Similarly, using alternative dark matter density profiles instead of NFW(e.g., Einasto or Hernquist) also offers robust and consistent constraints on the scaling relation parameters.

On the data side, we tested the impact of photometric redshift uncertainties on the inference of the mass-richness scaling relation by using outputs from two photometric redshift algorithms applied to cosmoDC2 galaxy magnitudes: BPZ (Benítez 2011) and FlexZBoost (Izbicki & Lee 2017). To accommodate the performance of each algorithm, we adapted the source galaxy selection accordingly. We found that FlexZBoost yields results that are fully consistent with the true-redshift case, indicating no significant bias. In contrast, the use of BPZ introduces a bias greater than 1σ in the normalization of the scaling relation. However, this bias can be mitigated by introducing a multiplicative correction factor, fitted jointly with the standard scaling relation parameters. In addition, we investigated the effect of shear-richness covariance, which arises from selection biases in the cluster-finding process. We found that within the radial range considered for fitting the cluster lensing profiles, the impact of covariance alone can shift μₘ by 0.5σ. This is comparable to the shift from the different c(M) models or using BPZ redshifts.

We complemented the summary plot by quantifying the level of tension between the results of various analyses and the fiducial constraints, as well as with respect to the baseline analysis. To this end, we employed the Gaussian tension metric (Raveri & Doux 2021; Leizerovich et al. 2024), which expresses the tension between a Gaussian posterior, P, and a reference posterior, P_{fid}, for a parameter θ in units of standard deviations. The tension is defined as

$$n_{\sigma}(\theta) = \frac{\langle\theta\rangle_{P} - \langle\theta\rangle_{P_{\text{fid}}}}{\sqrt{\sigma^2(\theta)_{P} + \sigma^2(\theta)_{P_{\text{fid}}}}}, \quad (39)$$

where ⟨θ⟩P (resp. ⟨θ⟩P_{fid}) is the mean of the parameter posterior P(θ) (resp. P_{fid}(θ)), and σ²(θ)P (resp. σ²(θ)P_{fid}) is the associated variance. For the different study cases, the level of tension with respect to the fiducial constraints (resp. the baseline case: weak lensing profiles combined with counts) is shown as circular dots (resp. star symbols) in Figure E.1. This approach allows us to quantitatively assess the Bayesian distance between the outcomes of various analysis setups and the fiducial model.

We find that combining probes (counts and lensing) does not significantly increase the level of tension, while it does enhance the statistical precision of the constraints. This demonstrates that our counts+lensing pipeline yields a consistent and reliable description of the observables. Furthermore, across all tested analysis configurations, the tension remains below the 2σ threshold. This highlights the robustness of our baseline analysis: variations in modeling assumptions or data choices result in posterior shifts that consistently stay within 2σ, underscoring the stability of the inferred scaling relation parameters.

## 6.2. Discussions and conclusions

To fully harness the cosmological constraining power of galaxy clusters, achieving a high level of precision on the cluster scaling relation is essential. This can be accomplished by jointly analyzing cluster counts and cluster weak gravitational lensing, enabling the inference of the mass-richness relation and its intrinsic scatter. For such analyses to be reliable, constraints from both lensing and counts must be robust against systematic effects. In this paper, we focus specifically on systematics related to cluster weak lensing.

In light of the first LSST data, we investigate the mass-richness relation of redMaPPer-detected galaxy clusters within the redshift range of 0.2 < z < 1 and richness range of 20 < λ < 200, using the DESC DC2 simulated dataset. We inferred
a six-parameter log-normal cluster scaling relation, accounting for mass and redshift dependence in both the mean and scatter, via a joint analysis of stacked redMaPPer weak lensing profiles and cluster counts, thereby enhancing the overall constraining power. A series of robustness tests were performed, exploring the impact of cluster lensing profile modeling choices, photometric redshift uncertainties, and shear-richness covariance effects. We find that our constraints remain stable when lensing profiles are fitted in the one-halo regime ($1 < R < 3.5$ Mpc), assuming either true or photometric redshifts, and perfectly known source shapes. Furthermore, our results yield unbiased estimates of the fiducial redMaPPer mass-richness relation. While the constraining power in this analysis is limited compared to LSST expectations-to be $\sqrt{\Omega_{LSST}/\Omega_{DC2}} \approx 6.4$ times larger – this work provides a solid foundation for scaling relation inference with future LSST cluster samples.

Moreover, the code work allows us to demonstrate the use of some of the software tools developed within DESC in concrete cluster analysis, representing a significant step in the calibration, validation, and methodological development. This paper thoroughly addresses all major systematics, except for shape measurement, which is discussed in detail below. The framework employed in this analysis is well integrated into the CLMM and CCL pipelines, paving the way for the development of cluster-based analyses in other official DESC pipeline codes such as TJPCov<fnref n="41" /> (official DESC covariance calculator interface), Firecrown<fnref n="42" /> (providing the DESC framework for implementing likelihoods and sampling parameters), and TXPipe<fnref n="43" /> (DESC library for summary statistics measurement, see Prat et al. 2023).

Due to the lack of ray-tracing resolution that motivated our conservative cut, we excluded the innermost regions, where the cluster lensing signal is the strongest. These scales are, however, crucial to enhancing the precision of weak lensing cluster masses. Analyzing such radial ranges requires a particular attention to modeling the dilution of the cluster lensing signal by member galaxies (contamination; see, e.g., Varga et al. 2019) and miscentering due to uncertainties in cluster finding methods (Zhang et al. 2019; Köhlinger et al. 2015), which can be significant when referring to stacks of galaxy clusters.

Our analysis could benefit from further improvement to infer the redMaPPer cluster mass-richness relation. First, we considered the fiducial halo mass function of Despali et al. (2015), but propagating calibration uncertainties of this relation in our pipeline could have a non-negligible impact (Artis et al. 2021; Kugel et al. 2025). Due to the noise level, we neglected the off-diagonal terms in the covariances of the stacked lensing profiles. However, rather close radial bins are known to be somewhat correlated due to the intrinsic variation of the halo density profile and large-scale structures (Wu et al. 2019), which could degrade the precision of our results, even though they are much smaller than the galaxy shot and shape noise contribution.

Moreover, for the mean mass (two-step) approach, the link between the “true” mean mass of a stack and the inferred mean weak lensing mass may not be trivial (Becker & Kravtsov 2011), contrary to what we used in Eq. (14), namely, a standard average mass over the redshift-richness cluster distribution. In reality, the excess surface density profile does not simply scale with the underlying cluster mass. For instance, Melchior et al. (2017) proposed a simple model where the excess surface density follows $\Delta\Sigma(R|M) \propto M^{\Gamma(R|M)}$ where $\Gamma$ is the logarithmic slope of

the excess surface density<fnref n="44" />; as a result, the mean weak lensing mass $\langle M_{\text{WL}}^{\text{stack}} \rangle$ has to be modeled by $\sim \langle m^\Gamma \rangle^{1/\Gamma}$. Such nonlinear dependence can be accounted for in the mean mass modeling in Eq. (14). Beyond this rather simple correction to link the mean lensing mass to the true average mass of the stack, McClintock et al. (2019) used simulations to reconstruct the stack lensing signal around clusters with masses drawn from the known $P(M|\lambda)$ relation (as well as for concentrations and miscentering offsets drawn from specific distributions). They found that the bias between the true average mass and mean weak lensing mass may reach 10% (see their Figure 9) and has a strong richness and redshift dependence, indicating there a strong impact coming from the intrinsic scatter (in richness, concentration, miscentering, etc.).

Even though it is not sensitive to miscentering in our analysis (due to our restrictive low radial cut), we considered a fixed concentration-mass relation. However, it is known that cluster concentrations are scattered around a mean $c(M)$ relation (see, e.g., Bullock et al. 2001; Darragh-Ford et al. 2023). This effect is particularly important in stacked analysis, especially at small radii. Accounting for this effect should degrade the precision of our obtained constraints by lowering the sensitivity of our model.

Finally, noise in galaxy shapes is one of the current main sources of scatter in weak lensing cluster mass estimation. First, as already mentioned earlier in this work, it originates from the intrinsic scatter of galaxy shapes, which is more and more mitigated by the increasing statistics of ongoing and future optical surveys. Second, it originates from the galaxy shape measurement algorithms, returning dispersed and biased estimations of intrinsic galaxy ellipticities. Shape reconstruction methods generally need to be calibrated from simulation or internally. Dispersed measured shapes as well as calibration uncertainties in these algorithms inevitably propagate to the cluster mass estimates, downgrading their precision. For that purpose, the DC2 object catalog (Abolfathi et al. 2021) was built using the Rubin LSST science pipeline by identifying galaxies on cosmoDC2-based realistic images of the sky, including several realistic observational effects, such as atmospheric turbulence, telescope optics, and some detector effects. Several codes are used to measure shapes of detected galaxies, such as HSM (Hirata & Seljak 2003; Mandelbaum et al. 2005) used on the HSC data (Mandelbaum et al. 2018), and METACALIBRATION (Sheldon & Huff 2017), used in the DES-Y3 cosmic shear analysis in Abbott et al. (2022). Ongoing efforts in DESC are aimed at studying the impact of galaxy detection and shape measurement in cluster fields on the accuracy and precision of cluster mass estimates (Ramel et al. in prep). These results will be essential to handling the full error budget on the mass-richness relation, in preparation for LSST cluster cosmology.

*Acknowledgements.* The authors thank the anonymous reviewer for their insightful comments and suggestions. This paper has undergone internal review by the LSST Dark Energy Science Collaboration. The authors thank the internal reviewers Tomomi Sunayama and Shenming Fu for their valuable comments. The authors thank Fabien Lacasa for useful discussions and for their help in using the PySSC code. The authors thank Dominique Boutigny for his help in using Qserv to access the cosmoDC2 simulated datasets at CC-IN2P3. CP conceptualized the project, developed the code, led the analysis, contributed to the paper's structure and text, and oversaw the overall writing process. ZZ contributed to the shear-richness covariance analysis, the early development of the analysis framework and code, and the paper's structure and text. MA developed the DESC matching library CIEvaR and contributed to the paper's structure and text. CC assisted with project conceptualization and advising, partici-

<footnote n="41">⁴¹ https://github.com/LSSTDESC/tjpcov</footnote>

<footnote n="42">⁴² https://github.com/LSSTDESC/firecrown</footnote>

<footnote n="43">⁴³ https://github.com/LSSTDESC/txpipe</footnote>

<footnote n="44">⁴⁴ Melchior et al. (2017) found a typical value of $\Gamma = 0.74$ for an NFW profile.</footnote>
pated in the early development and review of the analysis framework and code, and contributed to the paper's structure and text. TG provided the parametrization of the redMaPPer selection function. MR prepared the matched cosmoDC2-redMaPPer catalog, reviewed various aspects of the analysis, and contributed to the paper's structure and text. The remaining authors generated the data used in this work, contributed to the project's infrastructure, and participated in discussions interpreting the results. The DESC acknowledges ongoing support from the Institut National de Physique Nucléaire et de Physique des Particules in France; the Science & Technology Facilities Council in the United Kingdom; and the Department of Energy and the LSST Discovery Alliance in the United States. DESC uses the resources of the IN2P3 Computing Center (CC-IN2P3-Lyon/Villeurbanne - France) funded by the Centre National de la Recherche Scientifique; the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231; STFC DIRAC HPC Facilities, funded by UK BEIS National E-infrastructure capital grants; and the UK particle physics grid, supported by the GridPP Collaboration. This work was performed in part under DOE Contract DE-AC02-76SF00515. MA is supported by the PRIN 2022 project EMC2 - Euclid Mission Cluster Cosmology: unlock the full cosmological utility of the Euclid photometric cluster catalog (code no. J53D23001620006). CA acknowledges support from DOE grant DE-SC0019193 and the Leinweber Center for Theoretical Physics. AF acknowledges support from the National Science Foundation under Cooperative Agreement 2421782 and the Simons Foundation award MPS-AI-00010515. We thank the developers and maintainers of the following software tools used in this work: NumPy (van der Walt et al. 2011), SciPy (Jones et al. 2001), Matplotlib (Hunter 2007), GetDist (Lewis 2019), emcee (Foreman-Mackey et al. 2013), Astropy (Astropy Collaboration 2013), Jupyter (Kluyver et al. 2016).

## References

Abbott, T., Aguena, M., Alarcon, A., et al. 2022, Phys. Rev. D, 105, 023520

Abdullah, M. H., Klypin, A., & Wilson, G. 2020, ApJ, 901, 90

Abell, P. A., Allison, J., Anderson, S. F., et al. 2009, ArXiv e-prints [arXiv:0912.0201]

Abolfathi, B., Alonso, D., Armstrong, R., et al. 2021, ApJS, 253, 31

Aguena, M., & Lima, M. 2018, Phys. Rev. D, 98, 123529

Aguena, M., Benoist, C., da Costa, L. N., et al. 2021a, MNRAS, 502, 4435

Aguena, M., Avestruz, C., Combet, C., et al. 2021b, MNRAS, 508, 6092

Aguena, M., Alves, O., Annis, J., et al. 2023, ArXiv e-prints [arXiv:2309.06593]

Aihara, H., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S8

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

Anbajagane, D., Evrard, A. E., Farahi, A., et al. 2020, MNRAS, 495, 686

Artis, E., Melin, J.-B., Bartlett, J. G., & Murray, C. 2021, A&A, 649, A47

Astropy Collaboration 2013, A&A, 558, A33

Bartlett, J. G. 1997, ASP Conf. Ser., 126, 365

Baxter, E. J., Rozo, E., Jain, B., Rykoff, E., & Wechsler, R. H. 2016, MNRAS, 463, 205

Becker, M. R., & Kravtsov, A. V. 2011, ApJ, 740, 25

Bellagamba, F., Roncarelli, M., Maturi, M., & Moscardini, L. 2018, MNRAS, 473, 5221

Benítez, N. 2011, Astrophysics Source Code Library [record ascl:1108.011]

Benson, A. J. 2012, New Astron., 17, 175

Bernstein, G. M., & Nakajima, R. 2009, ApJ, 693, 1508

Bhattacharya, S., Habib, S., Heitmann, K., & Vikhlinin, A. 2013, ApJ, 766, 32

Binggeli, B. 1982, A&A, 107, 338

Bocquet, S., Grandis, S., Bleem, L. E., et al. 2024, Phys. Rev. D, 110, 083510

Bulbul, E., Liu, A., Kluge, M., et al. 2024, A&A, 685, A106

Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, MNRAS, 321, 559

Buote, D. A., & Lewis, A. D. 2004, ApJ, 604, 116

Chang, C., Jarvis, M., Jain, B., et al. 2013, MNRAS, 434, 2121

Chen, K. F., Chiu, I. N., Oguri, M., et al. 2025, Open J. Astrophys., 8, 2

Chisari, N. E., Alonso, D., Krause, E., et al. 2019, ApJS, 242, 2

Chiu, I. N., Chen, K. F., Oguri, M., et al. 2024, Open J. Astrophys., 7, 90

Corless, V. L., & King, L. J. 2009, MNRAS, 396, 315

Costanzi, M., Rozo, E., Rykoff, E. S., et al. 2018, MNRAS, 482, 490

Costanzi, M., Rozo, E., Simet, M., et al. 2019, MNRAS, 488, 4779

Cromer, D., Battaglia, N., Miyatake, H., & Simet, M. 2022, JCAP, 2022, 034

Darragh-Ford, E., Mantz, A. B., Rasia, E., et al. 2023, MNRAS, 521, 790

DES Collaboration 2020, Phys. Rev. D, 102, 023509

DES Collaboration 2025, ArXiv e-prints [arXiv:2503.13632]

Despali, G., Giocoli, C., & Tormen, G. 2014, MNRAS, 443, 3208

Despali, G., Giocoli, C., Angulo, R. E., et al. 2015, MNRAS, 456, 2486

Diemer, B., & Joyce, M. 2019, ApJ, 871, 168

Diemer, B., & Kravtsov, A. V. 2014, ApJ, 789, 1

Duffy, A. R., Schaye, J., Kay, S. T., & Dalla Vecchia, C. 2008, MNRAS, 390, L64

Einasto, J. 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87

Escoffier, S., Cousinou, M. C., Tilquin, A., et al. 2016, ArXiv e-prints [arXiv:1606.00233]

Euclid Collaboration 2019, A&A, 627, A23

Evrard, A. E., Arnault, P., Huterer, D., & Farahi, A. 2014, MNRAS, 441, 3562

Farahi, A., Evrard, A. E., Rozo, E., Rykoff, E. S., & Wechsler, R. H. 2016, MNRAS, 460, 3900

Farahi, A., Evrard, A. E., McCarthy, I., Barnes, D. J., & Kay, S. T. 2018, MNRAS, 478, 2618

Farahi, A., Anbajagane, D., & Evrard, A. E. 2022, ApJ, 931, 166

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Fumagalli, A., Saro, A., Borgani, S., et al. 2021, A&A, 652, A21

Fumagalli, A., Costanzi, M., Saro, A., Castro, T., & Borgani, S. 2024, A&A, 682, A148

Gao, L., Navarro, J. F., Cole, S., et al. 2008, MNRAS, 387, 536

Ghirardini, V., Bulbul, E., Artis, E., et al. 2024, A&A, 689, A298

Giocoli, C., Marulli, F., Moscardini, L., et al. 2021, A&A, 653, A19

Gouyou Beauchamps, S., Lacasa, F., Tutusaus, I., et al. 2022, A&A, 659, A128

Graham, M. L., Connolly, A. J., Ivezić, Z., et al. 2017, AJ, 155, 1

Grandis, S., Bocquet, S., Mohr, J. J., Klein, M., & Dolag, K. 2021, MNRAS, 507, 5671

Grandis, S., Ghirardini, V., Bocquet, S., et al. 2024, A&A, 687, A178

Grishin, K., Mei, S., & Ilić, S. 2023, A&A, 677, A101

Gruen, D., Seitz, S., Becker, M. R., Friedrich, O., & Mana, A. 2015, MNRAS, 449, 4264

Hasselfield, M., Hilton, M., Marriage, T. A., et al. 2013, JCAP, 2013, 008

Hearin, A., Korytov, D., Kovacs, E., et al. 2020, MNRAS, 495, 5040

Heitmann, K., Finkel, H., Pope, A., et al. 2019, ApJS, 245, 16

Hernández-Martín, B., Schrabback, T., Hoekstra, H., et al. 2020, A&A, 640, A117

Hernquist, L. 1990, ApJ, 356, 359

Hetterscheidt, M., Erben, T., Schneider, P., et al. 2005, A&A, 442, 43

Hirata, C., & Seljak, U. 2003, MNRAS, 343, 459

Hoekstra, H. 2003, MNRAS, 339, 1155

Hogg, D. W. 1999, ArXiv e-prints [arXiv:astro-ph/9905116]

Hu, W., & Kravtsov, A. V. 2003, ApJ, 584, 702

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

Izbicki, R., & Lee, A. B. 2017, ArXiv e-prints [arXiv:1704.08095]

James, F., & Roos, M. 1975, Comput. Phys. Commun., 10, 343

Jimenó, P., Diego, J. M., Broadhurst, T., De Martino, I., & Lazkoz, R. 2018, MNRAS, 478, 638

Jing, Y. P., & Suto, Y. 2002, ApJ, 574, 538

Jones, E., Oliphant, T., & Peterson, P. 2001, SciPy: Open Source Scientific Tools for Python

Kaiser, N. 1984, ApJ, 284, L9

Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, Jupyter Notebooks-a publishing format for reproducible computational workflows (IOS Press), 87

Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340

Köhlinger, F., Hoekstra, H., & Eriksen, M. 2015, MNRAS, 453, 3107

Komatsu, E., Dunkley, J., Nolta, M. R., et al. 2009, ApJS, 180, 330

Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18

Korytov, D., Hearin, A., Kovacs, E., et al. 2019, ApJS, 245, 26

Kovacs, E., Mao, Y.-Y., Aguena, M., et al. 2022, OJAp, 5, 1

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

Kugel, R., Schaye, J., Schaller, M., Forouhar Moreno, V. J., & McGibbon, R. J. 2025, MNRAS, 537, 2179

Lacasa, F., Lima, M., & Aguena, M. 2018, A&A, 611, A83

Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, ArXiv e-prints [arXiv:1110.3193]

Lee, B. E., Le Brun, A. M. C., Haq, M. E., et al. 2018, MNRAS, 479, 890

Lee, A., Wu, H. Y., Salcedo, A. N., et al. 2024, Phys. Rev. D, 111, 063502

Leizerovich, M., Landau, S. J., & Scóccola, C. G. 2024, Phys. Lett. B, 855, 138844

Lesci, G. F., Marulli, F., Moscardini, L., et al. 2022, A&A, 659, A88

Lewis, A. 2019, ArXiv e-prints [arXiv:1910.13970]

LSST Dark Energy Science Collaboration 2012, ArXiv e-prints [arXiv:1211.0310]

Mandelbaum, R., Hirata, C. M., Seljak, U., et al. 2005, MNRAS, 361, 1287

Mandelbaum, R., Lanusse, F., Leauthaud, A., et al. 2018, MNRAS, 481, 3170

Mantz, A. B. 2019, MNRAS, 485, 4863

Mantz, A., Allen, S. W., Ebeling, H., & Rapetti, D. 2008, MNRAS, 387, 1179

Mantz, A. B., von der Linden, A., Allen, S. W., et al. 2015, MNRAS, 446, 2205

McClintock, T., Varga, T. N., Gruen, D., et al. 2019, MNRAS, 482, 1352

Melchior, P., Gruen, D., McClintock, T., et al. 2017, MNRAS, 469, 4899
Mistele, T., & Durakovic, A. 2024, Open J. Astrophys., 7, 120

Mulroy, S. L., Farahi, A., Evrard, A. E., et al. 2019, MNRAS, 484, 60

Murata, R., Oguri, M., Nishimichi, T., et al. 2019, PASJ, 71, 107

Murray, C., Bartlett, J. G., Artis, E., & Melin, J.-B. 2022, MNRAS, 512, 4785

Myles, J., Gruen, D., Mantz, A. B., et al. 2021, MNRAS, 505, 33

Myles, J., Gruen, D., Jeltema, T., et al. 2025, ArXiv e-prints [arXiv:2506.06249]

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

Nord, B., Stanek, R., Rasia, E., & Evrard, A. E. 2008, MNRAS, 383, L10

Oguri, M., & Takada, M. 2011, Phys. Rev., D, 83

Oguri, M., Lee, J., & Suto, Y. 2003, ApJ, 599, 7

Oguri, M., Takada, M., Okabe, N., & Smith, G. P. 2010, MNRAS, 402, 307

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

Park, Y., Sunayama, T., Takada, M., et al. 2023, MNRAS, 518, 5171

Parroni, C., Mei, S., Erben, T., et al. 2017, ApJ, 848, 114

Payerne, C., Murray, C., Combet, C., et al. 2023, MNRAS, 520, 6223

Payerne, C., Murray, C., Combet, C., & Penna-Lima, M. 2024, MNRAS, 532, 381

Phriksee, A., Jullo, E., Limousin, M., et al. 2020, MNRAS, 491, 1643

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

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

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

Poisson, S. D. 1837, Bachelier, Imprimeur-Librairie

Prada, F., Klypin, A. A., Cuesta, A. J., Betancort-Rijo, J. E., & Primack, J. 2012, MNRAS, 423, 3018

Pranjal, R. S., Krause, E., Huang, H.-J., et al. 2023, MNRAS, 524, 3324

Prat, J., Zuntz, J., Chang, C., et al. 2023, Open J. Astrophys., 6, 13

Pratt, G. W., Arnaud, M., Biviano, A., et al. 2019, Space Sci. Rev., 215, 25

Raveri, M., & Doux, C. 2021, Phys. Rev. D, 104, 043504

Rozo, E., Evrard, A. E., Rykoff, E. S., & Bartlett, J. G. 2014, MNRAS, 438, 62

Rykoff, E. S., Rozo, E., Busha, M. T., et al. 2014, ApJ, 785, 104

Rykoff, E. S., Rozo, E., Hollowood, D., et al. 2016, ApJS, 224, 1

Sadibekova, T., Pierre, M., Clerc, N., et al. 2014, A&A, 571, A87

Sanderson, A. J. R., & Ponman, T. J. 2009, MNRAS, 402, 65

Saro, A., Bocquet, S., Rozo, E., et al. 2015, MNRAS, 454, 2305

Schmidt, S. J., Malz, A. I., Soo, J. Y. H., et al. 2020, MNRAS, 499, 1587

Schneider, P., Ehlers, J., & Falco, E. E. 1992, Gravitational Lenses

Schneider, M. D., Frenk, C. S., & Cole, S. 2012, JCAP, 2012, 030

Sereno, M., Fedeli, C., & Moscardini, L. 2016, JCAP, 2016, 042

Sheldon, E. S., & Huff, E. M. 2017, ApJ, 841, 24

Sheldon, E. S., Johnston, D. E., Frieman, J. A., et al. 2004, AJ, 127, 2544

Sheth, R. K., Mo, H. J., & Tormen, G. 2001, MNRAS, 323, 1

Shirasaki, M., & Takada, M. 2018, MNRAS, 478, 4277

Simet, M., McClintock, T., Mandelbaum, R., et al. 2017, MNRAS, 466, 3103

Sommer, M. W., Schrabback, T., Applegate, D. E., et al. 2022, MNRAS, 509, 1127

Spergel, D., Gehrels, N., Baltay, C., et al. 2015, ArXiv e-prints [arXiv:1503.03757]

Sunayama, T., Park, Y., Takada, M., et al. 2020, MNRAS, 496, 4468

Sunayama, T., Miyatake, H., Sugiyama, S., et al. 2024, Phys. Rev. D, 110, 083511

Tinker, J. L., Robertson, B. E., Kravtsov, A. V., et al. 2010, ApJ, 724, 878

Umetsu, K. 2020, A&ARv, 28

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22

Varga, T. N., DeRose, J., Gruen, D., et al. 2019, MNRAS, 489, 2511

Vogt, S. M. L., Bocquet, S., Davies, C. T., et al. 2024, Phys. Rev. D, 111, 043519

Wang, J., Bose, S., Frenk, C. S., et al. 2020, Nature, 585, 39

Wright, A. H., Hildebrandt, H., van den Busch, J. L., & Heymans, C. 2020, A&A, 637, A100

Wu, H.-Y., Weinberg, D. H., Salcedo, A. N., Wibking, B. D., & Zu, Y. 2019, MNRAS, 490, 2606

Wu, H.-Y., Costanzi, M., To, C.-H., et al. 2022, MNRAS, 515, 4471

Zhang, Y., Jeltema, T., Hollowood, D. L., et al. 2019, MNRAS, 487, 2578

Zhang, Z., Wu, H.-Y., Zhang, Y., et al. 2023, MNRAS, 523, 1994

Zhang, Z., Farahi, A., Nagai, D., et al. 2024, MNRAS, [arXiv:2310.18266]

Zhou, C., Wu, H.-Y., Salcedo, A. N., et al. 2024, Phys. Rev. D, 110, 103508

Zwicky, F. 1937, ApJ, 86, 217
## Appendix A: Stability of the fiducial cluster scaling relation

The fiducial scaling relation used throughout this paper was derived from the matched catalog between redMaPPer clusters and cosmoDC2 halos, restricted to $\lambda > 5$ and $M_{200c} > 10^{13} M_{\odot}$. In this appendix, we detail the various selection cuts applied to the matched catalog before inferring the fiducial relation using the likelihood defined in Eq. (18). For illustration, the posterior distributions of the mean scaling relation parameters, obtained using the matched catalog before applying any richness or mass cuts, are shown as green unfilled contours in Figure A.1.

Initially, we assume that the matched catalog is unaffected by the purity of the redMaPPer cluster finder, as false detections have been excluded after matching. However, the completeness of redMaPPer still impacts the catalog, particularly at the low-mass end where some halos are not detected. Ideally, the likelihood in Eq. (18) should account for such missing halos. To partially mitigate this effect, we apply a lower mass cut of $M_{200c} > 4 \times 10^{13} M_{\odot}$, a regime where redMaPPer's completeness exceeds 50%. The resulting posterior, reflecting this correction, is shown as magenta unfilled contours in Figure A.1.

Additionally, our cluster scaling relation model in Eq. (15) does not account for Poisson noise in redMaPPer richness determination. This noise can be incorporated by revising Eq. (15):

$$P(\lambda|m, z) = \int_{0}^{+\infty} \mathcal{P}(\lambda|\lambda_{\text{int}})P(\lambda_{\text{int}}|m, z)d\lambda_{\text{int}}, \quad (A.1)$$

where $\mathcal{P}(\lambda|\lambda_{\text{int}})$ is the Poisson distribution, and $P(\lambda_{\text{int}}|m, z)$ represents the probability density function for the intrinsic richness, which can be expressed with the same equation as in Eq. (15), but with intrinsic variance only. This updated $P(\lambda|m, z)$ can be approximated by using the same $P(\lambda_{\text{int}}|m, z)$ while updating the variance by

$$\sigma_{\ln(\lambda|m, z)}^2 \approx [\sigma_{\ln(\lambda|m, z)}]_{\text{int}}^2 + \frac{\exp(\langle\ln \lambda|m, z\rangle) - 1}{\exp(2\langle\ln \lambda|m, z\rangle)}, \quad (A.2)$$

where $[\sigma_{\ln(\lambda|m, z)}]_{\text{int}}^2$ is provided in Eq. (17), now only describing the intrinsic variance of richness. The impact of the Poisson noise, represented by the second term (see, e.g., Zhang et al. 2023), diminishes for large richness. Consequently, we apply a low richness cut $\lambda > 10$ to align with the Poisson-free distribution in Eq. (15).

Combining this with the mass cut $M_{200c} > 4 \times 10^{13} M_{\odot}$, the constraints are identified as our fiducial choice. We observe that the three fiducial constraints agree within $< 2\sigma$ of the baseline count+weak lensing posteriors. We also explore a stricter cut, $\lambda > 20$, combined with $M_{200c} > 4 \times 10^{13} M_{\odot}$. This $\lambda > 20$ scenario results in larger error bars (still compatible with count+weak lensing constraints) but biases the mass slope $\mu_m$ lower than other choices. In fact, from the left panel in Figure 1, the $\lambda > 20$ cut over-weights high-richness clusters in the fiducial fit, even after correcting the distribution $P(\lambda|m, z)$ using a truncated Gaussian to account for the minimum richness. While we adopted the black dashed contours as our fiducial constraints, further investigations are warranted.

## Appendix B: Impact of the radial fitting range on the scaling relation parameters

In this appendix we test the impact of the fitting radial range $[1, R_{\text{max}}]$ of the stacked excess surface density profiles on the

<figure><img src="image_14.png" /><figcaption>Fig. A.1. Posterior distribution of the mean scaling relation parameters from our count+lensing baselines (filled contours) and the matched redMaPPer-cosmoDC2 cluster-halo match catalog, for different richness and mass cuts.</figcaption></figure>

scaling relation parameters. Our baseline adopts a conservative choice of $R_{\text{max}} = 3.5$ Mpc for fitting the stacked excess surface density profiles, ensuring sensitivity to the one-halo regime only. At larger radii, the influence of the two-halo term is expected to become more significant, and modeling the stacked profiles using only the one-halo regime may be insufficient. From Oguri & Takada (2011), the two-halo excess surface density is given by

$$\Delta\Sigma_{2h}(R) = \frac{\rho_m(z)b_h(M_{200c}, z)}{(1+z)^3 D_A(z)^2} \int \frac{ldl}{(2\pi)} P_m(k_l, z) J_2(l\theta), \quad (B.1)$$

where $k_l = l/[D_A(z)(1+z)]$, $J_2$ is the second order Bessel function of the first kind, $\theta = R/D_A(z)$ is the separation angle, $\rho_m(z)$ is the matter density at redshift $z$, $b_h(M, z)$ is the halo bias at mass $M$ and redshift $z$ and $P_m(k, z)$ is the matter power spectrum at the halo redshift $z$. Finally, the total excess surface density profile is given by

$$\Delta\Sigma(R) = \Delta\Sigma_{1h}(R) + \Delta\Sigma_{2h}(R). \quad (B.2)$$

We repeated the scaling relation parameter fit using $R_{\text{max}} = 5.5$ and $R_{\text{max}} = 10$ Mpc, respectively. The mean scaling relation parameters for these different fitting radial ranges are shown in Figure B.1. Both new cuts yield conclusions consistent with the baseline, with shifts of less than $< 1\sigma$ compared to the baseline constraints.

Interestingly, the results for $R_{\text{max}} = 10$ Mpc are close to the baseline constraints, despite the absence of modeling for the two-halo term in the stacked lensing profile. We attribute this to the lower S/N of measurements at large scales, which reduces the impact of the two-halo term. However, with larger datasets, the two-halo term must be accurately modeled beyond $R = 4$ Mpc (Murata et al. 2019; McClintock et al. 2019; Melchior et al. 2017).
<figure><img src="image_15.png" /><figcaption>Fig. B.1. Posterior distribution of the mean scaling relation parameters from the combination of count and lensing profiles, for three different values of Rmax, used in the stacked lensing profile fitting range [1, Rmax].</figcaption></figure>

# Appendix C: Some other observational systematics

## C.1. Miscentering

When the inferred cluster center differs from the one of its underlying halo’s gravitational potential – which happens when deriving the cluster center from the intra-cluster galaxy distribution – the lensing profile modeling must account for a miscentering term. For a single miscentered cluster with offset radius Rmis, its radially averaged projected mass density is given by

$$ \Sigma_{\text{mis}}(R, R_{\text{mis}}) = \frac{1}{2\pi} \int_{0}^{2\pi} d\theta \Sigma \left( \sqrt{R^2 + R_{\text{mis}}^2 + 2RR_{\text{mis}}\cos\theta} \right). \quad (C.1) $$

For a stack of clusters, the above equation can be averaged over the distribution of miscentering P(Rmis), such as

$$ \Sigma_{\text{mis}}^{\text{stack}}(R) = \int_{0}^{+\infty} dR_{\text{mis}} P(R_{\text{mis}}) \Sigma_{\text{mis}}(R, R_{\text{mis}}). \quad (C.2) $$

From this equation, we can obtain ΔΣmis by injecting it into Eq. (10). The total one-halo term depends on the fraction fmis of miscentered clusters, as shown in, for example, Giocoli et al. (2021), and is given by

$$ \Delta\Sigma_{\text{1h}}(R) = f_{\text{mis}}\Delta\Sigma_{\text{1h,mis}}(R) + (1 - f_{\text{mis}})\Delta\Sigma_{\text{1h,cen}}(R), \quad (C.3) $$

where ΔΣ1h,cen(R) denotes the lensing profile of a perfectly centered stack of clusters. To test the impact of miscentering of redMaPPer clusters, we use the matched catalog between redMaPPer clusters and cosmoDC2 dark matter halos obtained in Sect. 3.4. In each cluster-halo pair, we consider the halo-member galaxy flagged as is_central==True in the cosmDC2 extragalactic catalog to trace the center of the halo’s gravitational potential, from which we compute the corresponding projected offset radius Rmis of the matched redMaPPer cluster. We found that the fraction of perfectly centered clusters (for

which Rmis = 0) is 85% so fmis = 0.15. For pairs with Rmis ≠ 0, the distribution of offsets between cosmoDC2 halo and redMaPPer cluster centers is well described by a Gamma distribution (McClintock et al. 2019) given by

$$ P(R_{\text{mis}}) = \frac{R_{\text{mis}}}{\sigma_{\text{mis}}^2} \exp\{-\frac{R_{\text{mis}}}{\sigma_{\text{mis}}}\}, \quad (C.4) $$

with σmis = 0.12 Mpc. We use the cluster_toolkit package<fnref n="45" /> to compute the miscentered profile evoked in Sect. 2 over a wide range of mass and redshift, and we found that the bias compared to a perfectly centered profile is at most 1% at R = 1 Mpc, and 0.25% at R = 2 Mpc. For simplicity, since the computation of the miscentering term is computationally demanding, we chose to neglect this contribution.

## C.2. Contamination

The background source selection may cause contamination of the source sample by cluster member galaxies (Varga et al. 2019), which dilutes the signal in the innermost region. The observed contaminated lensing profile accounting for this dilution is given by

$$ \Delta\Sigma_{ij}^{\text{cont}}(R) = [1 - f_{\text{cl}}(R)]\Delta\Sigma_{ij}(R), \quad (C.5) $$

where ΔΣij(R) is given in Eq. (13), fcl(R) is the fraction of "source-selected" member galaxies at a distance R from the cluster center, and is a decreasing function of R. At fixed R, the fraction fcl increases with mass, since massive clusters have more member galaxies, and thus are more subject to this effect (Varga et al. 2019).

For each redMaPPer cluster, we consider the corresponding matched dark matter halo in the matched catalog presented in Sect. 3. We applied the photo-z source selection and then identified the remaining halo "member galaxies," as labeled in the simulation (which are randomly drawn at the HOD level in the DC2 workflow to populate halos, see Korytov et al. 2019). This procedure enables us to recover fcl(R) for a variety of masses and redshifts, and we found that the correcting factor [1 − fcl(R = 1 Mpc)] ~ 1 after combining the two cuts, and for both BPZ and FlexZBoost algorithms. We find that our source selection enables us to discriminate efficiently between true sources and member galaxies (the effect is less than 0.1%).

# Appendix D: Derivation of the selection bias

In this appendix we use the notation

$$ N(m, z) = \frac{dn(m, z)}{dm} \frac{d^2V(z)}{dzd\Omega}. \quad (D.1) $$

We can consider that the measured shear and richness are correlated, such that they can be modeled as correlated random variables via P(ΔΣ, λ|m). Then, the average excess surface density profile within the i-th redshift bin and the j-th richness bin is given by

$$ \langle \Delta\Sigma \rangle \propto \int_{z_i}^{z_{i+1}} dz \int_{m_{\text{lim}}}^{m_{\text{max}}} dm N(m, z) \int_{\lambda_j}^{\lambda_{j+1}} d\lambda P(\lambda|m, z) \langle \Delta\Sigma | m, \lambda, z \rangle. \quad (D.2) $$

<footnote n="45">⁴⁵ https://github.com/tmcclintock/cluster_toolkit</footnote>
From Wu et al. (2022), we consider $P(\Delta\Sigma, \lambda|m)$ as a multivariate Gaussian distribution, with means $\langle\Delta\Sigma|m, z\rangle$ and $\langle\ln \lambda|m, z\rangle$, and with covariance matrix

$$C = \begin{pmatrix} \sigma_{\ln \lambda}^2 & r \sigma_{\ln \lambda} \sigma_{\Delta \Sigma} \\ r \sigma_{\ln \lambda} \sigma_{\Delta \Sigma} & \sigma_{\Delta \Sigma}^2 \end{pmatrix}. \qquad (D.3)$$

From Wu et al. (2019) (see their Eq. 8), the conditional mean of $\Delta\Sigma$ is given by

$$\langle\Delta\Sigma|m, \lambda, z\rangle = \langle\Delta\Sigma|m, z\rangle + r \frac{\sigma_{\Delta\Sigma}}{\sigma_{\ln\lambda}} [\ln \lambda - \langle\ln \lambda|m, z\rangle]. \qquad (D.4)$$

From Farahi et al. (2022), we can approximate the halo mass function to the form

$$N(m, z) \approx A(z) \exp[-\beta_1 \ln m]. \qquad (D.5)$$

By using the simplified mean mass-richness relation,

$$\langle \ln \lambda | m, z \rangle = \pi(z) + \bar{\mu}_m \ln m, \qquad (D.6)$$

where $\bar{\mu}_m$ is linked to the $\mu_m$ in Eq. (16) by $\bar{\mu}_m = \mu_m / \ln(10)$, and the results in Evrard et al. (2014) (see their Eq. 3 and 4), we have

$$\ln \lambda - \langle \ln \lambda | m, z \rangle = \frac{\sigma_{\ln \lambda}^2}{\bar{\mu}_m} \left( \frac{\langle \ln M | \lambda, z \rangle}{\sigma_{\ln M}^2} + \beta_1 \right) - \bar{\mu}_m \ln M. \qquad (D.7)$$

Using that $\sigma_{\ln \lambda} = \bar{\mu}_m \sigma_{\ln M}$ (Evrard et al. 2014), we have

$$\ln \lambda - \langle \ln \lambda | m, z \rangle = \frac{\sigma_{\ln \lambda}^2}{\bar{\mu}_m} \beta_1 + \bar{\mu}_m \langle \ln M | \lambda, z \rangle - \bar{\mu}_m \ln M. \qquad (D.8)$$

So, by combining Eq. (D.4) and Eq. (D.8), we get

$$\langle\Delta\Sigma|m, \lambda, z\rangle = \langle\Delta\Sigma|m, z\rangle + \frac{\beta_1}{\bar{\mu}_m} \mathrm{Cov}(\Delta\Sigma, \ln \lambda|m, z) \qquad (D.9)$$

$$+ \frac{\bar{\mu}_m}{\sigma_{\ln \lambda}^2} \mathrm{Cov}(\Delta\Sigma, \ln \lambda|m, z) (\langle \ln M | \lambda, z \rangle - \ln M). \qquad (D.10)$$

From above, we find that the stacked shear decomposes as $\langle\Delta\Sigma\rangle = \langle\Delta\Sigma\rangle_1 + \langle\Delta\Sigma\rangle_2 + \langle\Delta\Sigma\rangle_3$. The first term is given by

$$\langle\Delta\Sigma\rangle_1 = \int_{z_i}^{z_{i+1}} dz \int_{m_{\mathrm{lim}}}^{m_{\mathrm{max}}} dm N(m, z) \int_{\lambda_j}^{\lambda_{j+1}} d\lambda P(\lambda|m, z) \langle\Delta\Sigma|m, z\rangle. \qquad (D.11)$$

The second term is given by

$$\langle\Delta\Sigma\rangle_2 = \frac{\beta_1}{\bar{\mu}_m} \int_{z_i}^{z_{i+1}} dz \int_{m_{\mathrm{lim}}}^{m_{\mathrm{max}}} dm N(m, z) \int_{\lambda_j}^{\lambda_{j+1}} d\lambda P(\lambda|m, z) \mathrm{Cov}(\Delta\Sigma, \ln \lambda|m, z) \qquad (D.12)$$

$$\approx \frac{\beta_1}{\bar{\mu}_m} \int_{z_i}^{z_{i+1}} dz \int_{\lambda_j}^{\lambda_{j+1}} d\lambda P(\lambda, z) \mathrm{Cov}(\Delta\Sigma, \ln \lambda |m, z), \qquad (D.13)$$

$$\approx \frac{\beta_1}{\bar{\mu}_m} \langle \mathrm{Cov}(\Delta\Sigma, \ln \lambda|m, z) \rangle_{ij} \qquad (D.14)$$

where the second equation is obtained by considering that $\mathrm{Cov}(\Delta\Sigma, \ln \lambda|m, z)$ evolves linearly with mass $m$. The term $\langle\Delta\Sigma\rangle_3$ is negligible and not considered in this analysis (see, e.g., Wu et al. 2022; Zhang et al. 2024). Then, from the above equation, we see that a possible shear-richness correlation (e.g., $r \neq 0$, or equivalently $\mathrm{Cov}(\Delta\Sigma, \ln \lambda|m, z) \neq 0$), arising from selection bias in cluster finder algorithms, may lead to a misinterpretation of the stacked cluster profile as a naive average as given in Eq. (D.11).
## Appendix E: Constraints of the scaling relation: Summary table

We list in Table E.1 the values of the scaling relation parameters obtained from cluster count and/or cluster lensing in Sect. 5.1 (baseline analysis), Sect. 5.2.1 (impact of $c(M)$ relation), Sect. 5.2.2 (impact of density profile), Sect. 5.3.1 (impact of source galaxy photometric redshifts), and Sect. 5.3.2 (impact of shear-richness covariance). The Figure E.1 displays the tension metric for the different study cases with respect to (i) the fiducial constraints and (ii) the count+weak lensing profiles analysis, labeled as the "baseline" in Sect. 5.1.

<table><caption>Table E.1. Best-fit parameters of the cluster scaling relation in Eq. (15).</caption><thead><tr><th>Parameters</th><th>ln λ<sub>0</sub></th><th>μ<sub>z</sub></th><th>μ<sub>m</sub></th><th>σ<sub>lnλ<sub>0</sub></sub></th><th>σ<sub>z</sub></th><th>σ<sub>m</sub></th></tr></thead><tbody><tr><td colspan="7" style="text-align: center; font-weight: bold;">Baseline analysis, Sect. 5.1</td></tr><tr><td>Abundance alone</td><td>3.27 ± 0.14</td><td>0.04 ± 0.19</td><td>2.22 ± 0.20</td><td>0.53 ± 0.11</td><td>0.06 ± 0.25</td><td>-0.11 ± 0.14</td></tr><tr><td>Lensing masses alone</td><td>3.37 ± 0.20</td><td>0.18 ± 0.47</td><td>2.34 ± 0.44</td><td>0.52 ± 0.10</td><td>-0.03 ± 0.32</td><td>-0.00 ± 0.15</td></tr><tr><td>Lensing profiles alone</td><td>3.31 ± 0.19</td><td>0.36 ± 0.54</td><td>2.41 ± 0.45</td><td>0.54 ± 0.08</td><td>-0.05 ± 0.33</td><td>-0.03 ± 0.15</td></tr><tr><td>Lensing masses + Abundance</td><td>3.33 ± 0.04</td><td>0.14 ± 0.07</td><td>2.14 ± 0.06</td><td>0.58 ± 0.02</td><td>0.06 ± 0.11</td><td>0.13 ± 0.07</td></tr><tr><td>Lensing profiles + Abundance</td><td>3.37 ± 0.03</td><td>0.08 ± 0.07</td><td>2.18 ± 0.07</td><td>0.53 ± 0.03</td><td>0.20 ± 0.11</td><td>0.14 ± 0.05</td></tr><tr><td>Lensing profiles + Abundance (z &lt; 0.8)</td><td>3.34 ± 0.04</td><td>-0.04 ± 0.13</td><td>2.16 ± 0.07</td><td>0.54 ± 0.02</td><td>0.17 ± 0.15</td><td>0.08 ± 0.08</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">Impact of the concentration-mass relation, Sect. 5.2.1</td></tr><tr><td>Duffy et al. 2008</td><td>3.33 ± 0.04</td><td>0.14 ± 0.07</td><td>2.14 ± 0.06</td><td>0.58 ± 0.02</td><td>0.06 ± 0.11</td><td>0.13 ± 0.07</td></tr><tr><td>Prada et al. 2012</td><td>3.40 ± 0.04</td><td>0.17 ± 0.06</td><td>2.22 ± 0.07</td><td>0.52 ± 0.03</td><td>0.02 ± 0.10</td><td>0.16 ± 0.05</td></tr><tr><td>Bhattacharya et al. 2013</td><td>3.37 ± 0.04</td><td>0.14 ± 0.08</td><td>2.19 ± 0.06</td><td>0.55 ± 0.03</td><td>0.07 ± 0.12</td><td>0.14 ± 0.05</td></tr><tr><td>Diemer Kravtsov 2014</td><td>3.38 ± 0.03</td><td>0.15 ± 0.07</td><td>2.20 ± 0.07</td><td>0.54 ± 0.03</td><td>0.06 ± 0.12</td><td>0.15 ± 0.05</td></tr><tr><td>free concentration</td><td>3.37 ± 0.04</td><td>0.14 ± 0.08</td><td>2.19 ± 0.08</td><td>0.54 ± 0.03</td><td>0.08 ± 0.14</td><td>0.15 ± 0.06</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">Impact of the dark matter density profile, Sect. 5.2.2</td></tr><tr><td>Navarro-Frenk-White (NFW)</td><td>3.37 ± 0.04</td><td>0.14 ± 0.08</td><td>2.19 ± 0.08</td><td>0.54 ± 0.03</td><td>0.08 ± 0.14</td><td>0.15 ± 0.06</td></tr><tr><td>Einasto</td><td>3.36 ± 0.04</td><td>0.13 ± 0.08</td><td>2.18 ± 0.07</td><td>0.55 ± 0.03</td><td>0.08 ± 0.14</td><td>0.14 ± 0.06</td></tr><tr><td>Hernquist</td><td>3.36 ± 0.05</td><td>0.14 ± 0.08</td><td>2.17 ± 0.08</td><td>0.55 ± 0.03</td><td>0.08 ± 0.13</td><td>0.13 ± 0.06</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">Impact of source photometric redshifts, Sect. 5.3.1</td></tr><tr><td>true source redshifts</td><td>3.37 ± 0.03</td><td>0.08 ± 0.07</td><td>2.18 ± 0.07</td><td>0.53 ± 0.03</td><td>0.20 ± 0.11</td><td>0.14 ± 0.05</td></tr><tr><td>BPZ</td><td>3.29 ± 0.03</td><td>-0.02 ± 0.08</td><td>2.07 ± 0.06</td><td>0.59 ± 0.02</td><td>0.32 ± 0.11</td><td>0.14 ± 0.06</td></tr><tr><td>BPZ (with factor 1+b)</td><td>3.33 ± 0.07</td><td>-0.04 ± 0.09</td><td>2.12 ± 0.11</td><td>0.55 ± 0.06</td><td>0.40 ± 0.15</td><td>0.13 ± 0.06</td></tr><tr><td>FlexZBoost</td><td>3.37 ± 0.03</td><td>0.09 ± 0.08</td><td>2.18 ± 0.06</td><td>0.54 ± 0.02</td><td>0.17 ± 0.12</td><td>0.14 ± 0.05</td></tr><tr><td>FlexZBoost (with factor 1+b)</td><td>3.33 ± 0.07</td><td>0.08 ± 0.07</td><td>2.12 ± 0.10</td><td>0.57 ± 0.06</td><td>0.17 ± 0.13</td><td>0.14 ± 0.06</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">Impact of shear-richness covariance, Sect. 5.3.2</td></tr><tr><td>Without Cov(ΔΣ, λ) correction</td><td>3.37 ± 0.03</td><td>0.08 ± 0.07</td><td>2.18 ± 0.07</td><td>0.53 ± 0.03</td><td>0.20 ± 0.12</td><td>0.14 ± 0.04</td></tr><tr><td>With Cov(ΔΣ, λ) correction</td><td>3.36 ± 0.03</td><td>0.09 ± 0.07</td><td>2.17 ± 0.06</td><td>0.54 ± 0.02</td><td>0.17 ± 0.12</td><td>0.14 ± 0.05</td></tr><tr><td>With Cov(ΔΣ, λ) correction - BPZ</td><td>3.28 ± 0.03</td><td>-0.02 ± 0.07</td><td>2.05 ± 0.07</td><td>0.61 ± 0.02</td><td>0.30 ± 0.10</td><td>0.14 ± 0.06</td></tr><tr><td>With Cov(ΔΣ, λ) correction - FlexZBoost</td><td>3.36 ± 0.03</td><td>0.09 ± 0.08</td><td>2.17 ± 0.06</td><td>0.55 ± 0.02</td><td>0.16 ± 0.13</td><td>0.15 ± 0.05</td></tr><tr><td>With Cov(ΔΣ, λ) correction - BPZ (1+b)</td><td>3.34 ± 0.07</td><td>-0.01 ± 0.08</td><td>2.13 ± 0.10</td><td>0.55 ± 0.06</td><td>0.36 ± 0.14</td><td>0.14 ± 0.06</td></tr><tr><td>With Cov(ΔΣ, λ) correction - FlexZBoost (1+b)</td><td>3.34 ± 0.07</td><td>0.10 ± 0.08</td><td>2.16 ± 0.09</td><td>0.56 ± 0.06</td><td>0.13 ± 0.14</td><td>0.13 ± 0.06</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>Fig. E.1. Tension in numbers of $\sigma$ between the different count and/or lensing analyses and the fiducial constraints (circle dots) for the parameters listed in Figure 13, and between the different analyses and the baseline weak lensing and count analysis (star dots).</figcaption></figure>