The XXL survey: XLVI. Forward cosmological analysis of the C1 cluster sample

We present the forward cosmological analysis of an $XMM$ selected sample of galaxy clusters out to a redshift of unity. Following our previous 2018 study based on the dn/dz quantity alone, we perform an upgraded cosmological analysis of the same XXL C1 cluster catalogue (178 objects), with a detailed account of the systematic errors. We follow the ASpiX methodology: the distribution of the observed X-ray properties of the cluster population is analysed in a 3D observable space (count rate, hardness ratio, redshift) and modelled as a function of cosmology. Compared to more traditional methods, ASpiX allows the inclusion of clusters down to a few tens of photons. We obtain an improvement by a factor of 2 compared to the previous analysis by letting the normalisation of the M-T relation and the evolution of the L-T relation free. Adding constraints from the XXL cluster 2-point correlation function and the BAO from various surveys decreases the uncertainties by 23 and 53 % respectively, and 62% when adding both. Switching to the scaling relations from the Subaru analysis, and letting free more parameters, our final constraints are $\sigma_8$ = $0.99^{+0.14}_{-0.23}$, $\Omega_m$ = 0.296 $\pm$ 0.034 ($S_8 = 0.98^{+0.11}_{-0.21}$) for the XXL sample alone. Finally, we combine XXL ASpiX, the XXL cluster 2-point correlation function and the BAO, with 11 free parameters, allowing for the cosmological dependence of the scaling relations in the fit. We find $\sigma_8$ = $0.793^{+0.063}_{-0.12}$, $\Omega_m$ = 0.364 $\pm$ 0.015 ($S_8 = 0.872^{+0.068}_{-0.12}$), but still compatible with Planck CMB at 2.2$\sigma$. The results obtained by the ASpiX method are promising; further improvement is expected from the final XXL cosmological analysis involving a cluster sample twice as large. Such a study paves the way for the analysis of the eROSITA and future Athena surveys.


Introduction
As the largest gravitationally collapsed objects in the universe, clusters of galaxies occupy a twofold-privileged position in astrophysical studies.The cluster number counts and spatial distribution of galaxy clusters as a function of mass and redshift is sensitive to both the growth of structure and geometry of the universe, hence constituting a powerful cosmological probe.While the purely gravitational aspect is theoretically well understood the interplay between the three cluster components, 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 2 each, the XXL survey was specifically designed to obtain robust cosmological constraints from the X-ray cluster population out to a redshift of unity.It 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 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 vs. apparent-size parameter space, we can ensure a sample better than 95% pure 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 σ 8 and Ω m with a 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 (CR), hardness ratios (HR) and redshifts of the selected cluster population (X-ray observable diagrams, hereafter XOD).This method (named ASpiX) has been initiated by Clerc et al. (2012a,b) and further validated on analytical and numerical simulations (Pierre et al. 2017;Valotti Based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA et al. 2018).ASpiX is intrinsically equivalent to the study of the mass-redshift distribution since the mass information is encoded in the CR-HR-z distribution, but it is of much simpler use and less subject to physics/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 free; the predicted XOD are convolved by a realistic measurement error model.
The paper is organised as follows.Sec. 2 briefly recalls the main properties of the cluster sample.In Sec. 3, we describe the steps involved in the XOD construction.Sec. 4 performs 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 Section 5, we actualise 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 Sec.6, with constraints from other probes.Conclusions are drawn in Sec. 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 Section 4).We use the standard notation M ∆ to denote the cluster mass enclosed within a sphere of radius r ∆ , within which the mean overdensity equals ∆ times the critical density of the universe at a particular redshift z.

Cluster sample
The present paper deals with a sample of 178 XXL C1 clusters detected in the 47.36 deg 2 XXL survey.It 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 in two classes, namely the C1 and C2 class.The C1 sub-sample consists of 191 sources, achieving a high degree of purity (fewer than 5% of point sources 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.

Measurements
Along with a precise mapping of the selection function, the AS-piX method requires robust measurements of the CR and HR quantities with realistic error estimates.We describe below the adopted procedure.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 what 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, thus, 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.

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.

ASpiX method
Starting from a theoretical mass function, ASpiX reconstructs the XOD of a given cosmological plus cluster physics model, in   * ) The selection function is applied in the CR-θ c plane, but the θ c distribution is not directly used to constrain the cosmological parameters.The XOD is then integrated over the θ c and convolved with the error model, prior to the fit.
order to match the observed XOD.The strength of this method is to rely only 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 of the cluster characteristic size: 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: 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.

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.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.

Scaling relations
In the first part of the paper, in order to allow a direct comparison of the different methodologies applied, we stick to the scaling relations of XXL paper XXV modelled, as usual, by power laws: where M 500,WL is the weak lensing mass estimate within r 500 , T 300kpc the cluster X-ray temperature measured inside 300 kpc, L XXL 500,WL the luminosity within r 500 in the [0.5-2] keV energy band and E(z) is the redshift evolution of the Hubble parameter, E(z) ≡ H(z)/H 0 .The mass calibration only relies on weak lensing measurements based on the 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 500 and M 500,WL -T 300kpc relations within the base model for the purposes of comparing to XXL paper XXV.Subsequently, we add a scatter of 0.1 in the r c − r 500 in Section 5 to correspond to the actualised scaling relations.We discuss the impact of larger values of the scatter in Section 6.3.

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 −16 erg s −1 cm −2 ; the contribution of point sources below 4 × 10 −15 erg s −1 cm −2 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. We stress that the selection function is mapped back into the intrinsic CR-θ c space (the probability to detect a cluster that has these parameters -not the pipeline ones measured at those values); there is therefore no inconsistency in using for the cosmological analysis, CR values that were measured using the pyproffit package.

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 thus a key step of 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: 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) 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.

Likelihood
The log-likelihood model to infer the cosmological parameters is given, for each redshift bin, by 1 : where n j is the number of predicted clusters in the (CR j ,HR j ) 2D bin (and n the number of predicted clusters in the redshift bin i) and N j is the number of observed clusters.b j is the mean galaxy cluster bias for the 2D bin j and b the mean bias of the survey (see Eq. B.9 and B.10 of Appendix B).To calculate these quantities, we use 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: 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 Markov Chain Monte-Carlo procedure by an affine-invariant ensemble (Goodman & Weare 2010), following the EMCEE algorithm (Foreman-Mackey et al. 2013).In order to optimise computational time while having enough statistics to control the chain convergence, we ran five independent chains in parallel, each chain had 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.

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.

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 , σ 8 , Ω 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 ).
Notes.For the base model, we do not quote constraints on h since this parameter is poorly constrained by cluster counts and the posterior distributions are driven by the hard prior.S 8 is defined to be S 8 ≡ σ 8 (Ω m /0.3) 0.5 .The last column indicates the priors used in this analysis.N(µ, σ2 ) corresponds to a Gaussian prior with mean µ and variance σ 2 and U(A, B) a uniform prior within the range [A,B].We do not quote constraints on the nuisance parameter used in the analysis.by a factor of 5 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.
While keeping the same sample and scaling relations as in XXL paper XXV, the ASpiX method allows us to improve constraints in the Ω m − σ 8 plane by of a factor ∼ 2. This improvement 2 was not unexpected since the [CR, HR, z] combination is comparable to the mass information, which did not enter the first XXL cosmological analysis, based only on 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 3 in the Ω m − σ 8 plane between XXL ASpiX and XXL paper XXV) our central values now show a better agreement with the Planck results.We find Ω m = 0.342 +0.038 −0.046 versus 0.3165 ± 0.0085 for Planck, σ 8 = 0.829 ± 0.048 versus 0.8119 ± 0.0074 and S 8 = 0.882 ± 0.046 versus 0.834 ± 0.016 for Planck, leading to a 0.4σ posterior agreement in the Ω m − σ 8 plane.The results are summarised in Table 4.The Ω m − σ 8 contours are shown in Fig. 4 along with recent constraints from other cosmological probes.

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 previously, and use the XXL 2PCF mean and covarior loss in constraining power in the Ω m − σ 8 plane.We also quantify this improvement using a figure of merit (FoM) defined as FoM Ωm−σ 8 = Cov −1 Ωm−σ 8 where Cov Ωm−σ 8 denotes the covariance matrix. 3To compute the agreement between two different probes, we rely on the following process.We first draw a representative sample for each posterior of interest from the two probes.We compute the distance between each pair of points of these samples.We build, from this distance sample, the probability distribution using kernel density estimate.We then estimate the probability to exceed (PTE) by integrating the probability distribution over the interval [0-P(0)], with P(0) the probability of a distance equal to 0. The same formalism is used in Bocquet et al. (2019).The corresponding significance level is computed assuming Gaussian statistics.To insure that the results are not impacted by randomisation effects, we repeat this process one hundred times and present the mean significance level.ance results as additional priors for all the parameters that are left free during the analysis, namely : {Ω m , σ 8 , Ω 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 e f f (see Marulli et al. 2018, for the detailed procedure).
We model the (redshift dependent) effective sample bias following Matarrese et al. (1997) : and we define the averaged effective bias of the sample as: where N(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), the dark matter halos bias computed using the Tinker et al. ( 2010) model for M 500 .
While b e f f is free during the XXL clustering analysis, it is an output of the cluster counts (or ASpiX) analysis.It 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 on the scaling relations.We therefore implement the results for b e f f from the XXL clustering analysis as an additional Gaussian term in the likelihood from equation 7.
The results are shown in Table 4.The Ω m − σ 8 contours are shown in Fig. 5.The joint XXL ASpiX + 2PCF analysis reduces the uncertainties on Ω m and σ 8 by 23% (FoM increased by a factor of 1.3) compared to ASpiX alone.The Ω m result is slightly lower, Ω m = 0.314 ± 0.031, and σ 8 slightly higher, σ 8 = 0.840 ± 0.044 and with S 8 = 0.857 +0.042 −0.050 .Nevertheless, the results are still in good agreement with Planck CMB.

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 : where z d is the redshift at the drag epoch, z eq is the matterradiation 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 : where r f id s 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 equation 7 : where the inverse covariance matrix, C −1 WiggleZ , 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, σ 8 = 0.845 +0.035 −0.042 and S 8 = 0.861 +0.033 −0.042 ).The results are shown in Table 4.The Ω m − σ 8 contours are shown in Fig. 6.XXL ASpiX (base) refers to the results presented in section 4.1 (2 free scaling parameters); XXL-HSC ASpiX contours are the results following the methodology presented in section 5 (6 free, cosmology-dependent, scaling parameters); simple dark blue contours, same as XXL-HSC but the prior 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 ) due to the fact that the scaling relation priors do not depend on the cosmology and then introduce a bias in the analysis.
measurement) we reduce Planck uncertainties by 30% (FoM increased by a factor of 1.4) on Ω m and σ 8 and find a 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 comparable constraints to the Planck-2018 + Planck lensing combination.

Cosmological modelling with actualised scaling relations
Two independent mass-observable studies (Eckert et al. 2016, hereafter XXL paper XIII, andUmetsu 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. 2017McCarthy et al. , 2018)).They estimated that, at the typical mass for the XXL sample (M 500 = 7 × 10 13 h −1 M ), the bias is approximately -11%.We therefore apply a correction for the weak lensing mass bias by assuming a constant of -11%:  16is the intrinsic scatter of weak-lensing mass at fixed true cluster mass, M 500,True .
The L − T relation is given by: and we keep equation 5 for the relation between r c and r 500 .We fit the coefficients of the T − M and L − T relations (namely 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 andUmetsu 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 4 of weak lensing mass measurements and luminosities is expected to be small given the parameter range considered and the statistical/systematics errors inherent to our cluster sample.Nevertheless, to ensure better consistency, we model -a posteriori -the effect of cosmology on the scaling relations as follows: • We use an analytical approximation (Sereno 2015) to account for the dependence of the lensing mass on cosmology: where D l , D s , 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 : 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 • Masses from the M-T relations are used to rescale r 500 (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 L /d f id L ) 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 , σ 8 , Ω b , n s , h} plus 6 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 500 relation.
Resulting constraints on Ω m − σ 8 (we will refer to this model as XXL-HSC ASpiX) are shown in Fig. 8 and compared with the results of Sec.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.We find, for XXL-HSC ASpiX, Ω m and σ 8 results slightly higher and with larger error bars, Ω m = 0.378 0.068 0.13 , σ 8 = 0.89 0.12 0.28 (S 8 = 0.970 +0.067 −0.21 ).Nevertheless, the results are compatible at the 1σ level with our base model and the Planck CMB.Of course, since 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 σ 8 = 0.99 +0.14 −0.23 (S 8 = 0.98 +0.11 −0.21 ) letting us fully consistent with 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, σ 8 = 0.793 0.063 0.12 , S 8 = 0.872 +0.068 −0.12 ) and we find results in a better agreement with the Planck S-Z cluster sample while being consistent with Planck CMB at 2.2σ.Notes.S 8 is defined to be S 8 = σ 8 (Ω m /0.3) 0.5 .We do not quote constraints on the nuisance parameters used in this analysis.Furthermore, the mean and covariance of the Gaussian priors for the 6 free scaling relation parameters (namely {X 0,T −M , X 0,L−T , α T −M , α L−T , γ T −M , γ L−T }) are not shown here since they are rescaled as a function of cosmology (described in Section 5).

Result summary
Fig. 4 shows that, as expected, the constraints on Ω m − σ 8 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 et al. 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 Planck's; currently, the two data set 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 have rerun the ASpiX analysis by implementing the actualised scaling relations along with a modelling 0.80 0.85 0.90 0.95 S 8 = σ 8 (Ω m /0.3) 0.5
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 8 = σ 8 (Ω m /0.3) 0.5 product.

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 slightly increases the uncertainty on S 8 without any drastic change in the mean result.Nevertheless, the general pattern seems to indicate that assuming excessive measurement errors tends to decrease S 8 .The agreement between models is shown in Table 8.

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 500 relation.Then we study the effect of adding scatter in the M − T relation.
From now on, we will refer to the 0.5 scatter in the r c − r 500 relation as the σ 0.5,r c −r 500 model.In Fig. 10, we can see that adding a 50% scatter in the r c − r 500 relation favours a slightly higher S 8 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 4th dimension in the XOD and used it in the cosmological inference, the error bars would probably have been 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 will 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 − σ 8 − Ω b − n s space of 0.1σ and 0.3σ in the Ω m − σ 8 plane, see Table 8 and Fig. 10.

Impact on Ω b 's 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 will refer to this model from now on as Ω b prior × 4.
The resulting constraint on S 8 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 − σ 8 − Ω b − n s space and 0.04σ posteriors agreement in Ω m − σ 8 only, see Table 8).To conclude, because Ω b prior × 4 becomes computationally expensive, we determine that it is relevant to keep the base model prior for Ω b .

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 illdetermined 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: As easily understandable, in this way, fainter and smaller clusters are more affected.We will refer to this modelling from now on as Selfunc * Gauss.
The Selfunc * Gauss result on S 8 is shown in Fig. 10 and the agreement between models is shown in Table 8.While increasing uncertainties on S 8 , the blurred selection function also lowers the mean S 8 value.

Remaining sources of uncertainty
In addition to the sources of systematic uncertainties 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 (Sec.3.2) for this study.We aim to examine the impact of these assumptions in the subsequent and final XXL analysisconsisting of a larger number of clusters -to determine the most accurate, unbiased cosmological estimates from the XXL sample.

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 Figure 9.Nevertheless, the number of degrees of freedom left in the analysis reflect the accuracy of the recovered cosmological parameters, e,g.by 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.It 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 (Sec.4.1) on the 178 C1 clusters includes four free cosmological parameters plus two scaling relation coefficients as nuisance parameters; the REFLEX analysis let only the slope of the M-L relation free and 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 as better for σ 8 ; both parameter sets being compatible within the error bars.
Another cosmological analysis of RASS clusters has been conducted by 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 could be subsequently derived for 94 of them.Independent mass calibration was achieved by weak gravitational lensing for 27.This enabled the derivation of uniquely well-defined scaling relations and, subsequently, yielded a precision of the order of 5% on σ 8 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 ∼ 200 per cluster.We can anticipate that devoting very large amounts of X-ray follow-up time to the XXL clusters would outperform 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 i.e. 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, since 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 thus 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 ∼ 1.5 − 2 on the present constraints.
Using the same sample of 178 clusters, our immediate 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 2-point correlation function, while having little effect on the current study limited to Ω m and σ 8 , 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 shall be in a position to assess more precisely systematic uncertainties.
Because photometric redshifts are almost as efficient as spectroscopic redshifts in ASpiX (Clerc et al. 2012a), the application of the method to the up-coming eROSITA cluster sample should readily reveal most of the eROSITA sample's cosmological potential.To estimate the impact of the sample variance, we consider the covariance of the continuous number counts ni .We have where ξ 12 is the two-point correlation of halos of mass M 1 and M 2 .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 d k e ik (χ−χ 0 )+ik ⊥ •D(θ −θ) P(k, z 0 ), (B.16)where χ 0 is the comoving radial distance to the median redshift z 0 of the bin and P(k, z 0 ) is the matter density power spectrum at redshift z 0 .For redshift bins that are not too shallow, ∆χ Dθ s , the integral over χ along the line of sight suppresses the contributions from parallel wave numbers k > 1/(∆χ), so that ξ is dominated by transverse wave numbers k ⊥ ∼ 1/(Dθ 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 Introducing the 2D Fourier-space circular window W 2 (k ⊥ Dθ s ), where J 1 is the Bessel function of first order and first type, we obtain Appendix 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 (θ α ; N j ) in the redshift bin z i defined by Here we used the fact that the Poisson probabilities P N j , defined in Eq.(B.1), are governed by the continuous number counts n j , which are characterised by their means n j and the fluctuating part δ from Eq.(B.14).The means n j , the bias b j and the variance σ 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 which also reads as 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 δ, where we wrote δ = δN/N and b j ∼ b.For large survey sizes, with N 1, we typically expect δN ∼ √ N so that the approximation (B.23) is valid.Then, we can perform the Gaussian integration in Eq.(B.22), which gives For the estimation of the cosmological and cluster parameters θ α with the likelihood method (Cash 1979), we compare the logarithm L = − ln L obtained for different sets of parameters.The estimated parameters θ obs α are those that minimise L and the variation of L with θ α provides the confidence intervals, following a χ 2 law.Thus, we consider where we used j n j = n, within a given redshift bin.We discarded as usual the term ln( N j !), because it does not depend on the parameters θ α and cancels out in the difference L(θ α ) − L(θ α ).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 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 In practice, the parameters θ α should not be far from those derived from previous experiments, such as Planck (for the cosmological parameters).Then, as in Fisher matrix analysis where we neglect the cosmological dependence of the covariance matrix, we can neglect the dependence on θ α of the sample-variance quantities {σ 2 δ , b j , b}, which we compute for a reference cosmology labeled by the subscript (0), {σ 2 δ(0) , b j(0) , b(0) }.Then, we test the cosmological scenario through its predictions for the means n j .This implies that we can discard the second factor in Eq.(B.26), as it does not depend on θ α , and write Let us now investigate the behaviour of L z i as a function of the theoretical means n j .This will provide us some insight into the response of the cosmological parameters θ α to the measurements N j , through the associated means n j .Thus, the cosmology selected by the measurement corresponds to the set {n j } that maximises the likelihood L z i , i.e. that minimises the negative logarithm L z i .For N 2D bins at redshift z i , this gives the N equations where we used n = j n j and we defined First, we note that if the measurements are equal to the reference predictions we recover the reference cosmology: if N j = n j(0) then n j = n j(0) , (B.31) where we used the sum rules (B.10).Second, we note that the N equations (B.29) admit the solution where α is a solution of the single equation where we introduced 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 n j /n are equal to the measured ratios N j / N .Thus, this likelihood method selects cosmologies that predict the observed distribution profile in the 2D space (CR, HR), up to a uniform rescaling α.
Next, the quadratic equation (B.33) has two solutions, with the asymptotic behaviours for σ 2 → 0, The physical solution is α − , which goes to unity when the sample variance is negligible and we recover the shot-noise likelihood, where the inferred cosmological values are n j = N j .
The second solution α + is not physical and is due to the approximations in our treatment, such as (B.23).Indeed, it would correspond to a large uniform density fluctuation δ −1, where the approximation (B.23) is no longer valid.In practice α + does not appear and does not impair the likelihood algorithm because we restrict the search in the cosmological parameter space to a small realistic region, not too far from the Planck values (i.e., we do not consider cosmologies that would predict ten times more clusters or further than the Planck concordance cosmology).
We can see from Eq.(B.33) that α = 1 is a solution if N(b) = N.This leads to a generalisation of the solution (B.31), if N j = βn j(0) then n j = N j , (B.37) which applies for any β > 0. Therefore, if the measurements N j follow the same 2D profile as the reference cosmology, up to a uniform multiplicative factor β, there is no rescaling and the likelihood method selects cosmologies that predict the same number counts n j as those that are measured.
If N(b) > N, we can see from Eq.(B.36) that α < 1.Therefore, the likelihood selects cosmologies that predict mean counts n j that are lower than the measured values N j , and the required increase up to N j is explained by a local positive fluctuation δ > 0 of the density field, arising from a sample-variance effect.This can be understood from the fact that N(b) > N means that higher-bias pixels have a greater count than expected.This points towards a positive density fluctuations δ > 0. This is similar to the well-known Kaiser derivation of the bias of rare objects like clusters, with respect to the underlying dark matter density field.There, using for instance the Press-Schechter mass function or the peak formalism, it is noticed that positive largescale matter density fluctuations δ > 0 enhance the formation of rare massive objects, and the more extreme the object (i.e. a larger mass) the greater the enhancement.This means a larger bias for more massive halos, due to the increased sensitivity of the large-mass tail of the mass function.Reversing this picture, we can see that enhanced number counts of rare massive halos, i.e. of high-bias objects, arise from positive fluctuations of the underlying matter density field.Therefore, in our case (and more generally), N(b) > N signals an enhancement of high-bias objects and hence a positive underlying density fluctuation δ > 0. To accommodate this amplification with the observed values N j , the means n j must then be somewhat smaller than the targets N j .In  where n j is the typical number count in a 2D cell.Then, we expect the interval of confidence on the cosmological parameters to increase by a factor of the order of 1 + σ 2 n j /2, when we include the effect of the sample variance.Here n j should correspond to a binning that is well adapted to the survey, that is, which corresponds to the amount of information that can be drawn from the observations.By choosing increasingly small bins one decreases n j and the apparent magnitude of δλ in Eq.(B.43), but this is compensated by the larger size of the matrix H j and the greater number of constraints, which are mostly degenerate.
For the likelihood (B.28) to be meaningful, the Hessian H should be positive definite in the neighbourhood of the reference point (0), so that the solutions of Eq.(B.29) correspond indeed to minima of L, and not to local maxima or a saddle points.This requires the determinant (B.40) to be strictly positive.Substituting the expression (B.30), we find that det(H) > 0 if The weights p j are positive and sum to unity.Therefore, they can be interpreted as a probability distribution and we obtain b2 j(0) − b j(0) 2 ≥ 0. Thus, the determinant (B.40) is actually strictly positive at the reference point {n j(0) } for any value of σ δ(0) .By continuity on σ δ(0) , this also implies that the Hessian matrix H is always positive definite at this reference point.If the search for the cosmological parameters does not go too far from this reference, the Hessian matrix always remains positive definite.This ensures that the likelihood (B.28) is well behaved.The Ω m − σ 8 constraints from the Log-likelihood of Eq. (B.28) and a simple Poisson Log-Likelihood without sample variance contribution, Fig. 1: The 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.

Fig. 2 :
Fig.2: The 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.

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

Fig. 8 :
Fig.8: Impact of thawing parameters in the scaling relations.XXL ASpiX (base) refers to the results presented in section 4.1 (2 free scaling parameters); XXL-HSC ASpiX contours are the results following the methodology presented in section 5 (6 free, cosmology-dependent, scaling parameters); simple dark blue contours, same as XXL-HSC but the prior 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 ) due to the fact that the scaling relation priors do not depend on the cosmology and then introduce a bias in the analysis.
4 on Ω m and h, for a flat ΛCDM model 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 0.1 step 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] points of the grid to obtain the corresponding mean values and covariance matrices.

Fig. 9 :
Fig. 9: Same 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.

Notes.
Left of the diagonal: the posterior agreement in the Ω m − σ 8 plane.Right of the diagonal (blue cells): the posterior agreement in the 4D Ω m − σ 8 − Ω b − n s plane.
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.
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 theTinker et al. (2010) bias model during this analysis.Then, neglecting finite size effects associated with the borders of the survey volume, we can write the covariance matrix asC i j = δ z i ,z j ni bi n j b j ξi , (B.7)where δ z i ,z j is the Kronecker symbol with respect to the redshift bins i and j, ξi the mean correlation in the redshift bin i,Θ i [M, z]b(M, z).(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 N of clusters in the survey and its continuous counterpart n.For nonoverlapping bins j we haveN = j N j , n = j n j , n = j n j , nb = j n j b j .(B.10)Let us define the fluctuations δ j and δ of the continuous number counts n j and n in the bin j and in the full 2D volume (CR, HR),n j = (1 + δ j )n j , n = (1 + δ)n, (B.11)with means and covariances δ j = 0, δ j δ = b j b ξ, δ = 0, δ j δ = b j b ξthe fluctuations { δ j , δ} are linear functions of each other, and we obtainδ j = b j b δ. (B.14)This is a consequence of the factorisation (B.6) of the cluster correlation function.We also note σ 2 δ the variance of the total number density contrast in the redshift bin, windows and large enough redshift bins, it is possible to simplify the computation of the mean correlation ξ defined in Eq.(B.8).In the flat-sky approximation, for circular windows of angular radius θ s , it reads as ξ = dχ ∆χ dθdθ (πθ 2 s ) 2 28) Appendix B.4: Behaviour of the likelihood L Appendix B.4.1: Cosmology selected by the data Fig. B.1: The X-ray observable diagram (XOD) of the XXL C1 sample, integrated over the redshift range [0.05-1], used in this study, together with the 1D CR, HR distributions in blue.In red, the theoretical diagram predicted by our ΛCDM best fit parameters given in section 4.1.We can see where the CR-HR space is not dominated by the shot noise ( N j > 2-4 objects), we obtain, in average, that n j is smaller than N j .Error bars only account for shot noise.
B.40) which are decreased by the sample-variance term proportional to σ 2 .From Eqs.(B.39)-(B.40),we can estimate the change δλ j of the eigenvalues λ j of the Hessian H i j due to the sample-variance term.With N j ∼ n j , λ j ∼ λ + δλ, we write Tr(H + δH) ∼ Nλ 1 both the trace and the determinant, the comparison with Eqs.(B.39)-(B.40)gives the order-of-magnitude estimate δλ λ ∼ −σ 2 n j , (B.43) when σ δ(0) is small.Moreover, at the reference point, with N j = n j = n j(0) , this reads Fig. B.2: Ω m − σ 8 contours using the log-likelihood from equation B.28 (Base Log-Likelihood) and the Poisson log-likelihood from equation B.47 (Poisson Likelihood).The constraints when using the base log-likelihood compared to the Poisson one are increased by a factor of 10%.

Table 1 :
Sampling of the X-ray parameter distribution in the XOD.

Table 2 :
Cluster scaling laws used in the first part of the paper.

Table 3 :
The values of the {a x }, {b x } and {c x } coefficients from equation 6.
Notes.The functions, %err x (CR, θ c ) = a x CR bx θ cx c , are fitted with the data using the Levenberg-Marquardt algorithm.

Table 5 :
BAO data used in this analysis.D are not shown here since the three WiggleZ measurements are correlated and we take the full covariance matrix in the analysis.

Table 6 :
Umetsu et al. (2020) used in section 5.Notes.We disperse over the luminosity, temperature and core radius distributions in this case.Uncertainties on parameters indicate that these parameters are left free during the analysis.The values shown in this table are calculated assuming Ω m = 0.28 and h = 0.7.During the analysis, the parameter means and covariances are rescaled as a function of cosmology as described in section 5. with log 10 (M 500,WL ) =log 10 (M 500,T rue )+ log 10 (1 + b WL ) ± σ log 10 M WL (16)and we assume a Gaussian prior on log 10 (1 + b WL ) of log 10 (1 + b WL ) = log 10 (1 − 0.11) ± 5%/ln10 to marginalise over the mass calibration uncertainty of ±5%, seeUmetsu et al. (2020).Here σ log 10 M WL in Equation

Table 7 :
ASpiX cosmological constraints for the HSC-XXL ASpiX model and the joint analysis (flat ΛCDM).

Table 8 :
Posterior agreement between the various cases studied in Sec. 6.