<figure><img src="image_1.png" /></figure>

HAL
open science

# Impact of property covariance on cluster weak lensing scaling relations

Zhuowen Zhang, Arya Farahi, Daisuke Nagai, Erwin T Lau, Joshua Frieman, Marina Ricci, Anja von der Linden, Hao-Yi Wu, Anja von der Linden, Tamas Varga, et al.

<figure><img src="image_2.png" /></figure>

## ► To cite this version:

Zhuowen Zhang, Arya Farahi, Daisuke Nagai, Erwin T Lau, Joshua Frieman, et al. Impact of property covariance on cluster weak lensing scaling relations. Monthly Notices of the Royal Astronomical Society, 2024, 530 (3), pp.3127 - 3149. 10.1093/mnras/stae999. hal-04587419

HAL Id: hal-04587419

https://u-paris.hal.science/hal-04587419

Submitted on 24 May 2024

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

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Impact of property covariance on cluster weak lensing scaling relations

Zhuowen Zhang, Arya Farahi, Daisuke Nagai, Erwin T Lau, Joshua Frieman, Marina Ricci, Anja von der Linden, Hao-Yi Wu, Anja von der Linden, Tamas Varga, et al.

► To cite this version:

Zhuowen Zhang, Arya Farahi, Daisuke Nagai, Erwin T Lau, Joshua Frieman, et al.. Impact of property covariance on cluster weak lensing scaling relations. Monthly Notices of the Royal Astronomical Society, 2024, 530 (3), pp.3127 - 3149. 10.1093/mnras/stae999. hal-04587419

HAL Id: hal-04587419

https://u-paris.hal.science/hal-04587419v1

Submitted on 24 May 2024

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Impact of property covariance on cluster weak lensing scaling relations

Zhuowen Zhang<sup>1,1*</sup> Arya Farahi,<sup>2</sup> Daisuke Nagai<sup>1,3</sup> Erwin T. Lau,<sup>4</sup> Joshua Frieman,<sup>1</sup> Marina Ricci,<sup>5</sup> Anja von der Linden,<sup>6</sup> Hao-Yi Wu<sup>1,7</sup> and the LSST Dark Energy Science Collaboration

Accepted 2024 March 14. Received 2024 March 14; in original form 2023 October 27

## ABSTRACT

We present an investigation into a hitherto unexplored systematic that affects the accuracy of galaxy cluster mass estimates with weak gravitational lensing. Specifically, we study the covariance between the weak lensing signal, $\Delta\Sigma$, and the ‘true’ cluster galaxy number count, $N_{gal}$, as measured within a spherical volume that is void of projection effects. By quantifying the impact of this covariance on mass calibration, this work reveals a significant source of systematic uncertainty. Using the MDPL2 simulation with galaxies traced by the SAGE semi-analytic model, we measure the intrinsic property covariance between these observables within the three-dimensional vicinity of the cluster, spanning a range of dynamical mass and redshift values relevant for optical cluster surveys. Our results reveal a negative covariance at small radial scales ($R \lesssim R_{200c}$) and a null covariance at large scales ($R \gtrsim R_{200c}$) across most mass and redshift bins. We also find that this covariance results in a 2–3 per cent bias in the halo mass estimates in most bins. Furthermore, by modelling $N_{gal}$ and $\Delta\Sigma$ as multi-(log)-linear equations of secondary halo properties, we provide a quantitative explanation for the physical origin of the negative covariance at small scales. Specifically, we demonstrate that the $N_{gal}-\Delta\Sigma$ covariance can be explained by the secondary properties of haloes that probe their formation history. We attribute the difference between our results and the positive bias seen in other works with (mock)-cluster finders to projection effects. These findings highlight the importance of accounting for the covariance between observables in cluster mass estimation, which is crucial for obtaining accurate constraints on cosmological parameters.

**Key words:** gravitational lensing: weak – galaxies: clusters: general – cosmology: observations.

## 1 INTRODUCTION

Cluster abundance and its evolution with redshift are linked to the constituents of the Universe through the growth of cosmic structure (Allen, Evrard & Mantz 2011 for a review). Cluster abundance measured in large-scale galaxy surveys offers power constraints on cosmological parameters (e.g. Vikhlinin et al. 2009; Mantz et al. 2015; Dark Energy Survey Collaboration 2016; de Haan et al. 2016; Mantz et al. 2016a; Pierre et al. 2016; Costanzi et al. 2021). These constraints are based on accurate cluster mass measurements, which are not directly observable and must be inferred. Cluster mass calibration has been identified as one of the leading systematic uncertainties in cosmological constraints using galaxy cluster abundance (Mantz et al. 2010; Rozo et al. 2010; Applegate et al. 2014; von der Linden et al. 2014; Dodelson et al. 2016; Murata et al. 2019; Costanzi et al. 2021). Accurate mappings between a population of massive clusters and their observables are thus critical and essential in cluster

cosmology. Considerable effort has been put into measuring the statistical relationships between masses and observable properties that reflect their baryon contents (see Giodini et al. 2013 for a review) and quantifying the sources of uncertainties.

The Dark Energy Survey (DES) cluster cosmology from the Year 1 data set (Abbott et al. 2020) reported tension in $\Omega_m$ – the mean matter density of the universe – with the DES 3 x 2 pt probe that utilizes three two-point functions from the DES galaxy survey (Abbott et al. 2018). The tension between these two probes that utilize the same underlying data set may be attributed to systematics that bias the weak lensing mass of clusters low at the low-mass end (Costanzi et al. 2021; To et al. 2021). A possible origin for this discrepancy is that cluster masses are biased low due to systematics in cluster mass calibration. On the other hand, the tension can also originate from new physics that extends the Standard Cosmological model. Thus, it is important to understand the systematics of cluster mass calibration. Cluster masses estimated from X-ray and SZ data are known to suffer from hydrostatic bias (Pratt et al. 2019). Conversely, cluster masses estimated from weak lensing have the potential to be more accurate compared to X-ray and SZ cluster masses. The

<footnote n="1"><sup>1</sup>Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA</footnote>

<footnote n="2"><sup>2</sup>Departments of Statistics and Data Science, University of Texas at Austin, Austin, TX 78757, USA</footnote>

<footnote n="3"><sup>3</sup>Department of Physics, Yale University, New Haven, CT 06520, USA</footnote>

<footnote n="4"><sup>4</sup>Center for Astrophysics | Harvard & Smithsonian, Cambridge, MA, 02138, USA</footnote>

<footnote n="5"><sup>5</sup>Université Paris Cité, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France</footnote>

<footnote n="6"><sup>6</sup>Department of Physics and Astronomy, Stony Brooks University, Stony Brook, NY, 11794, USA</footnote>

<footnote n="7"><sup>7</sup>Department of Physics, Boise State University, Boise, ID 83725, USA</footnote>

<footnote>*E-mail: zzhang13@uchicago.edu</footnote>
# Impact of property covariance on cluster weak lensing scaling relations

Zhuowen Zhang<sup>1,1*</sup> Arya Farahi,<sup>2</sup> Daisuke Nagai<sup>1,3</sup> Erwin T. Lau,<sup>4</sup> Joshua Frieman,<sup>1</sup> Marina Ricci,<sup>5</sup> Anja von der Linden,<sup>6</sup> Hao-Yi Wu<sup>1,7</sup> and the LSST Dark Energy Science Collaboration

<sup>1</sup>Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA

<sup>2</sup>Departments of Statistics and Data Science, University of Texas at Austin, Austin, TX 78757, USA

<sup>3</sup>Department of Physics, Yale University, New Haven, CT 06520, USA

<sup>4</sup>Center for Astrophysics | Harvard & Smithsonian, Cambridge, MA, 02138, USA

<sup>5</sup>Université Paris Cité, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France

<sup>6</sup>Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY, 11794, USA

<sup>7</sup>Department of Physics, Boise State University, Boise, ID 83725, USA

Accepted 2024 March 14. Received 2024 March 14; in original form 2023 October 27

## ABSTRACT

We present an investigation into a hitherto unexplored systematic that affects the accuracy of galaxy cluster mass estimates with weak gravitational lensing. Specifically, we study the covariance between the weak lensing signal, $\Delta\Sigma$, and the ‘true’ cluster galaxy number count, $N_{gal}$, as measured within a spherical volume that is void of projection effects. By quantifying the impact of this covariance on mass calibration, this work reveals a significant source of systematic uncertainty. Using the MDPL2 simulation with galaxies traced by the SAGE semi-analytic model, we measure the intrinsic property covariance between these observables within the three-dimensional vicinity of the cluster, spanning a range of dynamical mass and redshift values relevant for optical cluster surveys. Our results reveal a negative covariance at small radial scales ($R \lesssim R_{200c}$) and a null covariance at large scales ($R \gtrsim R_{200c}$) across most mass and redshift bins. We also find that this covariance results in a 2–3 per cent bias in the halo mass estimates in most bins. Furthermore, by modelling $N_{gal}$ and $\Delta\Sigma$ as multi-(log)-linear equations of secondary halo properties, we provide a quantitative explanation for the physical origin of the negative covariance at small scales. Specifically, we demonstrate that the $N_{gal}-\Delta\Sigma$ covariance can be explained by the secondary properties of haloes that probe their formation history. We attribute the difference between our results and the positive bias seen in other works with (mock)-cluster finders to projection effects. These findings highlight the importance of accounting for the covariance between observables in cluster mass estimation, which is crucial for obtaining accurate constraints on cosmological parameters.

**Key words:** gravitational lensing: weak – galaxies: clusters: general – cosmology: observations.

## 1 INTRODUCTION

Cluster abundance and its evolution with redshift are linked to the constituents of the Universe through the growth of cosmic structure (Allen, Evrard & Mantz 2011 for a review). Cluster abundance measured in large-scale galaxy surveys offers power constraints on cosmological parameters (e.g. Vikhlinin et al. 2009; Mantz et al. 2015; Dark Energy Survey Collaboration 2016; de Haan et al. 2016; Mantz et al. 2016a; Pierre et al. 2016; Costanzi et al. 2021). These constraints are based on accurate cluster mass measurements, which are not directly observable and must be inferred. Cluster mass calibration has been identified as one of the leading systematic uncertainties in cosmological constraints using galaxy cluster abundance (Mantz et al. 2010; Rozo et al. 2010; Applegate et al. 2014; von der Linden et al. 2014; Dodelson et al. 2016; Murata et al. 2019; Costanzi et al. 2021). Accurate mappings between a population of massive clusters and their observables are thus critical and essential in cluster

cosmology. Considerable effort has been put into measuring the statistical relationships between masses and observable properties that reflect their baryon contents (see Giodini et al. 2013 for a review) and quantifying the sources of uncertainties.

The Dark Energy Survey (DES) cluster cosmology from the Year 1 data set (Abbott et al. 2020) reported tension in $\Omega_m$ – the mean matter density of the universe – with the DES 3 x 2 pt probe that utilizes three two-point functions from the DES galaxy survey (Abbott et al. 2018). The tension between these two probes that utilize the same underlying data set may be attributed to systematics that bias the weak lensing mass of clusters low at the low-mass end (Costanzi et al. 2021; To et al. 2021). A possible origin for this discrepancy is that cluster masses are biased low due to systematics in cluster mass calibration. On the other hand, the tension can also originate from new physics that extends the Standard Cosmological model. Thus, it is important to understand the systematics of cluster mass calibration. Cluster masses estimated from X-ray and SZ data are known to suffer from hydrostatic bias (Pratt et al. 2019). Conversely, cluster masses estimated from weak lensing have the potential to be more accurate compared to X-ray and SZ cluster masses. The

<footnote>*E-mail: zzhang13@uchicago.edu</footnote>
systematics in the weak lensing mass calibration has just started to be explored recently (Applegate et al. 2014; Schrabback et al. 2018; McClintock et al. 2019; Kiiveri et al. 2021; Wu et al. 2022).

A relatively unexplored category of cluster systematics is the covariance between different cluster properties, including cluster observables and mass proxies. In cluster mass calibration, it is often assumed that this property covariance is negligible. However, as initially pointed out by Nord et al. (2008) and later shown in Evrard et al. (2014) and Farahi et al. (2018), non-zero property covariances between cluster observables can induce non-trivial, additive bias in cluster mass. As property covariance is additive, the systematic uncertainties that it induces will not be mitigated with the reduction of statistical errors as the sample size of the cluster increases. To achieve accurate cosmological constraints with the next generation of large-scale cluster surveys, it is imperative that systematic uncertainties in the property covariance be accurately and precisely quantified (Rozo et al. 2014).

Although the property covariance linking mass to observable properties is becoming better understood and measured (Wu et al. 2015; Mantz et al. 2016a; Farahi et al. 2018, 2019; Sereno et al. 2020), studies that specifically investigate weak lensing property covariance are scarce, which poses a challenge for upcoming lensing surveys of galaxy clusters such as the Rubin (Ivezić et al. 2019) observatories. To achieve the percentage-level lensing mass calibration goals for the upcoming observations, the property covariance of weak lensing must be quantified.

The physical origins of property covariance in lensing signals of galaxy clusters can be attributed to the halo formation history of the cluster and baryonic physics (Xhakaj et al. 2022). Developing a first-principle physical model for the property covariance as a function of halo formation history and baryonic physics is a daunting task due to the highly non-linear and multiscale physics involved in cluster formation. To make progress, in this paper, we adopt a simulation-based, data-driven approach whereby we develop semi-analytical parametric models of property covariance, which we then calibrate with cosmological simulations. We then apply our model to quantify the bias induced due to a non-zero property covariance in the expected weak lensing signal and the mass-observable scaling relation.

As will be presented in Section 3, a key element of this analysis is the estimation of true cluster richness by encircling clusters within a three-dimensional (3D) radius within the physical vicinity of the halo centre, as opposed to a two-dimensional (2D) projected radius used by cluster finders as redMaPPer (Rykoff et al. 2014) by identifying galaxies within the red-sequence band in the colour–magnitude space – the major difference being the removal of projection effects, or the mis-identification of non-cluster galaxies in the 2D projected radius from the photometric redshift uncertainty of the red-sequence when estimating the true richness from a gravitationally bound region around the halo. Furthermore, as this simulation-based study does not introduce other observational systematics as shape noise of galaxies, point spread function, miscentring, among others, this study can be used to explore the intrinsic covariance between observables prior to the addition of extrinsic systematics as projection effects. Our results will not only provide insight into the physical origin of the covariance, the difference between the total covariance as measured by observations and the intrinsic covariance will provide estimates on the amplitude of the extrinsic component.

The goals of this work are to (i) develop an analytical model that accounts for and quantifies the effect of non-zero covariance on cluster mass calibration, (ii) quantify this property covariance utilizing cosmological simulations, (iii) update uncertainties on inferred cluster mass estimates, and (iv) explain the physical origin

of the covariance using secondary halo parameters. The rest of this paper is organized as follows. In Section 2, we present a population-based analytical framework. In Section 3, we describe the simulations and data-vector employed in this work. In Section 4, we present our measurements and the covariance model. In Section 5, we present the impact of the covariance on weak lensing mass calibration. In Section 6, we quantify the physical origin of the covariance by parametrizing it using secondary halo parameters. In Section 7, we compare our work with those that employ realistic cluster finders. We conclude in Section 8.

## 2 THEORETICAL FRAMEWORK

This section presents a theoretical framework that examines the impact of covariance on mass-observable scaling relations. In Section 2.1, we introduce the definitions of richness and weak lensing excess surface mass density and their scaling relations with cluster mass. We then describe the model of property covariance of richness and excess surface mass density in Section 2.2. In Section 2.3, we model the impact of covariance on stacked observable scaling relations. Finally, in Section 2.4, we develop a theoretical framework that explains the covariance based on a set of secondary halo parameters. A graphic representation of the outline of the paper is shown in Fig. 1. The notations used in this section to describe the covariance are listed in Table 1 and notations for scaling relations are listed in Table 2.

### 2.1 Observable-mass relations

#### 2.1.1 Excess surface mass density $\Delta\Sigma$ from weak lensing

In weak lensing measurements of galaxy clusters, the key observable is the excess surface mass density, denoted $\Delta\Sigma$. The excess surface mass density is defined as

$$ \Delta\Sigma(M, z, r_p) = \bar{\Sigma}(M, z, < r_p) - \Sigma(M, z, r_p), \quad (1) $$

where $\bar{\Sigma}(M, z, < r_p)$ denotes the average surface mass density within projected radius $r_p$, and $\Sigma(M, z, r_p)$ represents the average of the surface mass density at $r_p$. We model the average surface mass density $\Sigma$ as

$$ \Sigma(r_p) = \rho_m \int_{-\infty}^{+\infty} \left( 1 + \xi_{hm} \left( r = \sqrt{r_p^2 + \chi^2} \right) \right) d\chi, \quad (2) $$

where $\rho_m$ is the mean matter density at the redshift of the cluster, $R$ is the projected radius in the plane of the sky, $\chi$ is the comoving distance along the line of sight centred around the cluster, and $\xi_{hm}(r)$ is the halo-matter correlation function that characterizes the total mass density within a halo. Under the halo model, the halo-matter correlation function consists of a ‘one-halo’ term:

$$ \xi_{1h}(r|M) = \frac{\rho_{\text{NFW}}(r|M)}{\rho_{m0}} - 1, \quad (3) $$

and a ‘two-halo’ term:

$$ \xi_{2h}(r|M) = b(M)\xi_{\text{lin}}(r), \quad (4) $$

where $\rho_{\text{NFW}}$ is the Navarro–Frenk–White (NFW) density profile (Navarro, Frenk & White 1997), and $\xi_{\text{lin}}$ is the linear matter correlation function, and $b$ is the halo bias parameter.

In weak lensing, the excess surface density $\Delta\Sigma$ is tied to the tangential shear $\gamma_t$ of the galaxies relative to the centre of each foreground halo by the relation

$$ \Sigma_{\text{crit}}\gamma_t = \bar{\Sigma}(<R) - \Sigma(R) \equiv \Delta\Sigma(R), \quad (5) $$
systematics in the weak lensing mass calibration has just started to be explored recently (Applegate et al. 2014; Schrabback et al. 2018; McClintock et al. 2019; Kiiveri et al. 2021; Wu et al. 2022).

A relatively unexplored category of cluster systematics is the covariance between different cluster properties, including cluster observables and mass proxies. In cluster mass calibration, it is often assumed that this property covariance is negligible. However, as initially pointed out by Nord et al. (2008) and later shown in Evrard et al. (2014) and Farahi et al. (2018), non-zero property covariances between cluster observables can induce non-trivial, additive bias in cluster mass. As property covariance is additive, the systematic uncertainties that it induces will not be mitigated with the reduction of statistical errors as the sample size of the cluster increases. To achieve accurate cosmological constraints with the next generation of large-scale cluster surveys, it is imperative that systematic uncertainties in the property covariance be accurately and precisely quantified (Rozo et al. 2014).

Although the property covariance linking mass to observable properties is becoming better understood and measured (Wu et al. 2015; Mantz et al. 2016a; Farahi et al. 2018, 2019; Sereno et al. 2020), studies that specifically investigate weak lensing property covariance are scarce, which poses a challenge for upcoming lensing surveys of galaxy clusters such as the Rubin (Ivezić et al. 2019) observatories. To achieve the percentage-level lensing mass calibration goals for the upcoming observations, the property covariance of weak lensing must be quantified.

The physical origins of property covariance in lensing signals of galaxy clusters can be attributed to the halo formation history of the cluster and baryonic physics (Xhakaj et al. 2022). Developing a first-principle physical model for the property covariance as a function of halo formation history and baryonic physics is a daunting task due to the highly non-linear and multiscale physics involved in cluster formation. To make progress, in this paper, we adopt a simulation-based, data-driven approach whereby we develop semi-analytical parametric models of property covariance, which we then calibrate with cosmological simulations. We then apply our model to quantify the bias induced due to a non-zero property covariance in the expected weak lensing signal and the mass-observable scaling relation.

As will be presented in Section 3, a key element of this analysis is the estimation of true cluster richness by encircling clusters within a three-dimensional (3D) radius within the physical vicinity of the halo centre, as opposed to a two-dimensional (2D) projected radius used by cluster finders as redMaPPer (Rykoff et al. 2014) by identifying galaxies within the red-sequence band in the colour–magnitude space – the major difference being the removal of projection effects, or the mis-identification of non-cluster galaxies in the 2D projected radius from the photometric redshift uncertainty of the red-sequence when estimating the true richness from a gravitationally bound region around the halo. Furthermore, as this simulation-based study does not introduce other observational systematics as shape noise of galaxies, point spread function, miscentring, among others, this study can be used to explore the intrinsic covariance between observables prior to the addition of extrinsic systematics as projection effects. Our results will not only provide insight into the physical origin of the covariance, the difference between the total covariance as measured by observations and the intrinsic covariance will provide estimates on the amplitude of the extrinsic component.

The goals of this work are to (i) develop an analytical model that accounts for and quantifies the effect of non-zero covariance on cluster mass calibration, (ii) quantify this property covariance utilizing cosmological simulations, (iii) update uncertainties on inferred cluster mass estimates, and (iv) explain the physical origin

of the covariance using secondary halo parameters. The rest of this paper is organized as follows. In Section 2, we present a population-based analytical framework. In Section 3, we describe the simulations and data-vector employed in this work. In Section 4, we present our measurements and the covariance model. In Section 5, we present the impact of the covariance on weak lensing mass calibration. In Section 6, we quantify the physical origin of the covariance by parametrizing it using secondary halo parameters. In Section 7, we compare our work with those that employ realistic cluster finders. We conclude in Section 8.

## 2 THEORETICAL FRAMEWORK

This section presents a theoretical framework that examines the impact of covariance on mass-observable scaling relations. In Section 2.1, we introduce the definitions of richness and weak lensing excess surface mass density and their scaling relations with cluster mass. We then describe the model of property covariance of richness and excess surface mass density in Section 2.2. In Section 2.3, we model the impact of covariance on stacked observable scaling relations. Finally, in Section 2.4, we develop a theoretical framework that explains the covariance based on a set of secondary halo parameters. A graphic representation of the outline of the paper is shown in Fig. 1. The notations used in this section to describe the covariance are listed in Table 1 and notations for scaling relations are listed in Table 2.

### 2.1 Observable-mass relations

#### 2.1.1 Excess surface mass density $\Delta\Sigma$ from weak lensing

In weak lensing measurements of galaxy clusters, the key observable is the excess surface mass density, denoted $\Delta\Sigma$. The excess surface mass density is defined as

$$ \Delta\Sigma(M, z, r_p) = \bar{\Sigma}(M, z, < r_p) - \Sigma(M, z, r_p), \quad (1) $$

where $\bar{\Sigma}(M, z, < r_p)$ denotes the average surface mass density within projected radius $r_p$, and $\Sigma(M, z, r_p)$ represents the average of the surface mass density at $r_p$. We model the average surface mass density $\Sigma$ as

$$ \Sigma(r_p) = \rho_m \int_{-\infty}^{+\infty} \left( 1 + \xi_{hm} \left( r = \sqrt{r_p^2 + \chi^2} \right) \right) d\chi, \quad (2) $$

where $\rho_m$ is the mean matter density at the redshift of the cluster, $R$ is the projected radius in the plane of the sky, $\chi$ is the comoving distance along the line of sight centred around the cluster, and $\xi_{hm}(r)$ is the halo-matter correlation function that characterizes the total mass density within a halo. Under the halo model, the halo-matter correlation function consists of a ‘one-halo’ term:

$$ \xi_{1h}(r|M) = \frac{\rho_{\text{NFW}}(r|M)}{\rho_{m0}} - 1, \quad (3) $$

and a ‘two-halo’ term:

$$ \xi_{2h}(r|M) = b(M)\xi_{\text{lin}}(r), \quad (4) $$

where $\rho_{\text{NFW}}$ is the Navarro–Frenk–White (NFW) density profile (Navarro, Frenk & White 1997), and $\xi_{\text{lin}}$ is the linear matter correlation function, and $b$ is the halo bias parameter.

In weak lensing, the excess surface density $\Delta\Sigma$ is tied to the tangential shear $\gamma_t$ of the galaxies relative to the centre of each foreground halo by the relation

$$ \Sigma_{\text{crit}}\gamma_t = \bar{\Sigma}(<R) - \Sigma(R) \equiv \Delta\Sigma(R), \quad (5) $$
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Graphic representation of the modelling of Cov(ΔΣ, ln N<sub>gal</sub>|M, z) – the covariance between the halo weak lensing signal ΔΣ(R) and log-richness ln N<sub>gal</sub> conditioned on mass and redshift – and its dependence with secondary halo parameters Π. The labels marked S. XX point to the location in the text. A full list of the notations used in this paper is introduced in Tables 1 and 2.</figcaption></figure>

<table><caption>**Table 1.** Notations employed in our framework for the covariance in Section 2.3.</caption><thead><tr><th>Parameter</th><th>Explanation</th></tr></thead><tbody><tr><td>ΔΣ</td><td>Weak lensing signal</td></tr><tr><td>M</td><td>Halo mass in M<sub>☉</sub> h<sup>-1</sup></td></tr><tr><td>N<sub>gal</sub></td><td>Optical richness enclosed inside 3D radius</td></tr><tr><td>z</td><td>Redshift</td></tr><tr><td>r<sub>p</sub></td><td>Projected and normalized radius</td></tr></tbody></table>

<table><caption>**Table 2.** Scaling relation conventions.</caption><thead><tr><th>Parameter</th><th>Explanation</th></tr></thead><tbody><tr><td>π<sub>a</sub></td><td>Normalization in scaling relation &lang;a|M&rang;</td></tr><tr><td>α<sub>a</sub></td><td>Slope in scaling relation &lang;a|M&rang;</td></tr><tr><td>σ<sub>a</sub></td><td>Scatter about &lang;a|M&rang;</td></tr><tr><td>r<sub>a, b</sub></td><td>Correlation between a and b at fixed M</td></tr><tr><td>π<sub>a|b</sub></td><td>Normalization in scaling relation &lang;a|b&rang;</td></tr><tr><td>α<sub>a|b</sub></td><td>Slope in scaling relation &lang;a|b&rang;</td></tr><tr><td>σ<sub>a|b</sub></td><td>Scatter about &lang;a|b&rang;</td></tr><tr><td>a, b</td><td>a, b ∈ {ΔΣ, ln N<sub>gal</sub>}</td></tr></tbody></table>

where the critical surface density Σ<sub>crit</sub> defined as

$$ \Delta\Sigma_{\text{crit}} = \frac{c^2}{4\pi G} \frac{D_s}{D_l D_{ls}}, \qquad (6) $$

and where D<sub>s</sub>, D<sub>l</sub>, and D<sub>ls</sub> refer to the angular diameter distances to the source, the lens, and between the lens and source, respectively.

In this work, for each halo of mass M at redshift z, we compute the corresponding ΔΣ profile. We compare these measurements with theoretical predictions – in the one-halo regime we model the cluster overdensity as NFW profiles with their concentration determined by concentration-mass models of Prada et al. (2012), Ludlow et al. (2016) and Diemer & Joyce (2019), whereas in the two-halo term we adopt the linear matter correlation ξ<sub>mm</sub> multiplied by halo biases using the Tinker et al. (2010), Pillepich, Porciani & Hahn (2010) and Bhattacharya et al. (2011) models to derive the halo-matter correlation ξ<sub>hm</sub>. At the transition radius between the one- and two-halo regimes, we follow SDSS (Zu et al. 2014) in setting the halo-matter correlation to the maximum value of the two terms, i.e.

$$ \xi_{hm}(r|M) = \max\{\xi_{1h}(r|M), \xi_{2h}(r|M)\}. \qquad (7) $$

In Fig. 2, the theoretical models described above are compared with our measurements of ΔΣ in cosmological simulations to validate our data product.

We model the mean (ΔΣ|M, z, r<sub>p</sub>) at fixed mass M, redshift z, and projected radius r<sub>p</sub> as a log-linear relation given by

$$ \langle \Delta\Sigma | M, z, r_p \rangle_1 = \pi_{\Delta\Sigma}(M, r_p, z) + \alpha_{\Delta\Sigma}(M, r_p, z) \ln M, \qquad (8) $$

where α<sub>ΔΣ</sub> is the power-law slope of the relation and π<sub>ΔΣ</sub> is a normalization that is a function of redshift and mass.

### 2.1.2 Optical richness N<sub>gal</sub>

Optical richness N<sub>gal</sub> is an observable measure of the abundance of galaxies within a galaxy cluster. It is often defined as the number of detected member galaxies brighter than a certain luminosity threshold within a given aperture or radius around the cluster centre.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Graphic representation of the modelling of Cov(ΔΣ, ln N<sub>gal</sub>|M, z) – the covariance between the halo weak lensing signal ΔΣ(R) and log-richness ln N<sub>gal</sub> conditioned on mass and redshift – and its dependence with secondary halo parameters Π. The labels marked S. XX point to the location in the text. A full list of the notations used in this paper is introduced in Tables 1 and 2.</figcaption></figure>

<table><caption>**Table 1.** Notations employed in our framework for the covariance in Section 2.3.</caption><thead><tr><th>Parameter</th><th>Explanation</th></tr></thead><tbody><tr><td>ΔΣ</td><td>Weak lensing signal</td></tr><tr><td>M</td><td>Halo mass in M<sub>☉</sub> h<sup>-1</sup></td></tr><tr><td>N<sub>gal</sub></td><td>Optical richness enclosed inside 3D radius</td></tr><tr><td>z</td><td>Redshift</td></tr><tr><td>r<sub>p</sub></td><td>Projected and normalized radius</td></tr></tbody></table>

<table><caption>**Table 2.** Scaling relation conventions.</caption><thead><tr><th>Parameter</th><th>Explanation</th></tr></thead><tbody><tr><td>π<sub>a</sub></td><td>Normalization in scaling relation &lang;a|M&rang;</td></tr><tr><td>α<sub>a</sub></td><td>Slope in scaling relation &lang;a|M&rang;</td></tr><tr><td>σ<sub>a</sub></td><td>Scatter about &lang;a|M&rang;</td></tr><tr><td>r<sub>a, b</sub></td><td>Correlation between a and b at fixed M</td></tr><tr><td>π<sub>a|b</sub></td><td>Normalization in scaling relation &lang;a|b&rang;</td></tr><tr><td>α<sub>a|b</sub></td><td>Slope in scaling relation &lang;a|b&rang;</td></tr><tr><td>σ<sub>a|b</sub></td><td>Scatter about &lang;a|b&rang;</td></tr><tr><td>a, b</td><td>a, b ∈ {ΔΣ, ln N<sub>gal</sub>}</td></tr></tbody></table>

where the critical surface density Σ<sub>crit</sub> defined as

$$ \Delta\Sigma_{\text{crit}} = \frac{c^2}{4\pi G} \frac{D_s}{D_l D_{ls}}, \qquad (6) $$

and where D<sub>s</sub>, D<sub>l</sub>, and D<sub>ls</sub> refer to the angular diameter distances to the source, the lens, and between the lens and source, respectively.

In this work, for each halo of mass M at redshift z, we compute the corresponding ΔΣ profile. We compare these measurements with theoretical predictions – in the one-halo regime we model the cluster overdensity as NFW profiles with their concentration determined by concentration-mass models of Prada et al. (2012), Ludlow et al. (2016) and Diemer & Joyce (2019), whereas in the two-halo term we adopt the linear matter correlation ξ<sub>mm</sub> multiplied by halo biases using the Tinker et al. (2010), Pillepich, Porciani & Hahn (2010) and Bhattacharya et al. (2011) models to derive the halo-matter correlation ξ<sub>hm</sub>. At the transition radius between the one- and two-halo regimes, we follow SDSS (Zu et al. 2014) in setting the halo-matter correlation to the maximum value of the two terms, i.e.

$$ \xi_{hm}(r|M) = \max\{\xi_{1h}(r|M), \xi_{2h}(r|M)\}. \qquad (7) $$

In Fig. 2, the theoretical models described above are compared with our measurements of ΔΣ in cosmological simulations to validate our data product.

We model the mean (ΔΣ|M, z, r<sub>p</sub>) at fixed mass M, redshift z, and projected radius r<sub>p</sub> as a log-linear relation given by

$$ \langle \Delta\Sigma | M, z, r_p \rangle_1 = \pi_{\Delta\Sigma}(M, r_p, z) + \alpha_{\Delta\Sigma}(M, r_p, z) \ln M, \qquad (8) $$

where α<sub>ΔΣ</sub> is the power-law slope of the relation and π<sub>ΔΣ</sub> is a normalization that is a function of redshift and mass.

### 2.1.2 Optical richness N<sub>gal</sub>

Optical richness N<sub>gal</sub> is an observable measure of the abundance of galaxies within a galaxy cluster. It is often defined as the number of detected member galaxies brighter than a certain luminosity threshold within a given aperture or radius around the cluster centre.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** The measured $\Delta\Sigma$ profiles using downsampled particles for every 10 particles and theoretical $\Delta\Sigma$ as computed from the NFW profile using different concentration-mass relations (LHS in legend) in the one-halo regime and different halo-bias models (RHS in legend) in the two halo regimes, with errors taken to be $1 - \sigma$ standard deviations; the measurements are consistent with theoretical predictions and the size of the errors is too large to distinguish between models. The same conclusion (not shown) holds for $\Delta\Sigma$ binned by $R_{\text{vir}}$.</figcaption></figure>

Richness is often used as a proxy for cluster mass, as more massive clusters are expected to have more member galaxies (e.g. Rozo et al. 2014; Rykoff et al. 2014). The richness-mass scaling relation relates the richness of a galaxy cluster to its mass. In this work, we consider the mean $N_{\text{gal}}-M_{\Delta}$ scaling relation expressed as

$$ \langle \ln N_{\text{gal}} | M_{\Delta}, z \rangle_1 = \pi_{N_{\text{gal}}} (M, z) + \alpha_{N_{\text{gal}}} (M, z) \ln M_{\Delta}. \quad (9) $$

where $M_{\Delta}$ is the mass of the halo within a radius where the mean density is $\Delta$ times the critical density of the universe, $\alpha_{N_{\text{gal}}}(M, z)$ is

the power-law slope of the relation, $\pi_{N_{\text{gal}}}(M, z)$ is a normalization that is a function of redshift and mass.

### 2.1.3 Halo mass and radius definitions

A common approach to defining a radial boundary of a galaxy cluster is such that the average matter density inside a given radius is the product of a reference overdensity $\Delta_{\text{ref}}$ times the critical ($\rho_c$) or mean
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The measured $\Delta\Sigma$ profiles using downsampled particles for every 10 particles and theoretical $\Delta\Sigma$ as computed from the NFW profile using different concentration-mass relations (LHS in legend) in the one-halo regime and different halo-bias models (RHS in legend) in the two halo regimes, with errors taken to be $1 - \sigma$ standard deviations; the measurements are consistent with theoretical predictions and the size of the errors is too large to distinguish between models. The same conclusion (not shown) holds for $\Delta\Sigma$ binned by $R_{vir}$.</figcaption></figure>

Richness is often used as a proxy for cluster mass, as more massive clusters are expected to have more member galaxies (e.g. Rozo et al. 2014; Rykoff et al. 2014). The richness-mass scaling relation relates the richness of a galaxy cluster to its mass. In this work, we consider the mean $N_{gal}-M_{\Delta}$ scaling relation expressed as

$$ \langle \ln N_{\text{gal}} | M_{\Delta}, z \rangle_1 = \pi_{N_{\text{gal}}} (M, z) + \alpha_{N_{\text{gal}}} (M, z) \ln M_{\Delta}. \quad (9) $$

where $M_{\Delta}$ is the mass of the halo within a radius where the mean density is $\Delta$ times the critical density of the universe, $\alpha_{N_{gal}}(M, z)$ is

the power-law slope of the relation, $\pi_{N_{gal}}(M, z)$ is a normalization that is a function of redshift and mass.

### 2.1.3 Halo mass and radius definitions

A common approach to defining a radial boundary of a galaxy cluster is such that the average matter density inside a given radius is the product of a reference overdensity $\Delta_{ref}$ times the critical ($\rho_c$) or mean
density ($\rho_m$) of the universe at that redshift. The critical density is defined as

$$ \rho_c = \frac{3H_0^2}{8\pi G} E(z), \quad (10) $$

where $E(z)^2 = \Omega_{m,0}(1+z)^3 + \Omega_{\Lambda,0}$, $\Omega_{m,0}$ is the present-day matter fraction of the universe, $\Omega_{\Lambda,0}$ is the dark energy fraction at the present age such that $\Omega_{m,0} + \Omega_{\Lambda,0} = 1$ for a flat universe ignoring the minimal contribution from the radiation fraction. The mean (background) density is defined as

$$ \rho_b = \frac{3H_0^2}{8\pi G} (\Omega_{m,0}(1+z)^3). \quad (11) $$

The overdensity $\Delta_c = 200$ is commonly chosen as the reference overdensity in optical weak lensing studies and is closely related to the virial radius. Another radius definition is the virial radius $R_{vir}$, with overdensity values calibrated from cosmological simulations (Bryan & Norman 1998). In this work, we use $R_{200c}$ and $R_{vir}$ to scale various observations, including the $\Delta\Sigma$ measurements and richness. Since the covariance is close to zero at the outskirts $R \gtrsim R_{200c}$ as shown in Section 4, we adopt $r_p = R/R_{200c}$ and $r_p = R/R_{vir}$ as our normalized radii, as the cluster properties are more self-similar with respect to $\rho_c(z)$ compared to $\rho_b(z)$ (Diemer & Kravtsov 2014; Lau et al. 2015). To test for the robustness of our covariance against different radii definitions, we also introduce a physical radius of a toy model of a constant $R = 1$ Mpc $h^{-1}$; here $h = 0.6777$ is the reduced Hubble constant used in this study.

## 2.2 Covariance between $\Delta\Sigma$ and $N_{gal}$

In optical surveys, we cannot expect the covariance between richness $N_{gal}$ and the excess surface mass density $\Delta\Sigma$ to be zero. Ignoring this covariance will lead to bias in cluster mass inferred from the excess surface mass density of the cluster selected based on richness. This work aims to quantify and analyse this covariance and its impact on the mass calibration relation. To achieve this objective, we must first specify the joint probability distribution of excess surface mass density and richness, $p(\Delta\Sigma, \ln N_{gal}|M, z, r_p)$. In this work, we assume that the joint distribution follows a multivariate normal distribution (Stanek et al. 2010; Evrard et al. 2014; Mulroy et al. 2019; Miyatake et al. 2022), which is fully specified with two components, the mean vector and the property covariance. We have checked the goodness of this assumption in Appendix D.

From the mean observable-mass scaling relations in equations (8) and (9), the scaling relation between these two observables can be modelled as a local linear relation given by

$$ \langle \Delta\Sigma | N_{gal}, z, r_p \rangle = \pi \Delta\Sigma |_{N_{gal}} (N_{gal}, z, r_p) + \alpha \Delta\Sigma |_{N_{gal}} (N_{gal}, z, r_p) \ln N_{gal}, \quad (12) $$

where $\pi$ and $\alpha$ are the normalization and slope of the model.

The property covariance matrix is a combination of scatter and correlation between the scatter of $\Delta\Sigma$ and $\ln N_{gal}$ at a fixed halo mass, redshift, and projected radius. We use $\sigma_{N_{gal}}(M, z)$ and $\sigma_{\Delta\Sigma}(M, z, r_p)$ to denote the scatter of the observable-mass relation for $\ln N_{gal}$ and $\Delta\Sigma$, respectively, and use $r_{N_{gal}, \Delta\Sigma}(M, z, r_p)$ to denote the correlation between these scatters. The covariance matrix is then given by

$$ \mathrm{Cov}_{i,j}(M, z, r_p) = r_{i,j}(M, z, r_p) \sigma_i(M, z, r_p) \sigma_j(M, z, r_p), \quad (13) $$

where $i$ and $j \in \{\Delta\Sigma, \ln N_{gal}\}$. Specifically, the covariance between $\Delta\Sigma$ and $N_{gal}$ can be expressed in terms of the residuals about the mean quantities

$$ \mathrm{Cov}_{\Delta\Sigma, N_{gal}}(M, z, r_p) = \mathrm{Cov}(\mathrm{res}_{\Delta\Sigma}(M, z, r_p), \mathrm{res}_{N_{gal}}(M, z)), \quad (14) $$

where the residuals of the $\Delta\Sigma$ and $N_{gal}$ are, respectively, given by

$$ \mathrm{res}_{\Delta\Sigma}(M, z, r_p) = \Delta\Sigma - \langle \Delta\Sigma | M, z, r_p \rangle, \quad (15) $$

$$ \mathrm{res}_{N_{gal}}(M, z) = \ln N_{gal} - \langle \ln N_{gal} | M, z \rangle. \quad (16) $$

To model the mass dependencies of the mean profiles of $\Delta\Sigma$ and $\ln N_{gal}$, we employ the Kernel Localized Linear Regression (KLLR; Farahi, Anbajagane & Evrard 2022a) method. This regression method fits a locally linear model while capturing globally non-linear trends in data and has shown to be effective in modelling halo mass dependencies in scaling relations (Farahi et al. 2018; Anbajagane, Evrard & Farahi 2022; Wu et al. 2022). By developing a local-linear model of $\Delta\Sigma - \ln N_{gal}$ with respect to the halo mass and computing the residuals about the mean relation, we remove the bias in the scatter due to the mass dependence and reduce the overall size of the scatter. As shown in Fig. 4 the $1 - \sigma$ of the covariance is determined by bootstrap resampling.

## 2.3 Corrections to the $\Delta\Sigma - N_{gal}$ relation due to covariance

The shape of the halo mass function plays an important role in evaluating the conditional mean value of $\langle \Delta\Sigma | N_{gal}, z, r_p \rangle$ where the scatter between two observables with a fixed halo mass is correlated. Ignoring the contribution from the correlated scatter, to the zeroth order, the expected $\Delta\Sigma$ evaluated at fixed richness is given by equation (12). This is the model that has been used in mass calibration with stacked weak lensing profiles (Johnston et al. 2007; Kettula et al. 2015; McClintock et al. 2019; Chiu et al. 2020; Lesci et al. 2022).

The first- and second-order approximations of the scaling relation are given by

$$ \langle \Delta\Sigma | N_{gal}, z \rangle_1 = \langle \Delta\Sigma | N_{gal}, z \rangle_{\mathrm{fid}} + \frac{\gamma_1}{\alpha_{N_{gal}}} \times \mathrm{Cov}(\Delta\Sigma, \ln N_{gal}), \quad (17) $$

and

$$ \langle \Delta\Sigma | N_{gal}, z \rangle_2 = \langle \Delta\Sigma | N_{gal}, z \rangle_{\mathrm{fid}} + \mathrm{Cov}(\Delta\Sigma, \ln N_{gal}) \times \left[ \frac{x_s}{\alpha_{N_{gal}}^2} (\alpha_{N_{gal}} \gamma_1 + \gamma_2(\ln N_{gal} - \pi_{N_{gal}})) \right], \quad (18) $$

where $\langle \Delta\Sigma | N_{gal}, z \rangle_{\mathrm{fid}}$ is the fiducial relation taking into account the curvature of the HMF but independent of the covariance; $x_s = (1 + \gamma_2 \sigma_{M|N_{gal},1})^{-1}$ is the compression factor due to curvature of the HMF, the subscript 1 denoting that the scatter is taken from the HMF expanded to first order; here, we omit the $(M, z)$ dependence of the covariance as a shorthand notation. These expansions around the pivot mass are for haloes centred around a narrow enough mass bin. We show explicitly in Fig. 9 that the first-order expansion converges using our binning method. The derivations for the first and second-order expansion terms can be found in Evrard et al. (2014) and Farahi et al. (2018) and the derivation for this particular expression of the second-order term is shown in Appendix C.

Here, $\gamma_1$ and $\gamma_2$ are the parameters for the first- and second-order approximations to the mass dependence of the halo mass function (e.g. Evrard et al. 2014):

$$ \frac{\mathrm{d}n_{\mathrm{hmf}}(M, z)}{\mathrm{d} \ln M} \approx A(z) \exp \left[ -\gamma_1(M, z) \ln M - \frac{1}{2} \gamma_2(M, z) (\ln M)^2 \right]. \quad (19) $$
density ($\rho_m$) of the universe at that redshift. The critical density is defined as

$$ \rho_c = \frac{3H_0^2}{8\pi G} E(z), \qquad (10) $$

where $E(z)^2 = \Omega_{m,0}(1+z)^3 + \Omega_{\Lambda,0}$, $\Omega_{m,0}$ is the present-day matter fraction of the universe, $\Omega_{\Lambda,0}$ is the dark energy fraction at the present age such that $\Omega_{m,0} + \Omega_{\Lambda,0} = 1$ for a flat universe ignoring the minimal contribution from the radiation fraction. The mean (background) density is defined as

$$ \rho_b = \frac{3H_0^2}{8\pi G} (\Omega_{m,0}(1+z)^3). \qquad (11) $$

The overdensity $\Delta_c = 200$ is commonly chosen as the reference overdensity in optical weak lensing studies and is closely related to the virial radius. Another radius definition is the virial radius $R_{vir}$, with overdensity values calibrated from cosmological simulations (Bryan & Norman 1998). In this work, we use $R_{200c}$ and $R_{vir}$ to scale various observations, including the $\Delta\Sigma$ measurements and richness. Since the covariance is close to zero at the outskirts $R \gtrsim R_{200c}$ as shown in Section 4, we adopt $r_p = R/R_{200c}$ and $r_p = R/R_{vir}$ as our normalized radii, as the cluster properties are more self-similar with respect to $\rho_c(z)$ compared to $\rho_b(z)$ (Diemer & Kravtsov 2014; Lau et al. 2015). To test for the robustness of our covariance against different radii definitions, we also introduce a physical radius of a toy model of a constant $R = 1$ Mpc $h^{-1}$; here $h = 0.6777$ is the reduced Hubble constant used in this study.

## 2.2 Covariance between $\Delta\Sigma$ and $N_{gal}$

In optical surveys, we cannot expect the covariance between richness $N_{gal}$ and the excess surface mass density $\Delta\Sigma$ to be zero. Ignoring this covariance will lead to bias in cluster mass inferred from the excess surface mass density of the cluster selected based on richness. This work aims to quantify and analyse this covariance and its impact on the mass calibration relation. To achieve this objective, we must first specify the joint probability distribution of excess surface mass density and richness, $p(\Delta\Sigma, \ln N_{gal}|M, z, r_p)$. In this work, we assume that the joint distribution follows a multivariate normal distribution (Stanek et al. 2010; Evrard et al. 2014; Mulroy et al. 2019; Miyatake et al. 2022), which is fully specified with two components, the mean vector and the property covariance. We have checked the goodness of this assumption in Appendix D.

From the mean observable-mass scaling relations in equations (8) and (9), the scaling relation between these two observables can be modelled as a local linear relation given by

$$ \langle \Delta\Sigma | N_{gal}, z, r_p \rangle = \pi \Delta\Sigma |_{N_{gal}}(N_{gal}, z, r_p) + \alpha \Delta\Sigma |_{N_{gal}}(N_{gal}, z, r_p) \ln N_{gal}, \qquad (12) $$

where $\pi$ and $\alpha$ are the normalization and slope of the model.

The property covariance matrix is a combination of scatter and correlation between the scatter of $\Delta\Sigma$ and $\ln N_{gal}$ at a fixed halo mass, redshift, and projected radius. We use $\sigma_{N_{gal}}(M, z)$ and $\sigma_{\Delta\Sigma}(M, z, r_p)$ to denote the scatter of the observable-mass relation for $\ln N_{gal}$ and $\Delta\Sigma$, respectively, and use $r_{N_{gal}, \Delta\Sigma}(M, z, r_p)$ to denote the correlation between these scatters. The covariance matrix is then given by

$$ \mathrm{Cov}_{i,j}(M, z, r_p) = r_{i,j}(M, z, r_p) \sigma_i(M, z, r_p) \sigma_j(M, z, r_p), \qquad (13) $$

where $i$ and $j \in \{\Delta\Sigma, \ln N_{gal}\}$. Specifically, the covariance between $\Delta\Sigma$ and $N_{gal}$ can be expressed in terms of the residuals about the mean quantities

$$ \mathrm{Cov}_{\Delta\Sigma, N_{gal}}(M, z, r_p) = \mathrm{Cov}(\mathrm{res}_{\Delta\Sigma}(M, z, r_p), \mathrm{res}_{N_{gal}}(M, z)), \qquad (14) $$

where the residuals of the $\Delta\Sigma$ and $N_{gal}$ are, respectively, given by

$$ \mathrm{res}_{\Delta\Sigma}(M, z, r_p) = \Delta\Sigma - \langle \Delta\Sigma | M, z, r_p \rangle, \qquad (15) $$

$$ \mathrm{res}_{N_{gal}}(M, z) = \ln N_{gal} - \langle \ln N_{gal} | M, z \rangle. \qquad (16) $$

To model the mass dependencies of the mean profiles of $\Delta\Sigma$ and $\ln N_{gal}$, we employ the Kernel Localized Linear Regression (KLLR; Farahi, Anbajagane & Evrard 2022a) method. This regression method fits a locally linear model while capturing globally non-linear trends in data and has shown to be effective in modelling halo mass dependencies in scaling relations (Farahi et al. 2018; Anbajagane, Evrard & Farahi 2022; Wu et al. 2022). By developing a local-linear model of $\Delta\Sigma - \ln N_{gal}$ with respect to the halo mass and computing the residuals about the mean relation, we remove the bias in the scatter due to the mass dependence and reduce the overall size of the scatter. As shown in Fig. 4 the $1 - \sigma$ of the covariance is determined by bootstrap resampling.

## 2.3 Corrections to the $\Delta\Sigma - N_{gal}$ relation due to covariance

The shape of the halo mass function plays an important role in evaluating the conditional mean value of $\langle \Delta\Sigma | N_{gal}, z, r_p \rangle$ where the scatter between two observables with a fixed halo mass is correlated. Ignoring the contribution from the correlated scatter, to the zeroth order, the expected $\Delta\Sigma$ evaluated at fixed richness is given by equation (12). This is the model that has been used in mass calibration with stacked weak lensing profiles (Johnston et al. 2007; Kettula et al. 2015; McClintock et al. 2019; Chiu et al. 2020; Lesci et al. 2022).

The first- and second-order approximations of the scaling relation are given by

$$ \langle \Delta\Sigma | N_{gal}, z \rangle_1 = \langle \Delta\Sigma | N_{gal}, z \rangle_{\mathrm{fid}} + \frac{\gamma_1}{\alpha_{N_{gal}}} \times \mathrm{Cov}(\Delta\Sigma, \ln N_{gal}), \qquad (17) $$

and

$$ \langle \Delta\Sigma | N_{gal}, z \rangle_2 = \langle \Delta\Sigma | N_{gal}, z \rangle_{\mathrm{fid}} + \mathrm{Cov}(\Delta\Sigma, \ln N_{gal}) \times \left[ \frac{x_s}{\alpha_{N_{gal}}^2} (\alpha_{N_{gal}} \gamma_1 + \gamma_2(\ln N_{gal} - \pi_{N_{gal}})) \right], \qquad (18) $$

where $\langle \Delta\Sigma | N_{gal}, z \rangle_{\mathrm{fid}}$ is the fiducial relation taking into account the curvature of the HMF but independent of the covariance; $x_s = (1 + \gamma_2 \sigma_{M|N_{gal},1})^{-1}$ is the compression factor due to curvature of the HMF, the subscript 1 denoting that the scatter is taken from the HMF expanded to first order; here, we omit the $(M, z)$ dependence of the covariance as a shorthand notation. These expansions around the pivot mass are for haloes centred around a narrow enough mass bin. We show explicitly in Fig. 9 that the first-order expansion converges using our binning method. The derivations for the first and second-order expansion terms can be found in Evrard et al. (2014) and Farahi et al. (2018) and the derivation for this particular expression of the second-order term is shown in Appendix C.

Here, $\gamma_1$ and $\gamma_2$ are the parameters for the first- and second-order approximations to the mass dependence of the halo mass function (e.g. Evrard et al. 2014):

$$ \frac{\mathrm{d}n_{\mathrm{hmf}}(M, z)}{\mathrm{d} \ln M} \approx A(z) \exp \left[ -\gamma_1(M, z) \ln M - \frac{1}{2} \gamma_2(M, z) (\ln M)^2 \right]. \qquad (19) $$
<table><caption>**Table 3.** Notations employed in exploring the secondary halo parameter dependence.</caption><thead><tr><th>Parameter Explanation</th><th></th></tr></thead><tbody><tr><td>Π</td><td>Set of secondary halo parameters</td></tr><tr><td>Γ<sub>inst</sub></td><td>Instantaneous mass accretion rate (MAR)</td></tr><tr><td>Γ<sub>100 Myr</sub></td><td>Mean MAR over the past 100 Myr</td></tr><tr><td>Γ<sub>dyn</sub></td><td>Mean MAR over virial dynamical time</td></tr><tr><td>Γ<sub>2dyn</sub></td><td>Mean MAR over two virial dynamical time</td></tr><tr><td>Γ<sub>peak</sub></td><td>Growth rate of peak mass from current z to z + 0.5</td></tr><tr><td>a<sub>1/2</sub></td><td>Half mass scale factor</td></tr><tr><td>c<sub>vir</sub></td><td>R<sub>vir</sub> concentration</td></tr><tr><td>T/[U]</td><td>Absolute value of the kinetic to potential energy ratio</td></tr><tr><td>X<sub>off</sub></td><td>Offset of density peak from mean particle position (kpc h<sup>-1</sup>)</td></tr></tbody></table>

where $A(z)$ is the normalization of the mass function due to the redshift alone. In deriving the above approximations, we have made use of the fact that $\text{Cov}(\Delta\Sigma, \ln N_{\text{gal}}|M, z) \equiv r_{\Delta\Sigma, N_{\text{gal}}}\sigma_{\Delta\Sigma}\sigma_{N_{\text{gal}}}$. The terms $\sigma_{\Delta\Sigma}$, $\sigma_{N_{\text{gal}}}$, $r_{N_{\text{gal}}, \Delta\Sigma}$, $\gamma_1$ and $\gamma_2$ are evaluated at the mass implied by $\langle \ln M | N_{\text{gal}} \rangle$.

These property covariance correction terms are absent in the current literature. A key feature of this approximation method is that the second-order solution has better than per cent-level accuracy when the halo mass function is known (Farahi et al. 2016). In Fig. 9, we demonstrate that the statistical uncertainties for the first-order correction in equation (17) is larger than the uncertainty in the halo mass function and the uncertainty due to the second-order halo mass approximation.

## 2.4 Secondary halo parameter dependence of Cov(ΔΣ, ln Ngal | M, z)

We elucidate the physical origin of the covariance between ΔΣ and ln Ngal by developing a phenomenological model based on the secondary halo parameters listed in Table 3. These parameters are computed from the ROCKSTAR halo finder (Behroozi, Wechsler & Wu 2013). They capture the halo's mass accretion history, which we hypothesise is the driving force behind the observed covariance. To incorporate these parameters into our model, we extend equations (8) and (9) by introducing multilinear terms that include the secondary halo parameters denoted by the vector Π:

$$ \begin{aligned} (\ln N_{\text{gal}} | \Pi, M, z) &= \langle \ln N_{\text{gal}} | M, z \rangle_1 + \vec{\beta}_{N_{\text{gal}}}^{\top}(M, z) \cdot \Pi + \epsilon_{N_{\text{gal}}} \\ &= \pi_{N_{\text{gal}}}(M, z) + \alpha_{N_{\text{gal}}}(M, z) \ln M \\ &\quad + \vec{\beta}_{N_{\text{gal}}}^{\top}(M, z) \cdot \Pi + \epsilon_{N_{\text{gal}}}. \end{aligned} \quad (20) $$

$$ \begin{aligned} (\Delta\Sigma | \Pi, M, z) &= \langle \Delta\Sigma | M, z \rangle_1 + \vec{\beta}_{\Delta\Sigma}^{\top}(M, z) \cdot \Pi + \epsilon_{\Delta\Sigma} \\ &= \pi_{\Delta\Sigma}(M, z) + \alpha_{\Delta\Sigma}(M, z) \ln M \\ &\quad + \vec{\beta}_{\Delta\Sigma}^{\top}(M, z) \cdot \Pi + \epsilon_{\Delta\Sigma}, \end{aligned} \quad (21) $$

where Π is a vector of secondary halo parameters of potential interest listed in Table 3, and $\epsilon_{\Delta\Sigma}$ and $\epsilon_{N_{\text{gal}}}$ are normally distributed intrinsic scatter terms with zero means and uncorrelated variances. In Appendix E, we show that the residual conditioned on secondary halo parameters can largely be assumed to be Gaussian. Additionally, we assume $\langle \epsilon_{N_{\text{gal}}} \epsilon_{\Delta\Sigma} \rangle = 0$, which implies that the scatter about the mean relations is uncorrelated after factoring in the secondary properties.

Due to the bilinearity and distributive properties of covariance, combining equations (20) and (21) yields

$$ \begin{aligned} \text{Cov}(\Delta\Sigma, \ln N_{\text{gal}} | M, z) ={}& \text{Cov}(\langle\Delta\Sigma\rangle_1, \langle N_{\text{gal}} \rangle_1) + \text{Cov}(\langle\Delta\Sigma\rangle_1, \vec{\beta}_{N_{\text{gal}}}^{\top} \cdot \Pi) \\ & + \text{Cov}(\langle\Delta\Sigma\rangle_1, \epsilon_{N_{\text{gal}}}) + \text{Cov}(\vec{\beta}_{\Delta\Sigma}^{\top} \cdot \Pi, \langle \ln N_{\text{gal}} \rangle_1) \\ & + \text{Cov}(\vec{\beta}_{\Delta\Sigma}^{\top} \cdot \Pi, \vec{\beta}_{N_{\text{gal}}}^{\top} \cdot \Pi) + \text{Cov}(\vec{\beta}_{\Delta\Sigma}^{\top} \cdot \Pi, \epsilon_{N_{\text{gal}}}) \\ & + \text{Cov}(\epsilon_{\Delta\Sigma}, \langle \ln N_{\text{gal}} \rangle_1) + \text{Cov}(\epsilon_{\Delta\Sigma}, \vec{\beta}_{N_{\text{gal}}}^{\top} \cdot \Pi) \\ & + \text{Cov}(\epsilon_{N_{\text{gal}}}, \epsilon_{\Delta\Sigma}), \end{aligned} \quad (22) $$

where we omit the explicit $(M, z)$ dependence in $\langle \Delta\Sigma | M, z \rangle_1$, $\langle \ln N_{\text{gal}} | M, z \rangle_1$, $\vec{\beta}_{N_{\text{gal}}}^{\top}(M, z)$, $\vec{\beta}_{N_{\text{gal}}}^{\top}(M, z)$, $\epsilon_{N_{\text{gal}}}(M, z)$ and $\epsilon_{\Delta\Sigma}(M, z)$ to simplify the notation. The KLLR method is utilized to estimate these mass-dependent parameters. All terms involving $\langle \Delta\Sigma \rangle$ and $\langle \ln N_{\text{gal}} \rangle$ vanish, as these terms are independent of Π by definition. Terms involving $\epsilon_{N_{\text{gal}}}$ and $\epsilon_{\Delta\Sigma}$ also go to zero, as they are uncorrelated Gaussian scatters. Only the term $\text{Cov}(\vec{\beta}_{\Delta\Sigma}^{\top} \cdot \Pi, \vec{\beta}_{N_{\text{gal}}}^{\top} \cdot \Pi)$ remains, and hence our final expression for the covariance is

$$ \begin{aligned} \text{Cov}(\Delta\Sigma, \ln N_{\text{gal}} | M, z) &= \text{Cov}(\vec{\beta}_{\Delta\Sigma}^{\top} \cdot \Pi, \vec{\beta}_{N_{\text{gal}}}^{\top} \cdot \Pi) \\ &= \vec{\beta}_{\Delta\Sigma}^{\top} \text{Cov}(\Pi, \Pi) \vec{\beta}_{N_{\text{gal}}}^{\top}. \end{aligned} \quad (23) $$

To estimate the error in the covariance due to each of the secondary halo parameters, we compute $\text{Cov}(\Delta\Sigma, \Pi_i | M, z)$ for each secondary halo parameter $i$ in each $(r_p, M, z)$ bin and take their standard deviations as the error measurement. Modelling the richness-mass relation as in equation (20) and using the same derivation as in equation (22), we arrive at the expression

$$ \text{Cov}(\Delta\Sigma, \ln N_{\text{gal}} | M, z) = \sum_i \beta_{N_{\text{gal},i}(M, z)} \text{Cov}(\Delta\Sigma, \Pi_i | M, z), \quad (24) $$

in which the error from each contributing term in Π is the standard deviation for $\text{Cov}(\Delta\Sigma, \Pi_i | M, z)$ multiplied by the partial richness slope $\beta_{N_{\text{gal},i}$. The total variance of $\text{Cov}(\Delta\Sigma, \ln N_{\text{gal}} | M, z)$ are the errors of each term added in quadrature.

We test the validity of this model by checking how well the secondary halo parameters can explain covariance between lensing and richness in Section 6. After subtracting the covariance from each of the $\Pi_i$ terms, the full covariance should be consistent with null, given the uncertainty. Our results confirm that the dependency of secondary halo parameters can indeed explain the covariance.

# 3 DATA SET AND MEASUREMENTS

In this section, we describe the measurements on the individual ingredients that make up the covariance – ΔΣ the lensing signal in Section 3.1 and ln Ngal the log-richness measurement in Section 3.2.

## 3.1 Measurements of ΔΣ

We employ the MultiDark Planck 2 (MDPL2) cosmological simulation (Klypin et al. 2016) to measure halo properties. The MDPL2 is a gravity-only N-body simulation, consisting of 3840³ particles in a periodic box with a side length of $L_{\text{box}} = 1 h^{-1} \text{Gpc}$, yielding a particle mass resolution of approximately $m_p \approx 1.51 \times 10^9 h^{-1} M_{\odot}$. The simulation was conducted with a flat Lambda old dark matter (ΛCDM) cosmology similar to Planck Collaboration XVI (2014), with the following parameters: $h = 0.6777$, $\Omega_m = 0.307115$, $\Omega_\Lambda = 0.692885$, $\sigma_8 = 0.829$, and $n_s = 0.96$. We use the surface overdensity of down-sampled dark matter particles to measure the weak lensing signal. We selected cluster-sized haloes using the
<table><caption>**Table 3.** Notations employed in exploring the secondary halo parameter dependence.</caption><thead><tr><th>Parameter Explanation</th><th></th></tr></thead><tbody><tr><td>Π</td><td>Set of secondary halo parameters</td></tr><tr><td>Γ<sub>inst</sub></td><td>Instantaneous mass accretion rate (MAR)</td></tr><tr><td>Γ<sub>100 Myr</sub></td><td>Mean MAR over the past 100 Myr</td></tr><tr><td>Γ<sub>dyn</sub></td><td>Mean MAR over virial dynamical time</td></tr><tr><td>Γ<sub>2dyn</sub></td><td>Mean MAR over two virial dynamical time</td></tr><tr><td>Γ<sub>peak</sub></td><td>Growth rate of peak mass from current z to z + 0.5</td></tr><tr><td>a<sub>1/2</sub></td><td>Half mass scale factor</td></tr><tr><td>c<sub>vir</sub></td><td>R<sub>vir</sub> concentration</td></tr><tr><td>T/[U]</td><td>Absolute value of the kinetic to potential energy ratio</td></tr><tr><td>X<sub>off</sub></td><td>Offset of density peak from mean particle position (kpc h<sup>-1</sup>)</td></tr></tbody></table>

where A(z) is the normalization of the mass function due to the redshift alone. In deriving the above approximations, we have made use of the fact that $Cov(ΔΣ, ln N_{gal}|M, z) ≡ r_{ΔΣ, N_{gal}}σ_{ΔΣ}σ_{N_{gal}}$. The terms $σ_{ΔΣ}$, $σ_{N_{gal}}$, $r_{N_{gal}, ΔΣ}$, $γ_1$ and $γ_2$ are evaluated at the mass implied by $⟨ln M|N_{gal}⟩$.

These property covariance correction terms are absent in the current literature. A key feature of this approximation method is that the second-order solution has better than per cent-level accuracy when the halo mass function is known (Farahi et al. 2016). In Fig. 9, we demonstrate that the statistical uncertainties for the first-order correction in equation (17) is larger than the uncertainty in the halo mass function and the uncertainty due to the second-order halo mass approximation.

## 2.4 Secondary halo parameter dependence of Cov(ΔΣ, ln N<sub>gal</sub>|M, z)

We elucidate the physical origin of the covariance between ΔΣ and ln N<sub>gal</sub> by developing a phenomenological model based on the secondary halo parameters listed in Table 3. These parameters are computed from the ROCKSTAR halo finder (Behroozi, Wechsler & Wu 2013). They capture the halo’s mass accretion history, which we hypothesise is the driving force behind the observed covariance. To incorporate these parameters into our model, we extend equations (8) and (9) by introducing multilinear terms that include the secondary halo parameters denoted by the vector Π:

$$ 
\begin{aligned}
(\ln N_{\text{gal}} | \Pi, M, z) &= \langle \ln N_{\text{gal}} | M, z \rangle_1 + \vec{\beta}_{N_{\text{gal}}}^{\top}(M, z) \cdot \Pi + \epsilon_{N_{\text{gal}}} \\
&= \pi_{N_{\text{gal}}}(M, z) + \alpha_{N_{\text{gal}}}(M, z) \ln M \\
&\quad + \vec{\beta}_{N_{\text{gal}}}^{\top}(M, z) \cdot \Pi + \epsilon_{N_{\text{gal}}}.
\end{aligned} 
\quad (20) $$

$$ 
\begin{aligned}
(\Delta\Sigma | \Pi, M, z) &= \langle \Delta\Sigma | M, z \rangle_1 + \vec{\beta}_{\Delta\Sigma}^{\top}(M, z) \cdot \Pi + \epsilon_{\Delta\Sigma} \\
&= \pi_{\Delta\Sigma}(M, z) + \alpha_{\Delta\Sigma}(M, z) \ln M \\
&\quad + \vec{\beta}_{\Delta\Sigma}^{\top}(M, z) \cdot \Pi + \epsilon_{\Delta\Sigma},
\end{aligned} 
\quad (21) $$

where Π is a vector of secondary halo parameters of potential interest listed in Table 3, and $ϵ_{ΔΣ}$ and $ϵ_{N_{gal}}$ are normally distributed intrinsic scatter terms with zero means and uncorrelated variances. In Appendix E, we show that the residual conditioned on secondary halo parameters can largely be assumed to be Gaussian. Additionally, we assume $⟨ϵ_{N_{gal}}ϵ_{ΔΣ}⟩ = 0$, which implies that the scatter about the mean relations is uncorrelated after factoring in the secondary properties.

Due to the bilinearity and distributive properties of covariance, combining equations (20) and (21) yields

$$ 
\begin{aligned}
\text{Cov}(\Delta\Sigma, \ln N_{\text{gal}} | M, z) &= \text{Cov}(\langle\Delta\Sigma\rangle_1, \langle N_{\text{gal}}\rangle_1) + \text{Cov}(\langle\Delta\Sigma\rangle_1, \vec{\beta}_{N_{\text{gal}}}^{\top} \cdot \Pi) \\
&\quad + \text{Cov}(\langle\Delta\Sigma\rangle_1, \epsilon_{N_{\text{gal}}}) + \text{Cov}(\vec{\beta}_{\Delta\Sigma}^{\top} \cdot \Pi, \langle\ln N_{\text{gal}}\rangle_1) \\
&\quad + \text{Cov}(\vec{\beta}_{\Delta\Sigma}^{\top} \cdot \Pi, \vec{\beta}_{N_{\text{gal}}}^{\top} \cdot \Pi) + \text{Cov}(\vec{\beta}_{\Delta\Sigma}^{\top} \cdot \Pi, \epsilon_{N_{\text{gal}}}) \\
&\quad + \text{Cov}(\epsilon_{\Delta\Sigma}, \langle\ln N_{\text{gal}}\rangle_1) + \text{Cov}(\epsilon_{\Delta\Sigma}, \vec{\beta}_{N_{\text{gal}}}^{\top} \cdot \Pi) \\
&\quad + \text{Cov}(\epsilon_{N_{\text{gal}}}, \epsilon_{\Delta\Sigma}),
\end{aligned} 
\quad (22) $$

where we omit the explicit (M, z) dependence in $⟨\DeltaΣ|M, z⟩_1$, $⟨\ln N_{gal}|M, z⟩_1$, $\vec{\beta}_{N_{gal}}^{\top}(M, z)$, $\vec{\beta}_{N_{gal}}^{\top}(M, z)$, $ϵ_{N_{gal}}(M, z)$ and $ϵ_{ΔΣ}(M, z)$ to simplify the notation. The KLLR method is utilized to estimate these mass-dependent parameters. All terms involving $⟨\DeltaΣ⟩$ and $⟨\ln N_{gal}⟩$ vanish, as these terms are independent of Π by definition. Terms involving $ϵ_{N_{gal}}$ and $ϵ_{ΔΣ}$ also go to zero, as they are uncorrelated Gaussian scatters. Only the term $Cov(\vec{\beta}_{\DeltaΣ}^{\top} \cdot \Pi, \vec{\beta}_{N_{gal}}^{\top} \cdot \Pi)$ remains, and hence our final expression for the covariance is

$$ 
\begin{aligned}
\text{Cov}(\Delta\Sigma, \ln N_{\text{gal}} | M, z) &= \text{Cov}(\vec{\beta}_{\Delta\Sigma}^{\top} \cdot \Pi, \vec{\beta}_{N_{\text{gal}}}^{\top} \cdot \Pi) \\
&= \vec{\beta}_{\Delta\Sigma}^{\top} \text{Cov}(\Pi, \Pi) \vec{\beta}_{N_{\text{gal}}}^{\top}.
\end{aligned} 
\quad (23) $$

To estimate the error in the covariance due to each of the secondary halo parameters, we compute $Cov(ΔΣ, Π_i|M, z)$ for each secondary halo parameter i in each $(r_p, M, z)$ bin and take their standard deviations as the error measurement. Modelling the richness–mass relation as in equation (20) and using the same derivation as in equation (22), we arrive at the expression

$$ \text{Cov}(\Delta\Sigma, \ln N_{\text{gal}} | M, z) = \sum_{i} \beta_{N_{\text{gal},i}(M, z)} \text{Cov}(\Delta\Sigma, \Pi_i | M, z), \quad (24) $$

in which the error from each contributing term in Π is the standard deviation for $Cov(ΔΣ, Π_i|M, z)$ multiplied by the partial richness slope $β_{N_{gal},i}$. The total variance of $Cov(ΔΣ, ln N_{gal}|M, z)$ are the errors of each term added in quadrature.

We test the validity of this model by checking how well the secondary halo parameters can explain covariance between lensing and richness in Section 6. After subtracting the covariance from each of the $Π_i$ terms, the full covariance should be consistent with null, given the uncertainty. Our results confirm that the dependency of secondary halo parameters can indeed explain the covariance.

# 3 DATA SET AND MEASUREMENTS

In this section, we describe the measurements on the individual ingredients that make up the covariance – ΔΣ the lensing signal in Section 3.1 and ln N<sub>gal</sub> the log-richness measurement in Section 3.2.

## 3.1 Measurements of ΔΣ

We employ the MultiDark Planck 2 (MDPL2) cosmological simulation (Klypin et al. 2016) to measure halo properties. The MDPL2 is a gravity-only N-body simulation, consisting of 3840³ particles in a periodic box with a side length of $L_{\text{box}} = 1 h^{-1} \text{Gpc}$, yielding a particle mass resolution of approximately $m_p \approx 1.51 \times 10^9 h^{-1} \text{M}_\odot$. The simulation was conducted with a flat Lambda old dark matter (ΛCDM) cosmology similar to Planck Collaboration XVI (2014), with the following parameters: $h = 0.6777$, $\Omega_m = 0.307115$, $\Omega_\Lambda = 0.692885$, $\sigma_8 = 0.829$, and $n_s = 0.96$. We use the surface overdensity of down-sampled dark matter particles to measure the weak lensing signal. We selected cluster-sized haloes using the
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Using different prescriptions of the richness count, we compare with the SPT-DES data (Costanzi et al. 2021). The richness estimator, with no stellar or colour–magnitude cut, shows a similar trend with the data. In Section 4, we show that the results are robust to changes in the definition of the number count estimator.</figcaption></figure>

ROCKSTAR (Behroozi et al. 2013) halo catalogue, which includes the primary halo property of mass and redshift and a set of secondary halo properties listed in Table 3 that we utilize in this analysis. To capture the contribution of both the one- and two-halo terms to $\xi_{lm}$, we use a projection depth of $D_p = 200$ Mpc $h^{-1}$ to calculate $\Delta\Sigma$ (e.g. Costanzi et al. 2019; Sunayama et al. 2020). The MDPL2 data products are publicly available through the MultiDark Database (Riebe et al. 2013) and can be downloaded from the CosmoSim website.<fnref n="1" />

The excess overdensity, $\Delta\Sigma$, is calculated by integrating the masses of the dark matter particles in annuli of increasing radius centred around the halo centre. However, since clusters do not have a well-defined boundary, we compare the results of two radial binning schemes. The first scheme uses 20 equally log-spaced ratios between 0.1 and 10 times $R_{vir}$, while the second scheme spans 0.1 to 10 times $R_{200c}$. We consider the measurements binned at $R_{200c}$ as our final results to be consistent with the weak lensing literature. Fig. 2 shows that our measurements are consistent with most models of the concentration-mass and halo-bias models at a $1\sigma$ level.

At a projection depth of $D_p = 200$ Mpc $h^{-1}$, the projection effects can be modelled as a multiplicative bias (Sunayama 2023). In Sunayama (2023) the projection effects on $\Delta\Sigma$ are modelled as $\Delta\Sigma_{obs} = (1 + \alpha)\Delta\Sigma_{true}$, where $\alpha = 18.4 \pm 8.6$ per cent. Although the multiplicative bias of projection effects may increase the amplitude of $\text{Cov}(\Delta\Sigma, \ln N_{gal})$ by a factor of $(1 + \alpha)$, we argue that it does not introduce an additive bias into our model for $\text{Cov}(\Delta\Sigma, \ln N_{gal})$. This is because, under the richness model and $\Delta\Sigma$ in equations (20) and (21), only terms in richness that are correlated with projection effects will contribute to the covariance. As demonstrated in Section 3.2, we enclose the halo within a 3D physical radius, so the number count $N_{gal}$ of galaxies should not include projection effects. Therefore, projection effects should not introduce an additive bias to our covariance.

To remove the 2D integrated background density, we first computed the background density of the universe ($\rho_b$) at the cluster redshift using the cosmological parameters of the MDPL2 simulation. The integrated 2D background density is given by $\Sigma_b = 2D_p\rho_b$, where factor 2 comes from the integration of the foreground and background densities.

## 3.2 Measurements of $N_{gal}$

### 3.2.1 Data set for $N_{gal}$ – SAGE galaxy catalogue

The Semi-Analytic GALAXY Evolution (SAGE) is a catalogue of galaxies within MDPL2, generated through a post-processing step that places galaxies onto N-body simulations. This approach, known as a semi-analytic model (SAM), is computationally efficient compared to hydrodynamical simulations with fully self-consistent baryonic physics. SAMs reduce the computational time required by two to three orders of magnitude, allowing us to populate the entire 1 (Gpc $h^{-1}$)$^3$ simulation volume with galaxies. SAGE’s statistical power enables us to conduct stacked weak lensing analyses.

The baryonic prescription of SAGE is based on the work of Croton et al. (2016), which includes updated physics in baryonic transfer processes such as gas infall, cooling, heating, and reionization. It also includes an intracluster star component for central galaxies and addresses the orphan galaxy problem by adding the stellar mass of disrupted satellite galaxies as intra-‘cluster’ mass. SAGE’s primary data constraint is the stellar mass function at $z = 0$. Secondary constraints include the star formation rate density history (Somerville, Primack & Faber 2001), the Baryonic Tully-Fisher relation (Stark, McGaugh & Swaters 2009), the mass metallicity relation of galaxies (Tremonti et al. 2004), and the black hole–bulge mass relation (Stark et al. 2009).

### 3.2.2 Model for $N_{gal}$

To determine the number of galaxies inside a cluster-sized halo, we utilise the SAGE SAM and compute the total number of galaxies within a 3D radius around the halo centre. We compare the true richness ($N_{gal}$) to $M_{200c}$ scaling relations between different models and the observed richness-mass relations from Costanzi et al. (2021) using data from the DES Year-1 catalogue and mass-observable-relation from the South Pole Telescope (SPT) cluster catalogue (see Fig. 3). The observed richness-mass relation is fitted as a log-linear model with $2 - \sigma$ error bars that trace the posterior of the best-fitting richness-mass model parameters. The $M_{500c}$ mass definition in the catalogue is converted to $M_{200c}$ using an NFW profile for the surface density of the cluster and adopting the Diemer & Joyce (2019) concentration-mass relation anchored at $z = 0.35$, which is roughly the median redshift of the cluster sample.

We use the KLLR method to determine the local linear fit for our $N_{gal}$-mass model, which relaxes the assumption of global log-linearity (Anbajagane et al. 2020). Realistic cluster finders, such as redMaPPer (Rykoff et al. 2014), impose a colour-magnitude cut or a stellar mass cut, which are highly dependent on the red-sequence model or the spectral energy density model. We found that imposing a stellar mass cut of $10.5\log(M/M_{\odot})$ would correspond roughly to the bottom 5 per cent percentile of SDSS detected galaxies (Maraston et al. 2013). However, this drastically decreases the number of galaxies in a halo, with most having $N_{gal}$ in the single digits. As we are interested in the *intrinsic* covariance from the physical properties of the halo, we do not impose additional magnitude or stellar mass cuts. We confirm that, as described in Croton et al. (2016), the galaxy stellar mass distribution at $z = 0$ is consistent with the best-fitting double Schechter function calibrated with low-redshift galaxies from the Galaxy and Mass Assembly (Baldry et al. 2012) down to stellar masses of $M > 10^{8.5} M_{\odot}$.

Fig. 3 illustrates that our $N_{gal}$-mass models, which count the number of galaxies within a physical 3D radius and impose no colour-magnitude cut as redMaPPer does, resemble the general behaviour

<footnote n="1">¹https://www.cosmosim.org</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Using different prescriptions of the richness count, we compare with the SPT-DES data (Costanzi et al. 2021). The richness estimator, with no stellar or colour–magnitude cut, shows a similar trend with the data. In Section 4, we show that the results are robust to changes in the definition of the number count estimator.</figcaption></figure>

ROCKSTAR (Behroozi et al. 2013) halo catalogue, which includes the primary halo property of mass and redshift and a set of secondary halo properties listed in Table 3 that we utilize in this analysis. To capture the contribution of both the one- and two-halo terms to $\xi_{lm}$, we use a projection depth of $D_p = 200$ Mpc $h^{-1}$ to calculate $\Delta\Sigma$ (e.g. Costanzi et al. 2019; Sunayama et al. 2020). The MDPL2 data products are publicly available through the MultiDark Database (Riebe et al. 2013) and can be downloaded from the CosmoSim website.<fnref n="1" />

The excess overdensity, $\Delta\Sigma$, is calculated by integrating the masses of the dark matter particles in annuli of increasing radius centred around the halo centre. However, since clusters do not have a well-defined boundary, we compare the results of two radial binning schemes. The first scheme uses 20 equally log-spaced ratios between 0.1 and 10 times $R_{vir}$, while the second scheme spans 0.1 to 10 times $R_{200c}$. We consider the measurements binned at $R_{200c}$ as our final results to be consistent with the weak lensing literature. Fig. 2 shows that our measurements are consistent with most models of the concentration-mass and halo-bias models at a $1\sigma$ level.

At a projection depth of $D_p = 200$ Mpc $h^{-1}$, the projection effects can be modelled as a multiplicative bias (Sunayama 2023). In Sunayama (2023) the projection effects on $\Delta\Sigma$ are modelled as $\Delta\Sigma_{obs} = (1 + \alpha)\Delta\Sigma_{true}$, where $\alpha = 18.4 \pm 8.6$ per cent. Although the multiplicative bias of projection effects may increase the amplitude of $\text{Cov}(\Delta\Sigma, \ln N_{gal})$ by a factor of $(1 + \alpha)$, we argue that it does not introduce an additive bias into our model for $\text{Cov}(\Delta\Sigma, \ln N_{gal})$. This is because, under the richness model and $\Delta\Sigma$ in equations (20) and (21), only terms in richness that are correlated with projection effects will contribute to the covariance. As demonstrated in Section 3.2, we enclose the halo within a 3D physical radius, so the number count $N_{gal}$ of galaxies should not include projection effects. Therefore, projection effects should not introduce an additive bias to our covariance.

To remove the 2D integrated background density, we first computed the background density of the universe ($\rho_b$) at the cluster redshift using the cosmological parameters of the MDPL2 simulation. The integrated 2D background density is given by $\Sigma_b = 2D_p\rho_b$, where factor 2 comes from the integration of the foreground and background densities.

## 3.2 Measurements of $N_{gal}$

### 3.2.1 Data set for $N_{gal}$ – SAGE galaxy catalogue

The Semi-Analytic GALAXY Evolution (SAGE) is a catalogue of galaxies within MDPL2, generated through a post-processing step that places galaxies onto N-body simulations. This approach, known as a semi-analytic model (SAM), is computationally efficient compared to hydrodynamical simulations with fully self-consistent baryonic physics. SAMs reduce the computational time required by two to three orders of magnitude, allowing us to populate the entire 1 (Gpc $h^{-1}$)$^3$ simulation volume with galaxies. SAGE’s statistical power enables us to conduct stacked weak lensing analyses.

The baryonic prescription of SAGE is based on the work of Croton et al. (2016), which includes updated physics in baryonic transfer processes such as gas infall, cooling, heating, and reionization. It also includes an intracluster star component for central galaxies and addresses the orphan galaxy problem by adding the stellar mass of disrupted satellite galaxies as intra-‘cluster’ mass. SAGE’s primary data constraint is the stellar mass function at $z = 0$. Secondary constraints include the star formation rate density history (Somerville, Primack & Faber 2001), the Baryonic Tully-Fisher relation (Stark, McGaugh & Swaters 2009), the mass metallicity relation of galaxies (Tremonti et al. 2004), and the black hole–bulge mass relation (Stark et al. 2009).

### 3.2.2 Model for $N_{gal}$

To determine the number of galaxies inside a cluster-sized halo, we utilise the SAGE SAM and compute the total number of galaxies within a 3D radius around the halo centre. We compare the true richness ($N_{gal}$) to $M_{200c}$ scaling relations between different models and the observed richness-mass relations from Costanzi et al. (2021) using data from the DES Year-1 catalogue and mass-observable-relation from the South Pole Telescope (SPT) cluster catalogue (see Fig. 3). The observed richness-mass relation is fitted as a log-linear model with $2 - \sigma$ error bars that trace the posterior of the best-fitting richness-mass model parameters. The $M_{500c}$ mass definition in the catalogue is converted to $M_{200c}$ using an NFW profile for the surface density of the cluster and adopting the Diemer & Joyce (2019) concentration-mass relation anchored at $z = 0.35$, which is roughly the median redshift of the cluster sample.

We use the KLLR method to determine the local linear fit for our $N_{gal}$-mass model, which relaxes the assumption of global log-linearity (Anbajagane et al. 2020). Realistic cluster finders, such as redMaPPer (Rykoff et al. 2014), impose a colour-magnitude cut or a stellar mass cut, which are highly dependent on the red-sequence model or the spectral energy density model. We found that imposing a stellar mass cut of $10.5\log(M/M_{\odot})$ would correspond roughly to the bottom 5 per cent percentile of SDSS detected galaxies (Maraston et al. 2013). However, this drastically decreases the number of galaxies in a halo, with most having $N_{gal}$ in the single digits. As we are interested in the *intrinsic* covariance from the physical properties of the halo, we do not impose additional magnitude or stellar mass cuts. We confirm that, as described in Croton et al. (2016), the galaxy stellar mass distribution at $z = 0$ is consistent with the best-fitting double Schechter function calibrated with low-redshift galaxies from the Galaxy and Mass Assembly (Baldry et al. 2012) down to stellar masses of $M > 10^{8.5} M_{\odot}$.

Fig. 3 illustrates that our $N_{gal}$-mass models, which count the number of galaxies within a physical 3D radius and impose no colour-magnitude cut as redMaPPer does, resemble the general behaviour

<footnote n="1">¹https://www.cosmosim.org</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Measured against the left hand side y-axis are measurements of Cov(ΔΣ, N<sub>gal</sub>|M, z) with 1σ errors and different functional forms using the full model. The functions are classes of ‘Sigmoid’ functions. In all bins, the error function outperforms other functional forms in their DIC parameters, providing good χ<sup>2</sup> values. For M<sub>200c</sub> ∈ [5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>) at z = 0.49 and M<sub>200c</sub> ∈ [2 × 10<sup>14</sup>, 5 × 10<sup>15</sup>) at z = 1.03, the posteriors of the full models do not converge as the size of the covariance is too small. Measured against the right-hand side y-axis are the correlation coefficients r<sub>ΔΣ,N<sub>gal</sub>|M,z</sub> with smoothed bands representing the 1 − σ error. The errors are measured by bootstrap resampling.</figcaption></figure>

of the observed richness-mass relations in terms of both slope and intercept. However, we acknowledge that redMaPPer may suffer from projection effects that artificially inflate the number count of red-sequence galaxies within its aperture because of line-of-sight structures. Additionally, the N<sub>gal</sub> count within the R<sub>vir</sub> radius exceeds that of R<sub>200c</sub> as R<sub>vir</sub> is greater than R<sub>200c</sub>. In the toy model scenario, where we use a constant 1 Mpc h<sup>-1</sup> radius, the slope of the mass-richness relation starts to decrease as the mass increases due to the increasing physical size of the clusters, as expected. The diversity of cluster radii and the resulting variation in the local slope and intercept of the N<sub>gal</sub>-mass relations demonstrate the robustness of our covariance model. In Section 4, we show that different radii/mass

definitions have little impact on the parameters of our covariance model, thus establishing its independence from different reference radii, the definitions of cluster edges, and the resulting richness-mass relations.

4 RESULTS: COVARIANCE SHAPE AND
EVOLUTION

In this section, we report the measurements for our covariance. In Fig. 4, we find an anti-correlation between N<sub>gal</sub> and ΔΣ at small scales across most redshift and mass bins spanned by our data set, which we fit with the best-fitting ‘Sigmoid’ functional form of the
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Measured against the left hand side y-axis are measurements of Cov(ΔΣ, N<sub>gal</sub>|M, z) with 1σ errors and different functional forms using the full model. The functions are classes of ‘Sigmoid’ functions. In all bins, the error function outperforms other functional forms in their DIC parameters, providing good χ<sup>2</sup> values. For M<sub>200c</sub> ∈ [5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>) at z = 0.49 and M<sub>200c</sub> ∈ [2 × 10<sup>14</sup>, 5 × 10<sup>15</sup>) at z = 1.03, the posteriors of the full models do not converge as the size of the covariance is too small. Measured against the right-hand side y-axis are the correlation coefficients r<sub>ΔΣ,N<sub>gal</sub>|M,z</sub> with smoothed bands representing the 1 − σ error. The errors are measured by bootstrap resampling.</figcaption></figure>

of the observed richness-mass relations in terms of both slope and intercept. However, we acknowledge that redMaPPer may suffer from projection effects that artificially inflate the number count of red-sequence galaxies within its aperture because of line-of-sight structures. Additionally, the N<sub>gal</sub> count within the R<sub>vir</sub> radius exceeds that of R<sub>200c</sub> as R<sub>vir</sub> is greater than R<sub>200c</sub>. In the toy model scenario, where we use a constant 1 Mpc h<sup>-1</sup> radius, the slope of the mass-richness relation starts to decrease as the mass increases due to the increasing physical size of the clusters, as expected. The diversity of cluster radii and the resulting variation in the local slope and intercept of the N<sub>gal</sub>-mass relations demonstrate the robustness of our covariance model. In Section 4, we show that different radii/mass

definitions have little impact on the parameters of our covariance model, thus establishing its independence from different reference radii, the definitions of cluster edges, and the resulting richness-mass relations.

4 RESULTS: COVARIANCE SHAPE AND
EVOLUTION

In this section, we report the measurements for our covariance. In Fig. 4, we find an anti-correlation between N<sub>gal</sub> and ΔΣ at small scales across most redshift and mass bins spanned by our data set, which we fit with the best-fitting ‘Sigmoid’ functional form of the
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Evolution of Cov(ΔΣ, N<sub>gal</sub>|M, z) shape parameters with respect to mass and radial binning schemes and N<sub>gal</sub> definition at fixed redshift at z = 0. ΔΣ is binned in equal log-space radial bins in R/R<sub>200c</sub> or R/R<sub>vir</sub>; for each radial binning, the number count of galaxies inside the cluster is given by a constant radius of 1 Mpc h<sup>-1</sup> or R<sub>200c</sub> when binned by R<sub>200c</sub> and R<sub>vir</sub> when binned by R<sub>vir</sub>. We find no strong evolution in the shape or scale of the covariance under different binning schemes or N<sub>gal</sub> definitions. The trend is consistent across different redshift bins and demonstrates the robustness of the covariance under different true richness definitions. The error bars indicate the 1 − σ distribution of the posteriors.</figcaption></figure>

expression

$$
\mathrm{Cov}(\bar{x}) = s\left(\mathrm{erf}\left(\frac{\sqrt{\pi}}{2}\bar{x}\right) + g\right), \qquad (25)
$$

with $x \equiv \log R/R_{200c}$ the log-radius and $\bar{x} \equiv (x - \gamma)/\tau$ the scaled
and offset log-radius. In Appendix A, we offer statistical verification
of the best-fitting functional form.

We first describe the evolution of the covariance in Section 4.1
by binning across the (M, z) bins. Next, in Section 4.2, we present
an alternative binning scheme based on halo peak height that can
provide insight into the dependence of the time formation history of
the covariance scale.

4.1 Binned in (M, z)

Our best-fitting parameters in Table D1 indicate that in 9 out of 12
(M, z) bins, the Cov(ΔΣ, ln N<sub>gal</sub>|M, z) rejects the null-correlation
hypothesis with high statistical significance (p-value < 0.01). How-

ever, in two bins, specifically M<sub>200c</sub> ∈ [5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>) M<sub>☉</sub>h<sup>-1</sup>
at z = 0.49 and M<sub>200c</sub> ∈ [2 × 10<sup>14</sup>, 5 × 10<sup>15</sup>) M<sub>☉</sub>h<sup>-1</sup> at z = 1.03, the
magnitude of the covariance is relatively small compared to the size
of their errors. Consequently, it becomes challenging to constrain the
shape parameters in these two bins, and the covariance is consistent
with the null hypothesis. Furthermore, we exclude the bin M<sub>200c</sub> ∈
[5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>) M<sub>☉</sub>h<sup>-1</sup>, z = 1.03 due to the limited number of
haloes it contains.

Our results suggest that the shape of the covariances can be
accurately described by the full error function. Additionally, for R ≥
R<sub>vir</sub> or R ≥ R<sub>200c</sub>, the covariance aligns with the null-correlation
hypothesis. This alignment is reflected in the fact that all nine
bins with constrained posterior shape have best-fitting g values
within 2σ of g = -1. Deviations from g = -1 can be interpreted
as evidence of disagreements with the Press-Schechter formalism
(Press & Schechter 1974) of spherical collapse haloes, which can be
originated from the presence of anisotropic or non-Gaussian matter
distribution around haloes at large scales (Lokken et al. 2022), or it
can be an indicator of an open-shell model of haloes that allows for
the bulk transfer of baryonic and dark matter in and out of the halo
potential well during the non-linear collapse.

With g = -1 fixed, the reduced error function marginally im-
proves the constraints in most bins. However, with the reduced
model, we can provide posterior constraints for M<sub>200c</sub> ∈ [5 × 10<sup>14</sup>,
1 × 10<sup>15</sup>)M<sub>☉</sub>h<sup>-1</sup> at z = 0.49 and M<sub>200c</sub> ∈ [2 × 10<sup>14</sup>, 5 × 10<sup>15</sup>)M<sub>☉</sub>h<sup>-1</sup>
at z = 1.03, which the full model failed to constrain but with very
loose posterior constraints. The estimated parameters for both the full
and reduced models are presented in Tables D1 and D2, respectively.

To assess the impact of varying the definition of the halo radius
on our measurements of the shape of the covariance, we considered
two factors: the scale dependence of ΔΣ discussed in Sections 2.1.3
and 3.1, and the alteration of the richness-mass relation as shown in
Fig. 3 in Section 3.2. Fig. 5 demonstrates that there is no apparent
evolution of the shape parameters θ ∈ {τ, γ, g} when altering the
scale dependence for ΔΣ or the true richness count. However, we
find marginal 3σ evidence of a difference in the amplitude parameter
of the covariance s when changing the scale normalisation from r<sub>p</sub> =
R/R<sub>200c</sub> to r<sub>p</sub> = R/R<sub>vir</sub> while using the same true richness count. As
haloes exhibit more self-similarity in the inner regions when scaled
by R<sub>200c</sub> (Diemer & Kravtsov 2014), we adopt this as our radius
normalisation and use the number of galaxies enclosed within R<sub>200c</sub>
as our true richness count.

Subsequently, we explored the evolution of the shape parameters
with respect to (M, z) and found no strong mass dependence. How-
ever, we observed a monotonically decreasing redshift dependence
of the amplitude parameter s, as illustrated in Fig. 6. To explain both
the halo mass and the redshift dependence, we used the peak height
of the halo, ν(M, z).

4.2 Binned by peak height

An alternative binning scheme that encapsulates both the halo mass
and redshift information is to bin haloes by the peak height parameter,
defined as

$$
\nu = \frac{\delta_c}{\sigma(R, a)}, \qquad (26)
$$

where δ_c(z) is the collapse overdensity at which gravitational col-
lapses enter the non-linear regime and σ(R, a) is the smoothing scale
seen in equation (B2) at the radius of the cluster. For an Einstein-de
Sitter universe (Ω_m = 1, Ω_Λ = 0) δ_c ≈ 1.686 at the epoch of collapse
and is weakly dependent on cosmology and redshift (Percival 2005).
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Evolution of Cov(ΔΣ, N<sub>gal</sub>|M, z) shape parameters with respect to mass and radial binning schemes and N<sub>gal</sub> definition at fixed redshift at z = 0. ΔΣ is binned in equal log-space radial bins in R/R<sub>200c</sub> or R/R<sub>vir</sub>; for each radial binning, the number count of galaxies inside the cluster is given by a constant radius of 1 Mpc h<sup>-1</sup> or R<sub>200c</sub> when binned by R<sub>200c</sub> and R<sub>vir</sub> when binned by R<sub>vir</sub>. We find no strong evolution in the shape or scale of the covariance under different binning schemes or N<sub>gal</sub> definitions. The trend is consistent across different redshift bins and demonstrates the robustness of the covariance under different true richness definitions. The error bars indicate the 1 − σ distribution of the posteriors.</figcaption></figure>

expression

$$
\mathrm{Cov}(\bar{x}) = s\left(\mathrm{erf}\left(\frac{\sqrt{\pi}}{2}\bar{x}\right) + g\right), \qquad (25)
$$

with $x \equiv \log R/R_{200c}$ the log-radius and $\bar{x} \equiv (x - \gamma)/\tau$ the scaled
and offset log-radius. In Appendix A, we offer statistical verification
of the best-fitting functional form.

We first describe the evolution of the covariance in Section 4.1
by binning across the (M, z) bins. Next, in Section 4.2, we present
an alternative binning scheme based on halo peak height that can
provide insight into the dependence of the time formation history of
the covariance scale.

4.1 Binned in (M, z)

Our best-fitting parameters in Table D1 indicate that in 9 out of 12
(M, z) bins, the Cov(ΔΣ, ln N<sub>gal</sub>|M, z) rejects the null-correlation
hypothesis with high statistical significance (p-value < 0.01). How-

ever, in two bins, specifically M<sub>200c</sub> &isin; [5 &times; 10<sup>14</sup>, 1 &times; 10<sup>15</sup>) M<sub>&odot;</sub>h<sup>-1</sup>
at z = 0.49 and M<sub>200c</sub> &isin; [2 &times; 10<sup>14</sup>, 5 &times; 10<sup>15</sup>) M<sub>&odot;</sub>h<sup>-1</sup> at z = 1.03, the
magnitude of the covariance is relatively small compared to the size
of their errors. Consequently, it becomes challenging to constrain the
shape parameters in these two bins, and the covariance is consistent
with the null hypothesis. Furthermore, we exclude the bin M<sub>200c</sub> &isin;
[5 &times; 10<sup>14</sup>, 1 &times; 10<sup>15</sup>) M<sub>&odot;</sub>h<sup>-1</sup>, z = 1.03 due to the limited number of
haloes it contains.

Our results suggest that the shape of the covariances can be
accurately described by the full error function. Additionally, for R ≥
R<sub>vir</sub> or R ≥ R<sub>200c</sub>, the covariance aligns with the null-correlation
hypothesis. This alignment is reflected in the fact that all nine
bins with constrained posterior shape have best-fitting g values
within 2σ of g = -1. Deviations from g = -1 can be interpreted
as evidence of disagreements with the Press-Schechter formalism
(Press & Schechter 1974) of spherical collapse haloes, which can be
originated from the presence of anisotropic or non-Gaussian matter
distribution around haloes at large scales (Lokken et al. 2022), or it
can be an indicator of an open-shell model of haloes that allows for
the bulk transfer of baryonic and dark matter in and out of the halo
potential well during the non-linear collapse.

With g = -1 fixed, the reduced error function marginally im-
proves the constraints in most bins. However, with the reduced
model, we can provide posterior constraints for M<sub>200c</sub> &isin; [5 &times; 10<sup>14</sup>,
1 &times; 10<sup>15</sup>)M<sub>&odot;</sub>h<sup>-1</sup> at z = 0.49 and M<sub>200c</sub> &isin; [2 &times; 10<sup>14</sup>, 5 &times; 10<sup>15</sup>)M<sub>&odot;</sub>h<sup>-1</sup>
at z = 1.03, which the full model failed to constrain but with very
loose posterior constraints. The estimated parameters for both the full
and reduced models are presented in Tables D1 and D2, respectively.

To assess the impact of varying the definition of the halo radius
on our measurements of the shape of the covariance, we considered
two factors: the scale dependence of ΔΣ discussed in Sections 2.1.3
and 3.1, and the alteration of the richness-mass relation as shown in
Fig. 3 in Section 3.2. Fig. 5 demonstrates that there is no apparent
evolution of the shape parameters θ ∈ {τ, γ, g} when altering the
scale dependence for ΔΣ or the true richness count. However, we
find marginal 3σ evidence of a difference in the amplitude parameter
of the covariance s when changing the scale normalisation from r<sub>p</sub> =
R/R<sub>200c</sub> to r<sub>p</sub> = R/R<sub>vir</sub> while using the same true richness count. As
haloes exhibit more self-similarity in the inner regions when scaled
by R<sub>200c</sub> (Diemer & Kravtsov 2014), we adopt this as our radius
normalisation and use the number of galaxies enclosed within R<sub>200c</sub>
as our true richness count.

Subsequently, we explored the evolution of the shape parameters
with respect to (M, z) and found no strong mass dependence. How-
ever, we observed a monotonically decreasing redshift dependence
of the amplitude parameter s, as illustrated in Fig. 6. To explain both
the halo mass and the redshift dependence, we used the peak height
of the halo, ν(M, z).

4.2 Binned by peak height

An alternative binning scheme that encapsulates both the halo mass
and redshift information is to bin haloes by the peak height parameter,
defined as

$$
\nu = \frac{\delta_c}{\sigma(R, a)}, \qquad (26)
$$

where δ_c(z) is the collapse overdensity at which gravitational col-
lapses enter the non-linear regime and σ(R, a) is the smoothing scale
seen in equation (B2) at the radius of the cluster. For an Einstein-de
Sitter universe (Ω_m = 1, Ω_Λ = 0) δ_c ≈ 1.686 at the epoch of collapse
and is weakly dependent on cosmology and redshift (Percival 2005).
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Evolution of Cov(ΔΣ, N<sub>gal</sub>|M, z) shape parameters of the error function with respect to mass and redshift, binned in units of R<sub>200c</sub> and with N<sub>gal</sub> taken to be the number of clusters inside the R<sub>200c</sub> radius of the cluster. There is no strong dependence of τ, γ, and g with respect to mass and redshift and a strong monotonically decreasing s with respect to redshift. At M<sub>200c</sub> ∈ [5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>), z = 0.49 and M<sub>200c</sub> ∈ [2 × 10<sup>14</sup>, 5 × 10<sup>15</sup>), z = 1.03 the covariance is consistent with null at p = 0.01 and p = 0.05 levels, respectively. The error bars indicate the 1 − σ distribution of the posteriors.</figcaption></figure>

σ(R, a) scales as σ(M, a) = σ(M, a = 1)D<sub>+0</sub>(a) at the linear collapse regime, where D<sub>+0</sub>(a) ≡ D<sub>+</sub>(a)/D<sub>+</sub>(a = 1). Here, D<sub>+</sub>(a) is the linear growth factor defined as

$$D_{+}(a) = \frac{5\Omega_{M}}{2} E(a) \int_{0}^{a} \frac{da'}{[a'E(a')]^{3}}, \qquad (27)$$

for a ΛCDM cosmology, where E(a) ≡ H(a)/H<sub>0</sub> is the normalized Hubble parameter. σ(R, z) depends strongly on redshift, and hence, the peak height v strongly depends on the halo radius and the redshift of non-linear collapse.

The peak height has been adopted to simplify the mass and redshift dependence in various halo properties, such as halo concentration (Prada et al. 2012) and halo triaxiality (Allgood et al. 2006). Here, we explore whether the peak height can serve as a universal parameter to explain the scale and shape of Cov(ΔΣ, ln N<sub>gal</sub>|M, z). We bin the haloes into deciles of v and set posterior constraints on the shape of the covariance using our erf model in the full model case. In the highest decile (90–100 per cent percentile), we reject the null-correlation hypothesis at the p = 0.01 level, but due to the size of

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Best-fitting ‘full’ error function model for Cov(ΔΣ, ln N<sub>gal</sub>|M, z) when binned in deciles of halo peak height. The first nine bins reject the null hypothesis at a p < 0.01 level, and the highest decile rejects the null hypothesis at a p < 0.05 level. We can provide posterior constraints for all bins in peak height except for the one with the highest peak height value.</figcaption></figure>

the error bars, the shape of the parameters τ, γ, and g and largely unconstrained and s = 10<sup>12</sup> × 0.16<sup>+0.34</sup><sub>-0.12</sub>. Due to the large degeneracy, we exclude the highest decile from our data set and limit the range of our model to v ∈ [1.57, 3.40), which spans 0–90 per cent of our sample set. The large error bars may be due to the fact that the halo abundance as a function of v falls precipitously around v ~ 4, so the highest decile spans a wide tail of high v ∈ [3.4, 4.6). The plots in Fig. 7 are the best-fitting templates when binned by peak height, and Fig. 8 shows the best-fitting parameters as a function of peak height. We do not see a strong dependence on the peak height for τ, γ, and g. For s, its dependence on v can be modelled as a log-linear relation of the form

$$\log_{10}(s) = C_s + \alpha v, \qquad (28)$$

with C<sub>s</sub> = 13.07<sup>+0.26</sup><sub>-0.26</sub> and α = −0.44<sup>+0.11</sup><sub>-0.11</sub>. At the highest decile, the s = 10<sup>12</sup> × 0.16<sup>+0.34</sup><sub>-0.12</sub> falls within the 1σ confidence band of the log-linear fit. Compared to the first nine deciles, the fit yields a χ<sup>2</sup> p-value of 0.73. The negative slope between s and v indicates that more massive haloes at the cosmic era of their formation exhibit a lesser anticorrelation between ΔΣ and ln N<sub>gal</sub>.

## 5 IMPACT OF COV(ΔΣ, LN N<sub>GAL</sub>|M, z) ON WEAK LENSING MEASUREMENTS

To assess the impact of Cov(ΔΣ, ln N<sub>gal</sub>|M, z) on the scaling relation ⟨ΔΣ|N<sub>gal</sub>, z⟩, we utilize equation (17) for the first-order correction and equation (18) for the second order. The mean mass of the haloes in each (M, z) bin is chosen as the pivot mass around which the HMF is Taylor expanded, and the intercept πN<sub>gal</sub> and slope αN<sub>gal</sub> for the richness-mass scaling relation shown in equation (9) are computed locally at the pivot point in each bin (M, z) bin. Our mock data, binned in R<sub>200c</sub>, yield results that are consistent with the global richness-mass relation found in the literature (Bocquet et al. 2016; Costanzi et al. 2021; To et al. 2021), as shown in Fig. 3.

For the correction terms, we adopt the Tinker et al. (2008) halo mass function as our nominal model and compute the numeric log-derivative for values of γ<sub>1</sub> and γ<sub>2</sub>, the log-slope and curvature of
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Evolution of Cov(ΔΣ, N<sub>gal</sub>|M, z) shape parameters of the error function with respect to mass and redshift, binned in units of R<sub>200c</sub> and with N<sub>gal</sub> taken to be the number of clusters inside the R<sub>200c</sub> radius of the cluster. There is no strong dependence of τ, γ, and g with respect to mass and redshift and a strong monotonically decreasing s with respect to redshift. At M<sub>200c</sub> ∈ [5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>), z = 0.49 and M<sub>200c</sub> ∈ [2 × 10<sup>14</sup>, 5 × 10<sup>15</sup>), z = 1.03 the covariance is consistent with null at p = 0.01 and p = 0.05 levels, respectively. The error bars indicate the 1 − σ distribution of the posteriors.</figcaption></figure>

σ(R, a) scales as σ(M, a) = σ(M, a = 1)D<sub>+0</sub>(a) at the linear collapse regime, where D<sub>+0</sub>(a) ≡ D<sub>+</sub>(a)/D<sub>+</sub>(a = 1). Here, D<sub>+</sub>(a) is the linear growth factor defined as

$$D_{+}(a) = \frac{5\Omega_{M}}{2} E(a) \int_{0}^{a} \frac{da'}{[a'E(a')]^{3}}, \qquad (27)$$

for a ΛCDM cosmology, where E(a) ≡ H(a)/H<sub>0</sub> is the normalized Hubble parameter. σ(R, z) depends strongly on redshift, and hence, the peak height v strongly depends on the halo radius and the redshift of non-linear collapse.

The peak height has been adopted to simplify the mass and redshift dependence in various halo properties, such as halo concentration (Prada et al. 2012) and halo triaxiality (Allgood et al. 2006). Here, we explore whether the peak height can serve as a universal parameter to explain the scale and shape of Cov(ΔΣ, ln N<sub>gal</sub>|M, z). We bin the haloes into deciles of v and set posterior constraints on the shape of the covariance using our erf model in the full model case. In the highest decile (90–100 per cent percentile), we reject the null-correlation hypothesis at the p = 0.01 level, but due to the size of

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Best-fitting ‘full’ error function model for Cov(ΔΣ, ln N<sub>gal</sub>|M, z) when binned in deciles of halo peak height. The first nine bins reject the null hypothesis at a p < 0.01 level, and the highest decile rejects the null hypothesis at a p < 0.05 level. We can provide posterior constraints for all bins in peak height except for the one with the highest peak height value.</figcaption></figure>

the error bars, the shape of the parameters τ, γ, and g and largely unconstrained and s = 10<sup>12</sup> × 0.16<sup>+0.34</sup><sub>-0.12</sub>. Due to the large degeneracy, we exclude the highest decile from our data set and limit the range of our model to v ∈ [1.57, 3.40), which spans 0–90 per cent of our sample set. The large error bars may be due to the fact that the halo abundance as a function of v falls precipitously around v ~ 4, so the highest decile spans a wide tail of high v ∈ [3.4, 4.6). The plots in Fig. 7 are the best-fitting templates when binned by peak height, and Fig. 8 shows the best-fitting parameters as a function of peak height. We do not see a strong dependence on the peak height for τ, γ, and g. For s, its dependence on v can be modelled as a log-linear relation of the form

$$\log_{10}(s) = C_s + \alpha v, \qquad (28)$$

with C<sub>s</sub> = 13.07<sup>+0.26</sup><sub>-0.26</sub> and α = −0.44<sup>+0.11</sup><sub>-0.11</sub>. At the highest decile, the s = 10<sup>12</sup> × 0.16<sup>+0.34</sup><sub>-0.12</sub> falls within the 1σ confidence band of the log-linear fit. Compared to the first nine deciles, the fit yields a χ<sup>2</sup> p-value of 0.73. The negative slope between s and v indicates that more massive haloes at the cosmic era of their formation exhibit a lesser anticorrelation between ΔΣ and ln N<sub>gal</sub>.

## 5 IMPACT OF COV(ΔΣ, LN N<sub>GAL</sub>|M, z) ON WEAK LENSING MEASUREMENTS

To assess the impact of Cov(ΔΣ, ln N<sub>gal</sub>|M, z) on the scaling relation ⟨ΔΣ|N<sub>gal</sub>, z⟩, we utilize equation (17) for the first-order correction and equation (18) for the second order. The mean mass of the haloes in each (M, z) bin is chosen as the pivot mass around which the HMF is Taylor expanded, and the intercept πN<sub>gal</sub> and slope αN<sub>gal</sub> for the richness-mass scaling relation shown in equation (9) are computed locally at the pivot point in each bin (M, z) bin. Our mock data, binned in R<sub>200c</sub>, yield results that are consistent with the global richness-mass relation found in the literature (Bocquet et al. 2016; Costanzi et al. 2021; To et al. 2021), as shown in Fig. 3.

For the correction terms, we adopt the Tinker et al. (2008) halo mass function as our nominal model and compute the numeric log-derivative for values of γ<sub>1</sub> and γ<sub>2</sub>, the log-slope and curvature of
<figure><img src="image_11.png" /><figcaption>**Figure 8.** The evolution of shape parameters for $Cov(\Delta\Sigma, \ln N_{\text{gal}}|M, z)$ binned in deciles of the peak height $v$, excluding the highest decile. The parameters $\tau, \gamma$, and $g$ show little dependency with $v$ while the amplitude $s$ exhibits a log-linear relationship with $v$ of the form shown in equation (28). The mean $g$ is consistent with $-1$. The horizontal teal bands fill the $1\sigma$ range around the mean, and the pink line is the best log-linear fit between $s$ and $v$ with $1\sigma$ confidence bands.</figcaption></figure>

the halo mass function around the pivot mass. We compare the Tinker mass function results with others, including Watson et al. (2013), Bocquet et al. (2016), Despali et al. (2016), and find that the difference is subdominant to the first-order correction, which is at a ~1 per cent level at small scales, as shown in Fig. 9.

To estimate the mass bias in each bin, we stack $\langle\Delta\Sigma|M, z\rangle$ and model the profiles as if they were individual haloes with a mean mass, redshift and concentration as described in equations (2)–(7). We assume NFW profile using the concentration–mass model of Diemer & Joyce (2019) in the one-halo regime. The two-halo regime should not be affected, as the covariance is consistent with zero at $R \gtrsim R_{200c}$. We convert the 3D overdensity of the modelled halo $\xi_{hm}$ to $\Delta\Sigma$ using equations (2) and (5), and then apply the first-order correction in equation (17). Using a Monte Carlo method we obtain the expected mass with and without this correction and report the change in the mean halo mass with this correction for each $(M, z)$ bin. As shown in Fig. 9, we find that adding the correction leads to an *upward correction* of the stacked halo mass of approximately $\delta M/M \sim 2-3$ per cent for most $(M, z)$ bins.

# 6 EXPLAINING THE COVARIANCE

## 6.1 Secondary halo parameter dependence of $\ln N_{\text{gal}}$

We employed a multivariable linear regression model to determine the best fit when incorporating secondary properties in the regression.

Initially, considering the full set of parameters listed in Table 3, we applied a backward modelling scheme to identify the relevant parameters of interest. Details of this process can be found in Appendix E, which led to the selection of the following secondary halo parameters for our model: $\Pi \subset \{\Gamma_{2\text{dyn}}, a_{1/2}, c_{\text{vir}}, T/|U|, X_{\text{off}}\}$. The resulting model demonstrated good explanatory power, as indicated by a high $R^2$ coefficient. Additionally, the model passed various tests, including variance inflation, global F-statistic, partial F-statistic, T-statistic, scatter heteroscedasticity, and scatter normality in most cases. Specifically, through a comparison of $F_{\text{partial}}$ values, we found that richness could be modelled by a multilinear equation involving all secondary halo parameters. Further information can be found in Table E1, where the F-statistic demonstrates that all parameters are statistically significant. Only when considered collectively can they accurately reflect the dependence of richness on halo formation history.

To establish informative priors for upcoming weak lensing surveys such as Hyper Suprime-Cam (HSC) and LSST, we examined whether the dependence of $N_{\text{gal}}$ on secondary halo properties, as inferred from the slope $\beta_{N_{\text{gal}}}$, aligns with arguments based on halo formation physics. We expected that $\beta_{N_{\text{gal}},c_{\text{vir}}}$ resulting from the formation of satellite galaxies (equivalent to $N_{\text{gal}} - 1$ in the presence of a central galaxy) within haloes would exhibit a negative relation, i.e. $\beta_{N_{\text{gal}},c_{\text{vir}}} < 0$. Simulation-based studies have suggested that early forming haloes possess higher concentrations (Wechsler et al. 2002), and correspondingly, high-concentration haloes (which form early) have fewer satellite galaxies due to galaxy mergers within the haloes (Zentner et al. 2005). This effect is known as galaxy assembly bias (Wechsler & Tinker 2018) – the change in galaxy properties inside a halo at fixed mass due to the halo formation history. There is marginal evidence of the existence of assembly bias from recent observations using galaxy clustering techniques (Zentner et al. 2019; Wang et al. 2022), as well as measurements of the magnitude gap between the brightest central galaxy (BCG) and a neighbouring galaxy as a proxy for formation time (Hearin et al. 2013; Golden-Marx & Miller 2018; Farahi, Ho & Trac 2020).

As noted in Table E1, the signs of $\beta_{N_{\text{gal}},i}$ for the remaining parameters $i \in \{a_{1/2}, T/|U|, \Gamma_{2\text{dyn}}, X_{\text{off}}\}$ align with our expectations of assembly bias in most bins – late-forming clusters undergo more rapid mass accretion (higher $\Gamma_{2\text{dyn}}$) and are less virialized (higher $T/|U|$), and because they also from the galaxy assembly bias mentioned above are richer in galaxy number counts when conditioned on the mass, we expect a positive partial slope $\beta_{N_{\text{gal}},\Gamma_{2\text{dyn}}}$ and $\beta_{N_{\text{gal}},T/|U|}$. The case for $a_{1/2}$ and $X_{\text{off}}$ is more complicated. Under the isolated formation of haloes $a_{1/2}$ and $X_{\text{off}}$ would be smaller for earlier forming haloes due to the monotonic mass accretion and relaxation of haloes over long time-scales. However, as haloes undergo mergers and tidal stripping the monotonicity of the parameters over time is not guaranteed. Therefore, we see a mixture of positive and negative partial slopes $\beta_{N_{\text{gal}},a_{1/2}}$ and $\beta_{N_{\text{gal}},X_{\text{off}}}$ in these cases. To describe the physical mechanisms on a case-by-case basis would require that we probe into the halo merger tree history of individual haloes.

In this paper, we take a closer look at the sign of $\beta_{N_{\text{gal}},c_{\text{vir}}}$ and observe that while the partial scope matches our expectations in most bins, in some mass bins at medium and high redshifts it changes signs from negative at lower redshifts to positive at higher ones. While we observe a diminishing impact of secondary halo properties on richness (indicated by a smaller absolute value for $\beta_{N_{\text{gal}},c_{\text{vir}}}$), the reversal of the coefficient's sign cannot be solely attributed to statistical fluctuations around zero, as some values are inconsistent with zero at levels exceeding $3\sigma$.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** The evolution of shape parameters for $Cov(\Delta\Sigma, \ln N_{\text{gal}}|M, z)$ binned in deciles of the peak height $v$, excluding the highest decile. The parameters $\tau, \gamma$, and $g$ show little dependency with $v$ while the amplitude $s$ exhibits a log-linear relationship with $v$ of the form shown in equation (28). The mean $g$ is consistent with $-1$. The horizontal teal bands fill the $1\sigma$ range around the mean, and the pink line is the best log-linear fit between $s$ and $v$ with $1\sigma$ confidence bands.</figcaption></figure>

the halo mass function around the pivot mass. We compare the Tinker mass function results with others, including Watson et al. (2013), Bocquet et al. (2016), Despali et al. (2016), and find that the difference is subdominant to the first-order correction, which is at a ~1 per cent level at small scales, as shown in Fig. 9.

To estimate the mass bias in each bin, we stack $\langle\Delta\Sigma|M, z\rangle$ and model the profiles as if they were individual haloes with a mean mass, redshift and concentration as described in equations (2)–(7). We assume NFW profile using the concentration–mass model of Diemer & Joyce (2019) in the one-halo regime. The two-halo regime should not be affected, as the covariance is consistent with zero at $R \gtrsim R_{200c}$. We convert the 3D overdensity of the modelled halo $\xi_{hm}$ to $\Delta\Sigma$ using equations (2) and (5), and then apply the first-order correction in equation (17). Using a Monte Carlo method we obtain the expected mass with and without this correction and report the change in the mean halo mass with this correction for each $(M, z)$ bin. As shown in Fig. 9, we find that adding the correction leads to an *upward correction* of the stacked halo mass of approximately $\delta M/M \sim 2-3$ per cent for most $(M, z)$ bins.

# 6 EXPLAINING THE COVARIANCE

## 6.1 Secondary halo parameter dependence of $\ln N_{\text{gal}}$

We employed a multivariable linear regression model to determine the best fit when incorporating secondary properties in the regression.

Initially, considering the full set of parameters listed in Table 3, we applied a backward modelling scheme to identify the relevant parameters of interest. Details of this process can be found in Appendix E, which led to the selection of the following secondary halo parameters for our model: $\Pi \subset \{\Gamma_{2\text{dyn}}, a_{1/2}, c_{\text{vir}}, T/|U|, X_{\text{off}}\}$. The resulting model demonstrated good explanatory power, as indicated by a high $R^2$ coefficient. Additionally, the model passed various tests, including variance inflation, global F-statistic, partial F-statistic, T-statistic, scatter heteroscedasticity, and scatter normality in most cases. Specifically, through a comparison of $F_{\text{partial}}$ values, we found that richness could be modelled by a multilinear equation involving all secondary halo parameters. Further information can be found in Table E1, where the F-statistic demonstrates that all parameters are statistically significant. Only when considered collectively can they accurately reflect the dependence of richness on halo formation history.

To establish informative priors for upcoming weak lensing surveys such as Hyper Suprime-Cam (HSC) and LSST, we examined whether the dependence of $N_{\text{gal}}$ on secondary halo properties, as inferred from the slope $\beta_{N_{\text{gal}}}$, aligns with arguments based on halo formation physics. We expected that $\beta_{N_{\text{gal}},c_{\text{vir}}}$ resulting from the formation of satellite galaxies (equivalent to $N_{\text{gal}} - 1$ in the presence of a central galaxy) within haloes would exhibit a negative relation, i.e. $\beta_{N_{\text{gal}},c_{\text{vir}}} < 0$. Simulation-based studies have suggested that early forming haloes possess higher concentrations (Wechsler et al. 2002), and correspondingly, high-concentration haloes (which form early) have fewer satellite galaxies due to galaxy mergers within the haloes (Zentner et al. 2005). This effect is known as galaxy assembly bias (Wechsler & Tinker 2018) – the change in galaxy properties inside a halo at fixed mass due to the halo formation history. There is marginal evidence of the existence of assembly bias from recent observations using galaxy clustering techniques (Zentner et al. 2019; Wang et al. 2022), as well as measurements of the magnitude gap between the brightest central galaxy (BCG) and a neighbouring galaxy as a proxy for formation time (Hearin et al. 2013; Golden-Marx & Miller 2018; Farahi, Ho & Trac 2020).

As noted in Table E1, the signs of $\beta_{N_{\text{gal}},i}$ for the remaining parameters $i \in \{a_{1/2}, T/|U|, \Gamma_{2\text{dyn}}, X_{\text{off}}\}$ align with our expectations of assembly bias in most bins – late-forming clusters undergo more rapid mass accretion (higher $\Gamma_{2\text{dyn}}$) and are less virialized (higher $T/|U|$), and because they also from the galaxy assembly bias mentioned above are richer in galaxy number counts when conditioned on the mass, we expect a positive partial slope $\beta_{N_{\text{gal}},\Gamma_{2\text{dyn}}}$ and $\beta_{N_{\text{gal}},T/|U|}$. The case for $a_{1/2}$ and $X_{\text{off}}$ is more complicated. Under the isolated formation of haloes $a_{1/2}$ and $X_{\text{off}}$ would be smaller for earlier forming haloes due to the monotonic mass accretion and relaxation of haloes over long time-scales. However, as haloes undergo mergers and tidal stripping the monotonicity of the parameters over time is not guaranteed. Therefore, we see a mixture of positive and negative partial slopes $\beta_{N_{\text{gal}},a_{1/2}}$ and $\beta_{N_{\text{gal}},X_{\text{off}}}$ in these cases. To describe the physical mechanisms on a case-by-case basis would require that we probe into the halo merger tree history of individual haloes.

In this paper, we take a closer look at the sign of $\beta_{N_{\text{gal}},c_{\text{vir}}}$ and observe that while the partial scope matches our expectations in most bins, in some mass bins at medium and high redshifts it changes signs from negative at lower redshifts to positive at higher ones. While we observe a diminishing impact of secondary halo properties on richness (indicated by a smaller absolute value for $\beta_{N_{\text{gal}},c_{\text{vir}}}$), the reversal of the coefficient's sign cannot be solely attributed to statistical fluctuations around zero, as some values are inconsistent with zero at levels exceeding $3\sigma$.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** The percent level change in stacked ΔΣ measurements after including the covariance terms in equations (17) and (18) as denoted by ΔΣcov and without applying corrections as denoted by ΔΣfd. The slope and curvature of the halo mass function are calculated numerically from the Tinker et al. (2008) halo mass function in our nominal correction. The errors are taken from bootstrapped errors of the covariance. We compare the results with first-order corrections from other halo mass functions using Watson et al. (2013); Bocquet et al. (2016); Despali et al. (2016). We find that the percentile difference in ΔΣ far exceeds the uncertainty in the choice of halo mass functions, and that second-order corrections are subdominant to the first-order correction itself, which is at a ~ 1 per cent level at small scales for ΔΣ and propagates into an *upward correction* of stacked halo mass of δM/M ~ 2 – 3 per cent for most bins *after* applying the correction.</figcaption></figure>

This issue can be attributed to the effect of major mergers on concentration. Recent studies (Ludlow et al. 2012; Wang et al. 2022; Lee et al. 2023) have shown that haloes, during major merger events, experience a transient fluctuation in concentration before returning to the mean relation over a time period slightly less than the dynamical time of the halo. The measured concentration spike during major

mergers, particularly prominent at higher redshifts, could explain a positive $\beta_{[N_{gal}, c_{vir}]}$.

To test this hypothesis, we employ a toy model that divides haloes in each (M, z) bin based on the median $\Gamma_{1dyn}$ into low-$\Gamma_{1dyn}$ and high-$\Gamma_{1dyn}$ subsamples. Given the time-scale of mergers to be roughly the dynamical time of the halo, we choose $\Gamma_{1dyn}$ as a good proxy
<figure><img src="image_10.png" /><figcaption>**Figure 9.** The percent level change in stacked $\Delta\Sigma$ measurements after including the covariance terms in equations (17) and (18) as denoted by $\Delta\Sigma_{cov}$ and without applying corrections as denoted by $\Delta\Sigma_{fid}$. The slope and curvature of the halo mass function are calculated numerically from the Tinker et al. (2008) halo mass function in our nominal correction. The errors are taken from bootstrapped errors of the covariance. We compare the results with first-order corrections from other halo mass functions using Watson et al. (2013); Bocquet et al. (2016); Despali et al. (2016). We find that the percentile difference in $\Delta\Sigma$ far exceeds the uncertainty in the choice of halo mass functions, and that second-order corrections are subdominant to the first-order correction itself, which is at a ~1 per cent level at small scales for $\Delta\Sigma$ and propagates into an *upward correction* of stacked halo mass of $\delta M/M \sim 2-3$ per cent for most bins *after* applying the correction.</figcaption></figure>

This issue can be attributed to the effect of major mergers on concentration. Recent studies (Ludlow et al. 2012; Wang et al. 2022; Lee et al. 2023) have shown that haloes, during major merger events, experience a transient fluctuation in concentration before returning to the mean relation over a time period slightly less than the dynamical time of the halo. The measured concentration spike during major

mergers, particularly prominent at higher redshifts, could explain a positive $\beta_{[N_{gal}, c_{vir}]}$.

To test this hypothesis, we employ a toy model that divides haloes in each $(M, z)$ bin based on the median $\Gamma_{1dyn}$ into low-$\Gamma_{1dyn}$ and high-$\Gamma_{1dyn}$ subsamples. Given the time-scale of mergers to be roughly the dynamical time of the halo, we choose $\Gamma_{1dyn}$ as a good proxy
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Residual log-richness versus concentration relation in subsets of halo mass accretion rate (MAR). The figure consists of three panels – left, middle, and right panels corresponding to $z = 0$, $z = 0.49$, and $z = 1.03$, respectively, all on a benchmark mass bin of $M_{200c} \in [5 \times 10^{13}, 1 \times 10^{14}) M_{\odot} h^{-1}$. The sample is split into low and high $\Gamma_{1dyn}$ based on their median values. The scatter plot illustrates the data points, while the shaded regions show the best-fitting linear fit with $1\sigma$ confidence interval for the main sample and each subsample. At $z = 0$, the richness–concentration relation exhibits a negative slope, consistent with our expectations of halo formation physics. The slopes for the low and high $\Gamma_{1dyn}$ subsamples diverge due to the negative correlation between concentration and MAR. However, at $z = 0.49$ and $z = 1.03$, the slopes for the entire sample and/or the high $\Gamma_{1dyn}$ subsample become positive, contrary to our observations of the richness-concentration relation. In contrast, the low $\Gamma_{1dyn}$ subsample still shows a negative slope. These findings suggest that at medium to high redshifts, a subset of unrelaxed and recently merged haloes with high MAR could elevate the concentration from its expected value at hydrostatic equilibrium.</figcaption></figure>

for potential merger events even though this parameter is excluded in the final linear regression model due to multicollinearity (see Appendix E).

Fig. 10 displays the halo concentration plotted against the richness residuals, separated by $\Gamma_{1dyn}$, at benchmark bins of $M_{200c} \in [5 \times 10^{13}, 1 \times 10^{14}) M_{\odot} h^{-1}$ at three different redshift snapshots of $z = 0, 0.49, 1.03$. At $z = 0$, we observe a negative slope as expected from halo formation physics for both low-$\Gamma_{1dyn}$ and high-$\Gamma_{1dyn}$ subsamples, as well as for the overall sample. Furthermore, we observe a change in the slope between the low-$\Gamma_{1dyn}$ and high-$\Gamma_{1dyn}$ subsamples, which can be explained by the gradual increase (or decrease) in concentration ($\Gamma_{1dyn}$) over time, even without major merger events (Wechsler et al. 2002; Zhao et al. 2003; Lu et al. 2006). At redshifts of $z = 0.49, 1.03$, we observe a positive slope in the overall and/or high-$\Gamma_{1dyn}$ samples, which contradicts the scaling relations between HOD and concentration in models that track their gradual evolution over $T \gg T_{1dyn}$. However, in the presence of major mergers, when $\Gamma_{1dyn}$ is significantly enhanced, the halo concentration may also experience a transient spike after the merger. The deviation from hydrostatic equilibrium provides a plausible explanation for a positive $\beta\{N_{gal}, c_{vir}\}$, which could be fully tested on MDPL2 through the reconstruction of halo merger trees, an analysis beyond the scope of this paper.

## 6.2 Secondary halo parameter dependence of $\Delta\Sigma$

In this section, we employ a multilinear regression approach to model the lensing signal, similar to the methodology described in Section 6.1. We extend this approach to the model $P(\Delta\Sigma|N_{gal}(\Pi), M, z)$ as a linear function of $\Pi$. Upon analysing different $(M, z, r_p)$ bins, we observe that the reduced parameters $\Pi \subset \{a_{1/2}, c_{vir}, T/U, \Gamma_{2dyn}, X_{off}\}$ pass the variance inflation factor (VIF) test for multicollinearity, or in other words we showed that the variance is not inflated and thereby made less reliable in the case that the secondary halo parameters in the full model are highly correlated. As with the

case for the lensing signal, this indicates that $\ln N_{gal}$ can be described without redundancy by a linear decomposition of these reduced parameters. Furthermore, most bins exhibit homoscedasticity, as confirmed by passing the Breusch–Pagan Lagrange multiplier test. This implies that the scatter terms $\sigma_{\Delta\Sigma}$ and $\sigma_{\Pi_i}$ remain constant within each bin, with a few exceptions. Lastly, the scatter $\sigma_{\Delta\Sigma|N_{gal}}$ in most bins (with a few exceptions) meets the criteria of the Shapiro–Wilk test for Gaussianity, suggesting that the distribution closely resembles a Gaussian distribution.

The multilinear regression is a good fit to the conditioned lensing signal if we assume that $P(\ln N_{gal}|M, z)$ and $P(\Delta\Sigma|M, z)$ can be modelled with a normal distribution (e.g. Anbajagane et al. 2020; Costanzi et al. 2021; To et al. 2021). In this case $P(\Delta\Sigma|N_{gal}, M, z)$ is a multilinear equation with respect to the secondary halo parameters with mean

$$ \langle \Delta\Sigma | N_{gal}, M, z \rangle = \langle \Delta\Sigma | N_{gal}, M, z \rangle + C_1 \sigma_{\Delta\Sigma} \left( \sum_i \frac{\beta_{N_{gal},i}}{\sigma_{\Pi_i}} \rho_{\Delta\Sigma-\Pi_i} \times (\Pi_i - (\Pi_i | M, z)) \right) \quad (29) $$

and is normally distributed around the mean with variance

$$ \sigma_{\Delta\Sigma|\ln N_{gal}}^2 = \sigma_0^2 + C_2 \sum_i \beta_{N_{gal},i}^2 \sigma_{\Pi_i}^2 (1 - \rho_{\Delta\Sigma-\Pi_i}^2) + C_3 \sum_{i,j} \rho_{\Pi_i-\Pi_j} \sigma_{\Pi_i} \sigma_{\Pi_j} \quad (30) $$

The parameters $C_1, C_2, C_3$, and $\sigma_0$ can be explicitly derived where $P(\ln N_{gal}|M, z)$ and $P(\Delta\Sigma|M, z)$ are known, but the exact values are not essential for this paper. We refer the reader to Appendix F for derivations of equations (29) and (30).

We note that only in bins of $R \lesssim R_{200c}$ do the multilinear regression models pass the global F-statistic test and the T-statistic test for each parameter. This result suggests that, at $R \gtrsim R_{200c}$, we find little correlation between $\Delta\Sigma$ and $\Pi_i$. Because the scatter still passes
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Residual log-richness versus concentration relation in subsets of halo mass accretion rate (MAR). The figure consists of three panels – left, middle, and right panels corresponding to $z = 0$, $z = 0.49$, and $z = 1.03$, respectively, all on a benchmark mass bin of $M_{200c} \in [5 \times 10^{13}, 1 \times 10^{14}) M_{\odot} h^{-1}$. The sample is split into low and high $\Gamma_{1dyn}$ based on their median values. The scatter plot illustrates the data points, while the shaded regions show the best-fitting linear fit with $1\sigma$ confidence interval for the main sample and each subsample. At $z = 0$, the richness-concentration relation exhibits a negative slope, consistent with our expectations of halo formation physics. The slopes for the low and high $\Gamma_{1dyn}$ subsamples diverge due to the negative correlation between concentration and MAR. However, at $z = 0.49$ and $z = 1.03$, the slopes for the entire sample and/or the high $\Gamma_{1dyn}$ subsample become positive, contrary to our observations of the richness-concentration relation. In contrast, the low $\Gamma_{1dyn}$ subsample still shows a negative slope. These findings suggest that at medium to high redshifts, a subset of unrelaxed and recently merged haloes with high MAR could elevate the concentration from its expected value at hydrostatic equilibrium.</figcaption></figure>

for potential merger events even though this parameter is excluded in the final linear regression model due to multicollinearity (see Appendix E).

Fig. 10 displays the halo concentration plotted against the richness residuals, separated by $\Gamma_{1dyn}$, at benchmark bins of $M_{200c} \in [5 \times 10^{13}, 1 \times 10^{14}) M_{\odot} h^{-1}$ at three different redshift snapshots of $z = 0, 0.49, 1.03$. At $z = 0$, we observe a negative slope as expected from halo formation physics for both low-$\Gamma_{1dyn}$ and high-$\Gamma_{1dyn}$ subsamples, as well as for the overall sample. Furthermore, we observe a change in the slope between the low-$\Gamma_{1dyn}$ and high-$\Gamma_{1dyn}$ subsamples, which can be explained by the gradual increase (or decrease) in concentration ($\Gamma_{1dyn}$) over time, even without major merger events (Wechsler et al. 2002; Zhao et al. 2003; Lu et al. 2006). At redshifts of $z = 0.49, 1.03$, we observe a positive slope in the overall and/or high-$\Gamma_{1dyn}$ samples, which contradicts the scaling relations between HOD and concentration in models that track their gradual evolution over $T \gg T_{1dyn}$. However, in the presence of major mergers, when $\Gamma_{1dyn}$ is significantly enhanced, the halo concentration may also experience a transient spike after the merger. The deviation from hydrostatic equilibrium provides a plausible explanation for a positive $\beta\{N_{gal}, c_{vir}\}$, which could be fully tested on MDPL2 through the reconstruction of halo merger trees, an analysis beyond the scope of this paper.

## 6.2 Secondary halo parameter dependence of $\Delta\Sigma$

In this section, we employ a multilinear regression approach to model the lensing signal, similar to the methodology described in Section 6.1. We extend this approach to the model $P(\Delta\Sigma|N_{gal}(\Pi), M, z)$ as a linear function of $\Pi$. Upon analysing different $(M, z, r_p)$ bins, we observe that the reduced parameters $\Pi \subset \{a_{1/2}, c_{vir}, T/U, \Gamma_{2dyn}, X_{off}\}$ pass the variance inflation factor (VIF) test for multicollinearity, or in other words we showed that the variance is not inflated and thereby made less reliable in the case that the secondary halo parameters in the full model are highly correlated. As with the

case for the lensing signal, this indicates that $\ln N_{gal}$ can be described without redundancy by a linear decomposition of these reduced parameters. Furthermore, most bins exhibit homoscedasticity, as confirmed by passing the Breusch-Pagan Lagrange multiplier test. This implies that the scatter terms $\sigma_{\Delta\Sigma}$ and $\sigma_{\Pi_i}$ remain constant within each bin, with a few exceptions. Lastly, the scatter $\sigma_{\Delta\Sigma|N_{gal}}$ in most bins (with a few exceptions) meets the criteria of the Shapiro-Wilk test for Gaussianity, suggesting that the distribution closely resembles a Gaussian distribution.

The multilinear regression is a good fit to the conditioned lensing signal if we assume that $P(\ln N_{gal}|M, z)$ and $P(\Delta\Sigma|M, z)$ can be modelled with a normal distribution (e.g. Anbajagane et al. 2020; Costanzi et al. 2021; To et al. 2021). In this case $P(\Delta\Sigma|N_{gal}, M, z)$ is a multilinear equation with respect to the secondary halo parameters with mean

$$ \langle \Delta\Sigma | N_{gal}, M, z \rangle = \langle \Delta\Sigma | N_{gal}, M, z \rangle + C_1 \sigma_{\Delta\Sigma} \left( \sum_i \frac{\beta_{N_{gal},i}}{\sigma_{\Pi_i}} \rho_{\Delta\Sigma-\Pi_i} \times (\Pi_i - (\Pi_i | M, z)) \right) \quad (29) $$

and is normally distributed around the mean with variance

$$ \sigma_{\Delta\Sigma|\ln N_{gal}}^2 = \sigma_0^2 + C_2 \sum_i \beta_{N_{gal},i}^2 \sigma_{\Pi_i}^2 (1 - \rho_{\Delta\Sigma-\Pi_i}^2) + C_3 \sum_{i,j} \rho_{\Pi_i-\Pi_j} \sigma_{\Pi_i} \sigma_{\Pi_j} \quad (30) $$

The parameters $C_1, C_2, C_3$, and $\sigma_0$ can be explicitly derived where $P(\ln N_{gal}|M, z)$ and $P(\Delta\Sigma|M, z)$ are known, but the exact values are not essential for this paper. We refer the reader to Appendix F for derivations of equations (29) and (30).

We note that only in bins of $R \lesssim R_{200c}$ do the multilinear regression models pass the global F-statistic test and the T-statistic test for each parameter. This result suggests that, at $R \gtrsim R_{200c}$, we find little correlation between $\Delta\Sigma$ and $\Pi_i$. Because the scatter still passes
<figure><img src="image_14.png" /><figcaption>**Figure 11.** The dependence of $\Delta\Sigma$ on accretion history parameters in $M_{200c} \in [2 \times 10^{14}, 5 \times 10^{14}) M_{\odot} h^{-2}$, $z = 0$. In each of the five panels is plotted the $\Delta\Sigma$ in different quintiles of the accretion history parameter $\Pi \in \{a_{1/2}, c_{vir}, T/|U|, \Gamma_{2dyn}, X_{off}\}$ compared to the mean $\Delta\Sigma$. For $\{a_{1/2}, T/|U|, \Gamma_{2dyn}, X_{off}\}$ there is a strong negative correlation at small scales at $R \lesssim R_{200c}$ and for $c_{vir}$ we find a strong positive correlation at $R \lesssim R_{200c}$. Comparing the width of $\Delta\Sigma$ binned at different quintiles of $\Pi$ with the standard deviation of the profiles, we find that accretion history parameters account for much of the variance at small scales and play a negligible role at large scales.</figcaption></figure>

the Shapiro-Walk test for Gaussanity, the conditional probability $P(\Delta\Sigma|N_{gal}, M, z)$ at large scales is still normally distributed, but with $\rho_{\Delta\Sigma-\Pi_i} = 0$. By setting $\rho_{\Delta\Sigma-\Pi_i} = 0$ the variance can be reduced to

$$ \sigma_{\Delta\Sigma|N_{gal}}^2 = \sigma_0^2 + \sum_i \beta_{N_{gal},i}^2 \sigma_{\Pi_i}^2 + \sum_{i,j} \rho_{\Pi_i-\Pi_j} \sigma_{\Pi_i} \sigma_{\Pi_j}. \quad (31) $$

We visualize the dependence of $\rho_{\Delta\Sigma-\Pi_i}$ on $R/R_{200c}$ in Fig. 11. By dividing $\Delta\Sigma$ into quintiles of $\Pi_i$ we find a strong correlation for all parameters at $R \lesssim R_{200c}$ and a null correlation at $R \gtrsim R_{200c}$. On small scales, our results show a positive correlation for concentration and a negative correlation for $\{a_{1/2}, T/U, \Gamma_{2dyn}, X_{off}\}$. We observe that this trend holds for all $(M, z)$ bins plotted for a benchmark bin of $M_{200c} \in [2 \times 10^{14}, 5 \times 10^{14}) M_{\odot} h^{-1}$ at $z = 0$.

The dependence of $\Delta\Sigma$ on secondary halo parameters qualitatively agrees with Xhakaj et al. (2022) wherein they targeted a narrow mass bin, with residual mass dependency inside the bin resampled so that mass follows the same distribution. In our work, we remove the mass dependency with the KLLR method (Farahi et al. 2022a), which achieves the same effect. We extend their results to mass and redshift bins probed by the optical surveys and quantitatively show that the dependence of $\Delta\Sigma$ on $\Pi$ can be modelled as a multilinear equation.

## 6.3 Results: secondary halo parameter dependence of Cov($\Delta\Sigma$, ln $N_{gal}|M, z$)

In Fig. 12, we observe that the total covariance $Cov(\Delta\Sigma, \ln N_{gal}|M, z)$, which remains after removing the contribution of each secondary halo parameter $\beta_{N_{gal},i}Cov(\Delta\Sigma, \Pi_i|M, z)$, is consistent with zero at a significance threshold of 0.05 in all bins. The errors on the total covariance and individual contributions are computed using bootstrapping, and the errors on the remaining term are determined by adding the errors of the total and individual terms in quadrature.

Based on our hypothesis in equation (24), we conclude that the set of secondary halo parameters $\Pi$, which are related to the formation time and the mass accretion history of the haloes, can fully explain the joint distribution of $\Delta\Sigma$ and $\ln N_{gal}$ given the precision allowed by current errors, limited by the resolution limit (see Appendix B for information on particle resolution and measurement errors).

Since the joint distribution of $\Delta\Sigma$ and $\ln N_{gal}$ follows a multivariate normal distribution, $P(\Delta\Sigma, \ln N_{gal}|M, z)$ is completely characterized by its mean relation and $Cov(\Delta\Sigma, \ln N_{gal}|M, z)$. It should be noted that the contribution of each individual parameter to the total covariance, $\beta_{N_{gal},i}Cov(\Delta\Sigma, \Pi_i|M, z)$, is determined by the richness dependency captured by the slope $\beta_{N_{gal},i}$ and the $\Delta\Sigma$ dependence represented by $Cov(\Delta\Sigma, \Pi_i|M, z)$. Qualitatively, individual contributions to total covariance maintain their sign when both $\Delta\Sigma$ and $N_{gal}$ contributions preserve their sign. Consistent with the arguments of halo formation, $\Pi$ correlates with $\Delta\Sigma$ at small scales, as demonstrated in Fig. 11 across all $(M, z)$ bins. In most cases, the dependence of the richness on secondary halo parameters also maintains its sign across the $(M, z)$ bins. In instances where we encounter a sign reversal in the $\ln N_{gal} - c_{vir}$ relation, we speculate that it is due to a transient increase in concentration following a major merger.

Furthermore, the total and individual contributions to the covariance tend to decrease in magnitude at smaller scales with increasing redshift. This decrease in covariance can be attributed to two factors: the decreasing explanatory power of $\Pi$ on richness, as indicated by the decreasing values of $R^2$ and $F_{partial}$ values in Table E1, and the decreasing absolute value of $Cov(\Delta\Sigma, \Pi_i|M, z)$. This trend aligns with the idea that as haloes have more time to form, the secondary halo properties related to the mass accretion history become more significant both in richness and in $\Delta\Sigma$. As discussed in Section 4, the dependence of the mass and redshift on covariance can be explained by the halo peak height, $\nu(M, z)$.
<figure><img src="image_12.png" /><figcaption>**Figure 11.** The dependence of $\Delta\Sigma$ on accretion history parameters in $M_{200c} \in [2 \times 10^{14}, 5 \times 10^{14}) M_{\odot} h^{-2}$, $z = 0$. In each of the five panels is plotted the $\Delta\Sigma$ in different quintiles of the accretion history parameter $\Pi \in \{a_{1/2}, c_{vir}, T/|U|, \Gamma_{2dyn}, X_{off}\}$ compared to the mean $\Delta\Sigma$. For $\{a_{1/2}, T/|U|, \Gamma_{2dyn}, X_{off}\}$ there is a strong negative correlation at small scales at $R \lesssim R_{200c}$ and for $c_{vir}$ we find a strong positive correlation at $R \lesssim R_{200c}$. Comparing the width of $\Delta\Sigma$ binned at different quintiles of $\Pi$ with the standard deviation of the profiles, we find that accretion history parameters account for much of the variance at small scales and play a negligible role at large scales.</figcaption></figure>

the Shapiro-Walk test for Gaussanity, the conditional probability $P(\Delta\Sigma|N_{gal}, M, z)$ at large scales is still normally distributed, but with $\rho_{\Delta\Sigma-\Pi_i} = 0$. By setting $\rho_{\Delta\Sigma-\Pi_i} = 0$ the variance can be reduced to

$$ \sigma_{\Delta\Sigma|N_{gal}}^2 = \sigma_0^2 + \sum_i \beta_{N_{gal},i}^2 \sigma_{\Pi_i}^2 + \sum_{i,j} \rho_{\Pi_i-\Pi_j} \sigma_{\Pi_i} \sigma_{\Pi_j}. \quad (31) $$

We visualize the dependence of $\rho_{\Delta\Sigma-\Pi_i}$ on $R/R_{200c}$ in Fig. 11. By dividing $\Delta\Sigma$ into quintiles of $\Pi_i$ we find a strong correlation for all parameters at $R \lesssim R_{200c}$ and a null correlation at $R \gtrsim R_{200c}$. On small scales, our results show a positive correlation for concentration and a negative correlation for $\{a_{1/2}, T/U, \Gamma_{2dyn}, X_{off}\}$. We observe that this trend holds for all $(M, z)$ bins plotted for a benchmark bin of $M_{200c} \in [2 \times 10^{14}, 5 \times 10^{14}) M_{\odot} h^{-1}$ at $z = 0$.

The dependence of $\Delta\Sigma$ on secondary halo parameters qualitatively agrees with Xhakaj et al. (2022) wherein they targeted a narrow mass bin, with residual mass dependency inside the bin resampled so that mass follows the same distribution. In our work, we remove the mass dependency with the KLLR method (Farahi et al. 2022a), which achieves the same effect. We extend their results to mass and redshift bins probed by the optical surveys and quantitatively show that the dependence of $\Delta\Sigma$ on $\Pi$ can be modelled as a multilinear equation.

## 6.3 Results: secondary halo parameter dependence of Cov($\Delta\Sigma$, $\ln N_{gal}|M, z$)

In Fig. 12, we observe that the total covariance $Cov(\Delta\Sigma, \ln N_{gal}|M, z)$, which remains after removing the contribution of each secondary halo parameter $\beta_{N_{gal},i}Cov(\Delta\Sigma, \Pi_i|M, z)$, is consistent with zero at a significance threshold of 0.05 in all bins. The errors on the total covariance and individual contributions are computed using bootstrapping, and the errors on the remaining term are determined by adding the errors of the total and individual terms in quadrature.

Based on our hypothesis in equation (24), we conclude that the set of secondary halo parameters $\Pi$, which are related to the formation time and the mass accretion history of the haloes, can fully explain the joint distribution of $\Delta\Sigma$ and $\ln N_{gal}$ given the precision allowed by current errors, limited by the resolution limit (see Appendix B for information on particle resolution and measurement errors).

Since the joint distribution of $\Delta\Sigma$ and $\ln N_{gal}$ follows a multivariate normal distribution, $P(\Delta\Sigma, \ln N_{gal}|M, z)$ is completely characterized by its mean relation and $Cov(\Delta\Sigma, \ln N_{gal}|M, z)$. It should be noted that the contribution of each individual parameter to the total covariance, $\beta_{N_{gal},i}Cov(\Delta\Sigma, \Pi_i|M, z)$, is determined by the richness dependency captured by the slope $\beta_{N_{gal},i}$ and the $\Delta\Sigma$ dependence represented by $Cov(\Delta\Sigma, \Pi_i|M, z)$. Qualitatively, individual contributions to total covariance maintain their sign when both $\Delta\Sigma$ and $N_{gal}$ contributions preserve their sign. Consistent with the arguments of halo formation, $\Pi$ correlates with $\Delta\Sigma$ at small scales, as demonstrated in Fig. 11 across all $(M, z)$ bins. In most cases, the dependence of the richness on secondary halo parameters also maintains its sign across the $(M, z)$ bins. In instances where we encounter a sign reversal in the $\ln N_{gal} - c_{vir}$ relation, we speculate that it is due to a transient increase in concentration following a major merger.

Furthermore, the total and individual contributions to the covariance tend to decrease in magnitude at smaller scales with increasing redshift. This decrease in covariance can be attributed to two factors: the decreasing explanatory power of $\Pi$ on richness, as indicated by the decreasing values of $R^2$ and $F_{partial}$ values in Table E1, and the decreasing absolute value of $Cov(\Delta\Sigma, \Pi_i|M, z)$. This trend aligns with the idea that as haloes have more time to form, the secondary halo properties related to the mass accretion history become more significant both in richness and in $\Delta\Sigma$. As discussed in Section 4, the dependence of the mass and redshift on covariance can be explained by the halo peak height, $\nu(M, z)$.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** The dependence of Cov(ΔΣ, ln Ngal) on secondary halo parameters Π. The solid blue line is the total covariance. The dashed lines of the lines represent the covariance contribution coming from each of the secondary halo parameters modelled and governed by equation (24), where Cov(ΔΣ, Π|M, z) comes from the dependence of ΔΣ and the slope βᵢ is the dependence of richness. The thick black dashed line is the remaining covariance after removing the contribution from each Πᵢ term; the errors are obtained by adding the total and individual errors in quadrature without considering the correlations between terms. In agreement with our hypothesis in equation (24), the remaining term is consistent with null at a *p* < 0.01 level for all bins.</figcaption></figure>

# 7 DISCUSSIONS

**Intrinsic versus Extrinsic Covariance.** Our study distinguishes between the intrinsic covariance investigated here and that observed in empirical cluster data sets. This distinction arises from systematic biases introduced by the cluster-finding algorithm (extrinsic

component) and the underlying physics governing halo formation (intrinsic component). Specifically, our analysis involves counting galaxies in 3D physical space. In contrast, a realistic cluster finder like redMaPPer (Rykoff et al. 2014) employs a probabilistic assignment of galaxies to haloes in 2D physical space, considering projected radii and redshift through colour matching on to the red sequence.
<figure><img src="image_13.png" /><figcaption>**Figure 12.** The dependence of Cov(ΔΣ, ln Ngal) on secondary halo parameters Π. The solid blue line is the total covariance. The dashed lines of the lines represent the covariance contribution coming from each of the secondary halo parameters modelled and governed by equation (24), where Cov(ΔΣ, Π|M, z) comes from the dependence of ΔΣ and the slope βᵢ is the dependence of richness. The thick black dashed line is the remaining covariance after removing the contribution from each Πᵢ term; the errors are obtained by adding the total and individual errors in quadrature without considering the correlations between terms. In agreement with our hypothesis in equation (24), the remaining term is consistent with null at a *p* < 0.01 level for all bins.</figcaption></figure>

# 7 DISCUSSIONS

**Intrinsic versus Extrinsic Covariance.** Our study distinguishes between the intrinsic covariance investigated here and that observed in empirical cluster data sets. This distinction arises from systematic biases introduced by the cluster-finding algorithm (extrinsic

component) and the underlying physics governing halo formation (intrinsic component). Specifically, our analysis involves counting galaxies in 3D physical space. In contrast, a realistic cluster finder like redMaPPer (Rykoff et al. 2014) employs a probabilistic assignment of galaxies to haloes in 2D physical space, considering projected radii and redshift through colour matching on to the red sequence.
Our study does not account for the observational systematics in redMaPPer associated with uncertainties in photometric redshifts and projection effects (Rozo et al. 2015; Farahi et al. 2016).

We find that the fractional amplitude of the bias and the scale dependence on the lensing signal observed in our results (Fig. 9) are comparable to those reported by Farahi, Nagai & Anbajagane (2022b), who measured the covariance between the dark matter density and the galaxy number count enclosed inside a halo after applying a realistic stellar mass cut. The enclosed mass within a 3D radius using the IllustrisTNG100 simulation is anchored at $z = 0.24$. By comparing our findings to those obtained using a realistic cluster finder such as redMaPPer, we can unravel intrinsic and extrinsic contributions to the covariance between weak lensing observables (Wu et al. 2022). This work provides more profound insights into the distinct effects originating from the underlying physics and the methodology employed in cluster-finding algorithms (Euclid Collaboration 2019).

**Projection effects.** A noteworthy distinction arises regarding the covariance observed in our study compared to others examining clusters in projected space. Projection effects can potentially introduce a sign flip in covariance, as they can positively bias both $\Delta\Sigma$ and richness (Costanzi et al. 2019; Wu et al. 2022; Zhang & Annis 2022; Zhang et al. 2023). Particularly, Zhang et al. (2023) showed that the lensing signal can be affected both at large and small scales from the preferential alignment of halo orientation with the underlying large-scale structure filament. Wu et al. (2022) detected a positive correlation between $\Delta\Sigma$ and $\ln N_{\text{gal}}$ employing Buzzard simulations, where $\Delta\Sigma$ were measured using dark matter particles and galaxy counts were performed within a cylindrical region of depth $60 \text{ Mpc } h^{-1}$. Their investigation revealed that the positive correlation primarily stems from galaxy number counts beyond the halo's virial radius and within $60 \text{ Mpc } h^{-1}$. On the other hand, using the Dark Quest emulator and HOD-based galaxy catalogues (Nishimichi et al. 2019), Sunayama et al. (2020) found negligible deviations from the mean relation at small scales and an overall reduction in selection bias at large scales, approximately halving the effect observed by Wu et al. (2022). Furthermore, Huang et al. (2022), using data from the Subaru HSC survey, observed that the selection bias is most prominent in the vicinity of the transition from the one-halo to the two-halo regime, as evidenced by the comparison between the outer stellar mass proxy and richness. In a study based on the IllustrisTNG300 simulation, Zhang & Annis (2022) discovered a net positive correlation between the fitted weak lensing mass and the projected 2D number count of the halo when conditioned on the halo mass.

These results suggest that projection effects can potentially introduce a positively correlated bias to both $\Delta\Sigma$ and $N_{\text{gal}}$. We can estimate the impact of projection effects by comparing the intrinsic covariance measured in our study with the total covariance observed in the projected space. Consequently, our results serve two essential purposes: (i) elucidating the physical origins of the negative covariance and (ii) discerning intrinsic and extrinsic components to determine the covariance attributable to projection effects accurately.

**Radial dependence.** There is a notable difference in the reported amplitude and scale dependence of covariance, which can be attributed to discrepancies in the employed halo occupation density models. Notably, a distinctive scale dependence discrepancy exists between simulation-based investigations of projection effects (Salcedo et al. 2020; Sunayama et al. 2020; Wu et al. 2022) and observational data from the HSC (Huang et al. 2022). In particular, the analysis of observational data reveals a prominent 1 Mpc bump,

which could be explained by uncertainties inherent in observations, such as miscentring effects. It is crucial to gain insight into the sensitivity of the covariance with respect to the model parameters and the influence of selection effects. Understanding these factors is essential to comprehensively interpret and account for the observed covariance in galaxy cluster survey studies.

**Accuracy versus precision.** The statistical power of current and future surveys enables us to determine the normalization and slope of the mass-observable relations at a few per cent levels (e.g. Farahi et al. 2016; Mantz et al. 2016b; Mulroy et al. 2019; To et al. 2021). However, these estimates are susceptible to known and unknown sources of systematic errors that inflate the uncertainties. These uncertainties introduce biases and degrade the accuracy of the results. Therefore, it is essential to carefully identify, quantify, and account for these systematic effects to ensure robust and reliable measurements. In this work, we focus on studying one of these sources of systematic uncertainty that was not considered previously.

## 8 SUMMARY

This work reveals insights into the scale-dependent covariance between weak lensing observables and the physical properties of the halo. Using the MDPL2 N-body simulation with galaxies painted using the SAGE SAM, we present several key findings:

(i) We observe that the intrinsic covariance between $\Delta\Sigma$ and $\ln N_{\text{gal}}$ enclosed within a 3D radius is negative at small scales and null at large scales in $(\ln M, z)$ ranges that cover optical surveys.

(ii) We model the shape of the covariance across all bins using an error function that is insensitive to the radius definition used to define halo boundaries.

(iii) We find that the magnitude of the covariance is relatively insensitive to mass and decreases considerably with increasing redshift. The $(M, z)$ dependence of the shape of the covariance can be encapsulated by the peak height parameter $v(M, z)$, which suggests that the scale of the covariance is related to the formation history of haloes.

(iv) We show that incorporating the covariance into $\langle\Delta\Sigma|N_{\text{gal}}, z, r_p\rangle$ using the first-order expansion of the halo mass function yields about > 1 per cent bias on $\langle\Delta\Sigma|N_{\text{gal}}, z, r_p\rangle$ at small scales, which implies a mass bias of > 2 per cent in the halo mass estimates in most bins.

(v) Our analysis reveals that the covariance between $\ln N_{\text{gal}}$ and $\Delta\Sigma$ can be fully explained by secondary halo parameters related to the history of the halo assembly. This finding provides strong evidence that the non-zero covariance results from the variation in the formation history of dark matter haloes.

This work underscores the importance of accounting for covariance in cluster mass calibration. Incorporating the covariance between richness and the weak lensing signal and its characterization should be an essential component of weak lensing cluster mass calibration in upcoming optical cluster surveys. The results of this work can be introduced as a simulation-based prior for the forward modelling of scaling relations used by cluster cosmological analysis pipelines. Within the LSST-DESC framework, this systematic bias would be implemented in the CLMM pipeline (Aguena et al. 2021) to update the stacked weak lensing mass in each richness bin. Considering this covariance paves the path toward per cent-level accuracy cosmological constraints, thereby enhancing the precision and reliability of our scientific conclusions. Moving forward, it is
Our study does not account for the observational systematics in redMaPPer associated with uncertainties in photometric redshifts and projection effects (Rozo et al. 2015; Farahi et al. 2016).

We find that the fractional amplitude of the bias and the scale dependence on the lensing signal observed in our results (Fig. 9) are comparable to those reported by Farahi, Nagai & Anbajagane (2022b), who measured the covariance between the dark matter density and the galaxy number count enclosed inside a halo after applying a realistic stellar mass cut. The enclosed mass within a 3D radius using the IllustrisTNG100 simulation is anchored at $z = 0.24$. By comparing our findings to those obtained using a realistic cluster finder such as redMaPPer, we can unravel intrinsic and extrinsic contributions to the covariance between weak lensing observables (Wu et al. 2022). This work provides more profound insights into the distinct effects originating from the underlying physics and the methodology employed in cluster-finding algorithms (Euclid Collaboration 2019).

**Projection effects.** A noteworthy distinction arises regarding the covariance observed in our study compared to others examining clusters in projected space. Projection effects can potentially introduce a sign flip in covariance, as they can positively bias both $\Delta\Sigma$ and richness (Costanzi et al. 2019; Wu et al. 2022; Zhang & Annis 2022; Zhang et al. 2023). Particularly, Zhang et al. (2023) showed that the lensing signal can be affected both at large and small scales from the preferential alignment of halo orientation with the underlying large-scale structure filament. Wu et al. (2022) detected a positive correlation between $\Delta\Sigma$ and $\ln N_{\text{gal}}$ employing Buzzard simulations, where $\Delta\Sigma$ were measured using dark matter particles and galaxy counts were performed within a cylindrical region of depth $60 \text{ Mpc } h^{-1}$. Their investigation revealed that the positive correlation primarily stems from galaxy number counts beyond the halo's virial radius and within $60 \text{ Mpc } h^{-1}$. On the other hand, using the Dark Quest emulator and HOD-based galaxy catalogues (Nishimichi et al. 2019), Sunayama et al. (2020) found negligible deviations from the mean relation at small scales and an overall reduction in selection bias at large scales, approximately halving the effect observed by Wu et al. (2022). Furthermore, Huang et al. (2022), using data from the Subaru HSC survey, observed that the selection bias is most prominent in the vicinity of the transition from the one-halo to the two-halo regime, as evidenced by the comparison between the outer stellar mass proxy and richness. In a study based on the IllustrisTNG300 simulation, Zhang & Annis (2022) discovered a net positive correlation between the fitted weak lensing mass and the projected 2D number count of the halo when conditioned on the halo mass.

These results suggest that projection effects can potentially introduce a positively correlated bias to both $\Delta\Sigma$ and $N_{\text{gal}}$. We can estimate the impact of projection effects by comparing the intrinsic covariance measured in our study with the total covariance observed in the projected space. Consequently, our results serve two essential purposes: (i) elucidating the physical origins of the negative covariance and (ii) discerning intrinsic and extrinsic components to determine the covariance attributable to projection effects accurately.

**Radial dependence.** There is a notable difference in the reported amplitude and scale dependence of covariance, which can be attributed to discrepancies in the employed halo occupation density models. Notably, a distinctive scale dependence discrepancy exists between simulation-based investigations of projection effects (Salcedo et al. 2020; Sunayama et al. 2020; Wu et al. 2022) and observational data from the HSC (Huang et al. 2022). In particular, the analysis of observational data reveals a prominent 1 Mpc bump,

which could be explained by uncertainties inherent in observations, such as miscentring effects. It is crucial to gain insight into the sensitivity of the covariance with respect to the model parameters and the influence of selection effects. Understanding these factors is essential to comprehensively interpret and account for the observed covariance in galaxy cluster survey studies.

**Accuracy versus precision.** The statistical power of current and future surveys enables us to determine the normalization and slope of the mass-observable relations at a few per cent levels (e.g. Farahi et al. 2016; Mantz et al. 2016b; Mulroy et al. 2019; To et al. 2021). However, these estimates are susceptible to known and unknown sources of systematic errors that inflate the uncertainties. These uncertainties introduce biases and degrade the accuracy of the results. Therefore, it is essential to carefully identify, quantify, and account for these systematic effects to ensure robust and reliable measurements. In this work, we focus on studying one of these sources of systematic uncertainty that was not considered previously.

## 8 SUMMARY

This work reveals insights into the scale-dependent covariance between weak lensing observables and the physical properties of the halo. Using the MDPL2 N-body simulation with galaxies painted using the SAGE SAM, we present several key findings:

(i) We observe that the intrinsic covariance between $\Delta\Sigma$ and $\ln N_{\text{gal}}$ enclosed within a 3D radius is negative at small scales and null at large scales in $(\ln M, z)$ ranges that cover optical surveys.

(ii) We model the shape of the covariance across all bins using an error function that is insensitive to the radius definition used to define halo boundaries.

(iii) We find that the magnitude of the covariance is relatively insensitive to mass and decreases considerably with increasing redshift. The $(M, z)$ dependence of the shape of the covariance can be encapsulated by the peak height parameter $v(M, z)$, which suggests that the scale of the covariance is related to the formation history of haloes.

(iv) We show that incorporating the covariance into $\langle\Delta\Sigma|N_{\text{gal}}, z, r_p\rangle$ using the first-order expansion of the halo mass function yields about > 1 per cent bias on $\langle\Delta\Sigma|N_{\text{gal}}, z, r_p\rangle$ at small scales, which implies a mass bias of > 2 per cent in the halo mass estimates in most bins.

(v) Our analysis reveals that the covariance between $\ln N_{\text{gal}}$ and $\Delta\Sigma$ can be fully explained by secondary halo parameters related to the history of the halo assembly. This finding provides strong evidence that the non-zero covariance results from the variation in the formation history of dark matter haloes.

This work underscores the importance of accounting for covariance in cluster mass calibration. Incorporating the covariance between richness and the weak lensing signal and its characterization should be an essential component of weak lensing cluster mass calibration in upcoming optical cluster surveys. The results of this work can be introduced as a simulation-based prior for the forward modelling of scaling relations used by cluster cosmological analysis pipelines. Within the LSST-DESC framework, this systematic bias would be implemented in the CLMM pipeline (Aguena et al. 2021) to update the stacked weak lensing mass in each richness bin. Considering this covariance paves the path toward per cent-level accuracy cosmological constraints, thereby enhancing the precision and reliability of our scientific conclusions. Moving forward, it is
imperative to integrate this understanding into the design and analysis of future optical cluster surveys.

## ACKNOWLEDGEMENTS

Author contributions are as follows:

Z. Zhang—conceptualisation, code development, methodology, writing, edits

A. Farahi—conceptualisation, code development, methodology, writing, edits, supervision

D. Nagai—methodology, writing, edits, supervision

E. Lau—methodology, code development, writing, edits, supervision

J. Frieman—edits, supervision

M. Ricci—edits, supervision, administrative organisation

A. Linden—internal reviewer

H. Wu—internal reviewer

We thank the anonymous referee for careful edits that greatly improved the quality of this manuscript. This paper has undergone an internal review in the LSST Dark Energy Science Collaboration. We thank Anja von der Linden, Tamas Varga and Hao-Yi Wu for serving as the LSST-DESC publication review committee. In addition, the authors thank Andrew Hearin, Heather Kelly, Johnny Esteves, Enia Xhakaj, and Conghao Zhou for their helpful discussions.

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, the National Science Foundation, and the LSST Corporation in the United States. DESC uses 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.

The CosmoSim data base used in this paper is a service by the Leibniz-Institute for Astrophysics Potsdam (AIP). The MultiDark data base was developed with the Spanish MultiDark Consolidier Project CSD2009-00064.

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) and the Partnership for Advanced Supercomputing in Europe (PRACE, www.prace-ri.eu) for funding the MultiDark simulation project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, www.lrz.de). The Bolshoi simulations have been performed within the Bolshoi project of the University of California High-Performance AstroComputing Center (UC-HiPACC) and were run at the NASA Ames Research Center.

This research used resources of 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 using NERSC award HEP-ERCAP0022779.

DN acknowledges support by NSF (AST-2206055 & 2307280) and NASA (80NSSC22K0821 and TM3-24007X) grants. AF acknowledges the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have partially contributed to the research results reported within this paper.

## DATA AVAILABILITY

Raw data from the MDPL2 simulation are publicly available at CosmoSim.<fnref n="2" /> Derived data products are available upon request. Our code is made publicly available at https://github.com/BaryonPasters/BP_CorrelatedScatter_MDPL2.

## REFERENCES

Abbott T. et al., 2018, Phys. Rev. D, 98, 043526

Abbott T. et al., 2020, Phys. Rev. D, 102, 023509

Aguena M. et al., 2021, MNRAS, 508, 6092

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

Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS, 367, 1781

Anbajagane D., Evrard A. E., Farahi A., Barnes D. J., Dolag K., McCarthy I. G., Nelson D., Pillepich A., 2020, MNRAS, 495, 686

Anbajagane D., Evrard A. E., Farahi A., 2022, MNRAS, 509, 3441

Applegate D. E. et al., 2014, MNRAS, 439, 48

Baldry I. K. et al., 2012, MNRAS, 421, 621

Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, ApJ, 762, 109

Bhattacharya S., Heitmann K., White M., Lukić Z., Wagner C., Habib S., 2011, ApJ, 732, 122

Bocquet S., Carter F. W., 2016, J. Open Source Softw., 1, 46

Bocquet S., Saro A., Dolag K., Mohr J. J., 2016, MNRAS, 456, 2361

Breusch T., Pagan A., 1979, Econometrica, 47, 1287

Bryan G. L., Norman M. L., 1998, ApJ, 495, 80

Chiu I. N., Umetsu K., Murata R., Medezinski E., Oguri M., 2020, MNRAS, 495, 428

Costanzi M. et al., 2019, MNRAS, 482, 490

Costanzi M. et al., 2021, Phys. Rev. D, 103, 043522

Croton D. J. et al., 2016, ApJS, 222, 22

Dark Energy Survey Collaboration et al., 2016, MNRAS, 460, 1270

Despali G., Giocoli C., Angulo R. E., Tormen G., Sheth R. K., Baso G., Moscardini L., 2016, MNRAS, 456, 2486

de Haan T. et al., 2016, ApJ, 832, 95

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

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

Dodelson S., Heitmann K., Hirata C., Honscheid K., Roodman A., Seljak U., Slosar A., Trodden M., 2016, preprint (arXiv:1604.07626)

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. et al., 2019, Nature Commun., 10, 2504

Farahi A., Ho M., Trac H., 2020, MNRAS, 493, 1361

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

Farahi A., Nagai D., Anbajagane D., 2022b, ApJ, 933, 48

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

Gelman A., Rubin D. B., 1992, Stat. Sci., 7, 457

Giodini S., Lovisari L., Pointecouteau E., Ettori S., Reiprich T. H., Hoekstra H., 2013, Space Sci. Rev., 177, 247

Golden-Marx J. B., Miller C. J., 2018, ApJ, 860, 2

Goodman J., Weare J., 2010, Commun. Appl. Math. Comput. Sci., 5, 65

Hearin A. P., Zentner A. R., Newman J. A., Berlind A. A., 2013, MNRAS, 430, 1238

Huang S. et al., 2022, MNRAS, 515, 4722

Ivezić Ž. et al., 2019, ApJ, 873, 111

Johnston D. E. et al., 2007, preprint (arXiv:0709.1159)

Kettula K. et al., 2015, MNRAS, 451, 1460

Kiiveri K. et al., 2021, MNRAS, 502, 1494

<footnote n="2">²https://www.cosmosim.org/</footnote>
imperative to integrate this understanding into the design and analysis of future optical cluster surveys.

## ACKNOWLEDGEMENTS

Author contributions are as follows:

Z. Zhang—conceptualisation, code development, methodology, writing, edits

A. Farahi—conceptualisation, code development, methodology, writing, edits, supervision

D. Nagai—methodology, writing, edits, supervision

E. Lau—methodology, code development, writing, edits, supervision

J. Frieman—edits, supervision

M. Ricci—edits, supervision, administrative organisation

A. Linden—internal reviewer

H. Wu—internal reviewer

We thank the anonymous referee for careful edits that greatly improved the quality of this manuscript. This paper has undergone an internal review in the LSST Dark Energy Science Collaboration. We thank Anja von der Linden, Tamas Varga and Hao-Yi Wu for serving as the LSST-DESC publication review committee. In addition, the authors thank Andrew Hearin, Heather Kelly, Johnny Esteves, Enia Xhakaj, and Conghao Zhou for their helpful discussions.

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, the National Science Foundation, and the LSST Corporation in the United States. DESC uses 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.

The CosmoSim data base used in this paper is a service by the Leibniz-Institute for Astrophysics Potsdam (AIP). The MultiDark data base was developed with the Spanish MultiDark Consolidator Project CSD2009-00064.

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) and the Partnership for Advanced Supercomputing in Europe (PRACE, www.prace-ri.eu) for funding the MultiDark simulation project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, www.lrz.de). The Bolshoi simulations have been performed within the Bolshoi project of the University of California High-Performance AstroComputing Center (UC-HiPACC) and were run at the NASA Ames Research Center.

This research used resources of 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 using NERSC award HEP-ERCAP0022779.

DN acknowledges support by NSF (AST-2206055 & 2307280) and NASA (80NSSC22K0821 and TM3-24007X) grants. AF acknowledges the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have partially contributed to the research results reported within this paper.

## DATA AVAILABILITY

Raw data from the MDPL2 simulation are publicly available at CosmoSim.<fnref n="2" /> Derived data products are available upon request. Our code is made publicly available at https://github.com/BaryonPasters/BP_CorrelatedScatter_MDPL2.

## REFERENCES

Abbott T. et al., 2018, Phys. Rev. D, 98, 043526

Abbott T. et al., 2020, Phys. Rev. D, 102, 023509

Aguena M. et al., 2021, MNRAS, 508, 6092

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

Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS, 367, 1781

Anbajagane D., Evrard A. E., Farahi A., Barnes D. J., Dolag K., McCarthy I. G., Nelson D., Pillepich A., 2020, MNRAS, 495, 686

Anbajagane D., Evrard A. E., Farahi A., 2022, MNRAS, 509, 3441

Applegate D. E. et al., 2014, MNRAS, 439, 48

Baldry I. K. et al., 2012, MNRAS, 421, 621

Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, ApJ, 762, 109

Bhattacharya S., Heitmann K., White M., Lukić Z., Wagner C., Habib S., 2011, ApJ, 732, 122

Bocquet S., Carter F. W., 2016, J. Open Source Softw., 1, 46

Bocquet S., Saro A., Dolag K., Mohr J. J., 2016, MNRAS, 456, 2361

Breusch T., Pagan A., 1979, Econometrica, 47, 1287

Bryan G. L., Norman M. L., 1998, ApJ, 495, 80

Chiu I. N., Umetsu K., Murata R., Medezinski E., Oguri M., 2020, MNRAS, 495, 428

Costanzi M. et al., 2019, MNRAS, 482, 490

Costanzi M. et al., 2021, Phys. Rev. D, 103, 043522

Croton D. J. et al., 2016, ApJS, 222, 22

Dark Energy Survey Collaboration et al., 2016, MNRAS, 460, 1270

Despali G., Giocoli C., Angulo R. E., Tormen G., Sheth R. K., Baso G., Moscardini L., 2016, MNRAS, 456, 2486

de Haan T. et al., 2016, ApJ, 832, 95

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

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

Dodelson S., Heitmann K., Hirata C., Honscheid K., Roodman A., Seljak U., Slosar A., Trodden M., 2016, preprint (arXiv:1604.07626)

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. et al., 2019, Nature Commun., 10, 2504

Farahi A., Ho M., Trac H., 2020, MNRAS, 493, 1361

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

Farahi A., Nagai D., Anbajagane D., 2022b, ApJ, 933, 48

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

Gelman A., Rubin D. B., 1992, Stat. Sci., 7, 457

Giodini S., Lovisari L., Pointecouteau E., Ettori S., Reiprich T. H., Hoekstra H., 2013, Space Sci. Rev., 177, 247

Golden-Marx J. B., Miller C. J., 2018, ApJ, 860, 2

Goodman J., Weare J., 2010, Commun. Appl. Math. Comput. Sci., 5, 65

Hearin A. P., Zentner A. R., Newman J. A., Berlind A. A., 2013, MNRAS, 430, 1238

Huang S. et al., 2022, MNRAS, 515, 4722

Ivezić Ž. et al., 2019, ApJ, 873, 111

Johnston D. E. et al., 2007, preprint (arXiv:0709.1159)

Kettula K. et al., 2015, MNRAS, 451, 1460

Kiiveri K. et al., 2021, MNRAS, 502, 1494

<footnote n="2">²https://www.cosmosim.org/</footnote>
Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340

Lau E. T., Nagai D., Avestruz C., Nelson K., Vikhlinin A., 2015, ApJ, 806, 68

Lee W. et al., 2023, ApJ, 945, 71

Lesci G. F. et al., 2022, A&A, 665, A100

Lokken M. et al., 2022, ApJ, 933, 134

Lu Y., Mo H. J., Katz N., Weinberg M. D., 2006, MNRAS, 368, 1931

Ludlow A. D., Navarro J. F., Li M., Angulo R. E., Boylan-Kolchin M., Bett P. E., 2012, MNRAS, 427, 1322

Ludlow A. D., Bose S., Angulo R. E., Wang L., Hellwing W. A., Navarro J. F., Cole S., Frenk C. S., 2016, MNRAS, 460, 1214

Mantz A., Allen S. W., Ebeling H., Rapetti D., Drlica-Wagner A., 2010, MNRAS, 406, 1773

Mantz A. B. et al., 2015, MNRAS, 446, 2205

Mantz A. B., Allen S. W., Morris R. G., Schmidt R. W., 2016a, MNRAS, 456, 4020

Mantz A. B. et al., 2016b, MNRAS, 463, 3582

Maraston C. et al., 2013, MNRAS, 435, 2764

McClintock T. et al., 2019, MNRAS, 482, 1352

Miyatake H. et al., 2022, Phys. Rev. D, 106, 083519

Mulroy S. L. et al., 2019, MNRAS, 484, 60

Murata R. et al., 2019, PASJ, 71, 107

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

Nishimichi T. et al., 2019, ApJ, 884, 29

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

Percival W. J., 2005, A&A, 443, 819

Pierre M. et al., 2016, A&A, 592, A1

Pillepich A., Porciani C., Hahn O., 2010, MNRAS, 402, 191

Planck CollaborationXVI, 2014, A&A, 571, A16

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

Pratt G. W., Arnaud M., Biviano A., Eckert D., Ettori S., Nagai D., Okabe N., Reiprich T. H., 2019, Space Sci. Rev., 215, 25

Press W. H., Schechter P., 1974, ApJ, 187, 425

Riebe K. et al., 2013, Astron. Nachr., 334, 691

Rozo E. et al., 2010, ApJ, 708, 645

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

Rozo E., Rykoff E. S., Becker M., Reddick R. M., Wechsler R. H., 2015, MNRAS, 453, 38

Rykoff E. S. et al., 2014, ApJ, 785, 104

Salcedo A. N., Wibking B. D., Weinberg D. H., Wu H.-Y., Ferrer D., Eisenstein D., Pinto P., 2020, MNRAS, 491, 3061

Schrabback T. et al., 2018, MNRAS, 474, 2635

Sereno M. et al., 2020, MNRAS, 492, 4528

Shapiro S. S., Wilk M. B., 1965, Biometrika, 52, 591

Shin T., Diemer B., 2023, MNRAS, 521, 5570

Somerville R. S., Primack J. R., Faber S. M., 2001, MNRAS, 320, 504

Stanek R., Rasia E., Evrard A. E., Pearce F., Gazzola L., 2010, ApJ, 715, 1508

Stark D. V., McGaugh S. S., Swaters R. A., 2009, AJ, 138, 392

Sunayama T., 2023, MNRAS, 521, 5064

Sunayama T. et al., 2020, MNRAS, 496, 4468

Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709

Tinker J. L., Robertson B. E., Kravtsov A. V., Klypin A., Warren M. S., Yepes G., Gottlöber S., 2010, ApJ, 724, 878

To C. et al., 2021, Phys. Rev. Lett., 126, 141301

Tremonti C. A. et al., 2004, ApJ, 613, 898

Vikhlinin A. et al., 2009, ApJ, 692, 1060

von der Linden A. et al., 2014, MNRAS, 439, 2

Wang K., Mao Y.-Y., Zentner A. R., Guo H., Lange J. U., van den Bosch F. C., Mezini L., 2022, MNRAS, 516, 4003

Watson W. A., Iliev I. T., D’Aloisio A., Knebe A., Shapiro P. R., Yepes G., 2013, MNRAS, 433, 1230

Wechsler R. H., Tinker J. L., 2018, ARA&A, 56, 435

Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., Dekel A., 2002, ApJ, 568, 52

Wu H.-Y., Evrard A. E., Hahn O., Martizzi D., Teyssier R., Wechsler R. H., 2015, MNRAS, 452, 1982

Wu H.-Y. et al., 2022, MNRAS, 515, 4471

Xhakaj E., Leauthaud A., Lange J., Hearin A., Diemer B., Dalal N., 2022, MNRAS, 514, 2876

Zentner A. R., Berlind A. A., Bullock J. S., Kravtsov A. V., Wechsler R. H., 2005, ApJ, 624, 505

Zentner A. R., Hearin A., van den Bosch F. C., Lange J. U., Villarreal A. S., 2019, MNRAS, 485, 1196

Zhang Y., Annis J., 2022, MNRAS, 511, L30

Zhang Z. et al., 2023, MNRAS, 523, 1994

Zhao D. H., Jing Y. P., Mo H. J., Brner G., 2003, ApJ, 597, L9

Zu Y., Weinberg D. H., Rozo E., Sheldon E. S., Tinker J. L., Becker M. R., 2014, MNRAS, 439, 1628

## APPENDIX A: FUNCTIONAL FORM

This section aims to characterize the shape of the covariance across mass and redshift bins by fitting a template curve. The process involves several transformations and adjustments. First, a logarithmic transformation is applied to the radial bins, denoted as $x = \log_{10}(R/R_{200c})$. Then, a horizontal offset is introduced using a parameter $\gamma$, and scaling is applied using a parameter $\tau$. This results in a transformed variable $\tilde{x} = (x - \gamma)/\tau$.

To analyse the transformed data vector $f(\tilde{x})$, we test a set of functional forms presented in Table A2. The normalization factors and coefficients associated with these functions are chosen in such a way that $f(\tilde{x})$ approaches 1 at large scales, $-1$ at small scales, $f(0) = 1$ and $f(0) = 0$.

A linear transformation of $f(\tilde{x})$ is then performed, given by $s(f(\tilde{x}) + g)$, where $g$ represents a vertical shift and $s$ represents a scaling factor. The magnitude of $s$ is comparable to the magnitude of Cov($\Delta\Sigma$, $\ln N_{gal}|M, z$), while the parameters $\gamma$, $\tau$, $g$ are of the order of unity. These parameters, along with $s$, form the set of parameters denoted as $\theta \in \{\tau, \gamma, g, s\}$, which define our best-fitting model. If $g = -1$, it implies a zero covariance at large scales. We fit two models: a full model with all parameters free, and a reduced model with $g = -1$, and the rest of the parameters free. The priors for the parameters are specified in Table A1.

For the full model, we choose the error function as our fiducial functional form. The estimated parameters for both the full and reduced models are presented in Tables D1 and D2, respectively. Appendix D provides robustness testing to determine the best-fitting model for our covariance. In Table D1, we compare the model parameters, $\chi^2$ p-value, and the difference in Deviance Information Criterion (DIC) with our fiducial model using the candidate functions listed in Table A2. The error function generally outperforms other models when all parameters $\theta \in \{\tau, \gamma, g, s\}$ are allowed to vary. Table D2 shows that the DIC of the reduced error function ($g = -1$)

<table><caption>**Table A1.** Priors for the model. We introduce two sets of priors. In the 'full' models, the parameters are given physical (i.e. $\tau > 0$, $s > 0$) but non-informative uniform or log-uniform priors. In the 'reduced' case, assuming that Cov($\Delta\Sigma$, $\ln N_{gal}|M, z$) = 0 at large scales, we restrict $g = -1$ while assigning the same set of priors to all other parameters.</caption><thead><tr><th rowspan="2">Parameters</th><th colspan="2">Priors</th></tr><tr><th>Full</th><th>Reduced</th></tr></thead><tbody><tr><td>τ</td><td>Uniform (0, 10)</td><td>Uniform (0, 10)</td></tr><tr><td>γ</td><td>Uniform (-5, 5)</td><td>Uniform (-5, 5)</td></tr><tr><td>g</td><td>Uniform (-2, 1)</td><td>Fixed at g = -1</td></tr><tr><td>10<sup>-12</sup> × s</td><td>Log-uniform (0.01, 10)</td><td>Log-uniform (0.01, 10)</td></tr></tbody></table>
Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340

Lau E. T., Nagai D., Avestruz C., Nelson K., Vikhlinin A., 2015, ApJ, 806, 68

Lee W. et al., 2023, ApJ, 945, 71

Lesci G. F. et al., 2022, A&A, 665, A100

Lokken M. et al., 2022, ApJ, 933, 134

Lu Y., Mo H. J., Katz N., Weinberg M. D., 2006, MNRAS, 368, 1931

Ludlow A. D., Navarro J. F., Li M., Angulo R. E., Boylan-Kolchin M., Bett P. E., 2012, MNRAS, 427, 1322

Ludlow A. D., Bose S., Angulo R. E., Wang L., Hellwing W. A., Navarro J. F., Cole S., Frenk C. S., 2016, MNRAS, 460, 1214

Mantz A., Allen S. W., Ebeling H., Rapetti D., Drlica-Wagner A., 2010, MNRAS, 406, 1773

Mantz A. B. et al., 2015, MNRAS, 446, 2205

Mantz A. B., Allen S. W., Morris R. G., Schmidt R. W., 2016a, MNRAS, 456, 4020

Mantz A. B. et al., 2016b, MNRAS, 463, 3582

Maraston C. et al., 2013, MNRAS, 435, 2764

McClintock T. et al., 2019, MNRAS, 482, 1352

Miyatake H. et al., 2022, Phys. Rev. D, 106, 083519

Mulroy S. L. et al., 2019, MNRAS, 484, 60

Murata R. et al., 2019, PASJ, 71, 107

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

Nishimichi T. et al., 2019, ApJ, 884, 29

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

Percival W. J., 2005, A&A, 443, 819

Pierre M. et al., 2016, A&A, 592, A1

Pillepich A., Porciani C., Hahn O., 2010, MNRAS, 402, 191

Planck CollaborationXVI, 2014, A&A, 571, A16

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

Pratt G. W., Arnaud M., Biviano A., Eckert D., Ettori S., Nagai D., Okabe N., Reiprich T. H., 2019, Space Sci. Rev., 215, 25

Press W. H., Schechter P., 1974, ApJ, 187, 425

Riebe K. et al., 2013, Astron. Nachr., 334, 691

Rozo E. et al., 2010, ApJ, 708, 645

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

Rozo E., Rykoff E. S., Becker M., Reddick R. M., Wechsler R. H., 2015, MNRAS, 453, 38

Rykoff E. S. et al., 2014, ApJ, 785, 104

Salcedo A. N., Wibking B. D., Weinberg D. H., Wu H.-Y., Ferrer D., Eisenstein D., Pinto P., 2020, MNRAS, 491, 3061

Schrabback T. et al., 2018, MNRAS, 474, 2635

Sereno M. et al., 2020, MNRAS, 492, 4528

Shapiro S. S., Wilk M. B., 1965, Biometrika, 52, 591

Shin T., Diemer B., 2023, MNRAS, 521, 5570

Somerville R. S., Primack J. R., Faber S. M., 2001, MNRAS, 320, 504

Stanek R., Rasia E., Evrard A. E., Pearce F., Gazzola L., 2010, ApJ, 715, 1508

Stark D. V., McGaugh S. S., Swaters R. A., 2009, AJ, 138, 392

Sunayama T., 2023, MNRAS, 521, 5064

Sunayama T. et al., 2020, MNRAS, 496, 4468

Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709

Tinker J. L., Robertson B. E., Kravtsov A. V., Klypin A., Warren M. S., Yepes G., Gottlöber S., 2010, ApJ, 724, 878

To C. et al., 2021, Phys. Rev. Lett., 126, 141301

Tremonti C. A. et al., 2004, ApJ, 613, 898

Vikhlinin A. et al., 2009, ApJ, 692, 1060

von der Linden A. et al., 2014, MNRAS, 439, 2

Wang K., Mao Y.-Y., Zentner A. R., Guo H., Lange J. U., van den Bosch F. C., Mezini L., 2022, MNRAS, 516, 4003

Watson W. A., Iliev I. T., D’Aloisio A., Knebe A., Shapiro P. R., Yepes G., 2013, MNRAS, 433, 1230

Wechsler R. H., Tinker J. L., 2018, ARA&A, 56, 435

Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., Dekel A., 2002, ApJ, 568, 52

Wu H.-Y., Evrard A. E., Hahn O., Martizzi D., Teyssier R., Wechsler R. H., 2015, MNRAS, 452, 1982

Wu H.-Y. et al., 2022, MNRAS, 515, 4471

Xhakaj E., Leauthaud A., Lange J., Hearin A., Diemer B., Dalal N., 2022, MNRAS, 514, 2876

Zentner A. R., Berlind A. A., Bullock J. S., Kravtsov A. V., Wechsler R. H., 2005, ApJ, 624, 505

Zentner A. R., Hearin A., van den Bosch F. C., Lange J. U., Villarreal A. S., 2019, MNRAS, 485, 1196

Zhang Y., Annis J., 2022, MNRAS, 511, L30

Zhang Z. et al., 2023, MNRAS, 523, 1994

Zhao D. H., Jing Y. P., Mo H. J., Brner G., 2003, ApJ, 597, L9

Zu Y., Weinberg D. H., Rozo E., Sheldon E. S., Tinker J. L., Becker M. R., 2014, MNRAS, 439, 1628

## APPENDIX A: FUNCTIONAL FORM

This section aims to characterize the shape of the covariance across mass and redshift bins by fitting a template curve. The process involves several transformations and adjustments. First, a logarithmic transformation is applied to the radial bins, denoted as $x = \log_{10}(R/R_{200c})$. Then, a horizontal offset is introduced using a parameter $\gamma$, and scaling is applied using a parameter $\tau$. This results in a transformed variable $\tilde{x} = (x - \gamma)/\tau$.

To analyse the transformed data vector $f(\tilde{x})$, we test a set of functional forms presented in Table A2. The normalization factors and coefficients associated with these functions are chosen in such a way that $f(\tilde{x})$ approaches 1 at large scales, $-1$ at small scales, $f(0) = 1$ and $f(0) = 0$.

A linear transformation of $f(\tilde{x})$ is then performed, given by $s(f(\tilde{x}) + g)$, where $g$ represents a vertical shift and $s$ represents a scaling factor. The magnitude of $s$ is comparable to the magnitude of Cov($\Delta\Sigma$, $\ln N_{gal}|M, z$), while the parameters $\gamma$, $\tau$, $g$ are of the order of unity. These parameters, along with $s$, form the set of parameters denoted as $\theta \in \{\tau, \gamma, g, s\}$, which define our best-fitting model. If $g = -1$, it implies a zero covariance at large scales. We fit two models: a full model with all parameters free, and a reduced model with $g = -1$, and the rest of the parameters free. The priors for the parameters are specified in Table A1.

For the full model, we choose the error function as our fiducial functional form. The estimated parameters for both the full and reduced models are presented in Tables D1 and D2, respectively. Appendix D provides robustness testing to determine the best-fitting model for our covariance. In Table D1, we compare the model parameters, $\chi^2$ p-value, and the difference in Deviance Information Criterion (DIC) with our fiducial model using the candidate functions listed in Table A2. The error function generally outperforms other models when all parameters $\theta \in \{\tau, \gamma, g, s\}$ are allowed to vary. Table D2 shows that the DIC of the reduced error function ($g = -1$)

<table><caption>**Table A1.** Priors for the model. We introduce two sets of priors. In the 'full' models, the parameters are given physical (i.e. $\tau > 0$, $s > 0$) but non-informative uniform or log-uniform priors. In the 'reduced' case, assuming that Cov($\Delta\Sigma$, $\ln N_{gal}|M, z$) = 0 at large scales, we restrict $g = -1$ while assigning the same set of priors to all other parameters.</caption><thead><tr><th rowspan="2">Parameters</th><th colspan="2">Priors</th></tr><tr><th>Full</th><th>Reduced</th></tr></thead><tbody><tr><td>τ</td><td>Uniform (0, 10)</td><td>Uniform (0, 10)</td></tr><tr><td>γ</td><td>Uniform (-5, 5)</td><td>Uniform (-5, 5)</td></tr><tr><td>g</td><td>Uniform (-2, 1)</td><td>Fixed at g = -1</td></tr><tr><td>10<sup>-12</sup> × s</td><td>Log-uniform (0.01, 10)</td><td>Log-uniform (0.01, 10)</td></tr></tbody></table>