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

Pier-Stefano Corasaniti  LUTH, UMR 8102 CNRS, Observatoire de Paris, PSL
Research University, Université Paris Diderot, 5 Place Jules Janssen, 92190
Meudon, France  Sorbonne Université, CNRS, UMR 7095, Institut d’Astrophysique
de Paris, 98 bis bd Arago, 75014 Paris, France  Mauro Sereno  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, 40127 Bologna, Italy  Stefano Ettori  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, 40127
Bologna, Italy Pier-Stefano Corasaniti<EMAIL_ADDRESS>
###### 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 those 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 LC2-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 less extent the reduced
Hubble parameter $h$. Hence, in order to break the internal parameter
degeneracies we perform a combined likelihood analysis of cluster sparsities
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$.

clusters — cosmology: theory — cosmology: cosmological parameters

##  1 Introduction

There is a wide spread consensus that observations of galaxy clusters can
provide a wealth of cosmological information (see e.g. Allen, Evrard & Mantz,
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 datasets from
X-ray observations of the intra-cluster gas (see e.g. Ebeling et al., 2010;
Pierre et al., 2016; Böhringer et al., 2017), the detection of the Sunyaev-
Zeldovich (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, 2016a) 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, 2016b; Schellenberger & Reiprich, 2017;
Pacaud et al., 2018; Bocquet et al., 2019; DES Collaboration, 2020). Quite
remarkably, these studies have consistently found values of the amplitude of
linear matter density fluctuations on the $8$ 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, 2016c, 2018). As an example, the analyses of the SZ
cluster counts from Planck (Planck Collaboration, 2014, 2016b) 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, DES Collaboration, 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, 2016b; 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 which 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 expectation from hydrodynamical simulations of galaxy clusters
(see e.g. Rasia et al., 2012; Lau, Nagai & Nelson, 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, Hikage & Komatsu, 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, Douspis & Aghanim, 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 which host these
structures have assembled through a hierarchical process that depends on the
cosmic matter content, expansion rate, and the 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 formula (NFW,
Navarro, Frenk & White, 1997), such that for a halo of 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, Tormen & Sheth, 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. Firstly,
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). Secondly, selection effects may have a strong
impact on the cosmological parameter inference (Sereno, Giocoli, Ettori &
Moscardini, 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) have proposed that the sparsity, i.e. the
ratio of the halo mass within radii enclosing different overdensities, can
provide a non-parametric characterisation of the mass distribution in halos,
while retaining the cosmological information encoded in the halo mass profile.
More recently, Corasaniti et al. (2018) have shown that the halo sparsity
provide 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 sensitive probe of the screening mechanism in
Modified Gravity scenarios (Corasaniti, Giocoli & Baldi, 2020).

Here, we infer cosmological parameter constraints on $\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 (BAO) and
cluster gas mass fraction from archival data.

The paper is organised as the following. In Sec. 2 we describe the
methodology, in Sec. 3 the datasets used and in Sec. 4 we present the results.
In Sec. 5 we discuss the conclusions.

##  2 Methodolgy

###  2.1 Halo Sparsity

Originally introduced by Balmès et al. (2014), the halo sparsity provides a
non-parametric characterisation of the mass profile of halos in terms of the
ratio of masses within radii enclosing two different overdensities (in units
of the critical density or the background density), namely:

$s_{\Delta_{1},\Delta_{2}}=\frac{M_{\Delta_{1}}}{M_{\Delta_{2}}},$ (1)

where $M_{\Delta_{1}}$ and $M_{\Delta_{2}}$ are the halo masses at
overdensities $\Delta_{1}$ and $\Delta_{2}$ respectively, with
$\Delta_{1}<\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}}^{\rm min}}^{M_{\Delta_{2}}^{\rm
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}}^{\rm min}}^{\langle s_{\Delta_{1},\Delta_{2}}\rangle
M_{\Delta_{2}}^{\rm max}}\frac{dn}{dM_{\Delta_{1}}}d\ln{M_{\Delta_{1}}},$ (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 Eq. (2) has been
extensively tested using halo catalogs from high-resolution large volume
N-body simulations in Corasaniti et al. (2018); Corasaniti & Rasera (2019) and
we refer the 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 Eq. (2) holds valid to better than
percent level (see e.g. Table 1 in Corasaniti & Rasera, 2019). Hence, Eq. (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 Eq. (2). In the following, we assume a parametrisation of
the halo mass function as given by the Sheth-Tormen 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}}},$
(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),$
(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)],$ (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 cut-off of the halo mass function
at large masses, the upper limits in the integrals of Eq. (2) can be set to an
arbitrarily large value without affecting the evaluation. Thus, we have set
$M^{\rm max}_{500c}=10^{16}\,M_{\odot}\,h^{-1}$. In the case of the lower
limits, we have set $M^{\rm min}_{500c}=10^{13}\,M_{\odot}\,h^{-1}$
consistently with the lowest mass in our cluster catalogs. However, we have
found that the solution of Eq. (2) is insensitive to the specific choice of
$M^{\rm min}_{500c}$ 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 parametrisation of the
ST coefficients $a$, $p$ and $A$ given by Eq. (12) in Despali et al. (2016),
to which we refers 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 less extent on $h$.

The computation of the average halo sparsity has been performed using a
modified version of the numerical code Halo_Sparsity111The code Halo_Sparsity
is publicly available at https://github.com/pierste75/Halo_Sparsity.
(Corasaniti, 2021), which we have specifically developed for the work
presented here.

###  2.2 Sparsity Systematic Errors

####  2.2.1 Halo Mass Function Parametrisation

Using the N-body halo catalogs from the RayGalGroupSims simulation of a
$\Lambda$CDM model, Corasaniti et al. (2018) have shown that ST-Despali
parametrisation 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 $\sim 15\%$. In the case
of the Tinker et al. (2008) parametrisation of the halo mass function
deviations from the N-body simulation results are much larger and we do not
consider it here. Instead, a parametrisation of the ST coefficients on the
RayGalGroupSims halo catalogs (see Eq. (A4), (A5) in Corasaniti et al., 2018)
gives average sparsity estimates that reproduce the N-body simulation results
to sub-percent level at all redshifts. We will refer to this parametrisation
as ST-RayGal.

As discussed in Corasaniti et al. (2018), the systematic deviation from the
ST-Despali parametrisation is two-fold. On the one hand, the ST-Despali
parametrisation has been calibrated on halo catalogs from a suite of
cosmological simulations which covers a volume that is eight time 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 realised using halo catalogs which contain different halos at
different overdensities, rather than the same halos with masses estimated at
different overdensities, as expected from Eq. (2). Since modelling 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 of 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 unaltered the
external regions. 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 which differently from X-ray observations probe the outer regions
of the dark matter distribution in clusters. We estimate the level of bias
that such an effect induce 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
Fig. 2) for halos with $M_{200c}\geq 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. 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 $\sim 5\%$ at the low-mass end ($M_{200c}\sim 10^{14}$
$M_{\odot}$) to $\sim 3\%$ at 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, Frenk &
White, 1997) to tangential shear profile measurements. Deviations from
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
Table 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 these effects are smaller compared to the systematic error due to the
modelling of the mass function previously discussed, and much smaller than
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 order of $20-40\%$ per cent 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_{\rm gas}$, from X-ray observations of galaxy clusters can break the
internal parameter degeneracies. This is because $f_{\rm 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, Tozzi & Rosati, 2003; Allen et al., 2004; Ettori et al., 2009).

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

$f_{{\rm
gas},\Delta}=K_{\Delta}\,Y_{b,\Delta}\frac{\Omega_{b}}{\Omega_{m}}\left[\frac{D_{A}^{\rm
fid}(z)}{D_{A}(z)}\right]^{3/2}-f_{*,\Delta},$ (6)

where $D_{A}$ is the angular diameter distance, $D^{\rm fid}_{A}$ 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_{{\rm gas},\Delta}=M^{\rm
gas}_{\Delta}/M^{\rm tot}_{\Delta}$, where $M^{\rm tot}_{\Delta}=M^{\rm
HE}_{\Delta}/K_{\Delta}$ amd $M^{\rm HE}_{\Delta}$ is the cluster mass
obtained assuming hydrostatic equilibrium. Hence, differently from the halo
sparsity, $f_{{\rm gas},\Delta}$ also depends on uncertain quantities that
parametrise the baryonic content of clusters, such as the mass calibration
bias $K_{\Delta}\equiv 1-b_{\Delta}=M^{\rm HE}_{\Delta}/M^{\rm tot}_{\Delta}$,
the gas depletion factor $Y_{b,\Delta}$, and the stellar fraction
$f_{*,\Delta}$222Here, we assume that $f_{*,\Delta}$ is independent of the
cluster mass.. 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 which 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^{2}_{M}(z)\right]^{1/3},$ (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^{\prime}}{H(z^{\prime})}.$ (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}},$ (9)

with $\omega_{m,b}=\Omega_{m,b}h^{2}$.

##  3 Datasets

###  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 (LC2, 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 HSC-XXL sample.

The LC2-single cluster catalog consists of weak lensing cluster masses from
archival data; at the time of this analysis the latest version333We use the
version 3.8 publicly available at
http://pico.oabo.inaf.it/~sereno/CoMaLit/LC2/ 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 profile444For the
clusters analysed in Klein et al. (2019), we considered the masses from the
2-parameter fits reported among the others in LC2-all, rather than the
1-parameter fit results listed in LC2-single.. This is because mass estimates
from a one parameter fit, such a 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 dataset, we discard $6$ clusters which 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
LC2-single clusters containing $181$ entries. This includes mass estimates
from the PSZ2Lens catalog (Sereno et al., 2017), that is a statistically
complete and homogeneous sample of $35$ galaxy clusters at $z<0.7$ from the
second Planck Catalogue of Sunyaev-Zeldovich Sources (PSZ2Union, Planck
Collaboration, 2016a) with weak lensing data from the CFHTLenS (Canada France
Hawaii Telescope Lensing Survey, Heymans et al., 2012) and the RCSLenS (Red
Cluster Sequence Lensing Survey, Hildebrandt et al., 2016) surveys. 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.

In Fig. 1 we plot $M_{200c}$ vs $M_{500c}$ for the clusters in the Selected
LC2-single sample, the PSZ2Lens sub-sample, and the HSC-XXL dataset
respectively. Estimates of the correlation coefficient between the two mass
measurements are only available for the PSZ2Lens sub-sample, these are shown
in the inset plot. We see that, being the sparsity very weakly dependent on
selection effects, data from homogeneous samples such as PSZ2Lens or HSC-XXL
are in very good agreement with the heterogeneous Selected LC2-single
clusters. This appears more clearly in Fig. 2, where we have plotted the
average sparsity for the different samples in redshift bins of size
$\Delta{z}=0.1$ containing at least $2$ clusters per bin. We can see that the
different estimates are consistent with each other within 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 sample555It is worth remarking that
being a mass ratio, the error propagation of mass uncertainties on individual
sparsity measurements result in smaller errors the larger the correlation
between the mass estimates.. In the inset plot in Fig. 2 we show the average
sparsity in redshift bins of size $\Delta{z}=0.1$ for the combined sample
Selected LC2-single+HSC-XXL, which we have used in the cosmological analysis.
We have further tested the robustness of these estimates by testing the
validity of the average halo sparsity consistency relations (see Appendix D).

Figure 1: $M_{500c}$ vs $M_{200c}$ for the Selected LC2-single sample (red
circles) and the HSC-XXL clusters (empty triangles). The values of the
correlation coefficients of the PSZ2Lens cluster masses (yellow squares) are
shown in the inset plot.  Figure 2: Average cluster sparsity in redshift bins
of size $\Delta{z}=0.1$ for Selected LC2-single (red points), PSZ2Lens (yellow
squares), and HSC-XXL (empty triangles) clusters. The average sparsity
estimates for the combined sample Selected LC2-single and HSC-XXL (blue
circles) are shown in the inset plot, the number of clusters in each bin is
given below every data point. Figure 3: Averaged gas mass fraction estimates
in redshift bins of size $\Delta{z}=0.1$.

###  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 intra-cluster X-ray
emitting gas. Here, we focus on measurements of $f_{\rm gas}$ within
$R_{500c}$. Our dataset 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_{\odot}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 hydrostatic
equilibrium equation (Ettori, Grandi & Molendi, 2002). In addition, we have
used gas mass fraction estimates of 3 high-redshift X-ray clusters in the
redshift range $0.9\lesssim z\lesssim 1.1$ derived in Bartalucci et al.
(2018). The fiducial cosmology assumed in all these analyses is a flat
$\Lambda$CDM model with $\Omega_{m}=0.3$ and $h=0.7$. We bin these
measurements in redshift bins of size $\Delta{z}=0.1$ (see Fig. 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 dataset from
Ettori et al. (2010) exploited XMM data, whereas Ghirardini et al. (2017) is
based on Chandra observations. Gas mass measurements are mainly based on
surface brightness profiles, which are highly consistent among different X-ray
observatory, 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 1 we quote the
measurements we have used in our analysis. These consist of estimates from the
6dF Galaxy Survey (6dFGS, Beutler et al., 2011), the 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), 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.

Survey | $z$ | $D_{V}/r_{d}$
---|---|---
6dFGS | $0.106$ | $2.976\pm 0.133$
MGS | $0.15$ | $4.466\pm 0.168$
BOSS-DR 12 | $0.38$ | $9.994\pm 0.108$
BOSS-DR 12 | $0.51$ | $12.701\pm 0.128$
BOSS-DR 12 | $0.61$ | $14.481\pm 0.149$
eBOSS-QSO | $1.480$ | $26.51\pm 0.42$
Table 1: Spherically averaged cosmic distance measurements from the 6dF Galaxy
Survey (6dFGS, Beutler et al., 2011), the 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), SDSS-IV extended Baryon Acoustic Oscillation Spectroscopic Survey
Data Release 16 Quasar Sample (eBOSS-QSO, Neveux et al., 2020).

##  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 has been 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 which cover a wider range of redshifts, thus more
sensitive to the dependence on $D_{A}$, as well as the latest cosmic distance
estimates from BAO analyses.

###  4.1 Markov Chains & 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 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 $\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 results obtained under a 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_{\rm gas}$, we assume priors on the nuisance parameters
$K_{500c}$, $Y_{b,500c}$ and $f_{*,500c}$ within $R_{500c}$ (consistently 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
different prior on $K_{500c}$. In particular, we consider the following
Gaussian priors: $K^{\rm CLASH}_{500c}=0.78\pm 0.09$ consistent with the
estimate from the analysis of the CLASH sample (Sereno & Ettori, 2015);
$K^{\rm CCCP}_{500c}=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 (Hoekstra et al., 2015) (CCCP); $K^{\rm CMB}_{500c}=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 BAO in Salvati, Douspis & Aghanim (2018). We also consider the extreme
case with hard prior $K_{500c}=1$. As far as the baryon depletion factor is
concerned, we assume a Gaussian prior $Y^{\rm The300}_{b,500c}=0.938\pm 0.041$
consistently 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^{\rm
FABLE}_{b,500c}(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, Puchwein & Sijacki, 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 dataset (Eckert et al., 2019). Nevertheless, being a
subdominant term in Eq. (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.

####  4.2.1 Sparsity & BAO

In Fig. 4, we plot the $1$ 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 short-dash lines) and under
the HST-prior (blue thin solid lines) and the Planck-prior (green thin solid
lines) respectively. We also show the constraints derived from the combined
analysis of the average sparsity with the BAO (solid black lines). We can
clearly see that the constraints from the average sparsity data only are
highly degenerate along curves of approximate constant $S_{8}$ values (red
dotted lines).

In Table 2, we quote the mean and standard deviation of $S_{8}$ for the
Selected LC2-single clusters under different $h$-priors as well as the results
obtained for the PSZ2Lens subsample. In the former case find $S_{8}=0.75\pm
0.20$ for the uniform-$h$ prior, while assuming 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 $h$ associated to the HST and Planck priors
respectively. The average sparsity also carries information on the Hubble
parameter. This can be seen in Fig. 4, where the contours inferred from the
analysis of the sparsity data alone do shifts toward larger $S_{8}$ values for
the HST-$h$ and Planck-$h$ priors respectively. As it can be deduced from the
values quoted in Table 2, the analysis of the average sparsity from the
PSZ2Lens subsample gives results (bottom rows) for the different $h$-priors
which are largely consistent with those from the full cluster catalog, though
with slightly larger $1\sigma$ errors (of $\sim 10\%$ order) 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 only are consistent with those inferred from the sparsity analysis of
a sample of X-ray clusters presented in Corasaniti et al. (2018), where it was
found $S_{8}=0.87\pm 0.20$ for the HST-$h$ prior.

Figure 4: Marginalised $1$ and $2\sigma$ contours in the
$\Omega_{m}-\sigma_{8}$ from the analysis of the average sparsity data only
for the uniform-$h$ (dashed black lines), HST-$h$ (blue thin lines) and
Planck-$h$ (green thin lines) respectively, and in combination with BAO cosmic
distance measurements (solid black lines). The cross-point 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$) respectively. The red dotted
lines corresponds to constant $S_{8}$ value curves. For illustrative purposes
we show the marginalised 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 BBN prior, BAO data and CCCP prior on
the mass bias (Planck Collaboration, 2016b). Figure 5: Residuals of the
average sparsity (top panel) and BAO (bottom panel) data with respect to the
best-fit $\Lambda$CDM model with parameter values $\hat{\Omega}_{m}=0.274$,
$\hat{\sigma}_{8}=0.850$ and $\hat{h}=0.662$ respectively.

Including the BAO data666We 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$ allows
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 Fig. 4 we plot the $1$ and $2\sigma$ marginalised contours
from the Planck primary CMB cosmological analysis (red-yellow filled contours)
and the Planck-SZ cluster counts (blue-light blue filled contours). In Table
3, 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 marginalised 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 Fig. 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 $\sim 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-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 counts
baseline model777We 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, 2016b). with $\Omega_{m}=0.330\pm 0.030$,
$\sigma_{8}=0.760\pm 0.033$ and $S_{8}=0.796\pm 0.042$, that were derived
under similar assumptions (most notably the combination of the BBN prior and
BAO data).

The comparison indicates that Planck-SZ cluster counts provides tighter
constraints on $\sigma_{8}$, while the bounds on $\Omega_{m}$ have the same
level of statistical uncertainty. Limits on $h$ were derived in Planck
Collaboration (2016b), though the values of the marginal statistics were not
quoted. Here, we have analysed 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 dataset used in Planck Collaboration (2016b), as no constraints
would be inferred from the use of the cluster counts only (J. B. Melin,
private communication). This is not the case of the BAO data used here.

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

Figure 6: Marginalised $1$ and $2\sigma$ contours in the
$\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}^{\rm The300}$
and $K^{\rm CMB}_{500c}$ (dotted black lines), $Y_{b,500c}^{\rm The300}$ and
$K^{\rm CLASH}_{500c}$ (dash-dotted black lines), $Y_{b,500c}^{\rm FABLE}$ and
$K^{\rm CLASH}_{500c}$ (dash-dotted blue lines), $Y_{b,500c}^{\rm The300}$ and
$K^{\rm CCCP}_{500c}$ (dashed black lines), and $Y_{b,500c}^{\rm 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 4. As in Fig. 4,
we plot the contours from the Planck primary CMB and Planck-SZ cluster counts
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 corresponds to
the HST-$h$ and Planck-$h$ priors respectively.

####  4.2.2 Sparsity & Gas Mass Fraction

Cluster gas mass fraction data can be used to break the
$\Omega_{m}-\sigma_{8}$ degeneracy that characterises the constraints from the
average sparsity. However, differently from the BAO, the possibility to derive
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 because such
parameters appear as multiplicative factors in Eq. (6), thus resulting in the
$h^{3/2}$ dependence to be highly degenerate with $\Omega_{b}/\Omega_{m}$,
$K_{500c}$ and $Y_{b,500c}$.

In Fig. 6 we plot the $1$ 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 respectively. The marginal constraints on $\Omega_{m}$,
$\sigma_{8}$, $S_{8}$ and $h$ are summarised in Table 4. The red shaded area
in Fig. 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 4 (top
rows). The inset plots shows the marginal statistics of $h$ (left inset) and
$\Omega_{m}$ (right inset) as function of the $K_{500c}$ prior for the
$Y_{b,500c}^{\rm The300}$ case respectively. 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
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$, while the smaller the value of $\Omega_{m}$
necessary to fit the same gas mass fraction data. As a result, the gas mass
fraction data breaks 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 Fig. 6, where the contours shifts 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}^{\rm
FABLE}$ prior for the $K_{500c}^{\rm CLASH}$ case and find no statistically
significant differences with respect to the $Y_{b,500c}^{\rm The300}$ prior.
As it can be seen in Fig. 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 central and bottom rows in Table 4).
Nevertheless, the inconsistencies among the different gas mass fraction prior
model parameter assumptions remain unsolved. Most notably, in the case of the
$K_{500c}^{\rm CMB}$ prior, we find constraints on $\Omega_{m}$ and
$\sigma_{8}$ which lie at more than $\gtrsim 2\sigma$ from the Planck results.
This suggests that the low value of $K_{500c}^{\rm CMB}$ not only stands in
contrast with expectation from numerical simulations of clusters and direct
estimates of the mass calibration bias from galaxy cluster samples, but also
with gas mass fraction measurements we have considered here.

Given the large scatter in the prior values of the baryon depletion factor and
the 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 marginalise over large uniform priors. However, because of the
degeneracy with $h$, this may result into weak cosmological parameter
constraints. Hence, an external independent dataset 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.

####  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 marginalised 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 which include the Gaussian
priors previously discussed.

Figure 7: Marginalised 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. Figure 8: Marginalised $1$ and $2\sigma$ contours in the
$\Omega_{m}-\sigma_{8}$ from the combined analysis of the average cluster
sparsity, gas mass fraction and BAO data (dot black lines). As in Fig. 4, we
plot marginalised contours from the Planck primary CMB analysis (yellow and
red contours) and the Planck-SZ number counts (dark and light steal-blue
contours). The cross-point 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$) respectively.

In Fig. 7 we show the 1D marginalised posteriors, while in Fig. 8, we plot the
marginalised $1$ and $2\sigma$ contours in the $\Omega_{m}$-$\sigma_{8}$
plane. As in Fig. 4, we plot the contours from the Planck primary CMB and the
Planck-SZ cluster counts cosmological analyses respectively. In Table 5, we
quote the results of the marginal statistics of $\Omega_{m}$, $\sigma_{8}$,
$S_{8}$ and $h$.

Figure 9: Residuals of the average sparsity (top panel), gas mass fraction
(central 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.

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
$\hat{h}=0.690$. The data residuals with respect to the best-fit model are
shown in Fig. 9. Notice that the 1D marginalised 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 of 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 to the peak of the marginalised 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 $2$. The constraints on $\sigma_{8}$
and $h$ are improved with a gain on the $1\sigma$ errors at the $\sim 49\%$
and $\sim 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 Fig. 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}^{\rm CLASH}$
($K_{500c}^{\rm CCCP}$) does not noticeably improve the constraints on the
cosmological parameters. In contrast, it allows to infer a lower limit on
$Y_{b,500c}\gtrsim 0.89$ ($Y_{b,500c}\gtrsim 0.91$) at $1\sigma$.

In Fig. 10 and Fig. 11 we show summary plots of state-of-art estimates of
$S_{8}$ (see also Pratt et al., 2019) and $h$ from various probes.

Figure 10: Constraints on $S_{8}$ from the analysis of shear lensing data
(blue triangles) of KV450 (Hildebrandt et al., 2018), DES-Y1 (Troxel et al.,
2018), KV450+DES-Y1 (Joudaki et al., 2020), KV450+BOSS (Tröster et al., 2019);
cluster number counts (red squares) of Planck-SZ (Planck Collaboration,
2016b), HIFLUGCS (Schellenberger & Reiprich, 2017), SPT-SZ (Bocquet et al.,
2019), XXL (Pacaud et al., 2018) and DES-Y1 (DES Collaboration, 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$ and $2\sigma$ limits from the Planck primary CMB analysis
(TT,TE,EE+lowE+lensing, Planck Collaboration, 2018). 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), 8TDGL (Denzel et a., 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$ and $2\sigma$ limits from the Planck primary CMB
analysis (TT,TE,EE+lowE+lensing, Planck Collaboration, 2018).

##  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 source 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}$ which are lower than those
obtained from the analysis of the Planck anisotropy power spectra (see e.g.
Planck Collaboration, 2016b; Bocquet et al., 2019; DES Collaboration, 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
the effect 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) have shown 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 LC2-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 combination with gas mass
fraction and BAO data. We find the combination of these dataset to provide
competitive constraints on $\Omega_{m}$, $\sigma_{8}$ and $h$ which are
consistent with those from the Planck primary CMB and the 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 selections 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
($\sim 100$) of objects homogeneously analysed in their X-ray and lensing
signal out to $R_{500}$ and beyond, as the one that will be available from the
CHEX-MATE project888http://xmm-heritage.oas.inaf.it/ (CHEX-MATE Collaboration,
2020).

PSC 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.
SE and MS acknowledge financial contribution from contract ASI-INAF
n.2017-14-H.0 and INAF ‘Call per interventi aggiuntivi a sostegno della
ricerca di main stream di INAF’ n. 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 n. 279954. Data visualisations are prepared with the
matplotlib999http://matplotlib.org/ library (Hunter, 2007).

## References

  * Abbott et al. (2018) Abbott, T. M. C. et al., 2018, MNRAS, 480, 3879, doi: 10.1093/mnras/sty1939
  * Alam et al. (2017) Alam, S., et al., 2017, MNRAS, 470, 2617, doi: 10.1093/mnras/stx721
  * Allen et al. (2004) Allen, S. W., et al., MNRAS, 353, 457, doi: 10.1111/j.1365-2966.2004.08080.x
  * Allen et al. (2008) Allen, S. W., et al., 2008, MNRAS, 383, 879, doi: 10.1111/j.1365-2966.2007.12610.x
  * Allen, Evrard & Mantz (2011) Allen, S. W., Evrard, A. E. and Mantz, A. B., 2011, ARA&A, 49, 409, doi: 10.1146/annurev-astro-081710-102514
  * Aubourg et al. (2015) Aubourg, E., et al., 2015, Phys. Rev. D, 92,123516, doi: 10.1103/PhysRevD.92.123516
  * Balmès et al. (2014) Balmès, I., Rasera, Y., Corasaniti, P.-S., Alimi, J.-M., 2014, MNRAS, 437, 2328, doi: 10.1093/mnras/stt2050
  * Barnes et al. (2020) Barnes, D. J., et al., 2020, arXiv:2001.11508
  * Bartalucci et al. (2018) Bartalucci, I., Arnaud, M., Pratt, G. W., Le Brun, A. M. C., 2018, A&A, 617, A64, doi: 10.1051/0004-6361/201732458
  * Becker & Kravtsov (2011) Becker, M. R., Kravtsov, A. V., 2011, ApJ, 740, 25, doi: 10.1088/0004-637X/740/1/25
  * Beutler et al. (2011) Beutler, F., et al., MNRAS, 416, 3017, doi: 10.1111/j.1365-2966.2011.19250.x
  * Biffi et al. (2016) Biffi, V., et al., ApJ, 827, 112, doi: 10.3847/0004-637X/827/2/112
  * Birrer et al. (2020) Birrer, S., et al., 2020, arXiv:2007.02941
  * Bleem et al. (2015) Bleem, L. E., et al., ApJS, 216, 27, doi: 10.1088/0067-0049/216/2/27
  * Bocquet et al. (2019) Bocquet, S., et al., ApJ, 878, 55, doi: 10.3847/1538-4357/ab1f10
  * Böhringer et al. (2017) Böhringer, H., et al., 2017, AJ, 153, 220, doi: 10.3847/1538-3881/aa67ed
  * Bonafede et al. (2009) Bonafede, A., Giovannini, G., Feretti, L., Govoni, F., Murgia, M., 2009, A&A, 494, 429, doi: 10.1051/0004-6361:200810588
  * Bonvin et al. (2017) Bonvin, V., et al., 2017, MNRAS, 465, 491, doi: 10.1093/mnras/stw3006
  * Boschin, Barrena & Girardi (2010) Boschin, W., Barrena, R., Girardi, M., 2010, A&A, 521, 78, doi: 10.1051/0004-6361/201014839
  * Bullock et al. (2001) Bullock, J. S., et al., 2001, MNRAS, 321, 559, doi: 10.1046/j.1365-8711.2001.04068.x
  * Chen et al. (2020) Chen, K.-F., Oguri, M., Lin, Y.-T., Miyazaki, S., 2020, ApJ, 891, 139, doi: 10.3847/1538-4357/ab74d3
  * CHEX-MATE Collaboration (2020) CHEX-MATE Collaboration:, Arnaud, M., Ettori, S., et al., 2020, A&A, in press, arXiv:2010.11972
  * Clowe et al. (2000) Clowe, D., Luppino, G. A., Kaiser, N., Gioia, I. M., 2000, ApJ, 539, 540, doi: 10.1086/309242
  * Corasaniti et al. (2018) Corasaniti, P. S., et al., 2018, ApJ, 862, 40, doi: 10.3847/1538-4357/aaccdf
  * Corasaniti & Rasera (2019) Corasaniti, P. S., Rasera, Y., 2019, MNRAS, 487, 4382, doi: 10.1093/mnras/stz1579
  * Corasaniti, Giocoli & Baldi (2020) Corasaniti, P.S., Giocoli, C., Baldi, M., 2020, Phys. Rev. D, 102, 043501, doi: 10.1103/PhysRevD.102.043501
  * Corasaniti (2021) Corasaniti, P. S., 2021, Halo_Sparsity, v2.0, Zenodo, doi: 10.5281/zenodo.4558696
  * Cyburt et al. (2016) Cyburt, R. H., Fields, B. D., Olive, K. A., Yeh, T.-H., 2016, Rev. Mod. Phys., 410, 1911 doi: 10.1103/RevModPhys.88.015004
  * Cui et al. (2018) Cui, W., et al., 2018, MNRAS, 480, 2898, doi: 10.1093/mnras/sty2111
  * Denzel et a. (2020) Denzel, P., Coles, J. P., Saha, P., Williams, L. L. R., 2020, arXiv:2007.14398
  * DES Collaboration (2020) DES Collaboration, 2020, Phys. Rev. D, 102, 3509, doi: 10.1103/PhysRevD.102.023509
  * Despali et al. (2016) Despali, G., et al., 2016, MNRAS, 456, 2486, doi: 10.1093/mnras/stv2842
  * Dolag et al. (2004) Dolag, K., et al., 2004, A&A, 416, 853, doi: 10.1051/0004-6361:20031757
  * Duffy et al. (2008) Duffy, A. R., Schaye, J., Kay, S. T., Dalla Vecchia, C., 2008, MNRAS, 390, L64, doi: 10.1111/j.1745-3933.2008.00537.x10.1111/j.1745-3933.2008.00537.x
  * Ebeling et al. (2010) Ebeling, H., et al., 2010, MNRAS, 407, 83, doi: 10.1111/j.1365-2966.2010.16920.x
  * Eckert et al. (2019) Eckert, D., et al., 2019, A&A, 621, A40, doi: 10.1051/0004-6361/201833324
  * Eisenstein & Hu (1998) Eisenstein, D.J., Hu, W., 1998, ApJ, 496, 605, doi: 10.1086/305424
  * Ettori, Grandi & Molendi (2002) Ettori, S., De Grandi, S., Molendi, S., 2002, A&A, 391, 841, doi: 10.1051/0004-6361:20020905
  * Ettori, Tozzi & Rosati (2003) Ettori, S., Tozzi, P., Rosati, P., 2002, A&A, 398, 879, doi: 10.1051/0004-6361:20021706
  * Ettori et al. (2009) Ettori, S., et al., 2009, A&A, 501, 61, doi: 10.1051/0004-6361/200810878
  * Ettori et al. (2010) Ettori, S., et al., 2010, A&A, 524, A68, doi: 10.1051/0004-6361/201015271
  * Foëx, Chon & Böhringer (2017) Foëx, G. Chon, G., Böhringer, H., 2017, A&A, 601, 145, doi: 10.1051/0004-6361/201630086
  * Freedman et al. (2019) Freedman, W. L., et al., 2019, ApJ, 882, 34, doi: 10.3847/1538-4357/ab2f73
  * Gelman & Rubin (1992) Gelman, A., Rubin, D. B., 1992, Statist. Sci., 7, 457, doi: 10.1214/ss/1177011136
  * Ghirardini et al. (2017) Ghirardini, V., et al., 2017, A&A, 604, A100, doi: 10.1051/0004-6361/201630209
  * Giocoli, Tormen & Sheth (2012) Giocoli, C., Tormen, G., Sheth, R. K., 2012, MNRAS, 422, 185, doi: 10.1111/j.1365-2966.2012.20594.x
  * de Haan et al. (2016) de Haan, T., et al., 2016, ApJ, 832, 95, doi: 10.3847/0004-637X/832/1/95
  * Heymans et al. (2012) Heymans, C., et al., 2012, MNRAS, 427, 146, doi: 10.1111/j.1365-2966.2012.21952.x
  * Henden, Puchwein & Sijacki (2020) Henden, N. A., Puchwein, E., Sijacki, D., 2020, MNRAS, 498, 2114, doi: 10.1093/mnras/staa223
  * Henning et al. (2018) Henning, J. W., et al., 2018, ApJ, 852, 97, doi: 10.3847/1538-4357/aa9ff4
  * Herbonnet et al. (2020) Herbonnet, R. et al., 2020, MNRAS, 497, 4684, doi: 10.1093/mnras/staa2303
  * Hildebrandt et al. (2016) Hildebrandt, H., et al., 2016, MNRAS, 463, 635, doi: 10.1093/mnras/stw2013
  * Hildebrandt et al. (2017) Hildebrandt, H., et al., 2017, MNRAS, 465, 1454, doi: 10.1093/mnras/stw2805
  * Hildebrandt et al. (2018) Hildebrandt, H., et al., 2018, A&A, 633, A69.
  * Hoekstra et al. (2015) Hoekstra, H., et al., 2015, MNRAS, 449, 685, doi: 10.1093/mnras/stv275
  * Hunter (2007) Hunter, J. D., 2007, Comput. Sci. Eng., 9, 90, doi: 10.1109/MCSE.2007.55
  * Jee et al. (2011) Jee, M. J., et al., 2011, ApJ, 737, 59, doi: 10.1088/0004-637X/737/2/59
  * Joudaki et al. (2018) Joudaki, S., et al., 2018, MNRAS, 474, 4894, doi: 10.1093/mnras/stx2820
  * Joudaki et al. (2020) Joudaki, S., et al., 2019, A&A, 638, L1, doi: 10.1051/0004-6361/201936154
  * King & Mead (2011) King, L. J.,Mead, J. M. G., 2011, MNRAS, 416, 2539, doi: 10.1111/j.1365-2966.2011.19009.x
  * Kitayama & Suto (1996) Kitayama, T., Suto, Y., 1996, MNRAS, 280, 638, doi: 10.1093/mnras/280.3.638
  * Klein et al. (2019) Klein, M. et al., 2019, MNRAS, 488, 1704, doi: 10.1093/mnras/stz1491
  * Klypin et al. (2016) Klypin, A., Yepes, G., Gottlöber, S., Prada, F., Heß, S., 2016, MNRAS, 457, 4340, doi: 10.1093/mnras/stw248
  * Kravtsov & Borgani (2012) Kravtsov A. V., Borgani S., 2012, ARA&A, 50, 353, doi: 10.1146/annurev-astro-081811-125502
  * Lau, Nagai & Nelson (2013) Lau, E. T., Nagai, D., Nelson, K., 2013, ApJ, 777, 151, doi: 10.1088/0004-637X/777/2/151
  * Le Brun et al. (2014) Le Brun, A. M. C., McCarthy, I. G., Schaye, J., Ponman, T. J, MNRAS, 441, 1270, doi: 10.1093/mnras/stu608
  * Lee et al. (2018) Lee, B.E., et al., 2018, MNRAS, 479, 890, doi: 10.1093/mnras/sty1377
  * Makiya, Hikage & Komatsu (2020) Makiya, R., Hikage, C., Komatsu, E., 2020, PASJ, 72, 26, doi: 10.1093/pasj/psz147
  * Marriage et al. (2011) Marriage, T. A., et al., 2011, ApJ, 737, 61, doi: 10.1088/0004-637X/737/2/61
  * Marulli et al. (2018) Marulli, F., et al., 2018, A&A, 620, A1, doi: 10.1051/0004-6361/201833238
  * Maturi et al. (2019) Maturi, M., et al., 2019, MNRAS, 485, 498, doi: 10.1093/mnras/stz294
  * Mead et al. (2010) Mead, J. M. G., et al., 2010, MNRAS, 406, 434, doi: 10.1111/j.1365-2966.2010.16674.x
  * Navarro, Frenk & White (1997) Navarro, J. F., Frenk, C. S., White, S. D. M., 1997, ApJ, 490, 493, doi: 10.1086/304888
  * Neveux et al. (2020) Neveux, R., et al., 2020, MNRAS, 499, 210, doi: 10.1093/mnras/staa2780
  * Pacaud et al. (2018) Pacaud, F., et al., 2018, A&A, 620, A10, doi: 10.1051/0004-6361/201834022
  * Pierre et al. (2016) Pierre, M., et al., 2016, A&A, 592, 1, doi: 10.1051/0004-6361/201526766
  * Planck Collaboration (2014) Planck Collaboration, 2014, A&A, 571, A20, doi: 10.1051/0004-6361/201321521
  * Planck Collaboration (2014) Planck Collaboration, 2014, A&A, 571, A29, doi: 10.1051/0004-6361/201321523
  * Planck Collaboration (2016a) Planck Collaboration, 2016, A&A, 594, A27, doi: 10.1051/0004-6361/201525823
  * Planck Collaboration (2016b) Planck Collaboration, 2016, A&A, 594, A24, doi: 10.1051/0004-6361/201525833
  * Planck Collaboration (2016c) Planck Collaboration, 2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Planck Collaboration (2018) Planck Collaboration, 2018, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Postman et al. (2012) Postman, M., et al., 2012, ApJS, 199, 25, doi: 10.1088/0067-0049/199/2/25
  * Pratt et al. (2019) Pratt, G., et al., 2019, Space Sci. Rev., 215, 25, doi: 10.1007/s11214-019-0591-0
  * Ragagnin et al. (2021) Ragagnin, A., Saro, A., Singh, P., Dolag, K., 2020, MNRAS, 500, 5056, doi: 10.1093/mnras/staa3523
  * Rasia et al. (2012) Rasia, E., et al., 2012, New Journal of Physics, 14, 055018, doi: 10.1088/1367-2630/14/5/055018
  * Riess et al. (2016) Riess, A. G., et al., 2016, ApJ, 826, 56, doi: 10.3847/0004-637X/826/1/56
  * Riess et al. (2018) Riess, A. G., et al., 2018, ApJ, 861, 126, doi: 10.3847/1538-4357/aac82e
  * Ross et al. (2015) Ross, A. J.,Samushia, L., Howlett, C., Percival, W. J., Burden, A., Manera, M., 2015, MNRAS, 449, 835, doi: 10.1093/mnras/stv154
  * Rykoff et al. (2016) Rykoff, E. S., et al., 2016, ApJS, 224, 1, doi: 10.3847/0067-0049/224/1/1
  * Salvati, Douspis & Aghanim (2018) Salvati, L., Douspis, M., Aghanim, N., 2018, A&A, 614, 13, doi: 10.1051/0004-6361/201731990
  * Salvati et al. (2019) Salvati, L., et al., 2019, A&A, 626, 27, doi: 10.1051/0004-6361/201935041
  * Sereno, Giocoli, Ettori & Moscardini (2015) Sereno, M., Giocoli, C., Ettori, S., Moscardini, L., 2015, MNRAS, 449, 2024, doi: 10.1093/mnras/stv416
  * Sereno (2015) Sereno, M., 2015, MNRAS, 450, 3665, doi: 10.1093/mnras/stu2505
  * Sereno & Ettori (2015) Sereno, M., Ettori, S., 2015, MNRAS, 450, 3633, doi: 10.1093/mnras/stv810
  * Sereno et al. (2017) Sereno, M., et al., 2017, MNRAS, 472, 1946, doi: 10.1093/mnras/stx2085
  * Schellenberger & Reiprich (2017) Schellenberger, G., Reiprich, T. H., 2017, MNRAS, 471, 1370, doi: 10.1093/mnras/stx1583
  * Schombert et al. (2020) Schombert, J., McGaugh, S., Lelli, F., 2020, AJ, 160, 71, doi: 10.3847/1538-3881/ab9d88
  * Sheth & Tormen (1999) Sheth, R. K., Tormen, G., 1999, MNRAS, 308, 119, doi: 10.1046/j.1365-8711.1999.02692.x
  * Tinker et al. (2008) Tinker, J., et al., 2008, ApJ, 688, 709, doi: 10.1086/591439
  * Tröster et al. (2019) Tröster, T., et al., 2019, A&A, 633, L10, doi: 10.1051/0004-6361/201936772
  * Troxel et al. (2018) Troxel, M. A., et al., 2018, Phys. Rev. D, 98, 043528, doi: 10.1103/PhysRevD.98.043528
  * Umetsu et al. (2020) Umetsu, K., et al., 2020, ApJ, 890, 148, doi: 10.3847/1538-4357/ab6bca
  * White et al. (1993) White, S. D. M., Navarro, J. F., Evrard, A. E., Frenk, C.S., 1993, Nature, 366, 429, doi: 10.1038/366429a0
  * Zhao et al. (2003) Zhao, D. H, et al., 2003, ApJ, 597, 9, doi: 10.1086/379734
  * Zhao et al. (2009) Zhao, D. H, et al., 2009, ApJ, 707, 354, doi: 10.1088/0004-637X/707/1/354
  * Zubeldia & Challinor (2019) Zubeldia, I., Challinor, A., 2019, MNRAS, 489, 401, doi: 10.1093/mnras/stz2153

##  Appendix A Average Sparsity Redshift Model Correction

The validity of Eq. (2) in predicting the ensemble average halo sparsity has
been extensively tested using N-body halo catalogs from the RayGalGroupSims
simulation of a $\Lambda$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 then a few percent accuracy level.
On the other hand, as already mentioned in Section 2, assuming the ST-Despali
parametrisation, Corasaniti et al. (2018) found that the predicted average
sparsity is consistent with the RayGalGroupSims results within a few percent
level only in the redshift range $0.5\lesssim z\lesssim 1.2$, while
differences increase up to $\sim 10\%$ at lower and higher redshifts. Such
discrepancy is essentially due to the fact that the halo catalogs used in the
calibration of the ST-Despali parametrisation contain different halos at
different overdensities, rather than identical halos with masses estimated at
different overdensities as expected from Eq. (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 Fig. 12,
where we plot the difference between the value of the average sparsity from
the ST-Despali parametrisation and that measured from the N-body halos at the
redshift snapshots of the RayGalGroupSims (blue dots) and the MultiDark-
Planck2 (red triangles) simulations respectively. 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 and do not posses similar
characteristics101010The MultiDark-Planck2 simulation consists of a $(1\,{\rm
Gpc}\,h^{-1})^{3}$ volume with $3840^{3}$ N-body particles (corresponding to a
mass resolution of $m_{p}=1.51\cdot 10^{9}\,M_{\odot}h^{-1}$) of a flat
$\Lambda$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\,{\rm Gpc}\,h^{-1})^{3}$ volume with $4096^{3}$
N-body particles (corresponding to a mass resolution $m_{p}=1.88\cdot
10^{10}\,M_{\odot}h^{-1}$) of a flat $\Lambda$CDM model with
$\Omega_{m}=0.2573$, $\Omega_{b}=0.048356$, $h=0.72$, $n_{s}=0.963$ and
$\sigma_{8}=0.801$., thus highlighting the algorithmic nature of this
selection effect. As the RayGalGroupSims has a volume $\sim 8$ 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 Fig. 12, by introducing
a third-order polynomial function of redshift which well approximate the
numerical simulation points:

$\Delta\langle s_{200,500}\rangle_{\rm corr}^{fit}(z)=s_{0}+s_{1}\cdot
z+s_{2}\cdot z^{2}+s_{3}\cdot z^{3},$ (A1)

with $s_{0}=0.14073268$, $s_{1}=-0.37056373$, $s_{2}=0.32410084$ and
$s_{3}=-0.07044901$.

Figure 12: Difference between the prediction of the average sparsity at
$\Delta=200$ and $500$ using the ST-Despali parametrization and the average
sparsity from the N-body halo catalogs of the RayGalGroupSims (blue dots) and
MultiDark-Planck2 (red triangles) at the redshift snapshots of the
simulations.

##  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 estimate in
the literature fit the data with a 1-parameter profile such as the SIS profile
or the NFW (Navarro, Frenk & White, 1997) with 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 a SIS profile
results in a redshift and cosmology independent sparsity such that assuming
$\Delta=200$ and $500$ gives $s^{\rm SIS}_{200,500}=\sqrt{500/200}$ contrary
to the N-body simulation results. In Fig. 13, we plot the biased sparsity
values for a cluster sample with masses $M_{200c}$ equal to the LC2-single
clusters. The SIS estimates correspond to the strake 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 from N-body
simulations (see e.g. Ragagnin et al., 2021, for a recent dedicated analysis).
This can be seen in Fig. 13 where the blue dots corresponds 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
profiles. 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 dependence cancels out when taking mass ratios (Sereno, 2015).

Figure 13: Sparsity of clusters from mass measurements inferred fitting the
SIS profile (black dots), the one-parameter fit NFW-profile with the
concentration-mass relation from Duffy et al. (2008) (blue dots) and a two-
parameter fit NFW-profile (red triangles).

##  Appendix C Galaxy Cluster Sparsity Outliers

We have identified $6$ clusters in the LC2-single sample which exhibit
sparsity values that significantly differ from the distribution of the
sparsity of clusters in the same redshift bin as shown in Fig. 14. These are
Abell 2345, 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. Abell 2345 is a cluster at $z=0.176$, which
appears to be constituted by two merging sub-clusters (Bonafede et al., 2009;
Boschin, Barrena & Girardi, 2010). Similarly, RXC J0528.9-3927 ($z=0.284$)
(Foëx, Chon & Böhringer, 2017) and MS 1054.4-0321 ($z=0.83$) (Clowe et al.,
2000) appear to be not fully relaxed clusters. 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-centred 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 dynamical properties that might account for their large values
of the sparsity. In order to compare these outliers with the expected
distribution of halo sparsities from numerical simulations, we plot in Fig. 14
the $1$ and $2\sigma$ scatter of the halo sparsity from catalogs of the
RayGalGroupSims simulation presented in Corasaniti et al. (2018).

Figure 14: Distribution of the sparsity of the Selected LC2-single clusters in
redshift bins of size $\Delta{z}=0.2$ containing Abell 2345 (top left panel),
RXC J0528.9-3927 (top right panel), MS 1054.4-0321 and XMMU J1229.4+0151
(bottom left panel), XLSS J022303.0-043622 and ISCS J1429.3+3437 (bottom right
panel). The shaded area corresponds to the $1$ and $2\sigma$ scatter of the
halo sparsities estimated for the different redshift bins from the N-body halo
catalogs of the RayGalGroupSims simulation. Figure 15: Ensemble average
cluster sparsity for the combined sample Selected LC2-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.

##  Appendix D Average Sparsity Consistency Relations

Corasaniti & Rasera (2019) have shown that the validity of Eq. (2) also
implies the validity of the following relations:

$\langle
s_{200,500}\rangle\equiv\bigg{\langle}\frac{M_{200c}}{M_{500c}}\bigg{\rangle}\approx\frac{\langle
M_{200c}\rangle}{\langle M_{500c}\rangle}\approx\frac{\langle
1/M_{500c}\rangle}{\langle 1/M_{200c}\rangle}.$ (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 LC2-single+HSC-XXL sample in redshift bins of size
$\Delta{z}=0.1$. The results are shown in Fig. 15. We can see that the various
determinations are consistent with the ensemble average sparsity estimates
within the uncertainties due to mass measurement errors.

##  Appendix E MCMC Likelihood Analysis Marginal Statistics

We quote the mean and standard deviation of the different cosmological model
parameters inferred from the MCMC likelihood analysis of the average sparsity
data only in Table 2, the average sparsity in combination with BAO in Table 3
and in combination with the gas mass fraction in Table 4, while the results of
the joint analysis are given in Table 5.

Sparsity Only
Cluster Sample & $h$-prior $S_{8}$ $Selected$ LC${}^{2}-Single$ \+ Uniform
$0.75\pm 0.20$ $Selected$ LC${}^{2}-Single$ \+ HST $0.80\pm 0.18$ $Selected$
LC${}^{2}-Single$ \+ Planck $0.82\pm 0.16$ PSZ2Lens + Uniform $0.69\pm 0.22$
PSZ2Lens + HST $0.73\pm 0.21$ PSZ2Lens + Planck $0.75\pm 0.19$

Table 2: 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).

Sparsity + BAO
$h$-prior $\Omega_{m}$ $\sigma_{8}$ $S_{8}$ $h$ Uniform $0.277\pm 0.029$
$0.856\pm 0.100$ $0.818\pm 0.070$ $0.662\pm 0.023$ HST $0.318\pm 0.029$
$0.756\pm 0.089$ $0.775\pm 0.071$ - Planck $0.285\pm 0.023$ $0.833\pm 0.088$
$0.809\pm 0.069$ -

Table 3: Mean and standard deviation of $\Omega_{m}$, $\sigma_{8}$, $S_{8}$
and $h$ from the joint analysis of the average sparsity and BAO data.

Sparsity + Gas Mass Fraction
$f_{\rm gas}$ and $h$ priors $\Omega_{m}$ $\sigma_{8}$ $S_{8}$ $h$
$Y_{b,500c}^{\rm The300}$ \+ $K^{\rm CMB}_{500c}$ \+ Uniform $0.13\pm 0.01$
$1.35\pm 0.18$ $0.90\pm 0.09$ $1.10\pm 0.05$ $Y_{b,500c}^{\rm The300}$ \+
$K^{\rm CLASH}_{500c}$ \+ Uniform $0.19\pm 0.03$ $1.04\pm 0.14$ $0.82\pm 0.08$
$0.91\pm 0.06$ $Y_{b,500c}^{\rm FABLE}$ \+ $K^{\rm CLASH}_{500c}$ \+ Uniform
$0.19\pm 0.03$ $1.06\pm 0.15$ $0.83\pm 0.08$ $0.92\pm 0.07$ $Y_{b,500c}^{\rm
The300}$ \+ $K^{\rm CCCP}_{500c}$ \+ Uniform $0.20\pm 0.02$ $1.03\pm 0.12$
$0.83\pm 0.07$ $0.89\pm 0.05$ $Y_{b,500c}^{\rm The300}$ $(K_{500c}=1)$ \+
Uniform $0.25\pm 0.03$ $0.86\pm 0.10$ $0.79\pm 0.07$ $0.79\pm 0.05$
$Y_{b,500c}^{\rm The300}$ \+ $K^{\rm CMB}_{500c}$ \+ HST $0.16\pm 0.01$
$1.27\pm 0.04$ $0.92\pm 0.03$ - $Y_{b,500c}^{\rm The300}$ \+ $K^{\rm
CLASH}_{500c}$ \+ HST $0.26\pm 0.03$ $0.85\pm 0.10$ $0.79\pm 0.07$ -
$Y_{b,500c}^{\rm FABLE}$ \+ $K^{\rm CLASH}_{500c}$ \+ HST $0.26\pm 0.02$
$0.87\pm 0.07$ $0.81\pm 0.05$ - $Y_{b,500c}^{\rm The300}$ \+ $K^{\rm
CCCP}_{500c}$ \+ HST $0.24\pm 0.02$ $0.94\pm 0.11$ $0.83\pm 0.08$ -
$Y_{b,500c}^{\rm The300}$ $(K_{500c}=1)$ \+ HST $0.27\pm 0.03$ $0.83\pm 0.09$
$0.79\pm 0.07$ - $Y_{b,500c}^{\rm The300}$ \+ $K^{\rm CMB}_{500c}$ \+ Planck
$0.17\pm 0.01$ $1.25\pm 0.04$ $0.94\pm 0.03$ - $Y_{b,500c}^{\rm The300}$ \+
$K^{\rm CLASH}_{500c}$ \+ Planck $0.30\pm 0.03$ $0.78\pm 0.11$ $0.78\pm 0.06$
- $Y_{b,500c}^{\rm FABLE}$ \+ $K^{\rm CLASH}_{500c}$ \+ Planck $0.28\pm 0.01$
$0.85\pm 0.06$ $0.83\pm 0.01$ - $Y_{b,500c}^{\rm The300}$ \+ $K^{\rm
CCCP}_{500c}$ \+ Planck $0.23\pm 0.02$ $0.99\pm 0.01$ $0.87\pm 0.03$ -
$Y_{b,500c}^{\rm The300}$ $(K_{500c}=1)$ \+ Planck $0.31\pm 0.03$ $0.74\pm
0.10$ $0.75\pm 0.07$ -

Table 4: Mean and standard deviation of $\Omega_{m}$, $\sigma_{8}$ and $S_{8}$
and $h$ from the joint analysis of the average sparsity and gas mass fraction
data for different $Y_{b,500c}$, $K_{500,c}$ and $h$ priors.

Sparsity + Gas Mass Fraction + BAO
$h$-prior $\Omega_{m}$ $\sigma_{8}$ $S_{8}$ $h$ Uniform $0.316\pm 0.013$
$0.757\pm 0.067$ $0.776\pm 0.064$ $0.696\pm 0.017$ HST $0.325\pm 0.012$
$0.736\pm 0.066$ $0.766\pm 0.065$ - Planck $0.306\pm 0.010$ $0.781\pm 0.061$
$0.788\pm 0.059$ -

Table 5: Mean and standard deviation of $\Omega_{m}$, $\sigma_{8}$ and $S_{8}$
and $h$ from the joint analysis of the average sparsity, gas mass fraction and
BAO data.