CoMaLit-IV. Evolution and self-similarity of scaling relations with the galaxy cluster mass

The scaling of observable properties of galaxy clusters with mass evolves with time. Assessing the role of the evolution is crucial to study the formation and evolution of massive halos and to avoid biases in the calibration. We present a general method to infer the mass and the redshift dependence, and the time-evolving intrinsic scatter of the mass-observable relations. The procedure self-calibrates the redshift dependent completeness function of the sample. The intrinsic scatter in the mass estimates used to calibrate the relation is considered too. We apply the method to the scaling of mass M_Delta versus line of sight galaxy velocity dispersion sigma_v, optical richness, X-ray luminosity, L_X, and Sunyaev-Zel'dovich signal. Masses were calibrated with weak lensing measurements. The measured relations are in good agreement with time and mass dependencies predicted in the self-similar scenario of structure formation. The lone exception is the L_X-M_Delta relation whose time evolution is negative in agreement with formation scenarios with additional radiative cooling and uniform preheating at high redshift. The intrinsic scatter in the sigma_v-M_Delta relation is notably small, of order of 14 per cent. Robust predictions on the observed properties of the galaxy clusters in the CLASH sample are provided as cases of study. Catalogs and scripts are publicly available at http://pico.bo.astro.it/~sereno/CoMaLit/.


INTRODUCTION
Scaling relations among cluster global properties embody important clues on the formation and evolution of cosmic structures. They result from the main gravitational processes driving the cluster evolution (Kaiser 1986;Giodini et al. 2013;Battaglia et al. 2012;Ettori 2013). Accurate mass-observable relations are also needed to use the abundance of galaxy clusters to constrain cosmological parameters (Vikhlinin et al. 2009;Mantz et al. 2010bMantz et al. , 2014Rozo et al. 2010;Planck Collaboration et al. 2014b).
This paper is the fourth in a series titled 'CoMaLit' (COmparing MAsses in LITerature), which aims to assess our present capability to measure cluster masses, and to develop methods to measure scaling relations through Bayesian techniques. In the first paper , CoMaLit-I), we evaluated systematic differences in lensing and X-ray masses obtained from independent analyses and we quantified the overall level of bias and intrinsic scatter of these mass proxies. The second paper presented the formalism to calibrate an observable cluster property against the cluster mass and applied the methodology to the Sunyaev-Zel'dovich If gravity is the dominant process, the resulting self-similar model predicts scaling relations in form of scale-free power laws (Kaiser 1986;Giodini et al. 2013). Numerical simulations (Stanek et al. 2010;Fabjan et al. 2011) confirmed these scalings and showed that intrinsic scatters around the median relations approximately follow a log-normal distribution. This basic theoretical scheme is very successful in describing observed scaling relations in Xray and SZ (Ettori 2013(Ettori , 2015. Deviations from the self-similar scheme may indicate that non-gravitational processes, such as feedback and non-thermal processes, contribute significantly to the global energy budget in clusters (Maughan et al. 2012).
The precise measurement of the redshift evolution of the scaling relations is then crucial to understand how either gravitational or non-gravitational phenomena drive the formation and evolu-tion of clusters. Furthermore, if the time evolution is neglected or wrongly shaped, the estimated scaling with mass may be biased in cluster samples spanning a significant redshift range (Andreon 2014).
Any real time-dependence in the scaling of observables with mass and redshift has to be separated by other effects connected either to the evolution of the cluster mass function or to the redshift dependence of the selection function and of the completeness of the sample. Further complications are due to the fact that usual samples of clusters are not assembled according to well defined criteria but they may be just heterogeneous collections of systems with high quality data. In this case, the determination of the selection function is very problematic.
Here we develop a method that measures at the same time the evolution of the scaling relation and the completeness/selection function of the sample. This is a self-calibrating method which is intended to be optimised in large optical survey such as Euclid (Laureijs et al. 2011). If we calibrate a cluster observable against the cluster mass, this relation can be used to construct a mass proxy based on the observable. The optimal mass proxy is expected to be easy to measure, unbiased, and minimally scattered. A crucial aspect is that in the first step we cannot calibrate the observable against the true mass (which cannot be measured), but we have to rely on another mass proxy, such as the WL mass or the hydrostatic mass, which are scattered too (Rasia et al. 2012;CoMaLit-I). This scatter has to be considered to avoid biases (Mantz et al. 2014;CoMaLit-I;CoMaLit-II) We apply the method to calibrate mass proxies based on the line of sight velocity dispersion of cluster galaxies, which is supposedly the best mass proxy (Stanek et al. 2010;Saro et al. 2013), and three other observables, which are more scattered but optimised for large surveys, i.e., the optical richness, the X-ray luminosity, and the Sunyaev-Zel'dovich integrated Compton parameter.
The paper is organised as follows. Section 2 is devoted to general considerations on the redshift evolution of the scaling relations, of the intrinsic scatters, and of the selection/completeness function. Section 3 reviews the methodology employed to perform the regression and to recover at the same time the scaling relations and the completeness and selection functions. The cluster catalogs used in the analysis are introduced in Sec. 4. Results are presented in Sec. 5 whereas Sec. 6 is devoted to the comparison with theoretical predictions and previous works. Final considerations are contained in the Sec. 7. In App. A, we discuss how the masses of clusters in selected samples are usually distributed. Appendix B describe the format of the compiled catalogs of line-of-sight velocity dispersions.
Throughout the CoMaLit series of papers, we have been adopting the following conventions and notations. The framework cosmological model is the concordance flat ΛCDM universe with density parameter ΩM = 0.3, and Hubble constant H0 = 70 km s −1 Mpc −1 . H(z) is the redshift dependent Hubble parameter and Ez ≡ H(z)/H0. When H0 is not specified, h is the Hubble constant in units of 100 km s −1 Mpc −1 .
O∆ denotes a global property of the cluster measured within the radius which encloses a mean over-density of ∆ times the critical density at the cluster redshift, ρcr = 3H(z) 2 /(8πG). 'log' is the logarithm to base 10 and 'ln' is the natural logarithm.

REDSHIFT EVOLUTION
Scaling relations evolve with redshift. Numerical simulations (Stanek et al. 2010) and theoretical predictions (Giodini et al. 2013) agree that the relation between the mass M∆ and any observable quantity O can be summarised by the form Within this framework, the redshift evolution in the median scaling relation is accounted for by the factor E γ z whereas the slope β is redshift independent. In fact, in the self-similar scenario, the evolution does not depend on the mass scale and only the normalisation depends on cosmic time.
We tested this scheme in a number of cases. We considered observables connected either to the galaxy distribution, i.e., velocity dispersion of galaxies along the line of sight, σv, or optical richness, λ, or to the intra-cluster medium, i.e., the bolometric X-ray luminosity, LX, or the spherically integrated SZ Compton signal, YSZ. In the self-similar scenario, we expect that for clusters in equilibrium the scalings go like (Giodini et al. 2013;Ettori 2015) σv where DA is the angular diameter distance to the cluster. The above self-similar scaling relations evolve with redshift as E γss z , with γss = 1/3, 0, 7/3, or 2/3 for the galaxy velocity dispersion, the optical richness, the X-ray luminosity, or the spherical SZ signal, respectively. The scaling of the X-ray luminosity depends on the energy band. For the soft X-ray luminosity in the rest-frame energy band [0.1 − 2.4] keV, LX soft , the evolution can be expressed as (Ettori 2015), Together with the redshift dependence of the median relation, the intrinsic scatter of the relation and the scatter between the true mass and the mass proxy used to calibrate the relation may evolve as well. Furthermore, any apparent redshift evolution of the scaling may be degenerate with the evolution of either the mass or the selection function. We discuss these aspects in the following.

Intrinsic scatter
Broadly speaking, the intrinsic scatter of a scaling relation is related to the degree of regularity of the clusters. The larger the deviations from dynamical/hydrostatic equilibrium the larger the scatter (Fabjan et al. 2011;Saro et al. 2013). Scatter is then prominent in morphologically complex halos. Triaxiality is another major source of scatter, since clusters are usually studied under the simplifying assumption of spherical symmetry (Limousin et al. 2013;Sereno et al. 2013). Since high redshift clusters are more irregular and less spherical, the scatter is usually expected to increase with redshift.
Let us consider the evolution of scatter in a number of cases. Based on numerical simulations, Saro et al. (2013) showed that the scatter of dynamical mass estimates based on the line-of-sight velocity dispersion is approximately log-normal and that it increases with redshift as They argued that the dominant contributor to the scatter is the intrinsic triaxial structure of halos and that its evolution with redshift is also the dominant source of the increasing scatter of the 1D dynamical mass estimates with redshift. Fabjan et al. (2011) studied the scaling relations between the cluster mass and some proxies based on X-ray quantities with a set of cosmological hydrodynamical simulations. They found that the scatter distribution around the best-fitting relations is always close to a log-normal one and that the scatter increases with redshift.
The precise quantitative estimate of the scatter and of its evolution strongly depends on the details of the baryonic physics included in the simulations. We considered the results of Fabjan et al. (2011) for runs with non-radiative physics and standard viscosity. To study the time evolution, we fitted the values of the scatter obtained at different redshifts (z = 0.0, 0.25, 0.50, 0.80, and 1.0) under the assumption of self-similar scaling relation (Fabjan et al. 2011, tables 1, 2, and 3). The mass proxy MY X is based on YX, i.e., the product of the gas mass within r500 and the spectroscopic temperature outside the core (Kravtsov, Vikhlinin & Nagai 2006). We found that the scatter evolves as or, with an alternative form 1 , For the mass proxy based on the emission-weighted temperature, we found or The above results show that the scatter mildly increases with redshift. This suggest that the evolution of the scatter can be modelled as

Completeness
The completeness of a sample usually evolves with redshift. Very massive clusters are rare and difficult to be found in the local volume but they are still forming at high redshift. On the other hand, only clusters emitting very strong signals can be detected to very large distances. As detailed in App. A, the selection and the mass functions conjure to make the distribution of true masses in observed samples fairly unimodal. The evolution of the completeness of the sample can be characterised through the evolution of the peak and of the dispersion of this distribution. The mean (logarithmic) mass of the sample is connected to the observational threshold (see App. A), which may evolve with redshift, and to the scatter between the mass and the observable quantity used to select the clusters, which evolves too.
Let us first consider the evolution of the mass corresponding to a completeness limit. As a first example let us consider a fluxselected sample. The luminosity scales with mass as If we select only clusters above a limiting flux, f th , the corresponding luminosity evolves as L th (z) ∝ f th DL(z) 2 , where DL(z) is the luminosity distance. In absence of scatter, the corresponding limiting mass evolves as As a second example, let us consider a SZ-like signal, whose size increases with the projected physical surface covered by the cluster. In this case, the observable is proportional to DA(z) 2 θ 2 ∆ , where θ∆ is the angular extension of the cluster. The scaling can be written as The noise is proportional to the square root of the angular area, i.e., σY ∆ ∝ θ∆ = r∆/DA. In absence of scatter, if we select only clusters above a given signal-to-noise ratio (SNR), i.e., Y∆/σY ∆ > SNR th , the corresponding threshold mass evolves with redshift as The two above examples suggest that the evolution of the mass at a given completeness limit can be parametrized as where the factor E γ Ez z accounts for the evolution in both the mass threshold and the intrinsic scatter of the scaling relation. This modelling of the completeness limit was derived for samples selected with a cut on the detection observable but the functional form is flexible enough to address even more complicated cases. The choice of the angular diameter distance over the luminosity distance in Eq. (17) is irrelevant since the two differs for a factor (1 + z) 2 which can be approximately englobed in E γ Ez z for limited redshift baselines.
The evolution in the dispersion of the mass sample is mainly connected to the intrinsic scatter in the relation used to select the samples, see App. A. The redshift dependence can then be modelled as in Eq. (12).

REGRESSION SCHEME
When we calibrate a scaling relation we deal with: i) the true mass of the cluster M∆, which we cannot measure; ii) a scattered (and likely biased) proxy of the true mass, such as the weak lensing mass MWL,∆, which is the proxy we considered in following, or the hydrostatic mass MHE,∆ (see § 2 and app. A of CoMaLit-I); iii) an observable quantity O, which we assume to be on average related to the true mass with a power-law.
In logarithmic variables, the median scaling relation is approximatively linear and the scatter is Gaussian. As discussed in Sec. 2, the scaling can be expressed as Since we englobed the self-similar evolution in the left-hand of Eq. (18), values of the parameters γz which are different from zero denote deviations from the self-similar time dependence. In other words, the time evolution of the scaling relations γz is relative to that predicted by the self-similar model. Given a particular scaling law, there is negative, i.e., γz < 0 (positive, i.e., γz > 0) evolution if the normalisation at high redshift is lower (higher) than anticipated from the self-similar scaling.
In what follows, which is the general scheme we employed for Table 1. Parameters of the regression scheme and their description. The variables Z, X, and Y denote (the logarithm of) the true mass, the WL mass, and the self-similarly evolved observable, respectively.

Type
Meaning Symbol Intrinsic scatter of the WL mass Conditional scatter at z = 0 σ X|Z,0 Time evolution Mean of the mass function Dispersion of the mass function Dispersion at z = 0 σ Z,0 Time evolution γσ Z the regression analysis, we identify log M∆ with the variable Z, we identify the logarithm of the mass proxy, i.e., the weak lensing mass, with the random variable X, and we identify the logarithm of the self-similarly redshift evolved observable with the response Y . In this scheme, the mass is the covariate variable, as when using number counts of galaxy clusters to constrain cosmological parameters. The observed values are denoted with the lower case, i.e., x and y are the manifest measured estimates of the latent X and Y , respectively (Feigelson & Babu 2012). This notation is the convention adopted in the CoMaLit series. The conditional probability of X given Z is where N is the normal distribution. In Eq. (19), X is an unbiased proxy of Z. Any bias between X and Z would be degenerate with the estimated overall normalisation of the scaling between Y and Z. This bias cannot be determined with the data, which only constrain the relative bias between X and Y (see CoMaLit-I). As discussed in Section 2.1, the redshift evolution of the scatter is modelled as The mean observable for a given mass is linearly related to the (logarithm of the) mass and the relation evolves with redshift, the redshift z is deterministic and assumed to be known without measurement errors. Y is scattered and distributed according to the conditional probability with The distribution of the masses can be approximated with a Gaussian function P (Z) = N (µZ (z), σZ (z)).
The mass distribution resulting from usual selection procedures is fairly unimodal (see App. A) and can be approximated with the normal distribution of Eq. (24). The statistical improvement obtained considering more complex distributions, such as mixture of Gaussians with different means and variances, is usually negligible (Kelly 2007;CoMaLit-II). As discussed in Sec. 2.2, the evolution of the (mean of the) mass function can be modelled after Eq. (17) as The dispersion evolves as The completeness function at a given redshift can be computed by dividing the estimated mass function (Eq. 24) by the cosmological halo mass function. This approach requires the knowledge of the effective survey area of the sample, which may be difficult to estimate for heterogeneous samples. Alternatively, we can use the approximate formulae presented in App. A, which were derived under the assumptions that the completeness function can be approximated as a complementary error function and that the cosmological halo mass function can be approximated as a power-law.
If we assume that the uncertainty in the measurement process is Gaussian, the relation between the unknown Xi and Yi and the measured xi and yi is given by where N 2D and U are the bivariate Gaussian and the uniform distribution, respectively. In Eq. (27), Vi is the symmetric uncertainty covariance matrix whose diagonal elements are denoted as δ 2 x,i and δ 2 y,i , and whose off-diagonal elements are denoted as ρxyδx,iδy,i. The truncation, i.e., null probability for yi < y th,i , accounts for selection effects when only clusters above an observational limit (in the response variable) are included in the sample, i.e., the Malmquist bias (CoMaLit-II).
The treatment is complete once the priors on the parameters are made explicit. We choose non-informative priors as discussed in CoMaLit-I and CoMaLit-II. The priors on the intercept α Y |Z and on the meanμZ are taken to be flat, where is a small number. In our calculation we took = 10 −3 . A priori, the slopes follow the Student's t1 distribution with one degree of freedom, as suitable for uniformly distributed direction angles, The Student prior for the slopes is not informative. Negative time evolutions and scatters which decreases at early times are allowed.
For the variances, we adopted an inverse Gamma distribution, This regression scheme requires 12 parameters, i.e., three parameters characterising the scaling relation, two for the intrinsic scatter, two for the mass scatter, and five for the mass function, plus three variables for each cluster, i.e., the true weak lensing mass, the true mass, and the true observable. The parameters and their meanings are summarised in Table 1.
The relation in Eq. (21) expresses the conditional scaling relation, wherein YZ is the most likely value of the variable Y for a given Z. This is the relation to be used to predict the value of Y for a given Z. The relation between two random and scattered variables might be better described by the symmetric scaling relation, which goes along the direction where the probabilities of Z and Y are maximised at the same time (CoMaLit-II). In the above regression scheme, where Z and Y follow a bivariate normal distribution, the slope of the symmetric scaling can be expressed as (CoMaLit-II) where ρY Z is the correlation factor between Y and Z, and the variance in Y is related to the conditional scatter as The intercept of the symmetric relation can be expressed as The detailed regression scheme is simplified when we are interested in the scaling between the observable and the measured proxy mass, i.e., the WL or the X-ray mass, In this case, the adopted form for the scaling is which substitutes Eq. (21). The latent variable Z coincides now with the manifest one X and we do not have to model the conditional probability of X given Z, see Eqs. (19) and (20).

CLUSTER CATALOGS
There are different approaches to choose a sample of clusters to analyse. We may look for a statistical sample which is complete with respect to well-defined selection criteria. This sample would be ideal but most of the massive clusters with very good quality data might be excised. The alternative is to assemble samples as numerous as possible with the idea that variety and largeness can compensate for incompleteness and inhomogeneity. These two approaches are to some degree complementary and have been already discussed in CoMaLit-II and CoMaLit-III, which we refer to for further considerations. Here we are mainly interested in testing the regression algorithm and we focus on large samples. To this aim we assembled a catalog of clusters with measured velocity dispersions. As catalogs of weak lensing masses, optical richnesses, X-ray luminosities, and SZ effects, we used publicly available compilations. The subsample of weak lensing clusters also included in either the velocity dispersion, richness, X-ray, or SZ catalogs were used in the analysis presented in the next section. We briefly discuss the main properties of the catalogs and refer to the original references for further details.

Weak lensing masses
CoMaLit-III retrieved from literature 822 weak lensing analyses of clusters and groups with measured redshift and mass. Here, we consider the LC 2 -single, which contains 485 unique entries with reported coordinates, redshift, and WL masses to over-densities of 2500, 500, 200, and to the virial radius. 2 Duplicate entries from input references were carefully handled.
The cluster redshifts span a large interval, 0.02 < ∼ z < ∼ 1.5. The catalog is large and standardised but it is not statistically complete. We refer to CoMaLit-III for a detailed discussion of the catalog properties.

Velocity dispersions
We assembled some publicly available catalogs of clusters with measured velocity dispersions. We first review the source catalogs and then we introduce the merged compilation. Cava et al. (2009) presented the results from the spectroscopic survey WINGS (WIde-field Nearby Galaxy-cluster Survey)-SPE, which consists of 48 nearby clusters at 0.04 < ∼ z < ∼ 0.07 selected from three X-ray flux limited samples. They complemented the sample with 29 additional clusters not observed in the programme but for which literature data existed. The total sample contains 77 clusters over a broad range of richness, Bautz-Morgan class, and X-ray luminosity. Ebeling et al. (2007) presented the sample of the 12 most distant galaxy clusters detected at z > ∼ 0.5 by the Massive Cluster Survey (MACS). This catalog is statistically complete and comprehensive of measurements of radial velocity dispersions. Girardi & Mezzetti (2001) considered a sample of 51 distant galaxy clusters at 0.15 < ∼ z < ∼ 0.9, each cluster having at least 10 galaxies with available redshift in the literature. In some clusters, two peaks that are not clearly separable were identified in the velocity distribution. For these systems with uncertain internal dynamics, we considered the velocity dispersion measured by analysing the identified peaks together. We also discarded two systems with no major peak (CL J0023+0423 and CL J0949+44). Mazure et al. (1996) constructed a volume limited sample of 128 clusters out to z = 0.1 combining data from the ENACS (ESO Nearby Abell Clusters Survey) with pre-existing data from literature. They measured reliable velocity dispersions for a subset of 80 of them, based on at least 10 redshifts. They also analysed 26 additional clusters in the cone but with z > 0.1. The total catalog consists of 106 clusters. We discarded from our final catalog the secondary systems. Oegerle & Hill (2001) presented the spectroscopic study of a sample of 25 Abell clusters out to z = 0.1 containing a central cD galaxy. Redshifts measured with the MX Spectrometer were combined with those collected from the literature to obtain typically 50-150 observed velocities in each cluster. We used the estimates of the velocity dispersions within the smaller quoted aperture (∼ 1.2 Mpc/h at z = 0.1) Popesso et al. (2007) considered a sample of 137 optically selected and spectroscopically confirmed Abell clusters in the SDSS (Sloan Digital Sky Survey) database (Adelman-McCarthy et al. 2006). The clusters span the redshift range 0.04 < ∼ z < ∼ 0.17. 40 of the clusters were X-ray under-luminous, since they had a marginal X-ray detection or remained undetected in the ROSAT All Sky Survey. Rines & Diaferio (2006) studied the infall patterns of 72 nearby (z < 0.1) clusters from the Data Release (DR) 4 of the SDSS. The clusters were selected in X-ray flux from the ROSAT All-Sky Survey. Velocity dispersions were measured and masses were derived with the caustic method. Rines & Diaferio (2010) extended the approach to a sample of 16 groups with lower X-ray fluxes selected from the 400 deg 2 serendipitous survey of clusters. Spectroscopic data were taken from the SDSS DR5. Rines et al. (2013) selected 58 clusters by their X-ray flux and in the redshift interval 0.1 < z < 0.3 to build the Hectospec Cluster Survey (HeCS), the first systematic spectroscopic survey of cluster infall regions at z 0.1. For each cluster, high signal-to-noise spectra for ∼200 cluster members were acquired with MMT (Multi Mirror Telescope)/Hectospec. Ruel et al. (2014) presented optical spectroscopy of galaxies in clusters detected through the SZ effect with the South Pole Telescope (SPT). They reported measurements of 61 spectroscopic cluster redshifts, and 48 velocity dispersions each calculated with more than 15 member galaxies. After the inclusion of additional measurements of SPT-observed clusters previously reported in the literature, the final catalog presents 57 velocity dispersions. Being SZ selected, most of the clusters are at high redshift. The clusters span an interval 0.3 < ∼ z < ∼ 1.5. Sifón et al. (2013) presented the dynamical analysis of a sample of 16 SZ selected massive clusters detected with the Atacama Cosmology Telescope (ACT) over a 455 deg 2 area of the southern sky. 60 member galaxies on average per cluster were observed with deep multi-object spectroscopic observations. The sample spans the redshift range 0.3 < ∼ z < ∼ 1.1 with a median redshift z = 0.50. Zhang et al. (2011) presented a multi-wavelength analysis of 62 galaxy clusters in the HIFLUGCS (HIghest X-ray FLUx Galaxy Cluster Sample), an X-ray flux-limited sample. Velocity dispersions were computed thanks to 13439 cluster member galaxies with redshifts collected from literature. Most of the clusters (60 out of 62) are at z < 0.1.

Merged catalog
The catalogs listed before provides a total of 710 velocity dispersion estimates, comprehensive of multiple peaks and substructures which we did not consider in our final sample.
Cluster coordinates were taken from the original or from companion papers. When they were not reported, we used the coordinates listed by the NASA/IPAC Extragalactic Database (NED). 3 Not unique entries were identified by matching names and cluster coordinates. For clusters with multiple analysis, we preferred the study based on the larger number of identified cluster member galaxies with measured redshift, N members . The final catalog contains 564 unique clusters. The catalogs are publicly available at http://pico.bo.astro.it/~sereno/ CoMaLit/sigma/. Their format is detailed in App. B.
When original estimates were provided with asymmetric errors, we computed the mean value as suggested in D' Agostini (2004). To standardise the uncertainties, we followed Ruel et al. (2014). They found that the uncertainty in the velocity dispersion σv is well described by 3 http://ned.ipac.caltech.edu/.
when including the effect of membership selection. As we inferred from the matching with the LC 2 -single, WL masses are known for a subsample of 97 clusters. This size can be achieved only relying on a number of different source catalogs. 30 clusters are from Girardi & Mezzetti (2001), 23 from Rines et al. (2013), 13 from Zhang et al. (2011), 11 from Ebeling et al. (2007), 8 from Ruel et al. (2014), 5 from Popesso et al. (2007), 4 from Mazure et al. (1996), and 3 from Sifón et al. (2013). Rykoff et al. (2014) applied the redMaPPer (red-sequence Matched-filter Probabilistic Percolation), a red-sequence cluster finder designed for large photometric surveys, to ∼ 10000 deg 2 of SDSS DR8 data. The resulting catalog 4 contains ∼ 25000 candidate clusters over the redshift range 0.08 < ∼ z < ∼ 0.55. According to the catalog convention, the richness λ of a cluster is defined as the sum of the probabilities of the galaxies found near a cluster to be actually cluster members. The sum extends over all galaxies above a cut-off luminosity (0.2L * ) and below a radial cut which scales with richness. Clusters are included in the catalog if their richness exceeds 20 times the scale factor SRM (also provided in the catalog), which is a function of the photometric redshift of the cluster. This selection criterion approximately requires that every cluster has at least 20 galaxy counts above the flux limit of the survey or the cut-off luminosity at the cluster redshift, whichever is higher. 3) baselines. The sample was assembled from all publicly available Chandra data as of 2006 November. It consists of clusters at redshift greater than 0.1 listed in the NED which were the targets of observations made with the ACIS-I detector covering at least half of the area in the annulus [0.9-1.0] r500. The radius r500 was estimated assuming a M500-YX relation.

X-ray clusters
The sample was later reanalysed using updated softwares and calibration files in Maughan et al. (2012). Bolometric luminosities were measured either in the [0-1] r500 aperture, which we took as reference case to ease the comparison with theoretical predictions, or in the core-excised [0.15-1] r500 aperture, LX,ce. This large catalog is standardised in the measurement procedures but it is not statistically complete.
As an alternative we also looked at catalogs of X-ray luminosities measured in the [0.1-2.4] keV band, LX soft . We considered the MCXC (Meta-Catalogue of X-ray detected Clusters of galaxies, Piffaretti et al. 2011), which comprises 1743 unique Xray clusters collected from available ROSAT All Sky Survey-based and serendipitous cluster catalogues. X-ray luminosities were systematically homogenised and standardised to an over-density of ∆ = 500. Uncertainties are not provided in the catalog. For our tests, we fixed the statistical uncertainty to 10 per cent. As the LC 2 , the MCXC is not statistically complete.

Planck SZ catalog
The Planck SZ Catalogue (PSZ, Planck Collaboration et al. 2014a) contains 883 candidates identified with the Matched Multi-filter method MMF3 with detections above SNR = 4.5. The catalog spans a broad mass range from 0.1 to 16 × 10 14 M at a median redshift of z ∼ 0.22. The redshift determination is available for 664 candidates.
In CoMaLit-II, we computed the spherically integrated Y500 of the PSZ clusters within the weak-lensing determined r500. The measurements of MWL,500 and Y500 are then correlated. In our analysis, we used the full uncertainty covariance matrix. We refer to CoMaLit-II for a detailed discussion.

RESULTS
We analysed the scaling between mass and optical, X-ray, or SZ observables using the general regression scheme detailed in Section 3. In fact, the uncertainties on the redshifts were assumed to be negligible and the factors Ez were assumed to be known without errors. The Bayesian hierarchical model was implemented through JAGS. 5 According to the notation of Section 3, the X variable is the logarithm (to base 10) of the observed weak lensing mass (computed at an over-density of either 200 in case of scaling of optical observables or 500 for observables related to the gas); the Z variable is the logarithm of the unknown true mass; the observable Y is the logarithm of either the galaxy velocity dispersion, the optical richness, the X-ray luminosity, or the spherical SZ signal multiplied by E −γss z . Any deviation of the parameter γz from the null value implies a deviation from the self-similar evolution with redshift.
In the case of optical richness and SZ signal, we had to consider the correction for the Malmquist bias. The threshold value of the optical richness above which candidate clusters are included in the redMaPPer catalog was given by 20 times the scale factor at the cluster redshift (Rykoff et al. 2014). According to the notation of Section 4.3, the threshold for the i-th cluster is The limiting SZ flux of the Planck clusters was obtained by multiplying the minimum SNR(= 4.5) by the uncertainty on Y500 (CoMaLit-II). In this case, The results of the regression are summarised in Table 2 for the scaling relations and in Table 3 for the mass functions. Table 4 summarises the results of the scaling of the observables versus the WL mass. Parameter degeneracies are illustrated as bi-dimensional contour plots in Figs. 1, 2, 3, and 4. Figures 5, 6, 7, and 8 show the scaling relation and the evolution of the completeness function for σv-M200, λ-M200, LX-M500, and YSZ-M500, respectively.
To ease the comparison with theoretical predictions, we also computed the parameters of the the symmetric scaling relation (see Table 2). We did not require that βY -Z is redshift independent. However, the slopes turned out to be constant within the errors.
Slopes and intercepts of the symmetric relations in Tables 2 and 4 were computed at the median redshifts of the samples.
We obtained significant constraints on the evolution with redshift of the scaling relations and of the mass functions. On the other hand, the uncertainties on the evolution of the intrinsic scatters are too large to come to any conclusion.

Parameter degeneracy
Most of the regression parameters are uncorrelated (see Figs. 1-4). Since a significant percentage of the massive clusters is at high redshift, the time evolution can partially mimic the effects of the mass evolution (see the β Y |Z -γz panel). This degeneracy is most pronounced for the Planck selected clusters, see Fig. 4, whose mass completeness limits steadily increases with redshift, see Fig. 8. To obtain unbiased estimates of the scaling relations is then crucial to properly account for time-evolution and selection effects.
The parameters characterising the scaling, i.e., α Y |Z , β Y |Z , and γz, are not degenerate with the mass functions. Furthermore, as already noted in CoMaLit-I, since σ X|Z and σ Y |Z spread the observed relation in orthogonal directions, they are nearly uncorrelated too.
The only remaining significant degeneracy is among the parameters characterising the normalisation of the mean value of the mass function,μZ , and its evolution, parameterised in terms of γµ Z ,D and γµ Z . The evolution parameters of the (mean of the) mass function, γµ Z ,D and γµ Z , are degenerate too, see Eqs. (17) and (25). In fact, the redshift dependence in small intervals can be modelled either in terms of Ez or in terms of the distance.

Scaling with the weak-lensing mass
Results for the scalings with the WL mass are reported in Table 4. In this case, the general regression scheme simplifies as described in Sec. 3. Due to the intrinsic scatter of the WL mass with respect to the true mass, the conditional relation O-MWL,∆ is systematically flatter than the corresponding O-M∆, whereas the intrinsic scatter of the relation is larger (CoMaLit-I; CoMaLit-II). On the other hand, results for the symmetric scaling relations are consistent.
Parameters of the O-MWL,∆ are recovered with better accuracy but they cannot be used straight on in the conditional O-M∆ to make predictions based on cosmological functions. In absence of a full regression procedure modelling the intrinsic scatter between true mass and proxy mass, at least some corrections should be used (CoMaLit-I).
We tested that these general features of the LX-M∆ relation do not depend on the different ways in which the X-ray luminosity can be measured. Bolometric luminosities estimated in a region from which the core is excised are supposed to be slightly less sensitive to the baryonic physics. Considering the LX,ce from Maughan et al. (2012), we found βY -Z = 1.47 ± 0.18, γz = −1.39 ± 0.57, and σ Y |Z = 0.09 ± 0.05. As expected, the LX,ce-M∆ relation is more tight, with a smaller intrinsic dispersion, and it can be determined with a better statistical accuracy. Neverthe-less, the evolution in redshift is still negative and fully compatible with the result obtained for luminosities accounting for the core regions. The dependence with mass is slightly less pronounced but still steeper than the self-similar expectation.
We also evaluated the mass-X-ray luminosity relation in the soft band by considering the MCXC catalog, which comprises 193 clusters with measured WL mass. For the LX soft -M∆ relation we found βY -Z = 1.43 ± 0.15, which is steeper than the self-similar prediction βss = 1, and a negative redshift evolution which deviates from the expected γss = 2 (see e.g. Ettori (2015)   γz = −1.80 ± 0.59. As expected for the MCXC catalog, whose luminosities were standardised from heterogeneous data sets, the intrinsic scatter is over-estimated, σ Y |Z = 0.18 ± 0.07.

Intrinsic scatters
The intrinsic scatter of the weak lensing mass with respect to the true mass is estimated to be of the order of ∼ 25 ± 10 per cent. This is slightly larger but still compatible within errors with the scatter measured in Mantz et al. (2014)  CoMaLit-II, we noted that the scatter measured in heterogeneous samples, such as LC 2 -single, may be overestimated due to not coherent formulated statistical uncertainties on weak lensing masses.
In the case of the richness calibration, both estimated weak lensing scatter and related errors doubles, so that the difference in the central estimates of the scatter is not statistically significant. We confirm that mass proxies based on velocity dispersions are among the most accurate (scatter of ∼14 per cent). We verified that mass proxies based on X-ray luminosity are noisier (scatter of ∼30 per cent) than those based on the integrated SZ effect (∼25 per cent), in agreement with numerical simulations (Stanek et al.  2010). The richness, with an intrinsic scatter of ∼40 per cent, may compete with the X-ray luminosity. The scatter in proxies based on the X-ray luminosity strongly depend on the measurement/analysis process. Luminosities estimated in core-excised regions are less scattered (∼20 per cent), whereas inconsistencies in the measurement process and not uniform data-sets can boost the scatter up to ∼ 40 per cent. The catalog of velocity dispersion used to study the σv-M200 relation is heterogeneous which might bias the relation and inflate the estimated intrinsic scatter. Firstly, the central estimates of the velocity dispersion were estimated in the different source papers with different methods. However, these statistical differences are very small. Ruel et al. (2014) presented two independent measurements of σv, based either on the measured gapper scale or the biweight dispersion. The two estimates are very well consistent, with a distribution of relative differences whose centre deviates from zero by less than one per cent and whose scatter is less than 2 per cent.

Completeness function
The evolution of the completeness function is obtained through the analysis of the redshift dependence of the distribution of selected true masses. In Figs. 5, 6, 7, and 8, we plotted the 15, 50, and 85 per cent completeness limit as a function of the redshift. The completeness was approximated as a complementary error function, see App. A, whose scale and dispersion were derived from the parameters of the fitted mass function through Eqs. (A9) and (A10).
For the samples of clusters with weak lensing mass and velocity dispersion, SZ flux, or X-ray luminosity, the larger the redshift the larger the mass at a given completeness limit, see Figs. 5, 7, and 8. This is not the case for the WL clusters in the redMaPPer catalog, when the limits are nearly redshift independent, see Fig. 6. The apparent spike in some completeness functions at high redshift, see Table 2. Observed scaling relations, E −γss z O = 10 α M β ∆ E γz z . Conventions are as in Section 3 and Table 1: . γss is the self-similar evolution, which is equal to 1/3, 0, 7/3, 2/3 for σv-, λ-, L X -, and Y SZ -M ∆ , respectively. Units are 10 14 M for mass, km/s for σv, 10 44 ergs s −1 for the bolometric luminosity L X , and 10 −4 Mpc 2 for D 2 A Y 500 . Cols. 1-3: variables of the regression procedure. Col. 4: number of clusters in the sample (N cl ). Col. 5: median redshift of the sample. Cols. 6, 7, and 8: intercept, mass slope, and time evolution of the conditional scaling relation. Cols. 9-10: local scatter of the WL mass and its time-evolution.  intrinsic scatter of the scaling relation and its time-evolution.  intercept and slope of the symmetric scaling relation at the median redshift of the sample. Col. 15: self-similar prediction for the slope (βss).
Conditional scaling WL mass scatter Intrinsic scatter Symmetric scaling  Few high redshift clusters lie above the 85 per cent completeness limit. This is expected since high mass clusters are very rare at high redshifts and few of them have measured WL mass.
Even though the redMaPPer catalog of optical richnesses (Rykoff et al. 2014) and the Planck catalogs (Planck Collaboration et al. 2014a) are statistically complete, we do not expect to exactly recover the selection function of the full catalogs. In fact, the subsamples of clusters with measured weak lensing masses may be biased with respect to the full catalogs.
The derived completeness functions refer to the subsample of the cluster in the catalogs with measured WL mass. The catalog of WL masses is heterogeneous which makes the studied subsamples heterogeneous too, even in the case of parent catalogs constructed with well defined selection functions. However, despite the fact that we do not expect that the completeness function of the subsamples with WL masses strictly resembles the completeness of the parent catalog, some similarities are still in place. The redshift evolution of the derived completeness limits of the SZ flux-WL catalog follows that of the Planck clusters, see the upper panel of Fig. 8. The shift in normalisation reflects the fact that masses used in Planck Collaboration et al. (2014a) to derive the average limit were based on the Yz proxy which severely under-estimates the true masses (von der Linden et al. 2014;CoMaLit-II). Furthermore, the subsample with WL masses under-represents the clusters with low SNR just above the threshold (CoMaLit-II). The larger mean SNR of the subsample determines a larger average limit, as we observed.
The flatness of the completeness limits of the optical richness-WL subsample (see upper panel of Fig. 6) is also connected to the properties of the parent sample. The redMaPPer catalog is in fact nearly complete at z < ∼ 0.3 and λ > ∼ 30 (Rykoff et al. 2014).

Mass function
The modelling with a redshift-evolving Gaussian function, see Eqs. (24,25) and (26), is functional as far as the mass distribution is fairly unimodal and the redshift evolution is smooth. In the case of a cluster sample selected through a hard cut in the observable, these assumptions are well justified, see Sec. 2 and App. A. However, the scheme is flexible enough to work even with heterogeneous sample since detections and measurements are generally limited by observational thresholds. As far as the distribution is unimodal, a simple parameterisation of the distribution of the covariate variable in terms of a single Gaussian function can determine unbiased estimates of the scaling relation. Results are full consistent with more complex parameterisations adopting mixtures of Gaussians (Kelly 2007;CoMaLit-II).
The observed WL mass functions in different redshift bins are well reproduced by the regression model. Distributions for the velocity dispersion-, optical richness-, X-ray luminosity-, and SZ flux-WL samples are shown in Figs. 9, 10, 11, and 12, respectively. The regression model computes the distribution of the true masses. To compare them with the observed distribution of WL masses, we had to smooth the distribution firstly with a Gaussian function whose standard deviation is the intrinsic scatter between true and proxy mass, which is computed by the regression too, and then with a Gaussian whose dispersion is given by the observational uncertainty on the WL masses. We considered the median uncertainty in the redshift bin.

Predictions for the CLASH sample
We can now make use of the constraints on the investigated scaling relations to obtain robust predictions on the observed properties of the galaxy clusters that are part of the CLASH program (Cluster Lensing And Supernova survey with Hubble, Postman et al. 2012), which provides one of the samples best studied at different wavelengths.
The scaling relation σv-M200 can be used to predict the re- Table 4. Scaling relations as a function of the WL mass, E −γss z O = 10 α M β WL,∆ E γz z . Listed parameters refer to logarithmic variables. For the conditional scaling, the adopted form is Y X = α Y |X + β Y |X X + γz log Ez. Conventions and units are as in Sec. 3 and Tables 1 and 2.
For these predictions, we performed a new regression excising from the sample the eight CLASH clusters previously included. Since we are interested in predicting the velocity dispersion given the weak lensing mass, we considered the conditional σv-M WL 200 relation rather than the true mass-velocity dispersion σv-M200. The results of the regression were in full agreement with the regression of the full sample: α Y |X = 2.75 ± 0.05; β Y |X = 0.22 ± 0.05; γz = −0.02 ± 0.24; σ Y |X,0 = 0.07 ± 0.01; γσ Y |X = 0.07 ± 0.91. Alternatively, we might have used the σv-M200 scaling considering the additional source of error given by the intrinsic scatter between the known weak lensing masses and the unknown true masses. As a general remark, if we have weak lensing masses we cannot make predictions by plugging them in scaling relations which compare for example the observable to the hydrostatic mass.
We considered for the CLASH clusters the weak lensing masses reported in LC 2 -all (CoMaLit-III) and based on Umetsu et al. (2014), who performed a combined analysis of shear and magnification. The velocity dispersions based on the σv-MWL,200 relation are listed in Table 5. The main source of statistical uncertainty on the predictions is due to the intrinsic scatter of the relation, which abundantly tops the uncertainties due to the propagated error on the WL mass or due to the uncertainties in the scaling relation parameters.
The prediction for MACS J1206.2-0847 compares well with the first measurement from CLASH-VLT (Biviano et al. , σv = 1087 −55 km/s). The HeCS covered two clusters later on analysed in Umetsu et al. (2014), i.e., A2261 and RXJ2129. The prediction for RXJ2129 is in excellent agreement with the measurements (Rines et al. 2013, σv = 858 +71 −57 km/s). On the other hand, the prediction slightly exceeds the observed velocity dispersion of A2261 (Rines et al. 2013, σv = 780 +78 −60 km/s), even though the small discrepancy is fully covered by the uncertainty due to the intrinsic dispersion of the scaling. From top to bottom, the red, green, and blue lines show the 85, 50, and 15 per cent completeness levels, respectively. The shaded green region encloses the 68 per cent confidence region around the 50 per cent completeness level due to uncertainties on the mass function parameters. Note that the completeness is a function of the true mass, whereas the points refer to the weak lensing masses, which are scattered with respect to the true mass. Bottom panel: scaling between velocity dispersion and mass, M 200 . The black points mark the data (weak lensing mass and redshift evolved velocity dispersion), the blue (green) line represent the conditional (symmetric) scaling relation fitted to the data (true mass versus redshift evolved velocity dispersion). The dashed blue lines show the median scaling relation (full blue line) plus or minus the intrinsic scatter. The shaded blue region encloses the 68 per cent confidence region around the median relation due to uncertainties on the scaling parameters. The red line represents the theoretical prediction based on Munari et al. (2013) at the median redshift of the sample. Masses are in units of 10 14 M .

COMPARISONS
In this section, we compare our results to theoretical predictions or previous estimates.

Theoretical predictions
We first compare our results to predictions based either on theoretical models of structure formation or on numerical/hydro-dynamical simulations.

LX-M∆
Stanek et al. (2010) presented a computational study of the intrinsic covariance of cluster observables using the Millennium Gas Simulations. Two different physical treatments were proposed: shock heating driven by gravity only, or a second treatment with cooling and preheating. The predictions in Stanek et al. (2010) on the scaling between mass and bolometric X-ray luminosity are strongly dependent on the adopted scheme. Acceptable values of the slope are in the range 1.1 < ∼ β < ∼ 1.9, consistent with our estimate of βY -Z = 1.8 ± 0.3. Ettori et al. (2004a) noted a significant negative evolution in the LX-M∆ relation (γz ∼ −0.94) in a sample of local and highredshift galaxy clusters extracted from a large cosmological hydrodynamical simulation with cooling and preheating, supporting the first evidence of a negative evolution (with respect to the selfsimilar model) observed in Chandra data of high-z clusters in Ettori et al. (2004b). This negative evolution with z is in agreement with our findings (γz = −1.7 ± 0.8), even though the large statistical error make the estimate marginally compatible with selfsimilarity. We noted negative evolution in two different samples of X-ray clusters, i.e, the sample from Maughan et al. (2012) and the MCXC, whose luminosities were largely based on independent data-sets and procedures.
A exhaustive study of the joint effect of feedback from supernovae (SNe) and active galactic nuclei (AGNs) on the evolution of  the X-ray scaling laws presented in Short et al. (2010) predicts an opposite behaviour of the evolution of the LX-M∆ relation. They found that the energy output from SNe and AGNs as implemented through semi-analytic models of galaxy formation causes a positive evolution. On the other hand, simulations based on a pre-heating model where an entropy floor of 200 keV cm 2 is introduced at z = 4 confirmed the presence of a negative evolution (Short et al. 2010).
Positive evolution (i.e., higher luminosities at higher redshift, for a fixed mass, than the self-similar prediction) was also found in Pike et al. (2014) in a series of radiative hydrodynamical models, which suggested that radiative cooling is the main driver for departures from self-similarity.
The pre-heating scenario seems to be preferred from our results.

YSZ-M∆
The consensus from numerical simulations is that the YSZ-M∆ is approximately self-similar in mass (1.6 < ∼ β < ∼ 1.8) and it is characterised by a small intrinsic scatter (Stanek et al. 2010;Kay et al. 2012;Battaglia et al. 2012). Stanek et al. (2010) found that the evolution with redshift might be negative (γz ∼ −0.34) in presence of cooling and preheating. In agreement with simulations, we did not detect any departure from the self-similar scaling, with βY -Z = 1.50 ± 0.21. The uncertainty in the observed time-evolution is too large (γz = −0.8±0.8) to infer any statistically significant deviation from self-similarity. In fact, the evolution with mass is fully consistent with the findings of CoMaLit-II, where we found βY -Z = 1.37 ± 0.15 assuming a self-similar redshift evolution.

σv-M∆
Numerical simulations confirmed that the σv-M∆ relation is consistent with the self-similar scaling with mass (Evrard et al. 2008;Munari et al. 2013;Saro et al. 2013). Some differences may arise from the galaxy population used to estimate the velocity dispersion and from the impact of selection using galaxy colour, projected separation from the cluster centre, galaxy luminosity, and spectroscopic redshift (Saro et al. 2013). Whereas dark-matter particles in simulations trace a relation that is fully consistent with the theoretical expectations, sub-haloes and galaxies trace slightly steeper relations with β just above 1/3, and with slightly larger values of the normalisation (Munari et al. 2013). This is due to dynamical processes, namely dynamical friction and tidal disruption, which act on substructures and galaxies, but not on dark matter particles. The relevance of these effects depends on the halo mass and the effectiveness of baryon cooling, and may create a non-trivial dependence of the scaling relation on the tracer, the halo mass, and its redshift (Munari et al. 2013). A better statistical accuracy than that achieved with our results is needed to detect such effects. Saro et al. (2013) noted a substantial agreement between the time evolution of the σv-M∆ relation and the expected self-similar PDF Figure 11. Mass function of the clusters from the X-ray luminosity-WL sample in four redshift bins. Lines and conventions are as in Fig. 9. evolution, γz ∼ 0, which we confirm here within the statistical uncertainties.
The main sources of bias and scatter in velocity dispersion at fixed mass are the halo triaxiality, sampling noise, the presence of multiple kinematic populations within the cluster, and the effect of interlopers (Saro et al. 2013). Saro et al. (2013) found σ log(σv/M ∆ ) ∼ 0.05 locally, and that the intrinsic scatter increases with redshift, with velocity dispersions that are ∼25 per cent less accurate for estimating single cluster masses at z = 1 than at low redshift. Stanek et al. (2010) found a smaller scatter of σ log(σv/M ∆ ) ∼ 0.02 for dark matter particles. Our findings, i.e., σ log(σv/M ∆ ) ∼ 0.06±0.02, slightly exceed the theoretical predictions. The accuracy of our results is not good enough to appreciate any evolution of scatter with redshift.

Previous estimates
We now discuss our findings in relation to previous results. We limited the comparison mainly to scaling relations which employed direct mass measurements based on weak lensing or the assumption of hydrostatic equilibrium, whereas we mostly discarded works based on mass estimates based on external calibration. Previous analyses of the σv-M∆ were mostly based on mass measurements assuming the dynamical equilibrium and the virial theorem (Girardi & Mezzetti 2001;Rines et al. 2013). This mass measurement is direct too, but it is strongly correlated with the velocity dispersion, differently from the weak lensing masses we considered here.
We did not consider previous analyses of the conditional scaling relations in which the mass worked as the response variable. These relations are not easily inverted and cannot be compared straight on with our results.

LX-M∆
Observed slopes of the LX-M∆ relation from previous studies may be steeper than the self-similar prediction. Estimated values of β ranges from 1.3 to 1.9 (Vikhlinin et al. 2009;Pratt et al. 2009;Arnaud et al. 2010;Reichert et al. 2011;Ettori 2013). Source of disagreement may be various. The slope and normalisation of the relation depend on the energy band and method used for the flux extraction (Ettori 2015).
The intrinsic scatter of the LX-M∆ relation is ∼40 per cent (Giodini et al. 2013). It is the largest among the various X-ray scaling relations. We confirmed the large scatter in the LX-M∆ relation. The X-ray luminosity is heavily affected by non-gravitational processes, the presence of cool-cores, and the overall dynamical state of the halo (Giodini et al. 2013). Most of the scatter derives from the inner regions where cooling and merging effects are most pronounced. Our estimate of the scatter based on the the soft band luminosities provided by the MCXC is fully consistent with most of the previous results (∼ 40 per cent), whereas the estimate based on the core-excised bolometric luminosities from Maughan et al. (2012) is significantly smaller (∼ 20 per cent). This suggests that a careful choice of the energy band and of the methods for the flux measurement might significantly reduce the intrinsic scatter of the LX-M∆ relation.
A negative time-evolution of the relation was firstly noticed in Ettori et al. (2004b) and later confirmed by Reichert et al. (2011), which found γz = −1.3 ± 0.2 applying a tentative selection-bias correction. Our results confirm these previous findings.
From a multivariate analysis aimed to study X-ray luminosity, temperature, and gas mass fraction in a sample of clusters with well measured WL masses in the context of cosmological parameter determinations with cluster abundances, Mantz et al. (2014) estimated a slope of 1.71 ± 0.17 and an intrinsic scatter of 42 ± 5 per cent for the LX soft -M∆ relation assuming self-similar time evolution. Apart from the assumed redshift dependence, these results are directly comparable to ours, since Mantz et al. (2014) considered the scatter of the WL mass and the effects of the selection function. The estimated slope is in very good agreement with our result β Y |Z = 1.60±0.27. In principle, slopes obtained from a multivariate analysis may differ from the results of a single O-M∆ relation if the considered observables are strongly correlated (Ettori 2013). However, this is not the case of the X-ray properties considered in Mantz et al. (2014). Rozo et al. (2014a) developed a self-consistent method to derive scaling relations satisfying optical data from SDSS, X-ray data from ROSAT and Chandra, and SZ data from Planck. Assuming a self-similar time-evolution, they derived a slope of β Y |X = 1.55 ± 0.09 for the LX soft -M500 relation with a scatter of ∼ 39 ± 3 per cent. Slope and scatter from Rozo et al. (2014a) are consistent with our results for the LX soft -M∆ relation based on the MCXC (see Sec. 5) even though the analyses presents some major differences. In fact, Rozo et al. (2014a) used stacked data rather than measurements from single clusters and they assumed a self-similar time evolution.

YSZ-M∆
Observed slopes of the scaling relation between mass and SZ flux are discordant to some degree (CoMaLit-II). The Planck team determined the YSZ-M500 relation relying on masses estimates based on the YX proxy and assuming a self-similar time-evolution (Planck Collaboration et al. 2014b). Through the BCES-orthogonal regression, they found βY -X = 1.79 ± 0.06 and an intrinsic orthogonal scatter ∼ 15-20 per cent (Planck Collaboration et al. 2014b). A previous calibration based on 19 weak lensing clusters mainly from the LoCuSS sample (Local Cluster Substructure Survey, Okabe et al. 2010) gave βY -X = 1.7 ± 0.4 (Planck Collaboration et al. 2013). However, weak lensing masses of the LoCuSS clusters are biased low due to contamination effects and systematics in shape measurements (Okabe et al. 2013). The underestimate might be mass dependent and affect the estimated slope. With a self-consistent method, Rozo et al. (2014a) found a slope of β Y |X = 1.71 ± 0.08 with a scatter of ∼ 15 ± 2 per cent. These estimated slopes are steeper but consistent within the statistical uncertainty with our result of βY -X = 1.5 ± 0.2. Andreon (2014) claimed that the evolution of the YSZ -M500 relation is significantly inconsistent with the self-similar evolution. He found γz = 1.8 ± 0.4. The disagreement with our result, which is fully consistent with the self-similar prediction, might be due to the arbitrary choice in Andreon (2014) to use a pre-determined completeness function. An inappropriate modelling can bias the estimates of the scaling relations as well as neglecting time-evolution at all. Furthermore, Andreon (2014) neglected the intrinsic scatter in the mass estimate, and, following Planck Collaboration et al. (2014b), he employed a mass proxy based on YX, which is strongly correlated to the SZ signal. These choices likely explain the disagreement with our result.

λ-M∆
Most of the analyses correlating optical richness to mass were based on stacking techniques (Rozo et al. 2014a;Covone et al. 2014, and references therein). Here, we focus on not binned data. Wen, Han & Liu (2009) considered a compilation of clusters whose masses had been estimated by X-ray or weak-lensing methods to infer a slope of βY -Z = 1.17 ± 0.03.
The slower slope found in Andreon & Bergé (2012) might be due to their assumption that the masses, which were measured with the caustic method, were unscattered measurements of the true masses. However, masses based on this method are highly scattered (CoMaLit-II). Neglecting this effect makes relations flatter (CoMaLit-I; CoMaLit-II). The very large scatter of ∼ 58 ± 7 per cent measured in Andreon & Bergé (2012) is biased high for the same reason.

CONCLUSIONS
To assess the role of the evolution with time of the scaling relations between mass and observed properties it is crucial to avoid biases in the calibration. The evolution of cluster scaling relations is still debated. Small samples or biases due to the (unknown) selection function are two of the main problems. We developed a regression methodology that at the same time constrains the evolution and calibrates the completeness function of the studied sample.
The method is general and lets the data determine the timedependent scatter of the mass proxy or the intrinsic scatter between the observable and the mass. Selection functions and their timedependence are often not known a priori. In our approach, the selection function can be determined in the context of the regression procedure. The approach we implemented is Bayesian and can easily include any additional or a priori information on the completeness of the sample. This method is functional in the context of large photometric surveys such as Euclid (Laureijs et al. 2011), where selfcalibration of scaling relations is crucial to unbiased estimates of dark energy with study of cluster abundances (Majumdar & Mohr 2004) We tested the method with large heterogeneous samples to calibrate either optical properties, such as richness and velocity dispersion, or observables connected to the intra-cluster medium, such as X-ray luminosity and SZ flux. Masses were estimated with weak lensing analyses and intrinsic scatter in the mass estimate was considered. Weak lensing masses are reliable mass measurements up to high redshifts.
To our knowledge and not considering mass estimates based on the viral theorem, this is one of the first studies to compare galaxy velocity dispersions to direct estimates of masses of clusters (weak lensing masses in our case). This is the approach usually followed in numerical simulations (Evrard et al. 2008;Munari et al. 2013;Saro et al. 2013) to built mass proxies based on the velocity dispersion without assuming dynamical equilibrium and without exploiting the properties of the infall patterns. We found that observables scale self-similarly with respect to the mass and that they evolve self-similarly with cosmic time. The only exception is the LX-M∆ relation, which seems to show a negative evolution. A similar level of evolution can be obtained in hydrodynamical simulations of the intracluster medium by including additional radiative cooling and uniform preheating at high redshift as a simple model of non-gravitational heating from astrophysical sources (Ettori et al. 2004a;Short et al. 2010;Reichert et al. 2011).
The intrinsic scatter in the mass-velocity dispersion relation is notably small, which encourages the use of velocity dispersions as mass proxies. Our procedure accounts for time evolution of the intrinsic scatter. However, the large statistical uncertainties made the parameters characterising this dependence as de facto noise parameters which had to be marginalised over to avoid to bias low the uncertainties on the scaling relation. The determination of this evolution is out of reach in present samples of nearly one hundred of clusters (Mantz et al. 2010a), but is should be feasible in future large surveys that will be able to collect homogeneous multi-band information of an unprecedented number of clusters.  Figure A1. Distribution of true masses M 200 of a sample of clusters at z = 0.3 (grey histogram) whose observable proxy mass was larger than M th,200 = 0.5×10 14 M /h (vertical blue line). The assumed log-normal scatter is σ = 0.25. The empty histogram represents the full distribution of the true masses of all clusters before the selection; the grey histogram represent the distribution of selected clusters, i.e, clusters whose observed proxy mass is above the threshold (vertical blue line); the red line is the Gaussian approximation to the mass function.. The true masses were extracted from a cosmological halo mass function following Tinker et al. (2008).   Figure A2. Probability function of the log-mass µ, i.e., the logarithm of the over-density mass M ∆ . We assume that µ is intrinsically distributed as a local power-law and that clusters are selected if they have a proxy mass o > −1. The value of the slope of the halo function β 1 , and of the intrinsic scatter for the different cases are reported in the legend.