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

# Cosmological Constraints from Galaxy Cluster Sparsity, Cluster Gas Mass Fraction, and Baryon Acoustic Oscillation Data

Pier-Stefano Corasaniti¹,² Ⓓ, Mauro Sereno³,⁴ Ⓓ, and Stefano Ettori³,⁴ Ⓓ

¹ LUTH, UMR 8102 CNRS, Observatoire de Paris, PSL Research University, Université Paris Diderot, 5 Place Jules Janssen, F-92190 Meudon, France
Pier-Stefano.Corasaniti@obpm.fr

² Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, F-75014 Paris, France

³ INAF, Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via Piero Gobetti 93/3, I-40129 Bologna, Italy
⁴ INFN, Sezione di Bologna, viale Berti Pichat 6/2, I-40127 Bologna, Italy

Received 2020 May 20; revised 2021 February 22; accepted 2021 February 22; published 2021 April 19

## Abstract

In recent years, the availability of large, complete cluster samples has enabled numerous cosmological parameter inference analyses using cluster number counts. These have provided constraints on the cosmic matter density $\Omega_m$ and the amplitude of matter density fluctuations $\sigma_8$ alternative to that obtained from other standard probes. However, systematics uncertainties, such as the mass calibration bias and selection effects, may still significantly affect these data analyses. Hence, it is timely to explore other proxies of galaxy cluster cosmology that can provide cosmological constraints complementary to those obtained from cluster number counts. Here we use measurements of the cluster sparsity from weak-lensing mass estimates of the LC²-single and HSC-XXL cluster catalogs to infer constraints on a flat $\Lambda$CDM model. The cluster sparsity has the advantage of being insensitive to selection and mass calibration bias. On the other hand, it primarily constrains a degenerate combination of $\Omega_m$ and $\sigma_8$ (along approximately constant curves of $S_8 = \sigma_8\sqrt{\Omega_m}/0.3$ and, to a lesser extent, the reduced Hubble parameter $h$. Hence, in order to break the internal parameter degeneracies, we perform a combined likelihood analysis of the cluster sparsity estimates with cluster gas mass fraction measurements and BAO data. We find marginal constraints that are competitive with those from other standard cosmic probes: $\Omega_m = 0.316 \pm 0.013$, $\sigma_8 = 0.757 \pm 0.067$ (corresponding to $S_8 = 0.776 \pm 0.064$), and $h = 0.696 \pm 0.017$ at $1\sigma$. Moreover, assuming a conservative Gaussian prior on the mass bias of gas mass fraction data, we find a lower limit on the gas depletion factor $Y_{b,500c} \gtrsim 0.89$.

Unified Astronomy Thesaurus concepts: Cosmological parameters from large-scale structure (340); Cosmology (343); Galaxy clusters (584)

Supporting material: machine-readable tables

## 1. Introduction

There is a widespread consensus that observations of galaxy clusters can provide a wealth of cosmological information (see, e.g., Allen et al. 2011; Kravtsov & Borgani 2012, for a review). In recent years, the potential to probe cosmology with galaxy clusters has been explored thanks to numerous survey programs. These have provided increasingly large cluster data sets from X-ray observations of the intracluster gas (see, e.g., Ebeling et al. 2010; Pierre et al. 2016; Böhringer et al. 2017), the detection of the Sunyaev–Zel’dovich (SZ) effect in the cosmic microwave background (CMB) radiation (see, e.g., Marriage et al. 2011; Planck Collaboration 2014; Bleem et al. 2015; Planck Collaboration 2016), and measurements of galaxy overdensities (Rykoff et al. 2016; Maturi et al. 2019).

The availability of complete cluster samples has enabled cosmological parameter inference analyses from cluster number count measurements (de Haan et al. 2016; Planck Collaboration 2016; Schellenberger & Reiprich 2017; Pacaud et al. 2018; Bocquet et al. 2019; Abbott et al. 2020; Lesci et al. 2020). Quite remarkably, these studies have consistently found values of the amplitude of linear matter density fluctuations on the $8 \text{ Mpc } h^{-1}$ scale $\sigma_8$ (h being the reduced Hubble parameter) and the cosmic matter density $\Omega_m$, which differs from that of the Planck primary CMB analysis (Planck Collaboration 2016, 2018). As an example, the analyses of the SZ cluster counts from Planck (Planck Collaboration 2014, 2016)

and the South Pole Telescope (SPT; Bocquet et al. 2019), as well as the cluster counts from the Dark Energy Survey Year 1 (DES-Y1; Abbott et al. 2020), have resulted in lower values of $S_8 = \sigma_8\sqrt{\Omega_m}/0.3$ compared to those from the Planck primary CMB. Similar results have been obtained using measurements of galaxy clustering from gravitational-lensing shear data (see, e.g., Hildebrandt et al. 2017; Joudaki et al. 2018; Troxel et al. 2018; Joudaki et al. 2020). However, in the case of cluster number counts, the tension may result from systematic effects, since the statistical significance of the discrepancy depends on the amplitude of the cluster mass calibration bias (see, e.g., Planck Collaboration 2016; Salvati et al. 2019). This is because the masses of the Planck-SZ cluster have been estimated using scaling relations calibrated on X-ray mass estimates. The latter are derived under the hydrostatic equilibrium (HE) hypothesis; consequently, any departure from the HE can introduce a systematic bias that propagates into the cosmological parameter inference analysis. However, the level of mass bias necessary to reconcile the Planck-SZ cluster counts with the Planck primary CMB results contrasts with expectations from hydrodynamical simulations of galaxy clusters (see, e.g., Rasia et al. 2012; Lau et al. 2013; Biffi et al. 2016; Barnes et al. 2020), as well as the bias estimated from the analyses of clusters for which accurate lensing or X-ray data are available (Hoekstra et al. 2015; Sereno et al. 2017; Eckert et al. 2019; Zubeldia & Challinor 2019; Makiya et al. 2020). Departures from the standard cosmological scenario can also account for such
discrepancies. As an example, the presence of massive neutrinos has been shown to alleviate the tension (Salvati et al. 2018). Similarly, selection effects cannot be a priori excluded (see, e.g., Chen et al. 2020, for a study of the impact of mass bias in weak-lensing shear-selected cluster samples). Hence, it is timely to investigate other galaxy cluster observables that can provide model parameter constraints alternative to those inferred from cluster number counts.

Cosmological information is encoded in the internal mass distribution of clusters. This is because the massive dark matter halos that host these structures have assembled through a hierarchical process that depends on the cosmic matter content, expansion rate, and amplitude of initial matter density fluctuations. The analysis of N-body simulations has shown that the density profile of halos is described by the Navarro–Frenk–White (NFW) formula (Navarro et al. 1997), such that for a halo of a given mass M, its profile only depends on the concentration parameter c. Numerical studies have subsequently shown that the median halo concentration as function of the halo mass depends on the specificities of the simulated cosmological model (see, e.g., Bullock et al. 2001; Zhao et al. 2003; Dolag et al. 2004; Zhao et al. 2009; Giocoli et al. 2012). This has suggested that measurements of the halo concentration from observations of a sample of galaxy clusters can provide cosmological constraints. However, the use of the c–M relation as cosmological proxy suffers of several drawbacks. First, astrophysical processes affecting the baryon distribution in the inner region of clusters may alter the estimated concentration–mass relation, thus inducing a systematic error in the cosmological analysis (Duffy et al. 2008; Mead et al. 2010; King & Mead 2011). Second, selection effects may have a strong impact on the cosmological parameter inference (Sereno et al. 2015). Hence, using measurements of the concentration of galaxy clusters to test cosmology has proven to be very challenging (see, e.g., Ettori et al. 2010).

Alternatively, Balmès et al. (2014) proposed that the sparsity, i.e., the ratio of the halo mass within radii enclosing different overdensities, can provide a nonparametric characterization of the mass distribution in halos while retaining the cosmological information encoded in the halo mass profile. More recently, Corasaniti et al. (2018) showed that the halo sparsity provides a cosmological proxy that is insensitive to selection and mass calibration bias. Furthermore, as the halo sparsity quantifies the excess of mass within a spherical shell comprised between two radii relative to the mass enclosed in the inner one, it is a sensitive probe of the screening mechanism in modified gravity scenarios (Corasaniti et al. 2020).

Here we infer cosmological parameter constraints on the $\Omega_m$, $\sigma_8$, and h of a flat $\Lambda$CDM model using estimates of the average sparsity of clusters from weak-lensing cluster mass in combination with cosmic distance measurements from baryon acoustic oscillations (BAOs) and cluster gas mass fractions from archival data.

The paper is organized as follows. In Section 2, we describe the methodology; in Section 3, we describe the data sets used; and in Section 4, we present the results. In Section 5, we discuss the conclusions.

## 2. Methodology

### 2.1. Halo Sparsity

Originally introduced by Balmès et al. (2014), the halo sparsity provides a nonparametric characterization of the mass

profile of halos in terms of the ratio of masses within radii enclosing two different overdensities (in units of the critical or background density), namely,

$$s_{\Delta_i, \Delta_2} = \frac{M_{\Delta_i}}{M_{\Delta_2}}, \qquad (1)$$

where $M_{\Delta_i}$ and $M_{\Delta_2}$ are the halo masses at overdensities $\Delta_i$ and $\Delta_2$, respectively, with $\Delta_i < \Delta_2$. As shown in Balmès et al. (2014), this provides an observational proxy of the cosmological information encoded in the halo mass profile. This is because the mass distribution within halos is the result of the hierarchical halo assembly process that depends on the growth and initial amplitude of the matter density fluctuations, which, in turn, depend on the cosmological parameters.

The analysis of N-body simulations has shown that the halo sparsity varies weakly with halo mass. This has very important consequences, since it implies that the halo ensemble average value of the sparsity at a given redshift and for a given cosmology can be predicted from the halo mass function at the overdensities of interest as given by Balmès et al. (2014),

$$\int_{M_{\Delta_2}^{\min}}^{M_{\Delta_2}^{\max}} \frac{dn}{dM_{\Delta_2}} d \ln M_{\Delta_2} = \langle s_{\Delta_1, \Delta_2} \rangle \int_{\langle s_{\Delta_1, \Delta_2} \rangle M_{\Delta_2}^{\min}}^{\langle s_{\Delta_1, \Delta_2} \rangle M_{\Delta_2}^{\max}} \frac{dn}{dM_{\Delta_1}} d \ln M_{\Delta_1}; \qquad (2)$$

this can be solved numerically for $\langle s_{\Delta_1, \Delta_2} \rangle$ given the functional form of the halo mass functions $dn/dM_{\Delta_1}$ and $dn/dM_{\Delta_2}$, respectively. The validity of Equation (2) has been extensively tested using halo catalogs from high-resolution large-volume N-body simulations in Corasaniti et al. (2018) and Corasaniti & Rasera (2019), and we refer readers to these publications for further details. In particular, for the halo masses used in this work, namely, $M_{200c}$ and $M_{500c}$, which correspond to masses within radii enclosing the overdensity $\Delta = 200$ and 500, respectively (in units of the critical density $\rho_c$), the analysis of the N-body halo catalogs has shown that Equation (2) is valid to better than the percent level (see, e.g., Table 1 in Corasaniti & Rasera 2019). Hence, Equation (2) sets a quantitative framework to perform cosmological parameter inference analyses using measurements of the average galaxy cluster sparsity.

Here we predict the average halo sparsity for a given cosmology at a given redshift by solving Equation (2). In the following, we assume a parameterization of the halo mass function as given by the Sheth–Tormen (ST) formula (Sheth & Tormen 1999),

$$\frac{dn}{dM} = \frac{\rho_m}{M} \left( -\frac{1}{\sigma} \frac{d\sigma}{dM} \right) A \frac{\delta_c}{\sigma} \sqrt{\frac{2a}{\pi}} \left[ 1 + \left( a \frac{\delta_c^2}{\sigma^2} \right)^{-p} \right] e^{-\frac{a\delta_c^2}{2\sigma^2}}, \qquad (3)$$

where $\rho_m$ is the mean matter density; $\delta_c$ is the linearly extrapolated spherical collapse threshold, which we compute using the formula by Kitayama & Suto (1996); and

$$\sigma^2(M, z) = \frac{1}{2\pi} \int dk \, k^2 \tilde{W}^2[kR(M)] P(k, z) \qquad (4)$$

is the variance of the linear density field smoothed on a spherical volume of radius R enclosing the mass
$M = 4/3\pi\rho_m R^3$ with

$$ \tilde{W}[kR] = \frac{3}{(kR)^3} [\sin(kR) - kR \cos(kR)], \quad (5) $$

and $P(k, z)$ being the linear matter power spectrum at a given redshift $z$, which we compute assuming the linear transfer function from Eisenstein & Hu (1998). Notice that due to the exponential cutoff of the halo mass function at large masses, the upper limits in the integrals of Equation (2) can be set to an arbitrarily large value without affecting the evaluation. Thus, we have set $M_{500c}^{\text{max}} = 10^{16} M_{\odot} h^{-1}$. In the case of the lower limits, we have set $M_{500c}^{\text{min}} = 10^{13} M_{\odot} h^{-1}$ consistently with the lowest mass in our cluster catalogs. However, we have found that the solution of Equation (2) is insensitive to the specific choice of $M_{500c}^{\text{min}}$ due the weak dependence of the halo sparsity on the halo mass.

In order to account for the redshift and cosmology dependence of the mass function at the overdensities of interest, we adopt the parameterization of the ST coefficients $a$, $p$, and $A$ given by Equation (12) in Despali et al. (2016), which we refer to as ST-Despali. As shown in Corasaniti et al. (2018), the average sparsity depends primarily on $\sigma_8$ and $\Omega_m$, which are highly degenerate, and to a lesser extent on $h$.

The computation of the average halo sparsity has been performed using a modified version of the numerical code Halo_Sparsity<fnref n="5" /> (Corasaniti 2021), which we have specifically developed for the work presented here.

## 2.2. Sparsity Systematic Errors

### 2.2.1. Halo Mass Function Parameterization

Using the N-body halo catalogs from the RayGalGroupSims simulation of a ΛCDM model, Corasaniti et al. (2018) showed that ST-Despali parameterization reproduces the average halo sparsity to better than a few percent accuracy in the redshift range $0.5 \lesssim z \lesssim 1.2$, while at lower and higher redshifts, differences increase up to ~15%. In the case of the Tinker et al. (2008) parameterization of the halo mass function, deviations from the N-body simulation results are much larger, and we do not consider them here. Instead, a parameterization of the ST coefficients on the RayGalGroupSims halo catalogs (see Equations (A4) and (A5) in Corasaniti et al. 2018) gives average sparsity estimates that reproduce the N-body simulation results to a subpercent level at all redshifts. We will refer to this parameterization as ST-RayGal.

As discussed in Corasaniti et al. (2018), the systematic deviation from the ST-Despali parameterization is twofold. On the one hand, the ST-Despali parameterization has been calibrated on halo catalogs from a suite of cosmological simulations that covers a volume that is eight times smaller than that of the RayGalGroupSims, thus resulting in a less accurate determination of the mass function at the high-mass end. On the other hand, the calibration has been realized using halo catalogs that contain different halos at different overdensities, rather than the same halos with masses estimated at different overdensities, as expected from Equation (2). Since modeling errors on the average sparsity prediction due to the halo mass function can induce a systematic bias on the

cosmological parameter analysis (see Section 3.1 in Corasaniti et al. 2018), we correct the ST-Despali prediction using the results from the halo catalogs of the RayGalGroupSims simulation (see Appendix A).

### 2.2.2. Radial Dependent Halo Mass Bias

The halo sparsity, being a mass ratio, is by definition exempt from constant mass calibration bias. Furthermore, being nearly independent of the halo mass, it is also insensitive to selection effects, as shown in Corasaniti et al. (2018). In contrast, the presence of a radial dependent mass bias can introduce a systematic error. As an example, astrophysical processes can alter the inner region of the dark matter halo profile while leaving the external regions unaltered. This can induce a radial dependent bias on the mass estimation. In Corasaniti et al. (2018), this has been investigated in the case of hydrostatic masses. Here we focus instead on weak-lensing mass measurements that probe the outer regions of the dark matter distribution in clusters differently from X-ray observations. We estimate the level of bias that such an effect induces on the sparsity $s_{200,500}$ using the mass bias estimates at $M_{200c}$ and $M_{500c}$ obtained in Lee et al. (2018, see their Figure 2) for halos with $M_{200c} \ge 10^{14} M_{\odot}$ from the catalogs of the Cosmo-OWLS simulations (Le Brun et al. 2014). These are N-body/hydro simulations that account for astrophysical feedback from supernovae and active galactic nuclei (AGNs). In the case of an extreme astrophysical feedback model (AGN 8.7 in Le Brun et al. 2014), we find that the relative variation of the halo sparsity varies from ~5% at the low-mass end ($M_{200c} \sim 10^{14} M_{\odot}$) to ~3% at the high-mass end ($M_{200c} \sim 10^{15} M_{\odot}$). In the case of a more realistic feedback scenario, for which the profile of the simulated clusters reproduces that of the observed X-ray clusters (AGN 8.0 in Le Brun et al. 2014), we find a smaller effect with a relative variation limited to $\lesssim 3\%$ over the same mass interval.

Another source of radial dependent mass bias in weak-lensing observations results from fitting a spherically symmetric NFW profile (Navarro et al. 1997) to tangential shear profile measurements. Deviations from the sphericity of the halo mass distribution and projection effects may induce systematic errors on the estimated mass at different radii. This has been studied in Becker & Kravtsov (2011). Using the mass bias estimates at $\Delta = 200$ and $500$ (in units of the critical density) quoted in their Tables 3 and 4, we find a relative variation of the sparsity $\Delta s_{200,500}/s_{200,500} \lesssim 3\%$ at $z = 0.25$ and $\lesssim 2\%$ at $z = 0.5$.

All of these effects are smaller compared to the systematic error due to the modeling of the mass function previously discussed and much smaller than the current uncertainties on weak-lensing mass determinations from shear profile measurement errors (see Section 3.1). These sources of bias are also small when compared to the intrinsic scatter in weak-lensing mass estimates, which is of the order of 20%–40% per cluster (Becker & Kravtsov 2011; Sereno & Ettori 2015).

## 2.3. Gas Mass Fraction

The halo sparsity is primarily sensitive to a degenerate combination of $\sigma_8$ and $\Omega_m$, while it has a weaker dependence on $h$ (Corasaniti et al. 2018). Hence, estimates of the cluster gas mass fraction, $f_{\text{gas}}$, from X-ray observations of galaxy clusters can break the internal parameter degeneracies. This is because

<footnote n="5">⁵ The code Halo_Sparsity is publicly available at https://github.com/pierste75/Halo_Sparsity.</footnote>
$f_{\text{gas}}$ probes the cosmic baryon fraction, i.e., the ratio of the mean baryon density $\Omega_b$ to the mean matter density $\Omega_m$, $f_b = \Omega_b/\Omega_m$, as well as the angular diameter distance $D_A(z)$ (see, e.g., White et al. 1993; Ettori et al. 2002; Allen et al. 2004; Ettori et al. 2009).

Following Allen et al. (2008), we model the relation between the gas mass fraction $f_{\text{gas},\Delta}$ estimated within a radius enclosing a given overdensity $\Delta$ and the baryon fraction $f_{\text{bar}}$ at a given redshift as

$$f_{\text{gas},\Delta} = K_{\Delta} Y_{b,\Delta} \frac{\Omega_b}{\Omega_m} \left[ \frac{D_A^{\text{fid}}(z)}{D_A(z)} \right]^{3/2} - f_{*,\Delta}, \quad (6)$$

where $D_A$ is the angular diameter distance, and $D_A^{\text{fid}}$ is the angular diameter distance of the fiducial cosmological model assumed to infer the gas mass fraction measurement. This derives from the fact that we focus on X-ray-based estimates such that $f_{\text{gas},\Delta} = M_{\Delta}^{\text{gas}}/M_{\Delta}^{\text{tot}}$, where $M_{\Delta}^{\text{tot}} = M_{\Delta}^{\text{HE}}/K_{\Delta}$ and $M_{\Delta}^{\text{HE}}$ is the cluster mass obtained assuming HE. Hence, differently from the halo sparsity, $f_{\text{gas},\Delta}$ also depends on uncertain quantities that parameterize the baryonic content of clusters, such as the mass calibration bias $K_{\Delta} \equiv 1 - b_{\Delta} = M_{\Delta}^{\text{HE}}/M_{\Delta}^{\text{tot}}$, the gas depletion factor $Y_{b,\Delta}$, and the stellar fraction $f_{*,\Delta}^6$. These play the role of nuisance parameters. To be as conservative as possible, we have tested the stability of our results for different assumptions on these parameters.

## 2.4. BAO Cosmic Distance

Cosmic distance measurements from BAO analyses provide cosmological parameter constraints that are complementary to those inferred from estimates of the cluster sparsity and gas mass fraction. Here we focus on measurements of the spherically averaged cosmic distance relative to the sound horizon at the drag epoch at different redshifts, $D_V(z)/r_d$. The spherically averaged cosmic distance is defined as

$$D_V(z) = \left[ \frac{cz}{H(z)} D_M^2(z) \right]^{1/3}, \quad (7)$$

where c is the speed of light, $H(z)$ is the Hubble function, and, in a flat universe,

$$D_M(z) = \int_0^z \frac{cdz'}{H(z')}. \quad (8)$$

We approximate the cosmological dependence of the sound horizon at the drag epoch as (see, e.g., Aubourg et al. 2015)

$$r_d \approx \frac{55.154 e^{-72.3 \cdot 0.0006^2}}{\omega_m^{0.25351} \omega_b^{0.12807}}, \quad (9)$$

with $\omega_{m,b} = \Omega_{m,b}h^2$. In Table 3, we quote the measurements we have used in our analysis.

# 3. Data Sets

## 3.1. Weak-lensing Cluster Mass Measurements

We estimate the sparsity of galaxy clusters using lensing mass measurements from publicly available cluster catalogs at $M_{500c}$ and $M_{200c}$. In particular, we consider a selected sample of

clusters from the Literature Catalogs of Lensing Clusters (LC$^2$; Sereno 2015) that includes as subsets the PSZ2Lens (Sereno et al. 2017) and the Cluster Lensing and Supernova Hubble (CLASH) project (Postman et al. 2012) catalogs. In addition, we consider the Subaru Hyper Suprime-Cam (HSC) lensing mass determinations (Umetsu et al. 2020) of an X-ray-selected cluster sample from the XXL-XMM survey (Pierre et al. 2016), which we refer to as the HSC-XXL sample.

The LC$^2$-single cluster catalog consists of weak-lensing cluster masses from archival data; at the time of this analysis, the latest version$^7$ included 672 entries at redshifts $z < 1.7$. We specifically focus on a selected subset of 187 clusters for which cluster masses at $M_{200c}$ and $M_{500c}$ have been inferred from a two-parameter fit of the shear lensing profile.$^8$ This is because mass estimates from a one-parameter fit, such as those obtained assuming the singular isothermal sphere (SIS) or the NFW profile with a fixed concentration-mass relation, result in a biased determination of the halo sparsity (see Appendix B). Out of this data set, we discard six clusters that we found to be outliers in the distribution of the sparsity of clusters within the same redshift bin (see Appendix C). Hereafter, we will refer to this subsample as Selected LC$^2$-single clusters containing 181 entries. This includes mass estimates from the PSZ2Lens catalog (Sereno et al. 2017), which is a statistically complete and homogeneous sample of 35 galaxy clusters at $z < 0.7$ from the second Planck Catalogue of Sunyaev-Zel'dovich Sources (Planck Collaboration 2016) with weak-lensing data from the Canada-France-Hawaii Telescope Lensing Survey (Heymans et al. 2012) and the Red Cluster Sequence Lensing Survey (Hildebrandt et al. 2016). The HSC-XXL lensing mass catalog (Umetsu et al. 2020) consists of 136 X-ray-selected clusters in the redshift range $0.031 < z < 1.033$. Therefore, our total sample contains 317 clusters (see Table 1).

In Figure 1, we plot $M_{200c}$ versus $M_{500c}$ for the clusters in the Selected LC$^2$-single sample, the PSZ2Lens subsample, and the HSC-XXL data set. Estimates of the correlation coefficient between the two mass measurements are only available for the PSZ2Lens subsample; these are shown in the inset plot. We see that, because the sparsity is very weakly dependent on the selection effects, data from homogeneous samples such as PSZ2Lens or HSC-XXL are in very good agreement with the heterogeneous Selected LC$^2$-single clusters. This appears more clearly in Figure 2, where we have plotted the average sparsity for the different samples in redshift bins of size $\Delta z = 0.1$ containing at least two clusters per bin. We can see that the different estimates are consistent with each other within the statistical errors. The uncertainties on the sparsity of individual clusters have been evaluated through the error propagation of the cluster mass uncertainties, where we have conservatively assumed a $r = 0.9$ correlation, a value that is smaller than the correlation estimated in the PSZ2Lens sample.$^9$ In the inset plot in Figure 2, we show the average sparsity in redshift bins of size $\Delta z = 0.1$ for the combined sample Selected LC$^2$-single+HSC-XXL, which we have used in the cosmological analysis. We have further tested the robustness of these

<footnote>$^6$ Here we assume that $f_{*,\Delta}$ is independent of the cluster mass.</footnote>

<footnote>$^7$ We use version 3.8, publicly available at http://pico.oabo.inaf.it/~sereno/CoMaLit/LC2/.</footnote>

<footnote>$^8$ For the clusters analyzed in Klein et al. (2019), we considered the masses from the two-parameter fits reported among the others in LC$^2$-all, rather than the one-parameter fit results listed in LC$^2$-single.</footnote>

<footnote>$^9$ It is worth remarking that, being a mass ratio, the error propagation of mass uncertainties on individual sparsity measurements results in smaller errors the larger the correlation between the mass estimates.</footnote>
<table><caption>**Table 1**
Galaxy Clusters from the Selected LC²-Single and HSC-XXL Catalogs</caption><thead><tr><th>Name</th><th>z</th><th>M<sub>200c</sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>M<sub>500c</sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>s<sub>200,500</sub></th></tr></thead><tbody><tr><td>400d-0159 + 0030</td><td>0.390</td><td>4.89 ± 2.15</td><td>3.45 ± 1.52</td><td>1.42 ± 0.28</td></tr><tr><td>ABELL-291</td><td>0.196</td><td>8.51 ± 2.06</td><td>4.51 ± 0.99</td><td>1.89 ± 0.20</td></tr><tr><td>ACT-CL-J0204.8-0303</td><td>0.549</td><td>4.37 ± 1.61</td><td>2.99 ± 1.10</td><td>1.46 ± 0.24</td></tr><tr><td>ACT-CL-J0205.2-0439</td><td>0.968</td><td>3.61 ± 2.23</td><td>2.56 ± 1.58</td><td>1.41 ± 0.39</td></tr><tr><td>PSZ2-G167.98-59.95</td><td>0.140</td><td>2.65 ± 1.30</td><td>1.79 ± 0.81</td><td>1.48 ± 0.32</td></tr><tr><td>ACT-CL-J0215.3-0343</td><td>1.004</td><td>5.64 ± 3.25</td><td>4.05 ± 2.33</td><td>1.39 ± 0.36</td></tr><tr><td>RCS-J0220.9-0333</td><td>1.030</td><td>3.94 ± 2.09</td><td>2.89 ± 1.21</td><td>1.37 ± 0.32</td></tr><tr><td>RCS-J0221.6-0347</td><td>1.020</td><td>1.54 ± 0.61</td><td>1.24 ± 0.46</td><td>1.24 ± 0.21</td></tr><tr><td>ACT-CL-J0221.7-0346</td><td>0.432</td><td>6.05 ± 1.59</td><td>3.91 ± 1.03</td><td>1.55 ± 0.18</td></tr><tr><td>XLSS-J022303.0-043622</td><td>1.220</td><td>27.74 ± 27.91</td><td>9.67 ± 5.62</td><td>2.87 ± 1.56</td></tr></tbody></table>

**Note.** Columns from left to right specify the cluster's name, redshift, weak-lensing estimated masses $M_{200c}$ and $M_{500c}$, and sparsity. This table is published in its entirety in machine-readable format. A portion is shown here for guidance regarding its form and content.

(This table is available in its entirety in machine-readable form.)

<figure><img src="image_2.png" /><figcaption>**Figure 1.** The $M_{500c}$ vs. $M_{200c}$ for the Selected LC²-single sample (red circles) and the HSC-XXL clusters (open triangles). The values of the correlation coefficients of the PSZ2Lens cluster masses (yellow squares) are shown in the inset plot.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Average cluster sparsity in redshift bins of size $\Delta z = 0.1$ for the Selected LC²-single (red points), PSZ2Lens (yellow squares), and HSC-XXL (open triangles) clusters. The average sparsity estimates for the combined sample Selected LC²-single and HSC-XXL (blue circles) are shown in the inset plot, and the number of clusters in each bin is given below every data point.</figcaption></figure>

<table><caption>**Table 2**
Galaxy Cluster Gas Mass Fractions</caption><thead><tr><th>Name</th><th>z</th><th>f<sub>gas</sub></th></tr></thead><tbody><tr><td>A1644</td><td>0.0473</td><td>0.128 ± 0.008</td></tr><tr><td>A85</td><td>0.0555</td><td>0.150 ± 0.005</td></tr><tr><td>A2319</td><td>0.0557</td><td>0.189 ± 0.008</td></tr><tr><td>A3266</td><td>0.0589</td><td>0.132 ± 0.009</td></tr><tr><td>A3158</td><td>0.0590</td><td>0.145 ± 0.007</td></tr><tr><td>A1795</td><td>0.0622</td><td>0.139 ± 0.005</td></tr><tr><td>RXC-J1825</td><td>0.0650</td><td>0.133 ± 0.005</td></tr><tr><td>A644</td><td>0.0704</td><td>0.132 ± 0.012</td></tr><tr><td>A2029</td><td>0.0766</td><td>0.141 ± 0.005</td></tr><tr><td>ZwCl-1215</td><td>0.0766</td><td>0.106 ± 0.008</td></tr></tbody></table>

**Note.** Columns from left to right specify the cluster's name, redshift, and gas mass fraction. This table is published in its entirety in machine-readable format. A portion is shown here for guidance regarding its form and content.

(This table is available in its entirety in machine-readable form.)

estimates by testing the validity of the average halo sparsity consistency relations (see Appendix D).

## 3.2. Cluster Gas Mass Fraction Data

Estimates of the gas mass fraction in clusters can be obtained from measurements of the temperature and density of the intracluster X-ray-emitting gas. Here we focus on measurements of $f_{gas}$ within $R_{500c}$. Our data set consists of gas mass fraction measurements of 12 clusters at $z \lesssim 0.1$ from X-COP (Eckert et al. 2019), a sample of 44 clusters in the range $0.1 \lesssim z \lesssim 0.3$ investigated in Ettori et al. (2010), and a higher-redshift sample of 47 clusters from Ghirardini et al. (2017) in the redshift range $0.4 \lesssim z \lesssim 1.2$. The clusters from these samples are massive systems with $M_{200c} \gtrsim 10^{14} M_{☉} h^{-1}$. The gas mass fraction measurements have been inferred through a backward model analysis, which fits the measured temperature profile against that predicted by solving the HE equation (Ettori et al. 2002). In addition, we have used gas mass fraction estimates of three high-redshift X-ray clusters in the redshift range $0.9 \lesssim z \lesssim 1.1$ derived in Bartalucci et al. (2018). A compilation of the data set is given in Table 2. The fiducial cosmology assumed in all of these analyses is a flat $\Lambda$CDM
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Averaged gas mass fraction estimates in redshift bins of size $\Delta z = 0.1$.</figcaption></figure>

model with $\Omega_m = 0.3$ and $h = 0.7$. We bin these measurements in redshift bins of size $\Delta z = 0.1$ (see Figure 3).

The sparsity and gas mass fraction measurements cover a similar redshift range. It is worth remarking that the X-COP sample, as well as the data set from Ettori et al. (2010), exploited XMM data, whereas Ghirardini et al. (2017) was based on Chandra observations. Gas mass measurements are mainly based on surface brightness profiles, which are highly consistent among different X-ray observatories, whereas the dependence on temperature and metallicity is negligible.

### 3.3. BAO Data

We use cosmic distance measurements from BAO analyses to further reduce the effect of cosmological parameter degeneracies. In Table 3, we quote the measurements we have used in our analysis. These consist of estimates from the 6dF Galaxy Survey (6dFGS; Beutler et al. 2011), Sloan Digital Sky Survey (SDSS) Data Release 7 Main Galaxy Sample (MGS; Ross et al. 2015), SDSS-III Baryon Acoustic Oscillation Spectroscopic Survey Data Release 12 (BOSS-DR 12; Alam et al. 2017), and SDSS-IV extended Baryon Acoustic Oscillation Spectroscopic Survey Data Release 16 Quasar Sample (eBOSS-QSO; Neveux et al. 2020), which span a similar redshift interval as that of cluster sparsity and gas mass fraction measurements.

## 4. Cosmological Parameter Inference

A first attempt to perform a cosmological analysis with measurements of the internal structures of halos (using the concentration–mass relation) in combination with the gas mass fraction data was presented in Ettori et al. (2010). As already mentioned, the use of the halo concentration as cosmological proxy presents several pitfalls compared to the halo sparsity. Furthermore, here we take advantage of the availability of a larger sample of cluster mass fraction data that cover a wider range of redshifts and thus are more sensitive to the dependence on $D_A$, as well as the latest cosmic distance estimates from BAO analyses.

### 4.1. Markov Chains and Priors

We perform a likelihood Markov Chain Monte Carlo (MCMC) data analysis of the redshift distribution of average cluster sparsity estimates in combination with BAO and gas mass fraction data to infer constraints on the flat $\Lambda$CDM model specified by the following set of parameters: $\Omega_m, \sigma_8$, and $h$. In order to reduce the impact of parameter degeneracies, we assume a set of Gaussian priors on the remaining cosmological parameters. More specifically, we adopt a Big Bang nucleosynthesis (Cyburt et al. 2016) prior on the baryon density

<table><caption>**Table 3**
Spherically Averaged Cosmic Distance Measurements from the 6dFGS (Beutler et al. 2011), SDSS Data Release 7 MGS (Ross et al. 2015), BOSS-DR 12 (Alam et al. 2017), and eBOSS-QSO (Neveux et al. 2020)</caption><thead><tr><td>Survey</td><td>z</td><td>D<sub>v</sub>/r<sub>d</sub></td></tr></thead><tbody><tr><td>6dFGS</td><td>0.106</td><td>2.976 ± 0.133</td></tr><tr><td>MGS</td><td>0.15</td><td>4.466 ± 0.168</td></tr><tr><td>BOSS-DR 12</td><td>0.38</td><td>9.994 ± 0.108</td></tr><tr><td>BOSS-DR 12</td><td>0.51</td><td>12.701 ± 0.128</td></tr><tr><td>BOSS-DR 12</td><td>0.61</td><td>14.481 ± 0.149</td></tr><tr><td>eBOSS-QSO</td><td>1.480</td><td>26.51 ± 0.42</td></tr></tbody></table>

$\Omega_b h^2 = 0.022 \pm 0.002$ and a Planck prior (Planck Collaboration 2018) on the scalar spectral index $n_s = 0.965 \pm 0.004$.

We sample the target parameter space assuming uniform priors on $\Omega_m \sim U(0.1, 0.9)$, $\sigma_8 \sim U(0.1, 1.8)$, and $h \sim U(0.55, 1.20)$. Where specified, we quote the results obtained under an HST prior, $h = 0.7324 \pm 0.0174$ (Riess et al. 2016), and a Planck prior, $h = 0.674 \pm 0.005$ (Planck Collaboration 2018).

In the case of $f_{\text{gas}}$, we assume priors on the nuisance parameters $K_{500c}$, $Y_{b,500c}$, and $f_{*_{500c}}$ within $R_{500c}$ (consistent with the definition of the gas mass fraction data) as determined from various works in the literature. We test the results of the parameter inference assuming a different prior on $K_{500c}$. In particular, we consider the following Gaussian priors: $K_{500c}^{\text{CLASH}} = 0.78 \pm 0.09$, consistent with the estimate from the analysis of the CLASH sample (Sereno & Ettori 2015); $K_{500c}^{\text{CCCP}} = 0.84 \pm 0.04$, as given by the analysis of a sample of clusters in Herbonnet et al. (2020) from the Canadian Cluster Comparison Project (CCCP; Hoekstra et al. 2015); and $K_{500c}^{\text{CMB}} = 0.65 \pm 0.04$, consistent with the mass bias inferred from the joint analysis of the Planck primary CMB, the Planck-SZ number counts, the Planck-thermal SZ power spectrum, and the BAO in Salvati et al. (2018). We also consider the extreme case with a hard prior $K_{500c} = 1$. As far as the baryon depletion factor is concerned, we assume a Gaussian prior $Y_{b,500c}^{\text{The 300}} = 0.938 \pm 0.041$, consistent with the analysis of simulated clusters in Eckert et al. (2019) from the Three Hundred Project (Cui et al. 2018). In order to evaluate the impact of a redshift variation of the baryon depletion factor, we test the case of a Gaussian prior with mean $Y_{b,500c}^{\text{FABLE}}(z) = 0.931(1 + 0.017 \cdot z + 0.003 \cdot z^2)$ and scatter $\sigma_{Y_{b,500c}} = 0.04$, consistent with the results from the FABLE simulations (Henden et al. 2020). Finally, we assume a Gaussian prior on the stellar fraction $f_{*_{500c}} = 0.015 \pm 0.005$ consistent with estimates from a sample of clusters with masses in the same range of those of the gas mass fraction data set (Eckert et al. 2019). Nevertheless, being a subdominant term in Equation (6), our results are largely independent of such a prior.

### 4.2. Results

We use the Metropolis-Hastings algorithm to generate 15 independent random chains of $3 \times 10^5$ samples each. We evaluate the rejection rate every 100 steps and adjust the width of the random step of the parameters dynamically. We check the convergence of the chains using the Gelman-Rubin test (Gelman & Rubin 1992) with a threshold value of $R < 1.1$. We derive marginal constraints on $\Omega_m, \sigma_8, h$, and $S_8 = \sigma_8 \sqrt{\Omega_m/0.3}$ from the analysis of the MCMC chains. Summary tables are presented in Appendix E.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Marginalized $1\sigma$ and $2\sigma$ contours in $\Omega_m - \sigma_8$ from the analysis of the average sparsity data only for the uniform-$h$ (black dashed lines), HST-$h$ (blue solid lines), and Planck-$h$ (green solid lines) and in combination with BAO cosmic distance measurements (black solid lines). The plus sign corresponds to the best-fit $\Lambda$CDM model with parameter values $\hat{\Omega}_m = 0.27$ and $\hat{\sigma}_8 = 0.85$ (and $\hat{h} = 0.66$). The red dotted lines correspond to constant $S_8$ value curves. For illustrative purposes, we show the marginalized contours from the cosmological analysis of the primary Planck analysis of the CMB anisotropy power spectra (TT, TE, EE+lowE+lensing; Planck Collaboration 2018) and the Planck-SZ cluster analysis for their baseline model with a BBN prior, BAO data, and a CCCP prior on the mass bias (Planck Collaboration 2016).</figcaption></figure>

### 4.2.1. Sparsity and BAO

In Figure 4, we plot the $1\sigma$ and $2\sigma$ contours in the $\Omega_m - \sigma_8$ plane from the analysis of the average sparsity measurements alone with a uniform-$h$ prior (black dashed lines) and under the HST (blue solid lines) and Planck (green solid lines) priors. We also show the constraints derived from the combined analysis of the average sparsity with the BAO (black solid lines). We can clearly see that the constraints from the average sparsity data only are highly degenerate along curves of approximately constant $S_8$ values (red dotted lines).

In Table 4, we quote the mean and standard deviation of $S_8$ for the Selected LC$^2$-single clusters under different $h$ priors, as well as the results obtained for the PSZ2Lens subsample. In the former case, we find $S_8 = 0.75 \pm 0.20$ for the uniform-$h$ prior while assuming that the HST prior (Planck prior) gives $S_8 = 0.80 \pm 0.18$ ($S_8 = 0.82 \pm 0.16$). The mean value of $S_8$ increases for decreasing values of $h$ associated with the HST and Planck priors. The average sparsity also carries information on the Hubble parameter. This can be seen in Figure 4, where the contours inferred from the analysis of the sparsity data alone shift toward larger $S_8$ values for the HST-$h$ and Planck-$h$ priors. As can be deduced from the values quoted in Table 4, the analysis of the average sparsity from the PSZ2Lens subsample gives results (bottom rows) for the different $h$ priors that are largely consistent with those from the full cluster catalog, though with slightly larger $1\sigma$ errors (of order ~10%) due to the smaller size and reduced redshift interval of the subsample.

Overall, the marginal constraints on $S_8$ from the analysis of the cluster sparsity are consistent with those inferred from the sparsity analysis of a sample of X-ray clusters presented in

<caption>**Table 4**</caption>

Mean and Standard Deviation of $S_8$ from the Analysis of the Average Sparsity of the Selected LC$^2$-Single Catalog (Top Rows) and the PSZ2Lens Subsample (Bottom Rows)

<table><thead><tr><th colspan="2">Sparsity Only</th></tr><tr><th>Cluster Sample and h Prior</th><th>S<sub>8</sub></th></tr></thead><tbody><tr><td>Selected LC<sup>2</sup>-Single + uniform</td><td>0.75 ± 0.20</td></tr><tr><td>Selected LC<sup>2</sup>-Single + HST</td><td>0.80 ± 0.18</td></tr><tr><td>Selected LC<sup>2</sup>-Single + Planck</td><td>0.82 ± 0.16</td></tr><tr><td>PSZ2Lens + uniform</td><td>0.69 ± 0.22</td></tr><tr><td>PSZ2Lens + HST</td><td>0.73 ± 0.21</td></tr><tr><td>PSZ2Lens + Planck</td><td>0.75 ± 0.19</td></tr></tbody></table>

Corasaniti et al. (2018), where it was found that $S_8 = 0.87 \pm 0.20$ for the HST-$h$ prior.

Including the BAO data<fnref n="10" /> allows one to break the degeneracy along the $\Omega_m$ direction and infer close bounds on $\Omega_m$ and $\sigma_8$. Moreover, because of the dependence of the average sparsity on the Hubble parameter, we expect the joint analysis with the BAO data to provide closed bounds on $h$. In order to have an appreciation of the derived constraints with respect to those inferred from other probes, in Figure 4, we plot the $1\sigma$ and $2\sigma$ marginalized contours from the Planck primary CMB cosmological analysis (red and yellow filled contours) and the Planck-SZ cluster counts (blue and light blue filled contours). In Table 5, we quote the results of the marginal statistics from the joint sparsity and BAO analysis obtained under different $h$ priors.

We find the following marginal constraints (mean and standard deviation): $\Omega_m = 0.277 \pm 0.029$, $\sigma_8 = 0.856 \pm 0.100$ (corresponding to $S_8 = 0.818 \pm 0.070$), and $h = 0.662 \pm 0.023$. The marginalized posteriors are well approximated by Gaussian distributions; consequently, we find the best-fit model parameters to be close to the parameter average values. In particular, we have $\hat{\Omega}_m = 0.274$, $\hat{\sigma}_8 = 0.850$ (corresponding to $\hat{S}_8 = 0.813$), and $\hat{h} = 0.662$. The residuals are shown in Figure 5.

The joint analysis results in an improvement of the cosmological constraints with respect to those obtained from the BAO data alone of a factor of 3 on $\Omega_m$ and a factor of ~2 on $h$. On the other hand, imposing the HST-$h$ (Planck-$h$) prior does not significantly improve the uncertainties on $\Omega_m$ and $\sigma_8$.

These results are consistent to better than $1\sigma$ with those inferred from the Planck primary CMB cosmological analysis of the anisotropy power spectra (red and yellow shaded contours): $\Omega_m = 0.315 \pm 0.007$, $\sigma_8 = 0.811 \pm 0.006$, $S_8 = 0.832 \pm 0.013$, and $h = 0.6736 \pm 0.0054$ (TT, TE, EE+lowE+lensing; see Table 2 in Planck Collaboration 2018).

The constraints are also consistent with those of the Planck-SZ cluster count baseline model<fnref n="11" /> with $\Omega_m = 0.330 \pm 0.030$, $\sigma_8 = 0.760 \pm 0.033$, and $S_8 = 0.796 \pm 0.042$, which were derived under similar assumptions (most notably the combination of the BBN prior and BAO data).

The comparison indicates that Planck-SZ cluster counts provide tighter constraints on $\sigma_8$, while the bounds on $\Omega_m$ have the same level of statistical uncertainty. Limits on $h$ were derived in

<footnote n="10">¹⁰ We infer the following cosmological parameter constraints from the MCMC likelihood analysis of the BAO data alone for a flat $\Lambda$CDM model: $\Omega_m = 0.387 \pm 0.089$ and $h = 0.718 \pm 0.043$.</footnote>

<footnote n="11">¹¹ We are grateful to Richard Betty for providing us with the MCMC chains of the CCCP+BAO+BBN baseline model from the Planck-SZ data analysis (Planck Collaboration 2016).</footnote>
<table><caption>**Table 5**
Mean and Standard Deviation of $\Omega_m$, $\sigma_8$, $S_8$, and $h$ from the Joint Analysis of the Average Sparsity and BAO Data</caption><thead><tr><th rowspan="2">h Prior</th><th colspan="4">Sparsity + BAO</th></tr><tr><th>Ω<sub>m</sub></th><th>σ<sub>8</sub></th><th>S<sub>8</sub></th><th>h</th></tr></thead><tbody><tr><td>Uniform</td><td>0.277 ± 0.029</td><td>0.856 ± 0.100</td><td>0.818 ± 0.070</td><td>0.662 ± 0.023</td></tr><tr><td>HST</td><td>0.318 ± 0.029</td><td>0.756 ± 0.089</td><td>0.775 ± 0.071</td><td>...</td></tr><tr><td>Planck</td><td>0.285 ± 0.023</td><td>0.833 ± 0.088</td><td>0.809 ± 0.069</td><td>...</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Residuals of the average sparsity (top panel) and BAO (bottom panel) data with respect to the best-fit ΛCDM model with parameter values $\hat{\Omega}_m = 0.274$, $\hat{\sigma}_8 = 0.850$, and $\hat{h} = 0.662$.</figcaption></figure>

Planck Collaboration (2016), though the values of the marginal statistics were not quoted. Here we have analyzed the chains that were made available to us and obtained $h = 0.693 \pm 0.027$, which is consistent with the constraints on $h$ we have inferred from the joint sparsity and BAO analysis. Nonetheless, it is worth noticing that the constraints from the Planck-SZ analysis are dominated by the BAO data set used in Planck Collaboration (2016), as no constraints would be inferred from the use of the cluster counts only (J. B. Melin, 2021 private communication). This is not the case for the BAO data used here.

The comparison with other cluster number count analyses is less straightforward, since constraints on the cosmological parameters were derived under very different priors.

### 4.2.2. Sparsity and Gas Mass Fraction

Cluster gas mass fraction data can be used to break the $\Omega_m - \sigma_8$ degeneracy that characterizes the constraints from the average sparsity. However, different from the BAO, the possibility of deriving robust constraints on $\Omega_m$, $\sigma_8$, and $h$ depends on the sensitivity of the gas mass fraction data to the priors assumed on the baryon depletion factor and the mass calibration bias. This is

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Marginalized 1$\sigma$ and 2$\sigma$ contours in $\Omega_m - \sigma_8$ from the combined analysis of the average sparsity and gas mass fraction data for different sets of priors: $Y_{b,500c}^{\text{The300}}$ and $K_{500c}^{\text{CMB}}$ (dotted black lines), $Y_{b,500c}^{\text{The300}}$ and $K_{500c}^{\text{CLASH}}$ (dashed-dotted black lines), $Y_{b,500c}^{\text{FABLE}}$ and $K_{500c}^{\text{CLASH}}$ (dashed-dotted blue lines), $Y_{b,500c}^{\text{The300}}$ and $K_{500c}^{\text{CCCP}}$ (dashed black lines), and $Y_{b,500c}^{\text{The300}}$ with $K_{500c} = 1$ (solid black lines). The red shaded region corresponds to curves of constant $S_8 = 0.83 \pm 0.07$ values, that is, the mean and standard deviation of the different $S_8$ estimates quoted in Table 6. As in Figure 4, we plot the contours from the Planck primary CMB and Planck-SZ cluster count analyses, respectively. The inset plots show the mean and standard deviation of $h$ (left inset) and $\Omega_m$ (right inset) as function of the $K_{500c}$ prior. In the left inset, the light blue and blue shaded areas correspond to the HST-$h$ and Planck-$h$ priors, respectively.</figcaption></figure>

because such parameters appear as multiplicative factors in Equation (6), thus resulting in the $h^{3/2}$ dependence being highly degenerate with $\Omega_b/\Omega_m$, $K_{500c}$, and $Y_{b,500c}$.

In Figure 6, we plot the 1$\sigma$ and 2$\sigma$ contours in the $\Omega_m - \sigma_8$ plane from the joint analysis of the average sparsity and gas mass fraction data obtained under different $K_{500c}$ and $Y_{b,500c}$ priors. The marginal constraints on $\Omega_m$, $\sigma_8$, $S_8$, and $h$ are summarized in Table 6. The red shaded area in Figure 6 highlights the degeneracy of the average sparsity along curves of constant $S_8 = 0.83 \pm 0.07$ values, the mean and standard deviation of the $S_8$ estimates for the various prior configurations quoted in Table 6 (top rows). The inset plots show the marginal statistics of $h$ (left inset) and $\Omega_m$ (right inset) as a function of the $K_{500c}$ prior for the $Y_{b,500c}^{\text{The300}}$ case. We may notice that, depending on the $K_{500c}$ prior, the gas mass fraction data break the $\Omega_m - \sigma_8$ degeneracy along different locations of the $S_8$ curve. This is because, for a given $Y_{b,500c}$ and $\Omega_b$ prior, there
<table><caption>**Table 6**
Mean and Standard Deviation of $\Omega_m$, $\sigma_8$, $S_8$, and $h$ from the Joint Analysis of the Average Sparsity and Gas Mass Fraction Data for Different $Y_{b,500c}$, $K_{500c}$, and $h$ Priors</caption><thead><tr><th rowspan="2">f<sub>gas</sub> and h Priors</th><th colspan="4">Sparsity + Gas Mass Fraction</th></tr><tr><th>Ω<sub>m</sub></th><th>σ<sub>8</sub></th><th>S<sub>8</sub></th><th>h</th></tr></thead><tbody><tr><td>Y<sup>The 300</sup><sub>b,500c</sub> + K<sup>CMB</sup><sub>500c</sub> + uniform</td><td>0.13 ± 0.01</td><td>1.35 ± 0.18</td><td>0.90 ± 0.09</td><td>1.10 ± 0.05</td></tr><tr><td>Y<sup>The 300</sup><sub>b,500c</sub> + K<sup>CLASH</sup><sub>500c</sub> + uniform</td><td>0.19 ± 0.03</td><td>1.04 ± 0.14</td><td>0.82 ± 0.08</td><td>0.91 ± 0.06</td></tr><tr><td>Y<sup>FABLE</sup><sub>b,500c</sub> + K<sup>CLASH</sup><sub>500c</sub> + uniform</td><td>0.19 ± 0.03</td><td>1.06 ± 0.15</td><td>0.83 ± 0.08</td><td>0.92 ± 0.07</td></tr><tr><td>Y<sup>The 300</sup><sub>b,500c</sub> + K<sup>CCCP</sup><sub>500c</sub> + uniform</td><td>0.20 ± 0.02</td><td>1.03 ± 0.12</td><td>0.83 ± 0.07</td><td>0.89 ± 0.05</td></tr><tr><td>Y<sup>The 300</sup><sub>b,500c</sub> (K<sub>500c</sub> = 1) + uniform</td><td>0.25 ± 0.03</td><td>0.86 ± 0.10</td><td>0.79 ± 0.07</td><td>0.79 ± 0.05</td></tr><tr><td>Y<sup>The 300</sup><sub>b,500c</sub> + K<sup>CMB</sup><sub>500c</sub> + HST</td><td>0.16 ± 0.01</td><td>1.27 ± 0.04</td><td>0.92 ± 0.03</td><td>...</td></tr><tr><td>Y<sup>The 300</sup><sub>b,500c</sub> + K<sup>CLASH</sup><sub>500c</sub> + HST</td><td>0.26 ± 0.03</td><td>0.85 ± 0.10</td><td>0.79 ± 0.07</td><td>...</td></tr><tr><td>Y<sup>FABLE</sup><sub>b,500c</sub> + K<sup>CLASH</sup><sub>500c</sub> + HST</td><td>0.26 ± 0.02</td><td>0.87 ± 0.07</td><td>0.81 ± 0.05</td><td>...</td></tr><tr><td>Y<sup>The 300</sup><sub>b,500c</sub> + K<sup>CCCP</sup><sub>500c</sub> + HST</td><td>0.24 ± 0.02</td><td>0.94 ± 0.11</td><td>0.83 ± 0.08</td><td>...</td></tr><tr><td>Y<sup>The 300</sup><sub>b,500c</sub> (K<sub>500c</sub> = 1) + HST</td><td>0.27 ± 0.03</td><td>0.83 ± 0.09</td><td>0.79 ± 0.07</td><td>...</td></tr><tr><td>Y<sup>The 300</sup><sub>b,500c</sub> + K<sup>CMB</sup><sub>500c</sub> + Planck</td><td>0.17 ± 0.01</td><td>1.25 ± 0.04</td><td>0.94 ± 0.03</td><td>...</td></tr><tr><td>Y<sup>The 300</sup><sub>b,500c</sub> + K<sup>CLASH</sup><sub>500c</sub> + Planck</td><td>0.30 ± 0.03</td><td>0.78 ± 0.11</td><td>0.78 ± 0.06</td><td>...</td></tr><tr><td>Y<sup>FABLE</sup><sub>b,500c</sub> + K<sup>CLASH</sup><sub>500c</sub> + Planck</td><td>0.28 ± 0.01</td><td>0.85 ± 0.06</td><td>0.83 ± 0.01</td><td>...</td></tr><tr><td>Y<sup>The 300</sup><sub>b,500c</sub> + K<sup>CCCP</sup><sub>500c</sub> + Planck</td><td>0.23 ± 0.02</td><td>0.99 ± 0.01</td><td>0.87 ± 0.03</td><td>...</td></tr><tr><td>Y<sup>The 300</sup><sub>b,500c</sub> (K<sub>500c</sub> = 1) + Planck</td><td>0.31 ± 0.03</td><td>0.74 ± 0.10</td><td>0.75 ± 0.07</td><td>...</td></tr></tbody></table>

is a compensation between the fitting values of $h$ and $\Omega_m$ and the prior value of $K_{500c}$. In particular, the smaller the $K_{500c}$ prior, the larger the value of $h$ and the smaller the value of $\Omega_m$ necessary to fit the same gas mass fraction data. As a result, the gas mass fraction data break the $S_8$ degeneracy of the cluster sparsity at larger $\sigma_8$ values for smaller $K_{500c}$ priors. This is consistent with the trend shown in Figure 6, where the contours shift from the upper left to the lower right for increasing values of the $K_{500c}$ prior, while the average inferred value of $h(\Omega_m)$ shown in the left (right) inset plot decreases (increases). For comparison, we have tested the $Y_{b,500c}^{\text{FABLE}}$ prior for the $K_{500c}^{\text{CLASH}}$ case and find no statistically significant differences with respect to the $Y_{b,500c}^{\text{The300}}$ prior. As can be seen in Figure 6, only the constraints derived under the hard $K_{500c} = 1$ prior are marginally consistent with those from the Planck primary CMB analysis.

We find that imposing the HST and Planck priors tends to slightly shift the parameter constraints toward larger $\Omega_m$ and smaller $\sigma_8$ values (see values quoted in the middle and bottom rows in Table 6). Nevertheless, the inconsistencies among the different gas mass fraction prior model parameter assumptions remain unsolved. Most notably, in the case of the $K_{500c}^{\text{CMB}}$ prior, we find constraints on $\Omega_m$ and $\sigma_8$ that lie $\gtrsim 2\sigma$ from the Planck results. This suggests that the low value of $K_{500c}^{\text{CMB}}$ stands in contrast to not only expectations from numerical simulations of clusters and direct estimates of the mass calibration bias from galaxy cluster samples but also the gas mass fraction measurements we have considered here.

Given the large scatter in the prior values of the baryon depletion factor and mass calibration (and the impact that such priors have on the cosmological parameter inference), it is advisable to simply treat them as nuisance parameters and marginalize over large uniform priors. However, because of the degeneracy with $h$, this may result in weak cosmological parameter constraints. Hence, an external independent data set such as the BAO is required to isolate the $h$ dependence of gas mass fraction data and infer tighter bounds on the other cosmological parameters. This is what we investigate next.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Marginalized 1D posteriors on $\Omega_m$ (top left panel), $\sigma_8$ (top right panel), $h$ (bottom left panel), and $S_8$ (bottom right panel) from the joint analysis of cluster sparsity, gas mass fraction, and BAO data.</figcaption></figure>

### 4.2.3. Joint Analysis

Here we present the results of the joint likelihood analysis of average cluster sparsity, gas mass fraction, and BAO data. To be as conservative as possible, we infer marginalized bounds over uniform priors on the gas mass fraction nuisance parameters. More specifically, we uniformly sample the gas depletion factor $Y_{b,500c} \sim U(0.7, 1.0)$ and the mass calibration bias $K_{500c} \sim U(0.6, 1.0)$ over interval values that include the Gaussian priors previously discussed.

In Figure 7, we show the 1D marginalized posteriors, and in Figure 8, we plot the marginalized $1\sigma$ and $2\sigma$ contours in the $\Omega_m-\sigma_8$ plane. As in Figure 4, we plot the contours from the Planck primary CMB and the Planck-SZ cluster count cosmological analyses. In Table 7, we quote the results of the marginal statistics of $\Omega_m$, $\sigma_8$, $S_8$, and $h$.

We have $\Omega_m = 0.316 \pm 0.013$ and $\sigma_8 = 0.757 \pm 0.067$, which corresponds to $S_8 = 0.776 \pm 0.064$ and $h = 0.6958 \pm 0.0167$. We find the best-fit model parameters to be $\hat{\Omega}_m = 0.320$, $\hat{\sigma}_8 = 0.738$ (corresponding to $\hat{S}_8 = 0.763$), and
<table><caption>**Table 7**
Mean and Standard Deviation of $\Omega_m$, $\sigma_8$, $S_8$, and $h$ from the Joint Analysis of the Average Sparsity, Gas Mass Fraction, and BAO Data</caption><thead><tr><th rowspan="2">h Prior</th><th colspan="3">Sparsity + Gas Mass Fraction + BAO</th><th rowspan="2">h</th></tr><tr><th>Ω<sub>m</sub></th><th>σ<sub>8</sub></th><th>S<sub>8</sub></th></tr></thead><tbody><tr><td>Uniform</td><td>0.316 ± 0.013</td><td>0.757 ± 0.067</td><td>0.776 ± 0.064</td><td>0.696 ± 0.017</td></tr><tr><td>HST</td><td>0.325 ± 0.012</td><td>0.736 ± 0.066</td><td>0.766 ± 0.065</td><td>...</td></tr><tr><td>Planck</td><td>0.306 ± 0.010</td><td>0.781 ± 0.061</td><td>0.788 ± 0.059</td><td>...</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Marginalized $1\sigma$ and $2\sigma$ contours in the $\Omega_m-\sigma_8$ plane from the combined analysis of the average cluster sparsity, gas mass fraction, and BAO data (black lines). As in Figure 4, we plot marginalized contours from the Planck primary CMB analysis (yellow and red contours) and the Planck-SZ number counts (dark and light blue contours). The plus sign corresponds to the best-fit $\Lambda$CDM model with parameter values $\hat{\Omega}_m = 0.320$ and $\hat{\sigma}_8 = 0.738$ (and $\hat{h} = 0.690$).</figcaption></figure>

$\hat{h} = 0.690$. The data residuals with respect to the best-fit model are shown in Figure 9. Notice that the 1D marginalized posterior of $\sigma_8$ (and, consequently, $S_8$) has a slightly fatter tail than that of $\Omega_m$ and $h$, which are well approximated by Gaussian distributions. This is the reason for the small difference between the inferred average value of $\sigma_8$ (and $S_8$) and the best-fit value of $\hat{\sigma}_8$ (and $\hat{S}_8$) associated with the peak of the marginalized posterior. Compared to the constraints inferred from the joint analysis of the average sparsity with BAO, the addition of the gas mass fraction data improves the constraints on $\Omega_m$ by a factor of 2. The constraints on $\sigma_8$ and $h$ are improved with a gain on the $1\sigma$ errors at the ~49% and ~35% level, respectively. Quite remarkably, we infer $1\sigma$ constraints on $h$ that are competitive with those from the HST analysis (Riess et al. 2016). As we can see in Figure 9, the results of the joint analysis are consistent within $1\sigma$ with both the Planck primary CMB and the Planck-SZ cluster count analyses. Adding priors on $h$ does not significantly improve the determination of $\Omega_m$ and $\sigma_8$.

As expected, we find the gas mass fraction nuisance parameters to be unconstrained in the prior parameter interval due to the fact that they are perfectly degenerate. Hence, assuming a Gaussian prior $K_{500c}^{\text{CLASH}}$ ($K_{500c}^{\text{CCCP}}$) does not noticeably improve the constraints on the cosmological parameters. In contrast, it allows one to infer a lower limit on $Y_{b,500c} \gtrsim 0.89$ ($Y_{b,500c} \gtrsim 0.91$) at $1\sigma$.

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Residuals of the average sparsity (top panel), gas mass fraction (middle panel), and BAO (bottom panel) data with respect to the best-fit flat $\Lambda$CDM model with parameter values $\hat{\Omega}_m = 0.320$, $\hat{\sigma}_8 = 0.738$, and $\hat{h} = 0.690$, respectively.</figcaption></figure>

In Figures 10 and 11, we show summary plots of state-of-the-art estimates of $S_8$ (see also Pratt et al. 2019) and $h$ from various probes.

## 5. Conclusions

In recent years, the increased availability of large, complete samples has opened the way to probing cosmology with galaxy cluster observations. Cosmological parameter constraints have been primarily inferred from cluster number count data analyses. These have provided constraints complementary to those inferred from other standard probes, such as the CMB anisotropy power spectra. However, several sources of systematics may still affect the results of these studies. As an example, cluster number count measurements from SZ catalogs have resulted in values of $S_8$ that are lower than those obtained from the analysis of the Planck anisotropy power spectra (see, e.g., Planck Collaboration 2016; Bocquet et al. 2019; Abbott et al. 2020). Errors in the mass calibration of clusters, as well as selection effects, may be responsible for such discrepancies, although it cannot be a priori excluded that these are the effect
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Constraints on $S_8$ from the analysis of the shear lensing data (blue triangles) of KV450 (Hildebrandt et al. 2018), DES-Y1 (Troxel et al. 2018), KV450+DES-Y1 (Joudaki et al. 2020), and KV450+BOSS (Tröster et al. 2020) and the cluster number counts (red squares) of Planck-SZ (Planck Collaboration 2016), HIPLUGCS (Schellenberger & Reiprich 2017), SPT-SZ (Bocquet et al. 2019), XXL (Pacaud et al. 2018), and DES-Y1 (Abbott et al. 2020). The constraints inferred from the combined analysis of cluster sparsity and BAO (upper yellow circle) and together with the gas mass fraction (lower yellow circle) data are shown at the top of the plot. The shaded area corresponds to the $1\sigma$ and $2\sigma$ limits from the Planck primary CMB analysis (TT, TE, EE+lowE+lensing; Planck Collaboration 2018).</figcaption></figure>

of novel physics not in the standard $\Lambda$CDM model. For this reason, it is timely to investigate other probes of galaxy cluster cosmology that may provide independent constraints on the cosmological parameters. Besides cluster number counts, measurements of the spatial clustering of galaxy clusters (Marulli et al. 2018), as well as their internal mass distribution, carry cosmological information. As an example, Balmès et al. (2014) showed that measurements of the dark matter halo sparsity can retrieve the cosmological signal encoded in the mass profile of halos hosting galaxy clusters. Moreover, these measurements have been shown to be insensitive to selection effects and mass calibration errors (Corasaniti et al. 2018).

Here we have derived cosmological constraints from estimates of the average sparsity of galaxy clusters using lensing mass measurements. We have used lens mass measurements from a selected sample of clusters from the LC²-single catalog (Sereno 2015), together with HSC-XXL clusters (Umetsu et al. 2020). We have discussed different sources of systematic errors. In order to break cosmological parameter degeneracies, we have performed a combined MCMC likelihood analysis of average sparsity estimates in

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Constraints on $h$ from various probes: HST (Riess et al. 2016), DES-Y1+BAO+D/H (Abbott et al. 2018), SPTpol (Henning et al. 2018), H0LiCOW (Bonvin et al. 2017), SHOES (Riess et al. 2018), CCHP (Freedman et al. 2019), SPARC (Schombert et al. 2020), TDCOSMO-SLACS (Birrer et al. 2020), and 8TDGL (Denzel et al. 2021). The constraints inferred from the combined analysis of cluster sparsity and BAO (upper yellow circle) and together with the gas mass fraction (lower yellow circle) data are shown at the top of the plot. The shaded area corresponds to the $1\sigma$ and $2\sigma$ limits from the Planck primary CMB analysis (TT, TE, EE+lowE+lensing; Planck Collaboration 2018).</figcaption></figure>

combination with gas mass fraction and BAO data. We find that the combination of these data sets provides competitive constraints on $\Omega_m$, $\sigma_8$, and $h$ that are consistent with those from the Planck primary CMB and Planck-SZ cluster counts. In particular, we find $\Omega_m = 0.316 \pm 0.013$, $\sigma_8 = 0.757 \pm 0.067$ (resulting in $S_8 = 0.776 \pm 0.064$), and $h = 0.696 \pm 0.017$.

Compared to other cosmological proxies from galaxy cluster observations, the sparsity has the advantage of being less sensitive to systematic errors due to mass calibration bias and selection effects. In the future, cosmological parameter uncertainties can be further reduced thanks to the availability of larger galaxy cluster samples with improved mass estimates, as well as better gas mass fraction measurements. It is possible that the combined analysis of cluster sparsity and gas mass fraction data, together with estimates of the cluster abundance, may not only provide strong cosmological parameter constraints but also allow for an accurate determination of astrophysical model parameters, such as the gas depletion factor and the mass calibration bias, that are sources of systematic uncertainty associated with gas mass fraction and cluster number count data analyses. We leave this to a future investigation.
A first step in this direction will be the joint study of the sparsity and gas mass fraction in galaxy clusters from a well-selected sample of a large number (~100) of objects homogeneously analyzed in their X-ray and lensing signal out to $R_{500}$ and beyond, such as the one that will be available from the CHEX-MATE project<fnref n="12" /> (CHEX-MATE Collaboration et al. 2020).

P.S.C. is grateful to Iacopo Bartalucci for providing data from his X-ray cluster analysis, Richard Battye for providing the Planck-SZ chains, Jean-Baptiste Melin for discussions on the Planck-SZ analysis, and Florian Pacaud for providing the values of $S_8$ obtained from the XXL cluster number counts. S.E. and M.S. acknowledge financial contributions from contract ASI-INAF n.2017-14-H.0 and INAF “Call per interventi aggiuntivi a sostegno della ricerca di main stream di INAF” No. 1.05.01.86.10. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007–2013)/ERC grant agreement No. 279954. Data visualizations were prepared with the MATPLOTLIB<fnref n="13" /> library (Hunter 2007).

## Appendix A
Average Sparsity Redshift Model Correction

The validity of Equation (2) in predicting the ensemble average halo sparsity has been extensively tested using N-body halo catalogs from the RayGalGroupSims simulation of a ΛCDM model in Corasaniti et al. (2018) and the halo catalogs from the MultiDark-Planck2 simulation in Corasaniti & Rasera (2019). These studies have shown that the simulation calibrated mass functions can predict the average halo sparsity to better than a few percent accuracy. On the other hand, as already mentioned in Section 2, assuming the ST-Despali parameterization, Corasaniti et al. (2018) found that the predicted average sparsity is consistent with the RayGalGroupSims results within a few percent only in the redshift range $0.5 \lesssim z \lesssim 1.2$, while differences increase up to ~10% at lower and higher redshifts. Such a discrepancy is essentially due to the fact that the halo catalogs used in the calibration of the ST-Despali parameterization contain different halos at different overdensities, rather than identical halos with masses estimated at different overdensities, as expected from Equation (2). This induces a selection effect on the calibrated mass function at the two overdensities of interest that propagates into average sparsity prediction. This can be seen in Figure 12, where we plot the difference between the value of the average sparsity from the ST-Despali parameterization and that measured from the N-body halos at the redshift snapshots of the RayGalGroupSims (blue circles) and the MultiDark-Planck2 (red triangles) simulations. Quite importantly, we may notice that the systematic difference of the ST-Despali prediction with respect to these N-body results is similar, despite the fact that the two simulations do not share the same cosmology or possess similar characteristics,<fnref n="14" /> thus highlighting the algorithmic nature of this

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Difference between the prediction of the average sparsity at $\Delta = 200$ and 500 using the ST-Despali parameterization and the average sparsity from the N-body halo catalogs of the RayGalGroupSims (blue circles) and MultiDark-Planck2 (red triangles) at the redshift snapshots of the simulations.</figcaption></figure>

selection effect. As the RayGalGroupSims has a volume ~eight times larger than that of the MultiDark-Planck2 run, it allows us to precisely estimate the average sparsity at $z \gtrsim 1$, where the abundance of massive halos hosting galaxy groups and clusters drops rapidly. Henceforth, we can correct the systematic shift of the ST-Despali prediction shown in Figure 12 by introducing a third-order polynomial function of redshift that well approximate the numerical simulation points,

$$ \Delta \langle s_{200,500} \rangle_{\text{corr}}^{\text{fit}}(z) = s_0 + s_1 \cdot z + s_2 \cdot z^2 + s_3 \cdot z^3, \quad (A1) $$

with $s_0 = 0.14073268$, $s_1 = -0.37056373$, $s_2 = 0.32410084$, and $s_3 = -0.07044901$.

## Appendix B
Biased Sparsity Estimation from Shear Lensing Profile Fits

Weak-lensing masses are usually inferred by fitting the shear lensing measurements against a parametric shear profile. Several of these estimates in the literature fit the data with a one-parameter profile, such as the SIS profile or the NFW (Navarro et al. 1997), with a fixed concentration–mass relation. However, this results in biased estimates of the halo sparsity.

In the case of the SIS profile, this can be proven analytically, since assuming a radial density profile of the form $\rho(r) \propto 1/r^2$ gives a mass within a radius $r_\Delta$ enclosing an overdensity $\Delta$ that is $M_\Delta \propto 1/\sqrt{\Delta}$. Consequently, assuming an SIS profile results in a redshift and cosmology-independent sparsity, such that assuming $\Delta = 200$ and 500 gives $s_{200,500}^{\text{SIS}} = \sqrt{500/200}$, contrary to the N-body simulation results. In Figure 13, we plot the biased sparsity values for a cluster sample with masses $M_{200c}$ equal to the LC²-single clusters. The SIS estimates correspond to the streak of black points distributed along the straight line at $s_{200,500} \simeq 1.58$.

Similarly, sparsity estimates from mass measurements obtained assuming the NFW profile with a fixed concentration–mass relation, which leaves only one free parameter, cannot probe the full sparsity range either. Fixing the concentration–mass relation artificially reduces the scatter in the distribution of sparsity estimates compared to expectations

<footnote n="12">¹² http://xmm-heritage.oas.inaf.it/</footnote>

<footnote n="13">¹³ http://matplotlib.org/</footnote>

<footnote n="14">¹⁴ The MultiDark-Planck2 simulation consists of a (1 Gpc $h^{-1}$)³ volume with 3840³ N-body particles (corresponding to a mass resolution of $m_p = 1.51 \cdot 10^9 M_\odot h^{-1}$) of a flat ΛCDM model with $\Omega_m = 0.3071$, $\Omega_b = 0.048206$, $h=0.678$, $n_s = 0.96$, and $\sigma_8 = 0.823$ (Klypin et al. 2016); the RayGalGroupSims simulation has a (2.625 Gpc $h^{-1}$)³ volume with 4096³ N-body particles (corresponding to a mass resolution $m_p = 1.88 \cdot 10^{10} M_\odot h^{-1}$) of a flat ΛCDM model with $\Omega_m = 0.2573$, $\Omega_b = 0.048356$, $h = 0.72$, $n_s = 0.963$, and $\sigma_8 = 0.801$.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Sparsity of clusters from mass measurements inferred by fitting the SIS profile (black dots), the one-parameter fit NFW profile with the concentration–mass relation from Duffy et al. (2008; blue circles), and a two-parameter fit NFW profile (red triangles).</figcaption></figure>

from *N*-body simulations (see, e.g., Ragagnin et al. 2020, for a recent dedicated analysis). This can be seen in Figure 13, where the blue circles correspond to sparsity estimates from mass measurements obtained from a one-parameter fit NFW profile assuming the concentration–mass relation from Duffy et al. (2008). For comparison, the red triangles show the estimates from a two-parameter fit NFW profile. Hence, to avoid any bias in our analysis, we focus on lensing mass measurements inferred from two-parameter analyses. It is worth remarking that weak-lensing masses are measured assuming a cosmological framework; however, such a dependence cancels out when taking mass ratios (Sereno 2015).

## Appendix C
### Galaxy Cluster Sparsity Outliers

We have identified six clusters in the LC²-single sample that exhibit sparsity values that significantly differ from the distribution of the sparsity of clusters in the same redshift bin as shown in Figure 14. These are A2345, RXC J0528.9–3927, MS 1054.4–0321, XMMU J1229.4+0151, XLSS J022303.0–043622, and ISCS J1429.3+3437. As we can see, outliers have values of the sparsity $s_{200,500} > 2.2$. Such large values are a characteristic of unrelaxed/perturbed systems. The cluster A2345 is at $z=0.176$ and appears to be constituted of two merging subclusters (Bonafede et al. 2009; Boschin et al. 2010). Similarly, RXC J0528.9–3927 ($z=0.284$; Foëx et al. 2017) and MS 1054.4–0321 ($z=0.83$; Clowe et al. 2000) appear to be not fully relaxed clusters. The cluster XMMU J1229.4+0151 ($z=0.98$) may also be a perturbed system, since the lensing convergence map shows the presence of a second strong peak off-center with respect to the position of the peak of the X-ray emission and that of the cluster galaxies (Jee et al. 2011). Finally, XLSS J022303.0–043622 and ISCS J1429.3+3437 are among the most distant clusters detected to date at $z=1.22$ and 1.26, respectively. However, we found no information available in the literature on their structural and

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Distribution of the sparsity of the Selected LC²-single clusters in redshift bins of size $\Delta z = 0.2$ containing A2345 (top left panel), RXC J0528.9–3927 (top right panel), MS 1054.4–0321 and XMMU J1229.4+0151 (bottom left panel), and XLSS J022303.0–043622 and ISCS J1429.3+3437 (bottom right panel). The shaded area corresponds to the $1\sigma$ and $2\sigma$ scatter of the halo sparsities estimated for the different redshift bins from the *N*-body halo catalogs of the RayGalGroupSims simulation.</figcaption></figure>

dynamical properties that might account for their large values of sparsity. In order to compare these outliers with the expected distribution of halo sparsities from numerical simulations, we plot in Figure 14 the $1\sigma$ and $2\sigma$ scatter of the halo sparsity from catalogs of the RayGalGroupSims simulation presented in Corasaniti et al. (2018).
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Ensemble average cluster sparsity for the combined sample *Selected LC²-single* and HSC-XXL (red points) in redshift bins of size $\Delta z = 0.1$. The estimates given by $\langle M_{200c} \rangle/\langle M_{500c} \rangle$ (cyan stars) and $\langle 1/M_{500c} \rangle/\langle 1/M_{200c} \rangle$ (yellow squares) are consistent with the ensemble average values well within the uncertainties due to mass measurement errors.</figcaption></figure>

## Appendix D
Average Sparsity Consistency Relations

Corasaniti & Rasera (2019) showed that the validity of Equation (2) also implies the validity of the following relations:

$$ \langle s_{200,500} \rangle \equiv \left\langle \frac{M_{200c}}{M_{500c}} \right\rangle \approx \frac{\langle M_{200c} \rangle}{\langle M_{500c} \rangle} \approx \frac{\langle 1/M_{500c} \rangle}{\langle 1/M_{200c} \rangle}. \quad (D1) $$

These provide a set of consistency relations that can be used to test the consistency of sparsity measurements in galaxy cluster samples (Corasaniti & Rasera 2019). Here we compute the different average sparsity estimates for the combined *Selected LC²-single* + HSC-XXL sample in redshift bins of size $\Delta z = 0.1$. The results are shown in Figure 15. We can see that the various determinations are consistent with the ensemble average sparsity estimates within the uncertainties due to mass measurement errors.

## ORCID iDs

Pier-Stefano Corasaniti https://orcid.org/0000-0002-6386-7846

Mauro Sereno https://orcid.org/0000-0003-0302-0325

Stefano Ettori https://orcid.org/0000-0003-4117-8617

## References

* Abbott, T. M. C., Abdalla, F. B., Annis, J., et al. 2018, MNRAS, 480, 3879

* Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2020, PhRvD, 102, 3509

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

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

* Allen, S. W., Rapetti, D. A., Schmidt, R. W., et al. 2008, MNRAS, 383, 879

* Allen, S. W., Schmidt, R. W., Ebeling, H., Favian, A. C., & van Speybroeck, L. 2004, MNRAS, 353, 457

* Aubourg, E., Bailey, S., Bautista, J. E., et al. 2015, PhRvD, 92, 123516

* Balmès, I., Rasera, Y., Corasaniti, P.-S., & Alimi, J.-M. 2014, MNRAS, 437, 2328

* Barnes, D. J., Vogelsberger, M., Pearce, F. A., et al. 2020, arXiv:2001.11508

* Bartalucci, I., Arnaud, M., Pratt, G. W., & Le Brun, A. M. C. 2018, A&A, 617, A64

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

* Beutler, F., Blake, C., Colless, M., et al. 2011, MNRAS, 416, 3017

* Biffi, V., Borgani, S., Murante, G., et al. 2016, ApJ, 827, 112

* Birrer, S., Shajib, A. J., Galan, A., et al. 2020, A&A, 643, 165

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

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

* Böhringer, H., Chon, G., Retzlaff, J., et al. 2017, AJ, 153, 220

* Bonafede, A., Giovannini, G., Feretti, L., Govoni, F., & Murgia, M. 2009, A&A, 494, 429

* Bonvin, V., Courbin, F., Suyu, S. H., et al. 2017, MNRAS, 465, 491

* Boschin, W., Barrena, R., & Girardi, M. 2010, A&A, 521, 78

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

* Chen, K.-F., Oguri, M., Lin, Y.-T., & Miyazaki, S. 2020, ApJ, 891, 139

* CHEX-MATE Collaboration, Arnaud, M., Ettori, S., et al. 2020, arXiv:2010.11972

* Clowe, D., Luppino, G. A., Kaiser, N., & Gioia, I. M. 2000, ApJ, 539, 540

* Corasaniti, P. S. 2021, Halo_Sparsity, v2.0, Zenodo, doi:10.5281/zenodo.4558696

* Corasaniti, P. S., Ettori, S., Rasera, Y., et al. 2018, ApJ, 862, 40

* Corasaniti, P. S., Giocoli, C., & Baldi, M. 2020, PhRvD, 102, 043501

* Corasaniti, P. S., & Rasera, Y. 2019, MNRAS, 487, 4382

* Cui, W., Knebe, A., Yepes, G., et al. 2018, MNRAS, 480, 2898

* Cyburt, R. H., Fields, B. D., Olive, K. A., & Yeh, T.-H. 2016, RvMP, 410, 1911

* de Haan, T., Benson, B. A., Bleem, L.E., et al. 2016, ApJ, 832, 95

* Denzel, P., Coles, J. P., Saha, P., & Williams, L. L. R. 2021, MNRAS, 501, 784

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

* Dolag, K., Bertelmann, M., Perrotta, F., et al. 2004, A&A, 416, 853

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

* Ebeling, H., Edge, A. C., Mantz, A., et al. 2010, MNRAS, 407, 83

* Eckert, D., Ghirardini, V., Ettori, S., et al. 2019, A&A, 621, A40

* Eisenstein, D. J., & Hu, W. 1998, ApJ, 496, 605

* Ettori, S., De Grandi, S., & Molendi, S. 2002, A&A, 391, 841

* Ettori, S., Gastaldello, F., Leccardi, A., et al. 2010, A&A, 524, A68

* Ettori, S., Ghirardini, V., Ettori, S., et al. 2009, A&A, 501, 61

* Ettori, S., Tozzi, P., & Rosati, P. 2002, A&A, 398, 879

* Foëx, G., Chon, G., & Böhringer, H. 2017, A&A, 601, 145

* Freedman, W. L., Madore, B. F., Hatt, D., et al. 2019, ApJ, 882, 34

* Gelman, A., & Rubin, D. B. 1992, StaSc, 7, 457

* Ghirardini, V., Ettori, S., Amodeo, S., Capasso, R., & Sereno, M. 2017, A&A, 604, A100

* Giocoli, C., Tormen, G., & Sheth, R. K. 2012, MNRAS, 422, 185

* Henden, N. A., Puchwein, E., & Sijacki, D. 2020, MNRAS, 498, 2114

* Henning, J. W., Sayre, J. T., Reichardt, C. L., et al. 2018, ApJ, 852, 97

* Herbonnet, R., Sifón, C., Hoekstra, H., et al. 2020, MNRAS, 497, 4684

* Heymans, C., van Waerbeke, L., Miller, L., et al. 2012, MNRAS, 427, 146

* Hildebrandt, H., Choi, A., Blake, C., et al. 2016, MNRAS, 463, 635

* Hildebrandt, H., Köhlinger, F., van den Busch, J. L., et al. 2018, A&A, 633, A69

* Hildebrandt, H., Viola, M., Joudaki, S., et al. 2017, MNRAS, 465, 1454

* Hoekstra, H., Herbonnet, R., Muzzin, A., et al. 2015, MNRAS, 449, 685

* Hunter, J. D. 2007, CSE, 9, 90

* Jee, M. J., Dawson, K.S., Hoekstra, H., et al. 2011, ApJ, 737, 59

* Joudaki, S., Blake, C., Johnson, A., et al. 2018, MNRAS, 474, 4894

* Joudaki, S., Hildebrandt, H., Traykova, D., et al. 2020, A&A, 638, L1

* King, L. J., & Mead, J. M. G. 2011, MNRAS, 416, 2539

* Kitayama, T., & Suto, Y. 1996, MNRAS, 280, 638

* Klein, M., Israel, H., Nagarajan, A., et al. 2019, MNRAS, 488, 1704

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

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

* Lau, E. T., Nagai, D., & Nelson, K. 2013, ApJ, 777, 151

* Le Brun, A. M. C., McCarthy, I. G., Schaye, J., & Ponman, T. J. 2014, MNRAS, 441, 1270

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

* Lesci, G. F., Marulli, F., Moscardini, L., et al. 2020, arXiv:2012.12273

* Makiya, R., Hikage, C., & Komatsu, E. 2020, PASJ, 72, 26

* Marriage, T. A., Acquaviva, V., Ade, P. A. R., et al. 2011, ApJ, 737, 61

* Marulli, F., Veropalumbo, A., Sereno, M., et al. 2018, A&A, 620, A1
Maturi, M., Bellagamba, F., Radovich, M., et al. 2019, MNRAS, 485, 498

Mead, J. M. G., King, L. J., Sijacki, D., et al. 2010, MNRAS, 406, 434

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

Neveux, R., Burtin, E., de Mattia, A., et al. 2020, MNRAS, 499, 210

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

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

Planck Collaboration 2014, A&A, 571, A20

Planck Collaboration 2014, A&A, 571, A29

Planck Collaboration 2016, A&A, 594, A27

Planck Collaboration 2016, A&A, 594, A24

Planck Collaboration 2016, A&A, 594, A13

Planck Collaboration 2018, A&A, 641, A6

Postman, M., Coe, D., Benítez, N., et al. 2012, ApJS, 199, 25

Pratt, G., Arnaud, M., Biviano, A., et al. 2019, SSRv, 215, 25

Ragagnin, A., Saro, A., Singh, P., & Dolag, K. 2020, MNRAS, 500, 5056

Rasia, E., Meneghetti, M., Martino, R., et al. 2012, NJPh, 14, 055018

Riess, A. G., Casertano, S., Yuan, W., et al. 2018, ApJ, 861, 126

Riess, A. G., Macri, L., Hoffmann, S. L., et al. 2016, ApJ, 826, 56

Ross, A. J., Samushia, L., Howlett, C., et al. 2015, MNRAS, 449, 835

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

Salvati, L., Douspis, M., & Aghanim, N. 2018, A&A, 614, 13

Salvati, L., Douspis, M., Ritz, A., et al. 2019, A&A, 626, 27

Schellenberger, G., & Reiprich, T. H. 2017, MNRAS, 471, 1370

Schombert, J., McGaugh, S., & Lelli, F. 2020, AJ, 160, 71

Sereno, M. 2015, MNRAS, 450, 3665

Sereno, M., Covone, G., Izzo, L., et al. 2017, MNRAS, 472, 1946

Sereno, M., & Ettori, S. 2015, MNRAS, 450, 3633

Sereno, M., Giocoli, C., Ettori, S., & Moscardini, L. 2015, MNRAS, 449, 2024

Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119

Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709

Tröster, T., Sánchez, A. G., Asgari, M., et al. 2020, A&A, 633, L10

Troxel, M. A., MacCrann, N., Zuntz, J., et al. 2018, PhRvD, 98, 043528

Umetsu, K., Sereno, M., Lieu, M., et al. 2020, ApJ, 890, 148

White, S. D. M., Navarro, J. F., Evrard, A. E., & Frenk, C. S. 1993, Natur, 366, 429

Zhao, D. H., Jing, Y. P., Mo, H. J., & Börner, G. 2003, ApJ, 597, 9

Zhao, D. H., Jing, Y. P., Mo, H. J., & Börner, G. 2009, ApJ, 707, 354

Zubeldia, I., & Challinor, A. 2019, MNRAS, 489, 401