© C. Garrel et al. 2022

1. Introduction

As the largest gravitationally collapsed objects in the Universe, clusters of galaxies occupy a privileged position in astrophysical studies for two main reasons. The cluster number counts and spatial distribution of galaxy clusters as a function of mass and redshift are sensitive to both the growth of structure and the geometry of the Universe, and therefore these objects constitute powerful cosmological probes. While the purely gravitational aspect is theoretically well understood, the interplay between the three cluster components, namely galaxies (∼5%), gas (∼15%) and dark matter (∼80%), renders the physics of the intracluster medium (ICM) complex. A wide range of phenomena are involved: cooling through X-ray emission, enrichment and heating of gas through supernovae and AGN feedback, turbulence, and magnetic fields (see e.g., review by Allen et al. 2011). These processes make clusters interesting astrophysical laboratories and have motivated considerable computational efforts to reproduce their properties using hydrodynamic simulations (e.g., review by Borgani & Kravtsov 2011). The modelling of these properties is crucial in linking cluster observables like galaxy richness, velocity dispersion, gas mass, X-ray luminosity, and temperature (L_X, T_X) to the total cluster mass, a key component for cosmological studies.

In this context, a Very Large XMM programme was allocated in 2010: with its two spatially disconnected regions of 25 deg² each, the XXL survey was specifically designed to obtain robust cosmological constraints from the X-ray cluster population out to a redshift of unity. This survey is accompanied by an extensive multi-wavelength follow-up programme and has motivated the development of sophisticated detection and analysis procedures (Pierre et al. 2016, hereafter XXL Paper I). We refer the reader to this latter paper for a comprehensive bibliographical overview of cosmological X-ray cluster surveys. In the construction of the XXL cluster sample, two aspects were given special attention: (i) The cluster selection is solely described in terms of observed X-ray parameters: by selecting clusters in the two-dimensional count-rate versus apparent-size parameter space, we can ensure a sample purity of better than 95% and whose definition is independent of the cosmology. (ii) The cluster scaling relations entering the cosmological analysis are derived from the cluster sample data alone.

A first cosmological analysis of the brightest XXL clusters (the C1 sample containing 178 objects) was presented in Pacaud et al. (2018, hereafter XXL Paper XXV). This study, based on the modelling of the cluster redshift distribution (dn/dz), provided constraints on σ₈ and Ω_m with precision of the order of 10% and 20%, respectively. No cluster mass information was propagated in the analysis other than the resulting mass detection limit as a function of redshift and cosmology. A natural follow-up would be the subsequent analysis of the dn/dM/dz distribution, which is theoretically much more constraining than dn/dz. However, because the direct handling of the (cosmology-dependent) masses is difficult, we adopted a forward modelling based on the prediction of directly observable quantities; namely, the three-dimensional distribution of the count rates (CRs), hardness ratios (HRs), and redshifts of the selected cluster population (X-ray observable diagrams, hereafter XOD). This method (named ASpiX) was initiated by Clerc et al. (2012a,b) and further validated on analytical and numerical simulations (Pierre et al. 2017; Valotti et al. 2018). ASpiX is intrinsically equivalent to the study of the mass–redshift distribution because the mass information is encoded in the CR–HR–z distribution, but ASpiX is much simpler to use and is less affected by physics or cosmology degeneracies. The method consists in comparing the observed XOD with the predicted XODs as a function of cosmology and cluster evolutionary physics. In the present study, the comparison is performed adopting a Markov chain Monte-Carlo (MCMC) approach, in which selections of cosmological parameters and scaling relation coefficients are left free to vary; the predicted XOD are convolved by a realistic measurement error model.

The paper is organised as follows. Section 2 briefly recalls the main properties of the cluster sample. In Sect. 3, we describe the steps involved in the XOD construction. In Sect. 4 we present a first cosmological analysis under exactly the same hypotheses as in XXL Paper XXV; this allows a direct comparison of the two approaches. We further add constraints from the two-point correlation function from the same cluster sample. In Sect. 5, we update the study by using the revised scaling relations obtained from our recent lensing analysis of deep Hyper Suprime Camera images. The results, along with various sources of uncertainty, are discussed in Sect. 6 with constraints from other probes. Conclusions are drawn in Sect. 7. Appendix A describes the procedure used to measure the cluster quantities appearing in the XOD. Appendix B gives the details of the cosmological likelihood calculation in the CR–HR space, including the estimate of the sample variance.

Throughout the paper, we assume a spatially flat Λ cold dark matter (ΛCDM) model (see Sect. 4). We use the standard notation M_Δ to denote the cluster mass enclosed within a sphere of radius r_Δ, within which the mean overdensity is equal to Δ times the critical density of the Universe at a particular redshift z.

2. Cluster sample

The present paper deals with a sample of 178 XXL C1 clusters detected in the 47.36 deg² XXL survey. This is identical to the sample used in the previous XXL cosmological analysis (XXL Paper XXV). In this section, we recall the properties of this sample and describe the measurements of the cluster parameters used in the current study.

2.1. Sample

Adami et al. (2018), hereafter XXL Paper XX, published a sample of 365 clusters divided into two classes, namely the C1 and C2 class. The C1 subsample consists of 191 sources, achieving a high degree of purity (fewer than 5% of the point sources are misclassified as extended; Pacaud et al. 2006). We restricted the cosmological analysis to only the C1 spectroscopically confirmed clusters within the [0.05−1] redshift range. This led to the exclusion of 8 clusters that were outside of the redshift range and 5 without redshift estimates, resulting in a final sample of 178 clusters. In the present study, each cluster is characterised by four observable parameters: redshift (z), X-ray count rate (CR; defined as the number of X-ray counts received per second in the [0.5−2] keV energy band normalised to its on-axis value), hardness ratio (HR; defined as the ratio between the count rates in the [1−2] and [0.5−1] keV energy bands), and angular core radius (θ_c, assuming a β-profile with β = 2/3; Cavaliere & Fusco-Femiano 1976, 1978); CR and HR are equivalent to physical flux and colour. The CR–HR distribution of the XXL C1 sample is shown in Fig. 1.

Fig. 1. X-ray observable diagram (XOD) of the XXL C1 sample containing 178 objects, integrated over the [0.05−1] redshift range. The blue histograms show the 1D integrated CR and HR distributions. Error bars only account for shot noise.

2.2. Measurements

Along with a precise mapping of the selection function, the ASpiX method requires robust measurements of the CR and HR quantities with realistic error estimates. We describe the adopted procedure below.

The XXL detection pipeline (XAMIN; Pacaud et al. 2006) operates on the [0.5−2] keV images and provides a list of sources; a multi-PSF fit returns the source extent (θ_c) for the extended source model, and the resulting CR. The X-ray pipeline also provides the pixel segmentation mask for each source from the first detection pass. The fitted CR and θ_c values are ascribed measurement likelihoods by XAMIN, but errors are not provided for these quantities. In all that follows, the XAMIN output is used to deal with cluster selection issues only.

In order to obtain model-independent CR measurements along with associated errors, we apply a novel method based on Monte-Carlo sampling to fit the X-ray profile (pyproffit; Eckert et al. 2020) on the mosaic of overlapped XMM observations. Based on a multi-scale profile decomposition, this method allows robust CR and therefore HR measurements, together with an estimate of their uncertainties. Given that the mean number of collected photons per cluster is low (on average ∼200 counts for 10 ks exposures), we use a simplified minimisation algorithm for the θ_c measurements. The complete procedure is detailed in Appendix A.

3. Cosmological modelling

The main goal of the paper is to perform a forward cosmological analysis: the ASpiX method consists of fitting the CR − HR − z XOD. This approach was tested on simulations and described in a series of articles (Clerc et al. 2012a,b; Pierre et al. 2017; Valotti et al. 2018). In this section, we recall the principles and assumptions inherent to the method.

3.1. ASpiX method

Starting from a theoretical mass function, ASpiX reconstructs the XOD of a given cosmological-plus-cluster-physics model in order to match the observed XOD. The strength of this method is in the fact that it only relies on strictly observable X-ray parameters, which means that the cluster temperatures, luminosities, and masses are not explicitly computed.

We start from the differential mass function computed for a given cosmology, and expressed in terms of redshift (z) and sky area (Ω) folded with the XXL survey effective sky coverage. We compute the distribution in terms of M_Δ, z and cluster characteristic size:

$$\begin{array}{c}\hfill r_{\mathrm{\Delta }}={\left(\frac{3M_{\mathrm{\Delta }}}{4\pi \times \mathrm{\Delta }{\rho }_{\mathrm{c}}}\right)}^{1/3},\end{array}$$(1)

where r_Δ is the radius within which the average density is Δ times ρ_c, the critical density.

We use scaling relations linking mass and temperature [M_Δ–T], luminosity and temperature [L–T], and between r_Δ and the cluster core radius r_c, assuming a β model with β = 2/3, [r_c–r_Δ]. The relation between physical and apparent core radii reads:

$$\begin{array}{c}\hfill {\theta }_{\mathrm{c}}[\mathrm{arcsec}]=\frac{648000}{\pi }{r}_{\mathrm{c}}[\mathrm{Mpc}]/d_a(z)[\mathrm{Mpc}],\end{array}$$(2)

where d_a is the angular diameter distance.

We make use of the APEC model: the emission spectrum from collisionally ionised diffuse gas is calculated from the ATOMDB atomic database, (Smith et al. 2001) with a metallicity of 0.3 Z_⊙. Folding the spectrum with the EPIC (European Photon Imaging Camera) XMM response matrices provides us with count rates in the three energy bands of interest ([0.5−2] keV, [0.5−1] keV and [1−2] keV). From this, we subsequently construct the 4D z–CR–HR–θ_c diagram.

We stress that the mass information is implicitly encrypted in the θ_c, CR, and HR parameters via the scaling relations. We then apply the XXL survey selection function (f[CR,θ_c]) and finally convolve the XOD with the measurement-error model of each observable parameter except for z (spectroscopically measured and thus with negligible error). In the end, we integrate over θ_c to obtain the 3D z–CR–HR diagram expected for a given cosmology.

3.2. Assumptions and numerical inputs

For the purpose of this analysis, we use a Tinker mass function (Tinker et al. 2008) computed at an overdensity of Δ = 500. We disperse over the luminosity, temperature and core radius distributions by including a log-normal scatter around the mean scaling relations. The binning of the XOD is shown in Table 1.

3.2.1. Scaling relations

In the first part of the paper, in order to allow a direct comparison of the different methodologies applied, we continue to use the scaling relations of XXL Paper XXV modelled, as usual, by power laws:

$$\begin{array}{c}\hfill \frac{M_{500,\mathrm{WL}}}{X_{0,M-T}}={\left(\frac{T_{300\mathrm{kpc}}}{1\mathrm{keV}}\right)}^{{\alpha }_{M-T}}E{(z)}^{{\gamma }_{M-T}},\end{array}$$(3)

$$\begin{array}{c}\hfill \frac{L_{500,\mathrm{WL}}^{\mathrm{XXL}}}{X_{0,L-T}}={\left(\frac{T_{300\mathrm{kpc}}}{1\mathrm{keV}}\right)}^{{\alpha }_{L-T}}E{(z)}^{{\gamma }_{L-T}},\end{array}$$(4)

$$\begin{array}{c}\hfill \frac{r_{\mathrm{c}}}{X_{0,r_{\mathrm{c}}-r_{500}}}=r_{500},\end{array}$$(5)

where M_500, WL is the weak-lensing mass estimate within r₅₀₀, T_300 kpc is the cluster X-ray temperature measured inside 300 kpc, $L_{500,\mathrm{WL}}^{\mathrm{XXL}}$ is the luminosity within r₅₀₀ in the [0.5−2] keV energy band, and E(z) is the redshift evolution of the Hubble parameter, E(z)≡H(z)/H₀. The mass calibration only relies on weak lensing measurements based on the Canada-France-Hawaii Telescope (CFHT) lensing data (for a didactic review of cluster weak lensing, see Umetsu 2020). The scaling law parameters are summarised in Table 2.

During the analysis, two scaling relation parameters were introduced as nuisance parameters and marginalised. These parameters correspond to the ones indicated by uncertainties in Table 2. We do not include any scatter in the r_c − r₅₀₀ and M_500, WL–T_300 kpc relations within the base model so that comparisons can be made with XXL Paper XXV. Subsequently, we add a scatter of 0.1 in the r_c − r₅₀₀ in Sect. 5 to correspond to the updated scaling relations. We discuss the impact of larger values of the scatter in Sect. 6.3.

3.2.2. Selection function

Assuming a circular β = 2/3 model for extended sources, a cluster population with different count-rates and θ_c was generated for a range of XMM exposures. This process takes into account the instrumental characteristics (PSF distortion, vignetting, detector masks, background, and Poisson noise) for the three XMM detectors. Point sources are added at random over the XMM field of view, with a flux distribution following the log(N)−log(S) from Moretti et al. (2003) down to 5 × 10⁻¹⁶ erg s⁻¹ cm⁻²; the contribution of point sources below 4 × 10⁻¹⁵ erg s⁻¹ cm⁻² is included in the cosmic background component (Read & Ponman 2003). The XXL cluster-detection algorithm is then applied, allowing a statistical study to determine various levels of completeness and purity. The cluster selection is performed in the XAMIN output parameter space (EXT, EXT_STAT) and subsequently translated into the CR–θ_c plane. The details of the procedure are given in Pacaud et al. (2006). The XXL C1 selection function, matched to the XXL exposure and background maps, is shown in Fig. 2.

Fig. 2. XXL C1 selection function used in this analysis. From simulations of XMM cluster observations, the detection probability is expressed as a function of only observable quantities: the count rate and the apparent size of a β = 2/3 model, θc. The same selection function was used in XXL Paper XXV.

We stress that the selection function is mapped back into the intrinsic CR–θ_c space (the probability of detecting a cluster that has these parameters, and not the pipeline ones measured at those values); there is therefore no inconsistency in using CR values for the cosmological analysis that were measured using the pyproffit package.

3.2.3. Measurement errors

As shown in Clerc et al. (2012b), the inclusion of measurement errors changes the shape of the predicted X-ray observable diagrams. A precise estimate of the measurement errors is therefore a key step in the analysis.

The pyproffit package provides us with error estimates for each measurement. This allows us to subsequently model the relative measurement errors on CR, HR, and θ_c as a function of CR and θ_c using the following parametrisation:

$$\begin{array}{c}\hfill \%{\mathrm{err}}_x=a_x{\mathrm{CR}}^{b_x}{\theta }_{\mathrm{c}}^{c_x}(x=\mathrm{CR},\mathrm{HR},{\theta }_{\mathrm{c}}).\end{array}$$(6)

The choice of the CR–θ_c plane is a natural second-order approximation reflecting the fact that, physically, the brighter and more peaked a cluster, the better the measurement.

We perform a non-linear least square fit using the Levenberg-Marquardt algorithm (Levenberg 1944; Marquardt 1963) in order to constrain the {a_x}, {b_x}, and {c_x} coefficients. The resulting error models are shown in Fig. 3 and the coefficients are given in Table 3.

Fig. 3. Relative error model (in percent) for each observable used in this analysis. These models are computed as a function of CR and θc.

3.2.4. Likelihood

The log-likelihood model used to infer the cosmological parameters is given, for each redshift bin, by¹:

$$\begin{array}{cc}\hfill {\mathcal{L}}_{z_i}& =\overline{n}-{\displaystyle \sum _j{\widehat{N}}_{j}ln({\overline{n}}_j)}\hfill \\ \hfill & +\frac{1}{2}\ln (1+{\sigma }_{\delta }^2{\displaystyle \sum _j}{\widehat{N}}_j\frac{{\overline{b}}_j}{\overline{b}})\hfill \\ \hfill & -\frac{{\sigma }_{\delta }^2}{2}\times {{\displaystyle (\sum _j{\widehat{N}}_j\frac{{\overline{b}}_j}{\overline{b}}-\overline{n})}}^2{(1+{\sigma }_{\delta }^2{\displaystyle \sum _j}{\widehat{N}}_j\frac{{\overline{b}}_j^2}{{\overline{b}}^2})}^{-1},\hfill \end{array}$$(7)

where ${\overline{n}}_j$ is the number of predicted clusters in the (CR_j, HR_j) 2D bin (and $\overline{n}$ the number of predicted clusters in the redshift bin i) and ${\widehat{N}}_j$ is the number of observed clusters. Here, ${\overline{b}}_j$ is the mean galaxy cluster bias for the 2D bin j, and $\overline{b}$ is the mean bias of the survey (see Eqs. (B.9) and (B.10) of Appendix B). To calculate these quantities, we use the Tinker et al. (2010) bias model. ${\sigma }_{\delta }^2$ is the variance of the total number density contrast. Then, in Eq. (7), the first line is the usual shot-noise term and the second and third lines are the sample-variance terms. Finally,

$$\begin{array}{c}\hfill \mathcal{L}={\displaystyle \sum _i{\mathcal{L}}_{z_i}.}\end{array}$$(8)

Following the formalism presented in Valageas et al. (2011), we estimate that the sample variance value is ∼30% of the Poisson variance.

The cosmological parameters are constrained using a MCMC procedure by an affine-invariant ensemble (Goodman & Weare 2010), following the EMCEE algorithm (Foreman-Mackey et al. 2013). In order to optimise computation time while having enough statistics to control the chain convergence, we ran five independent chains in parallel, each with 2N walkers, with N specifying the number of free parameters. The chains are stopped when reaching the Gelman-Rubin convergence criterion of R − 1 < 0.03, after excluding a 20% burn-in phase.

4. Cosmological analysis with the scaling relations of XXL Paper XXV

We assume a flat ΛCDM model. We perform the Monte Carlo analysis and create contour plots by means of the GETDIST Python package (Lewis 2019). The displayed 1σ and 2σ confidence intervals show respectively the 68% and 95% limits.

4.1. XXL ASpiX alone

In this section, we present the cosmological constraints obtained from the XOD alone. We consider five free cosmological parameters within the ΛCDM framework: {Ω_m,  σ₈,  Ω_b,  n_s,  h}. Two scaling relation parameters are included as nuisance parameters: the M − T normalisation (X_{0, M − T}), and the L − T evolution index (γ_L − T) as summarised in Table 2; this already allows for significant freedom in the parametrisation of cluster physics unknowns and related cosmological dependence. These two parameters are marginalised during the Monte Carlo analysis.

We apply conservative Gaussian priors for the cosmological parameters which are not well constrained by cluster counts (namely Ω_b and n_s). These are centred on the Planck-2018 (Planck Collaboration VI 2020) values with errors multiplied by a factor of five in order not to force the agreement between our results and the Planck ones, i.e.: Ω_b = 0.0493 ± 0.0035, n_s = 0.9649 ± 0.022. The prior on the Hubble constant is chosen to be uniform within a [0.55−0.9] range.

We choose two different methods to estimate the improvement (or deterioration) in the parameter constraints in the Ω_m − σ₈ plane:

We first use Green’s theorem to compute the area, by line integral, inside the 1σ contours in the Ω_m − σ₈ plane for two different probes. Taking the square root of the ratio of these two areas gives us the estimated gain or loss in constraining power in the Ω_m − σ₈ plane.

We also quantify this improvement using a figure of merit (FoM) defined as ${\mathrm{FoM}}_{{\mathrm{\Omega }}_m-{\sigma }_8}=\sqrt{{\mathrm{Cov}}_{{\mathrm{\Omega }}_m-{\sigma }_8}^{-1}}$ where Cov_{Ωm − σ8} denotes the covariance matrix.

While keeping the same sample and scaling relations as in XXL Paper XXV, the ASpiX method allows us to improve constraints in the Ω_m − σ₈ plane by a factor of approximately two. This improvement was not unexpected because the [CR,  HR,  z] combination is comparable to the mass information, which did not enter the first XXL cosmological analysis, which was based on only dn/dz.

We note that, even though the XXL Paper XXV results, Planck constraints, and our new constraints are all compatible at the 2σ level (0.7σ posterior agreement² in the Ω_m − σ₈ plane between XXL ASpiX and XXL Paper XXV), our central values now show better agreement with the Planck results. We find ${\mathrm{\Omega }}_m=0.{342}_{-0.046}^{+0.038}$ versus 0.3165 ± 0.0085 for Planck, σ₈ = 0.829 ± 0.048 versus 0.8119 ± 0.0074, and S₈ = 0.882 ± 0.046 versus 0.834 ± 0.016 for Planck, leading to a 0.4σ posterior agreement in the Ω_m − σ₈ plane. The results are summarised in Table 4. The Ω_m − σ₈ contours are shown in Fig. 4 along with recent constraints from other cosmological probes.

Fig. 4. Ωm − σ8 cosmological constraints from XXL ASpiX (this work), Planck-2018 (Planck TT TE EE lowl lowE), Planck S-Z cluster counts (Planck15 Cluster), Planck lensing (CMB lensing potential analysis; Planck Collaboration VIII 2020), KIDS-450 (tomographic weak gravitational lensing of the 450 deg2 Kilo Degree Survey; Hildebrandt et al. 2017) and the previous XXL cosmological results, XXL Paper XXV (Pacaud et al. 2018).

4.2. XXL ASpiX + XXL clustering

Cosmological constraints from the 3D clustering analysis of the XXL cluster sample (two-point correlation function, 2PCF) were presented in Marulli et al. (2018). In this section, we combine the 2PCF and ASpiX results.

In order to perform the joint analysis, we run the MCMC procedure as before, and use the XXL 2PCF mean and covariance results as additional priors for all the parameters that are left free during the analysis, namely: {Ω_m,  σ₈,  Ω_b,  n_s,  h,  τ}. The 2PCF study was performed using seven free parameters, the six above-mentioned cosmological parameters and, in addition, the effective bias of the cluster sample, b_eff (see Marulli et al. 2018, for the detailed procedure).

We model the (redshift dependent) effective sample bias following Matarrese et al. (1997):

$$\begin{array}{c}\hfill b_{\mathrm{eff}}(z)={\mathcal{N}}^{-1}(z){\displaystyle {\int }_{\mathcal{M}}}\mathrm{d}ln{M}^{\prime }b(M^{\prime },z)\mathcal{N}(z,M^{\prime }),\end{array}$$(9)

and we define the averaged effective bias of the sample as:

$$\begin{array}{c}\hfill b_{\mathrm{eff}}=\sqrt{{\mathcal{N}}^{-2}{\int }_{\mathcal{Z}}\mathrm{d}{z}^{\prime }{b}_{\mathrm{eff}}^2(z^{\prime }){\mathcal{N}}^2(z^{\prime })},\end{array}$$(10)

where 𝒩(z, M) is the number of clusters with mass M and redshift z as predicted by a given cosmological scenario (including the selection effects), and b(M, z) is the dark matter halos bias computed using the Tinker et al. (2010) model for M₅₀₀.

While b_eff is free during the XXL clustering analysis, it is an output of the cluster counts (or ASpiX) analysis and depends on the selection function in the M − z space. The selection function in the M − z space depends, in turn, on the cosmology and the scaling relations. We therefore implement the results for b_eff from the XXL clustering analysis as an additional Gaussian term in the likelihood from Eq. (7).

The results are shown in Table 4. The Ω_m − σ₈ contours are shown in Fig. 5. The joint XXL ASpiX + 2PCF analysis reduces the uncertainties on Ω_m and σ₈ by 23% (figure of merit (FoM) increased by a factor of 1.3) compared to ASpiX alone. The Ω_m result is slightly lower, Ω_m = 0.314 ± 0.031, and σ₈ slightly higher, σ₈ = 0.840 ± 0.044, and with $S_8=0.{857}_{-0.050}^{+0.042}$. Nevertheless, the results are still in good agreement with Planck CMB.

Fig. 5. The Ωm − σ8 cosmological constraints of XXL ASpiX cluster counts alone (base), XXL 2PCF (clustering alone), and of the joint analysis.

4.3. XXL clusters + BAO joint analysis

In this section, we combine the ASpiX constraints with those obtained from baryon acoustic oscillation (BAO) measurements of clustering of galaxies. The BAO data used in this analysis are reported in Table 5. We describe the adopted methodology and present the results from the joint analysis.

We use two quantities to model the BAO distance measurements: (i) the spherically averaged distance,

$$\begin{array}{c}\hfill D_V(z)={\left[{(1+z)}^2{d}_a^2(z)\frac{\mathit{cz}}{H(z)}\right]}^{\frac{1}{3}},\end{array}$$(11)

and (ii) the sound horizon at the drag epoch (Eisenstein & Hu 1998),

$$\begin{array}{c}\hfill r_s(z_d)=\frac{2}{3k_{\mathrm{eq}}}\sqrt{\frac{6}{R(z_{\mathrm{eq}})}}ln\left[\frac{\sqrt{1+R(z_d)}+\sqrt{R(z_d)+R(z_{\mathrm{eq}})}}{1+\sqrt{R(z_{\mathrm{eq}})}}\right],\end{array}$$(12)

where z_d is the redshift at the drag epoch, z_eq is the matter-radiation equality redshift, k_eq is the scale of the particle horizon at the equality epoch, and R(z) is the ratio of the baryon to photon momentum density at redshift z. We model these four quantities following Eisenstein & Hu (1998) (with a CMB temperature of 2.725 K). The BAO distance is then given by

$$\begin{array}{c}\hfill {\mathcal{D}}_V(z)=\frac{D_V(z)}{r_s(z_d)}{r}_s^{\mathrm{fid}},\end{array}$$(13)

where $r_s^{\text{fid}}$ is the sound horizon computed for a chosen fiducial cosmology (see Table 5). The likelihood is therefore modified by adding a Gaussian log-likelihood term in Eq. (7) :

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\mathrm{BAO}}& =\frac{1}{2}{\displaystyle \sum _i{\left[\frac{{\mathcal{D}}^{\mathrm{th}}(z_i)-{\mathcal{D}}_i^{\mathrm{set}\mathrm{I}}}{{\sigma }_{{\mathcal{D}}_i^{\mathrm{set}\mathrm{I}}}}\right]}^2}\hfill \\ \hfill & +\frac{1}{2}[{\mathcal{D}}^{\mathrm{th}}-{\mathcal{D}}^{\mathrm{set}\mathrm{II}}]C_{\mathrm{WiggleZ}}^{-1}[{\mathcal{D}}^{\mathrm{th}}-{\mathcal{D}}^{\mathrm{set}\mathrm{II}}],\hfill \end{array}$$(14)

where the inverse covariance matrix, $C_{\mathrm{WiggleZ}}^{-1}$, comes from the fact that the three WiggleZ measurements are correlated (see Table 4 of Kazin et al. 2014).

Combining the BAO and the XXL 2PCF decreases the uncertainties by 62% (FoM increased by a factor of 2.6) and confirms the agreement with Planck (Ω_m = 0.317 ± 0.017, ${\sigma }_8=0.{845}_{-0.042}^{+0.035}$ and $S_8=0.{861}_{-0.042}^{+0.033}$). The results are shown in Table 4. The Ω_m − σ₈ contours are shown in Fig. 6.

Fig. 6. Ωm − σ8 constraints obtained by XXL ASpiX alone (base) and XXL ASpiX + XXL clustering + BAO analysis and Planck-2018 (Planck TT TE EE lowl lowE).

4.4. XXL ASpiX + Planck CMB

In this section, we combine our results with those of Planck-2018 (CMB anisotropy measurement) by means of importance sampling. The results are shown in Table 4. The Ω_m − σ₈ contours are shown in Fig. 7.

Fig. 7. Comparison between Planck CMB + XXL ASpiX and Planck CMB + Planck lens and Planck CMB alone in terms of the Ωm − σ8 constraints.

By combining our results with those of Planck-2018 (CMB anisotropies measurement) we reduce Planck uncertainties by 30% (FoM increased by a factor of 1.4) on Ω_m and σ₈ and find good agreement with the constraints from the combination of Planck-2018 and Planck lensing (lensing potential analysis of the temperature and polarisation data); see Fig. 6. We note that the XXL + Planck-2018 combination yields constraints comparable to those provided by the Planck-2018 + Planck lensing combination.

5. Cosmological modelling with updated scaling relations

Two independent mass-observable studies (Eckert et al. 2016, hereafter XXL Paper XIII, and Umetsu et al. 2020) suggest that cluster masses were overestimated in our first analysis based on the CFHT lensing data (XXL Paper IV). In this section, we rerun the cosmological analysis, assuming our newly determined scaling relations from the joint XXL-HSC (Hyper Suprime-Cam Survey) analysis by Umetsu et al. (2020) and Sereno et al. (2020).

For this purpose, we follow the formalism of Umetsu et al. (2020) for the scaling relations.

We consider the cluster true mass M_500, True as the fundamental property of galaxy clusters for the T − M relation and we use the weak lensing mass M_500, WL as a mass proxy.

Umetsu et al. (2020) characterised the weak lensing mass bias as a function of true cluster mass using cosmological N-body simulations (the dark-matter-only run from the BAHAMAS project; McCarthy et al. 2017, 2018). These authors estimated that, at the typical mass for the XXL sample (M₅₀₀ = 7 × 10¹³h⁻¹ M_⊙), the bias is approximately −11%. We therefore apply a correction for the weak lensing mass bias by assuming a constant of −11%:

$$\begin{array}{c}\hfill \frac{T_{300\mathrm{kpc}}}{X_{0,T-M}}={\left(\frac{M_{500,\mathrm{True}}}{7\times {10}^{13}{h}^{-1}{M}_{\odot }}\right)}^{{\alpha }_{T-M}}{\left(\frac{E(z)}{E(z=0.3)}\right)}^{{\gamma }_{T-M}},\end{array}$$(15)

with

$$\begin{array}{cc}\hfill {\log }_{10}(M_{500,WL})& ={\log }_{10}(M_{500,\mathrm{True}})\hfill \\ \hfill & +{\log }_{10}(1+b_{\mathrm{WL}})\pm {\sigma }_{{log}_{10}{M}_{\mathrm{WL}}},\hfill \end{array}$$(16)

and we assume a Gaussian prior on log₁₀(1 + b_WL) of log₁₀(1 + b_WL) = log₁₀(1 − 0.11)±5%/ln10 to marginalise over the mass calibration uncertainty of ±5%, see Umetsu et al. (2020). Here σ_log10MWL in Eq. (16) is the intrinsic scatter of weak-lensing mass at fixed true cluster mass, M_500, True.

The L − T relation is given by:

$$\begin{array}{c}\hfill \frac{L_{500,\mathrm{WL}}^{\mathrm{XXL}}}{X_{0,L-T}}={\left(\frac{T_{300\mathrm{kpc}}}{1\mathrm{keV}}\right)}^{{\alpha }_{L-T}}{\left(\frac{E(z)}{E(z=0.3)}\right)}^{{\gamma }_{L-T}},\end{array}$$(17)

and we keep Eq. (5) for the relation between r_c and r₅₀₀.

We fit the coefficients of the T − M and L − T relations (namely : {X_{0, T − M}, X_{0, L − T}, α_T − M, α_L − T, γ_T − M, γ_L − T, σ_T − M, σ_log10MWL, σ_L − T}, with σ_{T − M/L − T} the log-normal intrinsic scatters) using the publicly available LIRA package (Sereno 2016a,b) for the XXL C1 sample (using the procedure and measurements described in Sereno et al. 2020 and Umetsu et al. 2020). The results, computed for Ω_m = 0.28 and h = 0.7 in a flat ΛCDM universe, are shown in Table 6.

The effective impact of the cosmological dependence³ of weak lensing mass measurements and luminosities is expected to be small given the parameter range considered and the statistical and systematic errors inherent to our cluster sample. Nevertheless, to ensure better consistency, we model the effect of cosmology on the scaling relations a posteriori as follows:

• We use an analytical approximation (Sereno 2015) to account for the dependence of the lensing mass on cosmology:

  $$\begin{array}{c}\hfill M_{500,\mathrm{WL}}\propto D_{\mathrm{l}}^{-\frac{3{\delta }_{\gamma }}{2-{\delta }_{\gamma }}}{\left(\frac{D_{\mathrm{ls}}}{D_{\mathrm{s}}}\right)}^{-\frac{3}{2-{\delta }_{\gamma }}}H{(z)}^{-\frac{1+{\delta }_{\gamma }}{1-{\delta }_{\gamma }/2}},\end{array}$$(18)

  where D_l, D_s, and D_ls are the lens, the source, and the lens-source angular diameter distances respectively. In a first approximation, we assume a linear relation between the cluster redshift and mean redshift of the source galaxies:

  $$\begin{array}{c}\hfill ⟨z_{\mathrm{sources}}⟩=0.714z_{\mathrm{cluster}}+0.786.\end{array}$$(19)

  This results in a mean source–galaxy redshift of 1 for a cluster at 0.3 and of 1.5 for a cluster at redshift 1. We use δ_γ = 0.196, fitted on our C1 sample, to rescale the masses on a grid of Ω_m and h values. The ranges are defined to be Ω_m ∈ [0.1 − 0.8] and h ∈ [0.5 − 0.9], with a step of 0.1 for each parameter. This provides us with the masses for 40 combinations of Ω_m − h values. We then compute the T–M scaling relation for each [Ω_m − h] point of the grid to obtain the corresponding mean values and covariance matrices.

• Masses from the M–T relations are used to rescale r₅₀₀ (r_{500, rescale}); we then extrapolate the luminosities within r_{500, rescale} assuming a β-profile with a core radius r_c = 0.15 r_{500, rescale} and a slope β = 2/3. Finally, luminosities are normalised by the correction factor ${(d_{\mathrm{L}}/d_{\mathrm{L}}^{\mathrm{fid}})}^2$. This procedure provides us with rescaled luminosities for the 40 combinations of Ω_m − h values. We then compute the L–T scaling relation for each Ω_m − h point of the grid to obtain corresponding mean values and covariance matrices.

• Here, the cosmological analysis deals with five free cosmological parameters : {Ω_m,  σ₈,  Ω_b,  n_s,  h} plus six free scaling relation parameters: {X_{0, T − M}, X_{0, L − T}, α_T − M, α_L − T, γ_T − M,  γ_L − T}. In the MCMC, the values of the six scaling relation parameters are limited through adaptive Gaussian priors by interpolating the means and covariance matrices over the grid of 40 combinations of Ω_m − h values.

• We disperse temperatures and luminosities; these scatters are assumed to be independent of cosmology. We moreover introduce a log-normal scatter in the r_c − r₅₀₀ relation.

Resulting constraints on Ω_m − σ₈ (we refer to this model as XXL-HSC ASpiX) are shown in Fig. 8 and compared with the results of Sect. 4.1 (referring to this model as the base model). The constraints when adding the XXL 2PCF and BAO measurements are shown in Fig. 9. All the results are shown in Table 7.

Fig. 8. Impact of thawing parameters in the scaling relations. XXL ASpiX (base) refers to the results presented in Sect. 4.1 (two free scaling parameters); XXL-HSC ASpiX contours are the results following the methodology presented in Sect. 5 (six free cosmology-dependent scaling parameters); simple dark blue contours, same as XXL-HSC but the priors of the scaling coefficients are fixed to the indicated cosmology. We can see that, in this case, the results in Ωm − σ8 are shifted (lower Ωm and higher σ8) because of the fact that the scaling relation priors do not depend on the cosmology and then introduce a bias in the analysis.

Fig. 9. Same as Fig. 8 when adding constraints from XXL clustering and external BAO. The dashed line stands for the case where only two scaling coefficients are let free.

For XXL-HSC ASpiX, we find slightly higher Ω_m and σ₈ results and with larger error bars, ${\mathrm{\Omega }}_m=0.{378}_{0.13}^{0.068}$, ${\sigma }_8=0.{89}_{0.28}^{0.12}$ ($S_8=0.{970}_{-0.21}^{+0.067}$). Nevertheless, the results are compatible at the 1-σ level with our base model and the Planck CMB. Furthermore, because our uncertainties are now larger, we are compatible with the Planck S-Z cluster counts as well.

Adding the XXL clustering, we find a smaller Ω_m = 0.296 ± 0.034 and a higher ${\sigma }_8=0.{99}_{-0.23}^{+0.14}$ ($S_8=0.{98}_{-0.21}^{+0.11}$), which are fully consistent with the findings of Planck CMB at the 1-sigma level.

Combining XXL-HSC ASpiX with the XXL clustering and the BAO measurements, the results are shifted (Ω_m = 0.364 ± 0.015, ${\sigma }_8=0.{793}_{0.12}^{0.063}$, $S_8=0.{872}_{-0.12}^{+0.068}$) and we find results in better agreement with the Planck S-Z cluster sample while remaining consistent with Planck CMB at 2.2σ.

6. Discussion

6.1. Summary of results

Figure 4 shows that, as expected, the constraints on Ω_m − σ₈ have improved by more than a factor of two with respect to XXL Paper XXV under exactly the same hypotheses. This confirms, on real data, the power of the ASpiX forward modelling in terms of simplicity and accuracy. The size of the error bars is now comparable to that from the Planck S-Z cluster sample (Planck Collaboration XXIV 2016). At this point, it is important to recall that the Planck cluster sample contains almost three times as many clusters as XXL, that these clusters are much more massive, and that the XXL scaling relations do not rely on external cluster calibration samples or on hydrodynamical simulations, contrary to those of Planck; currently, the two data sets are still consistent at the 2-σ level, even though the favoured XXL cosmology is closer to the Planck CMB values. Combining the ASpiX constraints with the results from the cluster–cluster correlation function from the same sample improves the constraints by 23% (Fig. 5), and by 62 % (Fig. 6) when adding both cluster–cluster correlation function and BAO measurements.

In a second step, we reran the ASpiX analysis by implementing the updated scaling relations along with a modelling of their cosmological dependence, that is, by increasing the degrees of freedom on cluster physics from two to six (Fig. 8 and Table 6). As expected, this results in larger error bars: we now favour a higher σ₈ value, ${\sigma }_8=0.{89}_{-0.28}^{+0.12}$, but this is still compatible at the 1-σ with the results of the Planck CMB.

At this stage, we are close to having exhausted the cosmological information contained in the current data set relative to the XXL C1 cluster sample. It is instructive to review and discuss the various possible sources of systematic error that impinge upon these new results and, thus, assess their robustness.

A number of possible systematic uncertainties were already identified and discussed in the previous cosmological analysis of XXL Paper XXV. In the following section, we review further hypotheses and examine the impact on the initial base model. To quantify the robustness of our results, we analyse the posterior distribution of the S₈ = σ₈(Ω_m/0.3)^0.5 product.

6.2. Impact of error model

To study the impact of the error model, we arbitrarily modify the true error model by increasing errors by 20% and 50%, and we monitor the effect on the cosmological constraints (referring to them as Error + 20% and Error + 50% respectively).

We can see from Fig. 10 that increasing the relative measurement errors increases the uncertainty on S₈ slightly, without any drastic change in the mean result. Nevertheless, the general pattern seems to indicate that assuming excessive measurement errors tends to decrease S₈. The agreement between models is shown in Table 8.

Fig. 10. S8 = σ8(Ωm/0.3)0.5 results for all tests discussed in Sect. 6.

6.3. Scaling relation model

We now investigate how the results are impacted by different scaling relation models. First, we present the results when only adding a scatter of 0.5 in the r_c − r₅₀₀ relation. We then study the effect of adding scatter in the M − T relation.

Hereafter, we refer to the 0.5 scatter in the r_c − r₅₀₀ relation as the σ_{0.5, rc − r500} model. In Fig. 10, we can see that adding a 50% scatter in the r_c − r₅₀₀ relation favours a slightly higher S₈ value while increasing the uncertainties by only 3% compared to the base model. All in all, the results appear little affected (Table 8). If we had added θ_c as the fourth dimension in the XOD and used it in the cosmological inference, the error bars would probably be smaller, but more dependent on the scatter value (cf. Valotti et al. 2018).

In the base model, we did not implement a scatter in the M − T relation to keep the same configuration as in XXL Paper XXV (all scatter is supposed to be encapsulated in L − T). However, because HR directly depends on cluster temperature, it is logical to include a dispersion. We then include a 0.41 scatter in the M − T relation obtained from XXL Paper IV. We refer to this model from now on as σ_{0.41, M − T}.

This model is in good agreement with the base one with a significance level of posterior agreement in the Ω_m − σ₈ − Ω_b − n_s space of 0.1σ and 0.3σ in the Ω_m − σ₈ plane; see Table 8 and Fig. 10.

6.4. Impact on Ω_b priors

To ensure that the prior chosen for Ω_b is not too restrictive, we apply a Gaussian prior centred on the Planck 2018 values but with errors multiplied by a factor of 20 : Ω_b = 0.0493 ± 0.015 (i.e. our previous prior multiplied by 4). We refer to this model from now on as ‘Ω_b prior × 4’.

The resulting constraint on S₈ is shown in Fig. 10 and the agreement between models in Table 8.

We find that the results are fully consistent with the base model (0.02σ posteriors agreement in 4D Ω_m − σ₈ − Ω_b − n_s space and 0.04σ posteriors agreement in Ω_m − σ₈ only, see Table 8). In conclusion, because Ω_b prior × 4 becomes computationally expensive, we determine that it is relevant to keep the base model prior for Ω_b.

6.5. Error on the selection function

In this section, we study the impact of uncertainties on the selection function. This is a priori a key issue because an ill-determined selection directly biases the modelling of the cluster number counts. Currently, our selection is based on simulations assuming spherically symmetric and β = 2/3 profiles for the cluster emission. To quantify the impact of a poorly monitored cluster selection on the cosmological inference, we degrade the selection function; i.e. we blur the current function displayed in Fig. 2 by a 2D adaptive Gaussian filter characterised by:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{sel}.\mathrm{func}.}=0.05{\mathrm{CR}}^{-0.6}{\theta }_{\mathrm{c}}^{-0.4}.\end{array}$$(20)

As easily understandable, in this way, fainter and smaller clusters are more affected. We refer to this modelling from now on as Selfunc * Gauss.

The Selfunc * Gauss result on S₈ is shown in Fig. 10 and the agreement between models is shown in Table 8. While increasing uncertainties on S₈, the blurred selection function also lowers the mean S₈ value.

6.6. Remaining sources of uncertainty

In addition to the sources of systematic uncertainty reviewed above, we also note the main assumptions used in the course of the present study. Firstly, the covariance between the observable parameters (CR, HR, and θ_c) is neglected; the model has been slightly extrapolated in order to account for objects scattered out or in the measured domain. Furthermore, we do not consider the covariance between the scatters of the M–T and L–T scaling relations. In both cases, the scatter is assumed to be independent of the underlying cosmology. In the lensing analysis, we assume a linear relation between lens cluster redshift and the galaxy source photometric redshifts as stated in Eq. (19). Finally, we restrict our analysis to only one particular mass function (Sect. 3.2) for this study. We aim to examine the impact of these assumptions in the subsequent and final XXL analysis – consisting of a larger number of clusters – in order to determine the most accurate, unbiased cosmological estimates from the XXL sample.

7. Conclusions

Following Clerc et al. (2012b) and simulation case studies, we present the first application of the ASpiX cosmological forward modelling on real data with redshift information. The outcome confirms the flexibility and efficiency of the method. The constraints obtained from the 178 XXL C1 clusters, under various hypotheses, yield a precision comparable to that of the current BAO and Planck S-Z samples, as shown in Fig. 9. Nevertheless, the number of degrees of freedom left in the analysis reflects the accuracy of the recovered cosmological parameters, which is most easily seen when comparing the XXL base model alone (dashed blue contours) to the tightest constraints from this analysis (purple contours).

In short, the current results present an improvement by a factor of two compared to the preceding dn/dz analysis of the same sample. At this stage, we may recall the final cosmological modelling of the REFLEX survey number counts. This latter is based on the luminosity function of more than 800 clusters detected in the ROSAT All-Sky Survey (z < 0.4; Böhringer et al. 2014). Our base model analysis (Sect. 4.1) on the 178 C1 clusters includes four free cosmological parameters plus two scaling relation coefficients as nuisance parameters; in the REFLEX analysis, only the slope of the M–L relation was left free to vary free and it was assumed that the luminosity function does not evolve. Under these conditions, we find a precision on Ω_m comparable to that of REFLEX but almost twice better for σ₈; both parameter sets being compatible within the error bars.

Another cosmological analysis of RASS clusters was conducted as part of the ‘Weighing the Giants’ project. The 224 ‘Giants’ are massive clusters spanning the 0 < z < 0.3 redshift range. Gas masses from deep ROSAT and Chandra observations were subsequently derived for 94 of them. Independent mass calibration was achieved by weak gravitational lensing for 27 of these. This enabled the derivation of uniquely well-defined scaling relations and subsequently yielded a precision of the order of 5% on σ₈ and Ω_m (with standard priors on Ω_b, h, and n_s fixed; Mantz et al. 2015). Constraints are tighter than with the XXL clusters, but it is important to recall here that the only X-ray information used in the current study is the XMM 10ks survey data, which means a median number of photons of about 200 per cluster. We can anticipate that devoting very large amounts of X-ray follow-up time to the XXL clusters would lead to significant improvement on the WtG constraints thanks to the wider redshift range spanned by the XXL clusters.

Ultimately, the XXL XMM observation set will be reprocessed at full depth by running the detection algorithm on the mosaicked data (Faccioli et al. 2018, hereafter XXL Paper XXIV). This will not only increase the sensitivity but also the surveyed area, because the current cluster catalogue (XXL Paper XX) was extracted only from the single pointings, restricted to an off-axis distance of 13 arcminutes. It is therefore expected that the final XXL cosmological release will involve a sample twice as large as the current one, with a deeper C1 and C2 population. In the subsequent cosmological analysis, we shall add information from the third X-ray observable, the apparent core radius (Valotti et al. 2018). The final cosmological sample should bring an improvement of a factor of approximately 1.5 − 2 on the present constraints.

Using the same sample of 178 clusters, our next study will focus on the w parameter of the ΛCDM model. To this purpose, we shall make use of the HSC full depth information on the background galaxy photometric redshifts; the cosmological dependence of the cluster lensing masses will be rescaled as a function of w. The inclusion of the cluster two-point correlation function, while having little effect on the current study which is limited to Ω_m and σ₈, is expected to reduce the uncertainty on w by a factor of two (Pierre et al. 2011). Similarly, we shall allow for more flexibility in the determination of the cluster selection function: by considering a range of cluster ellipticities and quantifying the impact of cool cores or central AGN in the detection, we will be in a position to assess systematic uncertainties more precisely.

Because photometric redshifts are almost as efficient as spectroscopic redshifts in ASpiX (Clerc et al. 2012a), the application of the method to the upcoming eROSITA cluster sample should readily reveal most of the cosmological potential of the eROSITA sample.

Acknowledgments

XXL is an international project based around an XMM Very Large Programme surveying two 25 deg² extragalactic fields at a depth of ∼6  ×  10⁻¹⁵ erg cm⁻² s⁻¹ in the [0.5–2] keV band for point-like sources. The XXL website is http://irfu.cea.fr/xxl. The Saclay team acknowledges long term support from the Centre National d’Etudes Spatiales. This work was supported by the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. MS acknowledges 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’. LM acknowledges the grants PRIN-MIUR 2017 WSCC32 and ASI-INAF n. 2018-23-HH.0. KU acknowledges support from the Ministry of Science and Technology of Taiwan (grants MOST 106-2628-M-001-003-MY3 and MOST 109-2112-M-001-018-MY3) and from the Academia Sinica Investigator Award (grant no. AS-IA-107-M01). Finally, the authors would like to thanks the referee for the useful comments.

References

Appendix A: X-ray parameter measurements

The use of the XOD for cosmology requires accurate measurements along with realistic error estimates. In this section, we provide details on the various steps of the X-ray analysis for single clusters. Measurements are performed on the mosaicked co-added observations.

A.1. Step 1: Manually monitoring θ_bkg

We use the interactive fluxmes procedure (Clerc et al. 2012b) to determine the cluster-centric distance from which the background can be safely estimated (θ_bkg). In addition, this procedure allows us to check that all neighbouring sources have been correctly flagged by the pipeline; and possibly to re-adjust the corresponding masks. The fluxmes application is based on a curve of growth analysis and decides where to set the ‘no man’s land limit’ around the cluster. In the following, the particle background is measured using a very large number of ‘closed-filter’ observations, scaled to our XXL observations using the detector parts not exposed to the sky (corners); this results in two 25 deg² maps, corresponding to the mosaicked observations and exposure maps. The astrophysical background (passing through the XMM optics) is locally estimated in the [θ_bkg, 15 arcmin] annulus around each cluster in the course of the pyproffit procedure (see Sec. A.3).

A.2. Step 2: The θ_c measurements

Given that cluster angular sizes directly intervene in the selection (Fig. 2), we need to determine the apparent core radii along with associated uncertainties. This is performed under the hypothesis of a β = 2/3 model, as used to determine the selection function by means of simulations (Pacaud et al. 2006). The convolution by the PSF assumes a mean function over the entire mosaic (a cluster is generally seen at several off-axis positions on different observations):

$$\begin{array}{c}\hfill \mathrm{PSF}({\theta }_{on-axis})={[1+{\left(\frac{{\theta }_{on-axis}}{4.{765}^{″}}\right)}^2]}^{-1.505},\end{array}$$(A.1)

where θ_{on − axis} is the angular radius (on-axis). The fit is performed for each cluster on the [0.5-2] keV co-added XMM mosaic and includes the corresponding exposure and background maps; the neighbouring sources are removed from the data using the revised masks. We use iminuit (Dembinski et al. 2020) to fit the Cash log-likelihood (Cash 1979).

A.3. Step 3: Accurate count-rate measurements

The last step is to obtain best count-rate measurements along with realistic uncertainties for the three X-ray bands involved in the construction of the XOD (Sec. 2). To optimise the determination, we relax the hypothesis of a single β = 2/3 model and rather model the cluster emission by a linear combination of β-profiles following the pyproffit methodology (Eckert et al. 2020):

$$\begin{array}{c}\hfill \mathrm{SB}(\theta )={\displaystyle \sum _{i=1}^{N_{\mathcal{F}}}{\alpha }_i{\mathcal{F}}_i(\theta ),}\end{array}$$(A.2)

where N_ℱ is the total number of functions and α_i the model coefficients. {ℱ_i} is defined (Eckert et al. 2016) as

$$\begin{array}{c}\hfill {\mathcal{F}}_i(\theta )={[1+{\left(\frac{\theta }{{\theta }_{c,i}}\right)}^2]}^{-3{\beta }_i+1/2},\end{array}$$(A.3)

where θ is the distance to the cluster centre. The mean value of predicted counts (λ) in profile annulus a is given by

$$\begin{array}{c}\hfill {\lambda }_a=\mathrm{PSF}\ast \left(A_a{t}_a{\displaystyle \sum _{i=1}^{N_{\mathcal{F}}}}{\alpha }_i{\mathcal{F}}_i({\theta }_a)\right)+B_a.\end{array}$$(A.4)

Practically, we consider six β values within [0.6-3] and N = θ_bkg/5″ values for θ_c. The set of coefficients {α_i} in equation A.2 that maximises the log-likelihood given by equation A.5, namely

$$\begin{array}{c}\hfill -log(\mathcal{L})={\displaystyle \sum _{a=1}^{3N/2}{\mu }_a-N_{c,a}log({\lambda }_a),}\end{array}$$(A.5)

will describe the best-fit profile.

Eventually, we obtain the reconstructed profile by drawing 1000 posterior samples using No-U-Turn Sampler, NUTS (Hoffman & Gelman 2011). From this profile sample, we draw the count-rate posterior distribution. The various steps of the whole procedure are illustrated in Fig. A.1; for further technical details we refer the reader to (Eckert et al. 2020).

Fig. A.1. Example of cluster count-rate measurement with the pyproffit method. The displayed cluster is XLSSC 093 at a redshift of 0.429. (a) X-ray mosaic around the cluster; the image is 30 arcmin aside. (b) Particle background map. (c) Combined exposure map along with the masks hiding the neighbouring sources. (d) Extracted cluster profile (black crosses); the green line displays the particle background level extracted from map (b): this component is already subtracted from the displayed profile. (e) Overlaid on the extracted profile, the reconstructed (PSF-deconvolved) profile is shown in blue along with the 1-σ estimated uncertainty. (f) Count-rate posterior distribution (MOS1 normalised) of the reconstructed profile.

Appendix B: Likelihood

B.1. Probability distribution of galaxy clusters in the (z, CR, HR) space

We describe the distribution of X-ray clusters as a Poisson realisation of an underlying continuous field. Thus, the number of observed clusters ${\widehat{N}}_i$ in the bin i, in the 3D space (z_i, CR_i, HR_i) (or any other set of observables) follows the Poisson probability distribution

$$\begin{array}{c}\hfill P_{{\widehat{N}}_i}=\frac{{\widehat{n}}_i^{{\widehat{N}}_i}}{{\widehat{N}}_i!}{e}^{-{\widehat{n}}_i},\end{array}$$(B.1)

where ${\widehat{n}}_i$ is the continuous cluster density field. This continuous field ${\widehat{n}}_i$ is also a random variable because of the sample variance, following the fluctuations of the dark matter density field averaged over the survey volume. It is given by

$$\begin{array}{c}\hfill {\widehat{n}}_i={\displaystyle {\int }_{z_{i-}}^{z_{i+}}dz\frac{d\chi }{\mathit{dz}}{\mathcal{D}}^2\int d\mathbf{\Omega }\int dlnM\frac{d\widehat{n}}{dlnM}{\mathrm{\Theta }}_i[M,z],}\end{array}$$(B.2)

with

$$\begin{array}{cc}\hfill {\mathrm{\Theta }}_i[M,z]& ={\displaystyle {\int }_{-\infty }^{\infty }\frac{d{ϵ}_{\mathrm{CR}}}{\sqrt{2\pi }{\sigma }_{\mathrm{CR}}}{e}^{-{ϵ}_{\mathrm{CR}}^2/(2{\sigma }_{\mathrm{CR}}^2)}{\int }_{-\infty }^{\infty }\frac{d{ϵ}_{\mathrm{HR}}}{\sqrt{2\pi }{\sigma }_{\mathrm{HR}}}}\hfill \\ \hfill & \times e^{-{ϵ}_{\mathrm{HR}}^2/(2{\sigma }_{\mathrm{HR}}^2)}\hfill \\ \hfill & \times \mathrm{\Theta }({\mathrm{CR}}_{i-}<f_{\mathrm{CR}}(M,z)+{ϵ}_{\mathrm{CR}}<{\mathrm{CR}}_{i+})\hfill \\ \hfill & \times \mathrm{\Theta }({\mathrm{HR}}_{i-}<f_{\mathrm{HR}}(M,z)+{ϵ}_{\mathrm{HR}}<{\mathrm{HR}}_{i+}).\hfill \end{array}$$(B.3)

Here, f_CR(M, z) is the count rate associated with a cluster of mass M at redshift z, with a Gaussian scatter σ_CR. In practice, we can use a lognormal scatter by considering lnCR as our observable, or by replacing the Gaussian integral in (B.3) by a lognormal distribution. Then, the first factor Θ is a unit top-hat that is non-zero when the count rate falls in the bin i. Similar notations are used for the hardness ratio. Thus, Θ_i[M, z] is the probability that a cluster of mass M at redshift z falls in the 2D bin (CR_i, HR_i). In Eq.(B.2) we integrate the number of clusters over the redshift bin Δz_i = z_i+ − z_i−, the survey angular area ΔΩ, and the cluster mass M, where χ and 𝒟 are the radial and angular comoving distances and $\frac{d\widehat{n}}{dlnM}$ is the observed cluster mass function. Selection effects are included in the mass function $\frac{d\widehat{n}}{dlnM}$, which differs from the halo mass function and contains the response of the instrument.

At the level of the continuous field, the mean number of clusters in the bin i is

$$\begin{array}{c}\hfill {\overline{n}}_i\equiv ⟨{\widehat{n}}_i⟩=\mathrm{\Delta }\mathrm{\Omega }{\displaystyle {\int }_{z_{i-}}^{z_{i+}}dz\frac{d\chi }{\mathit{dz}}{\mathcal{D}}^2\int dlnM\frac{\mathit{dn}}{dlnM}{\mathrm{\Theta }}_i[M,z],}\end{array}$$(B.4)

where $\frac{\mathit{dn}}{dlnM}$ is the cluster mass function predicted by a given cosmological scenario (including the selection effects). If we neglect sample variance, we take ${\widehat{n}}_i={\overline{n}}_i$ without any scatter, and the Poisson distribution (B.1) has the fixed mean ${\overline{n}}_i$. This provides the shot-noise contribution to the measurement error bars that is associated with the discreteness of the cluster distribution.

To estimate the impact of the sample variance, we consider the covariance of the continuous number counts ${\widehat{n}}_i$. We have

$$\begin{array}{cc}\hfill C_{\mathit{ij}}& \equiv ⟨{\widehat{n}}_i{\widehat{n}}_j⟩-⟨{\widehat{n}}_i⟩⟨{\widehat{n}}_j⟩\hfill \\ \hfill & ={\displaystyle \int d{\chi }_{1}d{\mathbf{\Omega }}_{1}dln{M}_1\int d{\chi }_{2}d{\mathbf{\Omega }}_{2}dln{M}_2{\mathcal{D}}_1^2{\mathcal{D}}_2^2{\mathrm{\Theta }}_i[M_1,z_1]}\hfill \\ \hfill & \times {\mathrm{\Theta }}_j[M_2,z_2]\frac{\mathit{dn}}{dlnM}(M_1,z_1)\frac{\mathit{dn}}{dlnM}(M_2,z_2){\xi }_{12}({\mathit{x}}_1-{\mathit{x}}_2),\hfill \end{array}$$(B.5)

where ξ₁₂ is the two-point correlation of halos with masses of M₁ and M₂. We assume that the redshift bins are much larger than the correlation length of the clusters and that the correlation function can be factorised as

$$\begin{array}{c}\hfill {\xi }_{12}(x)=b_1{b}_2\xi (x),\end{array}$$(B.6)

where ξ(x) is the dark-matter correlation function and b(M, z) is the bias of clusters of mass M at redshift z. We use the Tinker et al. (2010) bias model during this analysis. Then, neglecting finite-sized effects associated with the borders of the survey volume, we can write the covariance matrix as

$$\begin{array}{c}\hfill C_{\mathit{ij}}={\delta }_{z_i,z_j}{\overline{n}}_i{\overline{b}}_i{\overline{n}}_j{\overline{b}}_j{\overline{\xi }}_i,\end{array}$$(B.7)

where δ_zi, zj is the Kronecker symbol with respect to the redshift bins i and j, ${\overline{\xi }}_i$ is the mean correlation in the redshift bin i, defined by

$$\begin{array}{c}\hfill {\overline{\xi }}_i={\displaystyle {\int }_{{\chi }_{i-}}^{{\chi }_{i+}}\frac{d\chi }{\mathrm{\Delta }{\chi }_i}\int \frac{d{\mathbf{\Omega }}_{1}d{\mathbf{\Omega }}_2}{{(\mathrm{\Delta }\mathrm{\Omega })}^2}\xi ({\mathit{x}}_1-{\mathit{x}}_2),}\end{array}$$(B.8)

and ${\overline{b}}_i$ is the mean bias defined by

$$\begin{array}{c}\hfill {\overline{n}}_i{\overline{b}}_i=\mathrm{\Delta }\mathrm{\Omega }{\displaystyle {\int }_{z_{i-}}^{z_{i+}}dz\frac{d\chi }{\mathit{dz}}{\mathcal{D}}^2\int dlnM\frac{\mathit{dn}}{dlnM}{\mathrm{\Theta }}_i[M,z]b(M,z).}\end{array}$$(B.9)

Because of the Kronecker redshift factor in the covariance matrix (B.7), different redshift bins are decoupled and we can analyse each redshift bin separately. Therefore, in the following we focus on a single redshift bin and the index j refers only to the 2D bins (CR_j, HR_j). It is also useful to consider the total number $\widehat{N}$ of clusters in the survey and its continuous counterpart $\widehat{n}$. For non-overlapping bins j we have

$$\begin{array}{c}\hfill \widehat{N}={\displaystyle \sum _j{\widehat{N}}_j,\widehat{n}=\sum _j{\widehat{n}}_j,\overline{n}=\sum _j{\overline{n}}_j,\overline{n}\overline{b}=\sum _j{\overline{n}}_j{\overline{b}}_j.}\end{array}$$(B.10)

Let us define the fluctuations ${\widehat{\delta }}_j$ and $\widehat{\delta }$ of the continuous number counts ${\widehat{n}}_j$ and $\widehat{n}$ in the bin j and in the full 2D volume (CR, HR),

$$\begin{array}{c}\hfill {\widehat{n}}_j=(1+{\widehat{\delta }}_j){\overline{n}}_j,\widehat{n}=(1+\widehat{\delta })\overline{n},\end{array}$$(B.11)

with means and covariances

$$\begin{array}{c}\hfill ⟨{\widehat{\delta }}_j⟩=0,⟨{\widehat{\delta }}_j{\widehat{\delta }}_{\ell }⟩={\overline{b}}_j{\overline{b}}_{\ell }\overline{\xi },⟨\widehat{\delta }⟩=0,⟨{\widehat{\delta }}_j\widehat{\delta }⟩={\overline{b}}_j\overline{b}\overline{\xi },\end{array}$$(B.12)

and unit correlation coefficients

$$\begin{array}{c}\hfill \frac{⟨{\widehat{\delta }}_j{\widehat{\delta }}_{\ell }⟩}{{⟨{\widehat{\delta }}_j^2⟩}^{1/2}{⟨{\widehat{\delta }}_{\ell }^2⟩}^{1/2}}=1,\frac{⟨{\widehat{\delta }}_j\widehat{\delta }⟩}{{⟨{\widehat{\delta }}_j^2⟩}^{1/2}{⟨{\widehat{\delta }}^2⟩}^{1/2}}=1.\end{array}$$(B.13)

This implies that the fluctuations $\{{\widehat{\delta }}_j,\widehat{\delta }\}$ are linear functions of each other, and we obtain

$$\begin{array}{c}\hfill {\widehat{\delta }}_j=\frac{{\overline{b}}_j}{\overline{b}}\widehat{\delta }.\end{array}$$(B.14)

This is a consequence of the factorisation (B.6) of the cluster correlation function. We also note ${\sigma }_{\delta }^2$ the variance of the total number density contrast in the redshift bin,

$$\begin{array}{c}\hfill {\sigma }_{\delta }^2\equiv ⟨{\widehat{\delta }}^2⟩={\overline{b}}^2\overline{\xi }.\end{array}$$(B.15)

B.2. Mean correlation $\overline{\xi }$

For small angular windows and large enough redshift bins, it is possible to simplify the computation of the mean correlation $\overline{\xi }$ defined in Eq.(B.8). In the flat-sky approximation, for circular windows of angular radius θ_s, this computation reads as

$$\begin{array}{c}\hfill \overline{\xi }={\displaystyle \int \frac{d\chi }{\mathrm{\Delta }\chi }\int \frac{d\mathit{\theta }d{\mathit{\theta }}^{\prime }}{{(\pi {\theta }_s^2)}^2}\int d\mathit{k}{e}^{i{k}_{\Vert }(\chi -{\chi }_0)+i{\mathit{k}}_{\perp }·\mathcal{D}({\mathit{\theta }}^{\prime }-\mathit{\theta })}P(k,z_0),}\end{array}$$(B.16)

where χ₀ is the comoving radial distance to the median redshift z₀ of the bin and P(k, z₀) is the matter density power spectrum at redshift z₀. For redshift bins that are not too shallow, Δχ ≫ 𝒟θ_s, the integral over χ along the line of sight suppresses the contributions from parallel wave numbers k_∥ > 1/(Δχ), so that $\overline{\xi }$ is dominated by transverse wave numbers k_⊥ ∼ 1/(𝒟θ_s) ≫ k_‖ and k ≃ k_⊥. This is Limber’s approximation in its Fourier form. Then, the integral over χ gives a Dirac factor 2πδ_D(k_∥), and the integration over k_∥ yields

$$\begin{array}{c}\hfill \overline{\xi }=\frac{2\pi }{\mathrm{\Delta }\chi }{\displaystyle \int \frac{d\mathit{\theta }d{\mathit{\theta }}^{\prime }}{{(\pi {\theta }_s^2)}^2}\int d{\mathit{k}}_{\perp }{e}^{i{\mathit{k}}_{\perp }·\mathcal{D}({\mathit{\theta }}^{\prime }-\mathit{\theta })}P(k_{\perp }).}\end{array}$$(B.17)

Introducing the 2D Fourier-space circular window W₂(k_⊥𝒟θ_s),

$$\begin{array}{c}\hfill W_2(k_{\perp }\mathcal{D}{\theta }_s)={\displaystyle \int \frac{d\mathit{\theta }}{\pi {\theta }_s^2}{e}^{i{\mathit{k}}_{\perp }·\mathcal{D}\mathit{\theta }}=\frac{2J_1(k_{\perp }\mathcal{D}{\theta }_s)}{k_{\perp }\mathcal{D}{\theta }_s},}\end{array}$$(B.18)

where J₁ is the Bessel function of first order and first type, and we obtain

$$\begin{array}{c}\hfill \overline{\xi }=\frac{4{\pi }^2}{\mathrm{\Delta }\chi }{\displaystyle {\int }_0^{\infty }kdkP(k)W_2{(k_{\perp }\mathcal{D}{\theta }_s)}^2.}\end{array}$$(B.19)

B.3. Likelihood for the number counts in (CR, HR)

We now extend the likelihood-ratio analysis of Cash (1979) to our case. Denoting θ_α the parameters of the model, such as the set of cosmological parameters and additional cluster parameters, we consider the likelihood $L_{z_i}({\theta }_{\alpha };{\widehat{N}}_j)$ in the redshift bin z_i defined by

$$\begin{array}{c}\hfill L_{z_i}({\theta }_{\alpha };{\widehat{N}}_j)=P({\widehat{N}}_j;{\theta }_{\alpha })={\displaystyle \int d\delta P(\delta )\prod _j{P}_{{\widehat{N}}_j}.}\end{array}$$(B.20)

Here we used the fact that the Poisson probabilities $P_{{\widehat{N}}_j}$ defined in Eq. (B.1) are governed by the continuous number counts n_j, which are characterised by their means ${\overline{n}}_j$ and the fluctuating part δ from Eq.(B.14). The means ${\overline{n}}_j$, the bias ${\overline{b}}_j$, and the variance ${\sigma }_{\delta }^2$ themselves depend on the cosmological parameters θ_α. At the level of the second-order moment for δ, assuming the survey is large enough so that the relative fluctuations δ of the total number of clusters in a redshift bin are small, we take P(δ) to be Gaussian so that its probability distribution is fully determined by its variance. This yields

$$\begin{array}{cc}\hfill L_{z_i}({\theta }_{\alpha };{\widehat{N}}_j)& ={\displaystyle {\int }_{-\infty }^{\infty }\frac{d\delta }{\sqrt{2\pi }{\sigma }_{\delta }}{e}^{-{\delta }^2/(2{\sigma }_{\delta }^2)}\prod _j\frac{{[(1+\frac{{\overline{b}}_j}{\overline{b}}\delta ){\overline{n}}_j]}^{{\widehat{N}}_j}}{{\widehat{N}}_j!}}\hfill \\ \hfill & \times e^{-(1+\frac{{\overline{b}}_j}{\overline{b}}\delta ){\overline{n}}_j},\hfill \end{array}$$(B.21)

which also reads as

$$\begin{array}{cc}\hfill L_{z_i}({\theta }_{\alpha };{\widehat{N}}_j)& ={\displaystyle \prod _j\frac{{\overline{n}}_j^{{\widehat{N}}_j}}{{\widehat{N}}_j!}{e}^{-{\overline{n}}_j}{\int }_{-\infty }^{\infty }\frac{d\delta }{\sqrt{2\pi }{\sigma }_{\delta }}{e}^{-{\delta }^2/(2{\sigma }_{\delta }^2)}}\hfill \\ \hfill & \times e^{{\sum }_j{\widehat{N}}_{j}ln(1+\frac{{\overline{b}}_j}{\overline{b}}\delta )-\overline{n}\delta }.\hfill \end{array}$$(B.22)

Here we used the last property in (B.10). The first product, independent of δ, is the usual shot-noise contribution, whereas the integral over δ of the second product gives the contribution from the sample variance. If the latter is negligible, σ_δ → 0, it goes to unity and we recover the shot-noise value. If the volume is large enough, the relative fluctuation σ_δ of the total number of clusters is small and we can expand the logarithm up to second order in δ,

$$\begin{array}{c}\hfill ln(1+\frac{{\overline{b}}_j}{\overline{b}}\delta )=\frac{{\overline{b}}_j}{\overline{b}}\delta -\frac{1}{2}\frac{{\overline{b}}_j^2}{{\overline{b}}^2}{\delta }^2+\cdots \end{array}$$(B.23)

This is valid if we have

$$\begin{array}{c}\hfill {\displaystyle \sum _j{\widehat{N}}_j\frac{1}{3}\frac{{\overline{b}}_j^3}{{\overline{b}}^3}{\delta }^3\sim \frac{{(\delta N)}^3}{3N^2}\ll 1,}\end{array}$$(B.24)

where we write δ = δN/N and ${\overline{b}}_j\sim \overline{b}$. For large survey sizes, with N ≫ 1, we typically expect $\delta N\sim \sqrt{N}$ so that the approximation (B.23) is valid. We can then perform the Gaussian integration in Eq.(B.22), which gives

$$\begin{array}{cc}& L_{z_i}({\theta }_{\alpha };{\widehat{N}}_j)={\displaystyle \prod _j\frac{{\overline{n}}_j^{{\widehat{N}}_j}}{{\widehat{N}}_j!}{e}^{-{\overline{n}}_j}{(1+{\sigma }_{\delta }^2\sum _j{\widehat{N}}_j\frac{{\overline{b}}_j^2}{{\overline{b}}^2})}^{-1/2}}\hfill \\ \hfill & \times exp\left[\frac{{\sigma }_{\delta }^2}{2}{{\displaystyle (\sum _j{\widehat{N}}_j\frac{{\overline{b}}_j}{\overline{b}}-\overline{n})}}^2{(1+{\sigma }_{\delta }^2{\displaystyle \sum _j}{\widehat{N}}_j\frac{{\overline{b}}_j^2}{{\overline{b}}^2})}^{-1}\right].\hfill \end{array}$$(B.25)

For the estimation of the cosmological and cluster parameters θ_α with the likelihood method (Cash 1979), we compare the logarithm ℒ = −lnL obtained for different sets of parameters. The estimated parameters ${\theta }_{\alpha }^{\text{obs}}$ are those that minimise ℒ and the variation of ℒ with θ_α provides the confidence intervals, following a χ² law. Thus, we consider

$$\begin{array}{cc}\hfill {\mathcal{L}}_{z_i}({\theta }_{\alpha };{\widehat{N}}_j)& =\overline{n}-{\displaystyle \sum _j{\widehat{N}}_{j}ln({\overline{n}}_j)}\hfill \\ \hfill & +\frac{1}{2}ln(1+{\sigma }_{\delta }^2{\displaystyle \sum _j}{\widehat{N}}_j\frac{{\overline{b}}_j^2}{{\overline{b}}^2})\hfill \\ \hfill & -\frac{{\sigma }_{\delta }^2}{2}{{\displaystyle (\sum _j{\widehat{N}}_j\frac{{\overline{b}}_j}{\overline{b}}-\overline{n})}}^2\times {(1+{\sigma }_{\delta }^2{\displaystyle \sum _j}{\widehat{N}}_j\frac{{\overline{b}}_j^2}{{\overline{b}}^2})}^{-1},\hfill \end{array}$$(B.26)

where we use $\sum _j{\overline{n}}_j=\overline{n}$, within a given redshift bin. Following common practice, we discard the term $ln({\widehat{N}}_j!)$ because it does not depend on the parameters θ_α and cancels out in the difference $\mathcal{L}({\theta }_{\alpha })-\mathcal{L}({\theta }_{{\alpha }^{\prime }})$. The first term is the usual shot-noise contribution while the other two terms are the sample-variance contribution, which vanishes for σ_δ → 0. The sums over j only need to run over the bins in the 2D space (CR, HR) that are not empty, as they come with a factor ${\widehat{N}}_j$. This ensures that the results are not affected if we enlarge the box in the 2D space (CR, HR) to a large volume far beyond the realistic domain, including regions that are always empty.

Going back to the 3D space (z, CR, HR), because the redshift bins are independent, we simply have for the full likelihood

$$\begin{array}{c}\hfill L={\displaystyle \prod _i{L}_{z_i},\mathcal{L}=\sum _i{\mathcal{L}}_{z_i}.}\end{array}$$(B.27)

In practice, the parameters θ_α should not be far from those derived from previous experiments, such as Planck (for the cosmological parameters). Then, as in the Fisher matrix analysis where we neglect the cosmological dependence of the covariance matrix, we can neglect the dependence on θ_α of the sample-variance quantities $\{{\sigma }_{\delta }^2,{\overline{b}}_j,\overline{b}\}$, which we compute for a reference cosmology labeled by the subscript (0), $\{{\sigma }_{\delta (0)}^2,{\overline{b}}_{j(0)},{\overline{b}}_{(0)}\}$. We then test the cosmological scenario through its predictions for the means ${\overline{n}}_j$. This implies that we can discard the second factor in Eq.(B.26), as it does not depend on θ_α, and write

$$\begin{array}{cc}\hfill {\mathcal{L}}_{z_i}({\theta }_{\alpha };{\widehat{N}}_j)& =\overline{n}-{\displaystyle \sum _j{\widehat{N}}_{j}ln({\overline{n}}_j)-\frac{{\sigma }_{\delta (0)}^2}{2}}\hfill \\ \hfill & \times {{\displaystyle (\sum _j{\widehat{N}}_j\frac{{\overline{b}}_{j(0)}}{{\overline{b}}_{(0)}}-\overline{n})}}^2{(1+{\sigma }_{\delta (0)}^2{\displaystyle \sum _j}{\widehat{N}}_j\frac{{\overline{b}}_{j(0)}^2}{{\overline{b}}_{(0)}^2})}^{-1}.\hfill \end{array}$$(B.28)

B.4. Behaviour of the likelihood ℒ

B.4.1. Cosmology selected by the data

Let us now investigate the behaviour of ℒ_zi as a function of the theoretical means ${\overline{n}}_j$. This will provide us some insight into the response of the cosmological parameters θ_α to the measurements ${\widehat{N}}_j$, through the associated means ${\overline{n}}_j$. Thus, the cosmology selected by the measurement corresponds to the set $\{{\overline{n}}_j\}$ that maximises the likelihood L_zi, i.e. that minimises the negative logarithm ℒ_zi. For 𝒩 2D bins at redshift z_i, this gives the 𝒩 equations $\frac{\partial \mathcal{L}}{\partial {\overline{n}}_j}=0$,

$$\begin{array}{c}\hfill 1\le j\le \mathcal{N}:1-\frac{{\widehat{N}}_j}{{\overline{n}}_j}+{\sigma }^2{\displaystyle \sum _{\ell }({\widehat{N}}_{\ell }\frac{{\overline{b}}_{\ell (0)}}{{\overline{b}}_{(0)}}-{\overline{n}}_{\ell })=0,}\end{array}$$(B.29)

where we use $\overline{n}=\sum _j{\overline{n}}_j$ and define

$$\begin{array}{c}\hfill {\sigma }^2={\sigma }_{\delta (0)}^2{(1+{\sigma }_{\delta (0)}^2{\displaystyle \sum _j}{\widehat{N}}_j\frac{{\overline{b}}_{j(0)}^2}{{\overline{b}}_{(0)}^2})}^{-1}>0.\end{array}$$(B.30)

First, we note that if the measurements are equal to the reference predictions, we recover the reference cosmology:

$$\begin{array}{c}\hfill \text{if}{\widehat{N}}_j={\overline{n}}_{j(0)}\text{then}{\overline{n}}_j={\overline{n}}_{j(0)},\end{array}$$(B.31)

where we use the sum rules (B.10).

Second, we note that the 𝒩 equations (B.29) admit the solution

$$\begin{array}{c}\hfill {\overline{n}}_j=\alpha {\widehat{N}}_j,\end{array}$$(B.32)

where α is a solution of the single equation

$$\begin{array}{c}\hfill 1-\frac{1}{\alpha }+{\sigma }^2({\widehat{N}}_{(b)}-\alpha \widehat{N})=0,\end{array}$$(B.33)

where we introduce

$$\begin{array}{c}\hfill {\widehat{N}}_{(b)}={\displaystyle \sum _j{\widehat{N}}_j\frac{{\overline{b}}_{j(0)}}{{\overline{b}}_{(0)}}.}\end{array}$$(B.34)

We can understand this from the fact that sample-variance fluctuations of the total number of clusters do not affect the relative counts in the different pixels (CR_i, HR_i), as seen in Eq.(B.14). Therefore, the inferred ratios ${\overline{n}}_j/{\overline{n}}_{\ell }$ are equal to the measured ratios ${\widehat{N}}_j/{\widehat{N}}_{\ell }$. Thus, this likelihood method selects cosmologies that predict the observed distribution profile in the 2D space (CR, HR), up to a uniform rescaling α.