Dark Energy Survey Year 1 Results: Methods for Cluster Cosmology and Application to the SDSS
Abstract
We perform the first blind analysis of cluster abundance data. Specifically, we derive cosmological constraints from the abundance and weaklensing signal of redMaPPer clusters of richness in the redshift range as measured in the Sloan Digital Sky Survey (SDSS). We simultaneously fit for cosmological parameters and the richness–mass relation of the clusters. For a flat CDM cosmological model with massive neutrinos, we find . This value is both consistent and competitive with that derived from cluster catalogues selected in different wavelengths. Our result is also consistent with the combined probes analyses by the Dark Energy Survey (DES) and the KiloDegree Survey (KiDS), and with the Cosmic Microwave Background (CMB) anisotropies as measured by Planck. We demonstrate that the cosmological posteriors are robust against variation of the richness–mass relation model and to systematics associated with the calibration of the selection function. In combination with Baryon Acoustic Oscillation (BAO) data and BigBang Nucleosynthesis (BBN) data (Cooke et al., 2016), we constrain the Hubble rate to be , independent of the CMB. Future work aimed at improving our understanding of the scatter of the richness–mass relation has the potential to significantly improve the precision of our cosmological posteriors. The methods described in this work were developed for use in the forthcoming analysis of cluster abundances in the DES. Our SDSS analysis constitutes the first part of a stagedunblinding analysis of the full DES data set.
keywords:
cosmology: cluster, cluster: general1 Introduction
Galaxy clusters form from the high density peaks of the initial matter distribution. As such, they bear the imprints of the statistical properties of the matter density field and its growth (for reviews, see e.g. Allen et al., 2011; Kravtsov & Borgani, 2012). The abundance of galaxy clusters has been used since the late 90’s to constrain the mean matter density of the Universe, , and the amplitude of the density fluctuations in terms of , the presentday rms of the linear density field in spheres of radii (e.g. Eke et al., 1998; Henry, 2000; Borgani et al., 2001; Pierpaoli et al., 2001; Reiprich & Böhringer, 2002; Henry, 2004, for early works). Constraints on are especially powerful in combination with measurements of the amplitude of the matter power spectrum at high redshift — e.g. from CMB data — enabling us to study the growth of density perturbations over cosmic time. These studies allow us to place constraints on parameters such as the total neutrino mass, the dark energy equation of state, and parameters governing modified gravity models.
Current studies using cluster catalogs selected in the Xray, optical, and millimeter wavelengths, provide consistent constraints on and (e.g. Vikhlinin et al., 2009; Rozo et al., 2010; Mantz et al., 2010; Mantz et al., 2015; Planck Collaboration et al., 2016b; de Haan et al., 2016). These data have been also used in combination with and Baryon Acoustic Oscillation (BAO) priors to place competitive constraints on the dark energy equation of state parameter, modification of gravity and neutrino masses (e.g. Burenin & Vikhlinin, 2012; Mantz et al., 2010; Mantz et al., 2015; Cataneo et al., 2015; Planck Collaboration et al., 2016b; de Haan et al., 2016). Ongoing — e.g. the Dark Energy Survey (DES)^{2}^{2}2https://www.darkenergysurvey.org, the Hyper SuprimeCam^{3}^{3}3http://hsc.mtk.nao.ac.jp/ssp/ — and forthcoming — the Large Synoptic Survey Telescope^{4}^{4}4https://www.lsst.org/, Euclid^{5}^{5}5http://sci.esa.int/euclid/, eRosita^{6}^{6}6http://www.mpe.mpg.de/eROSITA — widearea surveys aim to use clusters samples with an order of magnitude more systems than previous analyses in order to improve upon current constraints.
The most critical difficulty that cluster abundance studies must confront is the fact that cluster masses are not easily measured, forcing us to rely on observational proxies that correlate with mass. Specifically, while it is possible to predict the abundance of dark matter halos as a function of mass in an arbitrary cosmology with percent level accuracy (e.g. Sheth & Tormen, 1999; Tinker et al., 2008; Crocce et al., 2010; Bocquet et al., 2016; McClintock et al., 2018a), halo masses themselves are not directly observable. At present, cosmological constraints from cluster abundance analyses at all wavelengths are limited by the uncertainty in the calibration of the relation between the cluster mass and the observable property used as a mass proxy, be it richness (i.e. count of member galaxies), Xray luminosity, or the thermal SunyaevZeldovich signal.
Currently, weak gravitational lensing measurements provide the goldstandard technique for estimating cluster masses (see e.g. von der Linden et al., 2014b, for a discussion). The weak lensing signal, i.e. the tangential alignment of background galaxies around the foreground cluster due to gravitational lensing, is a wellunderstood effect, sensitive to both dark and baryonic matter. Moreover, in contrast to other techniques (e.g. velocity dispersion and hydrostatic mass measurements), weak lensing mass measurements do not rely on assumptions about the dynamical state of the cluster. Despite these advantages, many sources of systematic error do affect this type of measurement, e.g. shear and photometric redshift biases, halo triaxiality, and projection effects. These systematics represent a significant amount of the total error budget of many recent studies (e.g. von der Linden et al., 2014a; Hoekstra et al., 2015; Simet et al., 2017; Melchior et al., 2017; Medezinski et al., 2018; Miyatake et al., 2018). Not surprisingly, as the statistical uncertainty continues to decrease, these systematics have come to dominate the total error budget (e.g. McClintock et al., 2018b).
In this work we combine cluster abundances and stacked weak lensing mass measurements from the Sloan Digital Sky Survey data release 8 (SDSS DR8, Aihara et al., 2011) to simultaneously constrain cosmology and the richness–mass relation of galaxy clusters. Our cluster sample is selected using the red sequence Matchedfilter Probabilistic Percolation algorithm (redMaPPer; Rykoff et al., 2014), and the stacked weak lensing mass estimates are presented in Simet et al. (2017). The analysis is similar in spirit to that of Rozo et al. (2010) but with significant updates, particularly with regards to the modeling of the cluster selection function. Our observables are the number of clusters and the mean cluster mass in bins of richness — our mass proxy — and redshift. We explicitly account for the small cosmological dependence of the recovered weak lensing masses in our analysis. To avoid confirmation bias, the bulk of the analysis has been performed blind: the values of the cosmological parameters sampled by the Monte Carlo Markov Chains (MCMC) were randomly displaced by an amount that was unknown to us, and were shifted back only after a broad set of validation tests were passed. No changes were done to the analysis pipeline postunblinding. This is the first cosmology analysis to be performed blind in the cosmological parameters. The methods presented in this paper were developed for the forthcoming cosmological analysis of the DES Y1 redMaPPer cluster catalog (DES Collaboration in prep.).
This paper is organized as follows. In Section 2 we introduce the data used for this study. In Section 3 we present the model used to perform the cosmological analysis, including expectation values for the two observables, modeling of the systematics and likelihood model. We validate our model by means of synthetic data in Section 4. We detail our blinding procedure in Section 5.The results of our analysis are presented in Section 6 and 7. Finally we conclude in Section 8.
2 Data
2.1 Cluster and Weak Lensing Catalogs
Both the cluster and weak lensing shear catalogs used in this analysis are based on the Sloan Digital Sky Survey data release 8 (SDSS DR8, Aihara et al., 2011). A summary of the data employed in this analysis is presented in Table 1. Throughout the paper, all masses are given in units of , where , and refer to an overdensity of 200 with respect to the mean. We use "" and "" to refer to the logarithm with base and respectively.
Number counts  

We use photometrically selected galaxy clusters identified in the SDSS DR8 with the redMaPPer cluster finding algorithm (Rykoff et al., 2014). In brief, redMaPPer utilizes all five bands (ugriz) of the SDSS imaging to selfcalibrate a model for redsequence galaxies. This model is then used to identify galaxy clusters as red galaxy overdensities, while simultaneously estimating the probability that each red galaxy is a cluster member. The cluster richness is the sum of the membership probabilities over all the red galaxies within an empirically calibrated scale radius, . Due to the probabilistic approach adopted by redMaPPer an individual galaxy can contribute to the richness of more than one cluster. We use the cluster richness as an observational massproxy in our analysis, calibrating the relation between galaxy richness and cluster mass using weak lensing data. Of the covered by SDSS DR8, we restrict ourselves to the of high quality contiguous imaging defined by the Baryon Oscillation Spectroscopic Survey (BOSS) experiment (Dawson et al., 2013). Typical photometric redshift uncertainties for SDSS redMaPPer clusters are . We further restrict ourselves to the redshift range . The lower redshift limit ensures that galaxies are sufficiently dim that the SDSS photometric pipeline produces reliable galaxy magnitudes, while the upper redshift limit is set by requiring that the cluster catalog be volume limited. We also apply a richness threshold . In numerical simulations this richness threshold ensures that 99% of the redMaPPer galaxy clusters can be unambiguously mapped to individual dark matter halos (Farahi et al., 2016). In this work, we use v5.10 of the SDSS redMaPPer cluster catalog (Rozo et al., 2015).
The weak lensing mass estimates employed in this analysis are a slight update from those presented in Simet et al. (2017), and rely on the shear catalog presented in Reyes et al. (2012a). This catalog covers of the SDSS footprint and contains 39 million galaxies, corresponding to a source density of . Shear estimates were derived from the SDSS imaging using the reGaussianization algorithm of Hirata & Seljak (2003) and the appropriately calibrated responsivity to convert the measured shape distortions into shear estimates. The multiplicative shear bias appropriate for this catalog was characterized in Mandelbaum et al. (2012) and Mandelbaum et al. (2013). In a recent work, Mandelbaum et al. (2017) revised the multiplicative shear estimate of this source catalog due to a previously undetected bias arising from the impact of neighboring galaxies and unrecognized blends. Our error budget accounts for this additional source of error. The photometric redshifts for the sources in the shear catalog were obtained using the Zurich Extragalactic Bayesian Redshift Analyzer (ZEBRA, Feldmann et al., 2006), and the associated systematic uncertainties were calibrated in Nakajima et al. (2012).
2.2 Cluster Number Counts Data
Following the weak lensing analysis of Simet et al. (2017), we collect our galaxy clusters in five richness bins and a single redshift bin (see Table 1). The richness limits of our bins are . The two key observational systematics in our analysis are photometric redshift errors and cluster miscentering. Photometric redshift uncertainties are forwardmodeled as described in section 3.1.
Turning to cluster centering, in this work we assume that the correct center of a galaxy cluster is always coincident with a bright cluster galaxy (though not necessarily the brightest). This assumption is motivated by the fact that modern halo finders (e.g. Behroozi et al., 2013) define the center of dark matter halos as the position of the dominant dark matter substructure within the halo. It is expected that such a substructure will host a bright cluster galaxy, and therefore a proper comparison of theory to data should rely on identifying the correct central galaxy for each cluster. In redMaPPer, this is accomplished through an iteratively selftrained algorithm that combines luminosity and local galaxy density information to determine which cluster galaxy is the best central candidate in the cluster. This algorithm is demonstrably superior to centering clusters on the brightest clustermember galaxy (Hikage et al., 2017). Nevertheless, miscentering systematics still arise from the fact that redMaPPer does not always successfully identify the correct central galaxy in each cluster. As we shall see shortly, the uncertainty in the systematic corrections due to cluster miscentering are very nearly negligible, so rather than incorporating these corrections into our forwardmodel, we have opted for the simpler route of applying a correction to the observed datavector, and adding the corresponding systematic uncertainty in the abundances to the covariance matrix characterizing the uncertainty in our measurements.
Our miscentering correction is based on the analyses in Zhang et al. (in prep.) and von der Linden et al. (in prep.). To characterize the probability of a redMaPPer cluster being miscentered, we use subsamples of redMaPPer clusters with either Chandra or Swift Xray imaging. Specifically, we compared the redMaPPer central galaxies to the Xray derived cluster centers. Doing so, we demonstrated that of redMaPPer galaxy clusters are correctly centered. This value is somewhat lower than but consistent with the estimate of Hikage et al. (2017) on the basis of the weak lensing profile of redMaPPer clusters. In addition, we also characterized the radial offset distribution of the miscentered clusters using a 2dimensional Gaussian to describe the distribution of positional offsets for the miscentered clusters.
To determine the impact that cluster centering has on the recovered cluster richness, we measured the richness of all redMaPPer clusters at the location of the second most likely central galaxy, and characterized how cluster richness changes as a function of the centering offset. As expected, we find that miscentering systematically underestimates cluster richnesses, albeit with some scatter. Both the bias and scatter increase as the radial offset increases (Zhang et al., in prep.).
Having characterized (i) the fraction of miscentered clusters; (ii) the distribution of radial offsets of miscentered clusters; and (iii) how the richness of a cluster changes when it is miscentered by a given radial offset, we can readily estimate the impact of miscentering on the cluster abundance function. Specifically, we assigned richness values to halos in a numerical simulation using the model of Costanzi et al. (2018). The assigned halo richnesses are then scattered using our centering model. We compute the ratio between , the number of clusters in a bin in the absence of miscentering, and , the number of clusters in a bin including the impact of miscentering. We multiply the cluster abundance data vector by to correct for the impact of miscentering on the observed cluster abundances. From our Monte Carlo realizations, we are able to characterize both the mean correction and the corresponding covariance matrix. We apply the (see Table 1) systematic correction to the number of galaxy clusters we observe, and add the corresponding covariance matrix to the covariance matrix of the cluster counts. The uncertainty associated with cluster miscentering in the abundance function is . That is, cluster miscentering is subdominant to Poisson noise in all richness bins. Its impact on the cosmological posterior is easily subdominant to the uncertainty in our calibration of the richness–mass relation.
The shaded regions in the topleft panel of Figure 1 shows the miscenteringcorrected cluster abundances. The width of the regions along the axis is set by the diagonal entries of the corresponding covariance matrix, whereas their width along the axis reflects the width of the richness bin. Our bestfit model is shown as the blue points. The bottomleft panel shows the corresponding percent residuals.
2.3 Weak Lensing Cluster Masses
2.3.1 Measurement
We calculate the mean mass of clusters in a richness bin using the stacked weak lensing mass profiles of the clusters as described in Simet et al. (2017). These profiles are modeled using a Navarro, Frenk, and White profile (NFW, Navarro et al., 1997), accounting for the effects of cluster miscentering, halo triaxiality, and projection effects. The concentration of the bestfit NFW profiles is modeled using the massconcentration relation presented in Diemer & Kravtsov (2015) with a free parameter for the amplitude.
In this work, we follow the methodology of Simet et al. (2017) to estimate the mean cluster mass in each bin with three critical exceptions:

We correct for the multiplicative shear bias due to undetected blends that was first identified in Mandelbaum et al. (2017). The original multiplicative bias in Simet et al. (2017) had a systematic uncertainty of (tophat). To that we add a second bias term with a tophat prior of width . The sense of the correction is to increase the recovered weak lensing masses. The standard deviation of the resulting distribution is , i.e. roughly equivalent to a 3.4% Gaussian scatter.

We update the centering priors employed in the weak lensing analysis based on the work of Zhang et al. and von der Linden et al. (in preparation). In Simet et al. (2017), the fraction of correctly centered galaxy clusters was estimated as on the basis of published Xray data. The analyses in Zhang et al. and von der Linden et al. (in preparation) benefit from significantly improved statistics and custom Xray analyses specifically designed to calibrate this source of systematic uncertainty. The revised fraction of correctly centered clusters is . The radial offset of miscentered clusters is described by a two dimensional Gaussian of width , where is the cluster radius assigned by redMaPPer.

Rather than simultaneously modeling all galaxy clusters to derive a scaling relation that describes all richness bins, we estimate the mean mass of each individual richness bin. To account for the finite width of the bin, we fit a scaling relation to each individual bin, and then use the posterior of the scaling relation parameters to determine the mean cluster mass in that bin. It is these mean cluster masses that we use as the weak lensing observables in our cosmological analysis.
The recovered mean mass in any single bin is only weakly dependent on the scatter assumed in the mass–richness relation when performing the fits. The degeneracy of the weak lensing posterior takes the form . That is, varying over the range modifies the recovered masses by an amount ranging from 0 to 0.015. We adopt a fiducial correction of . Following the analysis described in Simet et al. (2017, see section 5.5), we further estimate the systematic uncertainty in our weak lensing masses due to modeling the lensing profile with a pure NFW halo, without accounting for a twohalo term, or due to deviations from the NFW profile (e.g. Diemer & Kravtsov, 2014). Briefly, we artificially populate dark matter halos in an body simulation by drawing from the richness–mass relation obtained by inverting our best fit mass–richness relation (Evrard et al., 2014). We then compute the corresponding cluster–mass correlation function, and project along the line of sight to obtain a predicted weak lensing profile, which we fit in the same way as we fit the data. The recovered biases for each of our richness bins varies from 2% on the low mass end to 3% on the high mass end. We apply these corrections to our data, and assign a systematic uncertainty on this correction equal to half the magnitude of the correction.
Our results are summarized in Table 1, where we collect the best weak lensing estimates for the mean mass of the galaxy clusters in each richness bin. The logarithm of the mean mass, , for each of our five richness bins is shown in the topright panel of Figure 1, along with the bestfit model from our cosmological analysis. The bottomright panel shows the corresponding residuals. We again emphasize that all aspects of this weak lensing analysis are as per Simet et al. (2017) except as detailed above; our results are a minor update to that work.
2.3.2 Weak Lensing Systematic Error Budget
Cluster cosmology has long been limited by systematic uncertainties in mass calibration. Here, we summarize the observational systematics we have explicitly accounted for in our analysis and their contributions to the error budget of the mass estimates. For further details on how these systematics were estimated, we refer the reader to Simet et al. (2017). Table 2 summarizes the sources of error detailed below.

Multiplicative shear and photo bias ( tophat). This bias reflects the total multiplicative systematic uncertainty arising from two distinct sources of error, each contributing roughly the same amount of uncertainty: 1) systematic uncertainties in shear calibration as evaluated from realistic image simulations (Reyes et al., 2012b; Mandelbaum et al., 2012; Mandelbaum et al., 2013), and 2) photometric redshift systematics (Nakajima et al., 2012). This uncertainty is on the amplitude of the lensing profile. The equivalent mass error is obtained by multiplying by , resulting in a tophat error in the mass.^{7}^{7}7The factor of 1.3 can be estimated by boosting a Navarro et al. (1997) profile by a constant factor, estimating the new aperture, and then verifying the ratio of the boosted mass to the original mass. The end result is that the boost to the mass is 1.3 times the boost to the profile. In practice, the mass measurements are marginalized over the tophat amplitude prior above.

Blending of lensing sources ( top hat). This refers to the systematic bias introduced by the blending of sources discovered by Mandelbaum et al. (2017), and characterized using deep imaging from the HyperSuprime Camera. Again, the 3% uncertainty is in the correction of the amplitude of the profile, leading to a error in the mass.

Cluster triaxiality and projection effects (3%, Gaussian). Cluster triaxiality and projection effects both introduce correlated cluster scatter between cluster richness and weak lensing masses. As detailed in Simet et al. (2017, section 5.3 and 5.4), we can readily compute the necessary corrections using the formalism laid out in Evrard et al. (2014). The estimate of the impact of triaxiality derived in this way is in excellent agreement with the simulationbased calibration of Dietrich et al. (2014). The quoted 3% error is the uncertainty in the mass.

Cluster centering (). Cluster centering is explicitly incorporated into the weak lensing mass estimates via forward modeling. The miscentering parameters have informative priors obtained from studies of the subsample of redMaPPer galaxy clusters with highresolution Xray data (Zhang et al. and von der Linden et al., in preparation). The weak lensing posterior fully accounts for the systematic uncertainty associated with cluster miscentering. The quoted 1% error is the uncertainty in the mass.

Modeling systematics (, Gaussian, richness dependent). Modeling systematics refers to biases in the inferred cluster mass due to our observational procedure. We explicitly calibrated this systematic using numerical simulations as detailed in the previous section. The quoted 2% error is the uncertainty in the mass.

Impact of the scatter on the recovered weak lensing masses (1.6%, Gaussian). This error refers to the sensitivity of the weak lensing masses to the assumed scatter in the mass–richness relation detailed in section 2.3.1. The error in corresponds to a 1.6% error on the mass.

Systematic uncertainty in the boost factor corrections (likely negligible). Misidentification of cluster galaxies as lensing sources dilutes the lensing signal. However, these galaxies are obviously clustered around the cluster. Consequently, measurements of the source–cluster correlation function allow one to directly measure the contamination rate of the source galaxy population as a function of radius, and thus estimate the boost factor needed to correct the lensing signal. We explicitly checked that the statistical uncertainty associated with this measurement is negligible compared to the systematic errors in the multiplicative shear bias and photoz corrections. However, smallscale features in the cluster or source selection function that are not adequately modeled by the associated random catalogs could lead to systematic uncertainties in the boost factor correction. Here, we assume these errors are negligible.
The total systematic error budget for our observational estimate of the mean mass of the galaxy clusters in a given richness bin is (Gaussian). Two additional, less wellunderstood potential systematics remain. These are:

Intrinsic alignments: Radial alignments of member galaxies that are mistakenly selected as lensing sources will bias the corresponding cluster mass. Simet et al. (2017) derives a 6% upper limit for the impact of intrinsic alignments, but does not correct for this effect nor enlarge the systematic error budget since this intrinsic alignment signal has yet to be detected in the lowluminosity cluster member galaxies that can be potentially misidentified as lensing sources. Note intrinsic alignments are unidirectional: explicitly modeling intrinsic alignments would increase the recovered cluster masses.

Baryonic impact on the recovered cluster masses. In Simet et al. (2017) it is argued that baryonic effects were unlikely to have a large impact on the recovered cluster masses for two reasons. First, baryonic effects are strongest in the cluster core, which we exclude from our fits. Second, we do not impose priors on the amplitude and slope of the concentration–mass relation, allowing for the concentration to absorb the effects that baryons have on the cluster lensing profile, so long as an NFW profile remains a reasonable description of the halo profiles (Schaller et al., 2015). This hypothesis has received additional support from the work of Henson et al. (2017), who found that while baryonic physics can clearly impact the lensing profile of galaxy clusters, the bias incurred from fitting said profiles is independent of baryonic physics, provided the concentration–mass relation of the clusters is allowed to vary. Based on Figure 11 in that work, it is difficult to imagine baryonic effects introducing biases larger than .
We do not increase our systematic error budget on the expectation that both of these effects will eventually be proven to be subdominant to the systematic errors quoted above. We are not aware of any unaccountedfor systematic effects impacting the recovered weak lensing masses at a level that exceeds that or is comparable to the systematics detailed above. A complete summary of the error budget associated with the weak lensing mass calibration of the SDSS redMaPPer clusters is detailed in Table 2.
Source  Associated Error 

Shear and Photo Bias  6.5% tophat (3.8% Gaussian equivalent) 
Source blending  3.9% tophat (2.3% Gaussian equivalent) 
Cluster triaxiality and projections  3% Gaussian 
Cluster centering  
Modeling systematics  2.0% Gaussian (richness dependent) 
Scatter corrections  1.6% Gaussian 
Total Systematic Error  6.0% Gaussian 
Statistical Error  4.8% Gaussian 
Total  7.7% 
2.3.3 Cosmology Dependence of the Recovered Masses
The weak lensing masses derived above depend on the angular diameter distance to the sources and the mean matter density at the redshift of the clusters. These in turn depend on cosmological parameters. Assuming a flat CDM model, and using mass estimates in units of , the only cosmological parameter that can appreciably impact the resulting weak lensing mass measurements is the matter density . Consequently, we repeat our analysis along a grid of values for ranging from to while setting . As in Simet et al. (2017), we find the recovered logmasses scale as a linear function in the matter density , i.e.
(1) 
The slopes derived from fitting Equation 1 to the data are listed in Table 1, and are used in our cosmological analysis to rescale at each step of the MCMC by the appropriate value.^{8}^{8}8Our posterior extends to matter densities below . We verified a posteriori that the linear matter density scaling extends to .
3 Theory and Methods
In the previous section, we described how we constructed the observational vectors (cluster abundances and weak lensing masses) employed in our analysis. We now turn to describe how we model the expectation values of these observables and derive their covariance in order to simultaneously constrain both cosmology and the richness–mass relation of the redMaPPer clusters. In what follows, all quantities labeled with “ob” denote quantities inferred from observation, while quantities labeled with “true” indicate intrinsic halo properties. denotes the conditional probability of given .
3.1 Expectation values
3.1.1 Base Model
Let denote richness bin , and denote the redshift bin . The expectation value of the number of galaxy clusters is given by
(2) 
where is the comoving volume element per unit redshift and solid angle, and is the effective survey area at redshift . The survey area depends on redshift because galaxy clusters are not pointlike: whether a cluster is formally within the survey area or not depends not just on the location of the galaxy cluster in the sky, but also on how the survey boundaries (including star holes and any other masked regions) intersect the projected area of the cluster in the sky. To estimate the survey area, we randomly place clusters in the sky, and compute the fraction of the galaxy cluster that does not fall within the survey footprint. The footprint of the cluster survey is defined by the collection of all points for which at least 80 per cent of the cluster falls within the photometric survey boundaries. This 80% criterion is the fiducial choice for all redMaPPer runs, and is chosen as a compromise between requiring clusters not be heavily masked, and losing a minimal amount of area due to masking. In principle, this masking criteria implies that the survey area depends on cluster richness (via the scale radius ), but we find this dependence to be negligible ( over the redshift range employed in this study).
The second integral of Eq. 2 accounts for the uncertainty in the photometric redshift estimate. We model — the probability of assigning to a cluster at redshift a photometric redshift — with a Gaussian distribution having mean and a redshift and richnessdependent variance. The variance is set by the reported photometric redshift uncertainty in the redMaPPer cluster catalog. Specifically, we fit a thirdorder polynomial to the redMaPPer photometric redshift errors as a function of redshift for galaxy clusters in each of our five richness bins (we find this is sufficient to fully describe our data). Photometric redshift uncertainties range from at to at , with richer clusters having somewhat smaller photometric redshift uncertainties than low richness clusters.
As shown in Figure 9 of Rykoff et al. (2014), the redMaPPer photometric redshifts are excellent: they are nearly unbiased, and the reported photometric redshift uncertainties are both small and they accurately describe the width of the photometric redshift offsets relative to the spectroscopic cluster redshifts (where available). Using the specific cluster sample employed in this work , we evaluate the systematic bias of the redMaPPer photometric redshift by fitting a Gaussian distribution to the redshift offset , where is our photometric cluster redshift estimate, and is the spectroscopic redshift of the central galaxy assigned to the cluster (when available). We find a mean bias of , i.e. . Likewise, from a Gaussian fit to the distribution of scaled errors we find that the reported photometric redshift uncertainties in redMaPPer should be boosted by a factor of 1.014. In both cases, the statistical uncertainties in the estimates are negligible. We have verified that the above systematic errors are completely negligible by running two versions of our analysis, one without applying these corrections, and one after applying these corrections. The impact on the posteriors is negligible. For specificity, from here on out we apply the above corrections, that is, we set and increase the photometric redshift errors by a factor of .
Finally, in Eq. 2 is the expected number density of halos in the richness bin . This quantity is given by
(3) 
where denotes the probability that a halo of mass at redshift is observed with a richness (see subsubsection 3.1.2) and is the halo mass function which is assumed to follow the form of Tinker et al. (2008):
(4) 
Several studies have explored how the mass function from body simulations should be extended in order to incorporate the effects of massive neutrinos (e.g. Brandbyge et al., 2010; VillaescusaNavarro et al., 2014; Castorina et al., 2014; Liu et al., 2017). A common finding of these studies is that massive neutrinos play a negligible role in the collapse of dark matter halos, while they suppress the growth of matter density fluctuations on scales smaller than the neutrino freestreaming length. Here, we account for these effects following the prescription of Costanzi et al. (2013): (i) we neglect the density neutrino component in the relation between mass and scale — i.e. — and (ii) use only the cold dark matter and baryon power spectrum components to compute the variance of the density field, .
We can use similar arguments to the ones used to derive Eq. 2 to compute the expectation value for the mean mass of galaxy clusters within a specific richness and redshift bin. This is given by:
(5) 
i.e. the ratio of the expected total mass inside the bin over the total number of clusters inside said bin (Eq. 2). Note that observationally, we stack , not , and since the two are not linearly related to each other — is proportional to the integrated mass density profile — it is not necessarily the case that the recovered weak lensing mass is identical to the mean mass of the clusters. However, as described in section 2.3.1, we calibrate the relation between the recovered weak lensing mass and the mean mass in a bin so as to ensure our measurement is an unbiased proxy of the mean mass of the clusters in a bin.
The total mass of all clusters in a bin is calculated via
(6) 
where the total mean mass per unit volume in the th richness bin is given by:
(7) 
In practice, the integrals over the observed redshift in the numerator and denominator of equation 5 are each weighted by the appropriate weak lensing weight per clusters , where is the mean weight applied to sources as a function of redshift — — times the number of sources per cluster. We have found including this weight changes the predicted masses by less than 1%. Nevertheless, our fiducial result includes this additional redshift weighting.
3.1.2 The Observed Richness–Mass Relation
Turning to the probability distribution , in addition to the stochastic nature of the relation between cluster richness and halo mass, the observed richness of a galaxy cluster is subject to projection effects. Indeed, there are now multiple lines in support of the existence of projection effects in the SDSS redMaPPer cluster catalog (Farahi et al., 2016; Zu et al., 2017; Busch & White, 2017; Sohn et al., 2017). Following Costanzi et al. (2018), we model as the convolution of two probability distributions:
(8) 
The first term inside the integral accounts for projection effects and observational noise in the richness estimates. The second term inside the integral accounts for the stochastic relation between halo mass and the intrinsic halo richness .
Below, we start describing the model adopted for the intrinsic richness–mass relation . The probability distribution was the focus of a detailed numerical study in a companion paper (Costanzi et al., 2018). A brief overview of that work is presented in the subsequent subsection.
The Intrinsic Richness–Mass Relation
Different parameterizations for have been proposed in the literature, typically assuming a lognormal distribution with the expectation value for the richness modeled as a power law (see e.g. Rozo et al., 2010; Mana et al., 2013; Murata et al., 2017). In this analysis we opt for a model which more closely resembles halo occupation distribution functions typically used to study galaxy clustering (Berlind & Weinberg, 2002; Bullock et al., 2002). A recent review of halo occupation modeling and other approaches to modeling the galaxy–halo connection can be found in Wechsler & Tinker (2018) (see also e.g. the discussion in section 4 of Jiang & van den Bosch 2016 and Appendix D of Reddick et al. 2013). The total richness of a halo of mass is given by , where is the number of central galaxies (either zero or unity), and is the number of satellite galaxies in the cluster (Kravtsov et al., 2004; Zheng et al., 2005). We model the expectation value of as a simple step function, for , and otherwise. While in practice these step functions have a finite width, we expect all clusters in our sample to have masses , so that the stepfunction approximation should be easily sufficient.
Turning to the satellite galaxy population, it has long been known that the scatter in the number of satellites is superPoissonian at large occupations numbers (BoylanKolchin et al., 2010), but is close to Poissonian otherwise. More recently, the number of satellites has been shown to be subPoissonian at very low occupancy (Mao et al., 2015; Jiang & van den Bosch, 2016). Since we are interested in galaxy clusters, we ignore this small deviation from Poisson statistics at low in our analysis. We add variance to a Poisson distribution by modeling as the convolution of a Poisson distribution with a Gaussian distribution. Operationally, this is equivalent to assuming the number of satellites galaxies in a halo of mass is a Poisson realization of some expectation value , where exhibits halotohalo scatter, e.g. due to formation history (Mao et al., 2015). We model the halotohalo scatter as a Gaussian with variance . This additional halotohalo scatter enables us to recover the superPoisson scatter in the halo occupation at large occupancy numbers.
In detail, the expectation value of the satellite contribution to is given by (Kravtsov et al., 2004; Zehavi et al., 2011)
(9) 
Here, is the minimum mass for a halo to form a central galaxy, while is the characteristic mass at which halos acquire one satellite galaxy. Our parameterization enforces when . Finally, the expectation value of the Gaussian component is set to zero, while the variance of the Gaussian term is set to .
The convolution of a Poisson distribution with a Gaussian distribution does not have an analytic closed form. However, we have found that a skewnormal distribution is an excellent approximation to the resulting distribution. Its model parameters — skewness and variance — depend on and , and are obtained by fitting the skewnormal model to the appropriate GaussianPoisson convolution (see Appendix B). By creating a lookup table for these parameters, we can avoid having to numerically compute the convolution of the Poisson and Gaussian distributions, significantly increasing the computational efficiency of our model. In Section 6.2 we discuss how our choice of parameterization impacts the cosmological constraints derived from the SDSS redMaPPer sample.
Modeling Observational Scatter in Richness Estimates
The scatter in the distribution is due to observational scatter, i.e. noise on the estimated richness values due to photometric noise, uncertainties in background subtraction, and projection effects. The latter refers to the boosting of the richness of a galaxy cluster due to additional galaxies along the line of sight that are mistakenly associated with the galaxy cluster. In Costanzi et al. (2018) we developed a formalism that quantitatively characterizes these effects, demonstrated the accuracy and precision of this formalism, and combined it with numerical simulations to calibrate the impact of projection effects and observational noise in the SDSS DR8 redMaPPer catalog. Here, we provide a succinct description of that work, and refer the reader to Costanzi et al. (2018) for details.
The scatter in due to background subtraction and photometric noise in the data is estimated as follows. First, we inject synthetic clusters of known richness and redshift at a random position in the survey footprint. We then measure the recovered richness, and repeat the procedure 10,000 times to characterize the associated scatter. The scatter is decomposed into a Gaussian component — the observational noise — and a decaying exponential component that characterizes the tail of high richness clusters due to projection effects. The biases in the recovered richness, as well as the standard deviation of the Gaussian component, are recorded. It is this Gaussian distribution that we use to model observational noise.
The tail due to projection effects is modeled explicitly as a sum over the true richnesses of galaxy clusters along the line of sight, weighted by a “projection kernel” that is calibrated directly from the data. The validity of this model is explicitly demonstrated by comparing the impact of projection effects around random points in the data to the impact of projection effects we predict when applying our model to body simulations (see also §4). Having validated the model, we measure the impact of projection effects at the location of massive halos in the body simulations, and use this calibrated measurement to characterize the tail due to projections arising from correlated large structures.
The end result of this procedure is a calibrated model for composed of a Gaussian observational noise, an exponentially decaying noise contribution due to projection effects, and a tail that extends to which accounts for incompleteness due to percolation effects: lowrichness clusters along the same line of sight of a richer objects can lose some of their galaxies to the richer object, leading to a severe underestimation for their richness. Figure 2 shows the distribution as a function of at redshift . The inset shows a crosssection of the distribution for . The three main components of the distribution are obvious by eye: a Gaussian kernel due to observational noise, a large tail to high richness due to projection effects, and a lowrichness tail due to percolation effects.
We note that the calibration of depends on the input cosmology and richness–mass relation parameters adopted to generate the synthetic cluster catalog. However, we verified in Costanzi et al. (2018) that this assumption has no impact on the posterior of the cosmological parameters. In section 6 we explicitly test the robustness of our cosmological constraints to different calibrations of .
3.2 MassFunction Systematics
The modeling of the cluster counts and mean cluster masses depends on the halo mass function. Here, we use the Tinker et al. (2008) halo mass function. Tinker et al. (2008) report their mass function formula to be accurate at the level, but we do not have a robust estimate of the associated systematic uncertainty as a function of mass. Moreover, a number of studies comparing different halo finders and fitting functions suggest a systematic uncertainty of the order of (e.g. Knebe et al., 2013; Hoffmann et al., 2015; Despali et al., 2016). To characterize the systematic uncertainty in the halo mass function in dark matter only simulations we introduce two nuisance parameters and relating the Tinker et al. (2008) mass function to the true mass function via
(10) 
where the pivot mass is set to . Note that if and , then the halo mass function is given by the Tinker et al. (2008) formula. We set the priors on and using the ensemble of simulations developed as part of the Aemulus project (DeRose et al., 2018). This is a set of 40 body simulations spanning a range of cosmologies in the redshift range . Each simulation box has a length and contains particles. The particle mass is cosmology dependent, and averages . The cosmologies sampled by the simulation spans the confidence interval from WMAP (Hinshaw et al., 2013) and Planck (Planck Collaboration et al., 2016a) in combination with Baryon Acoustic Oscillation data (Anderson et al., 2014) and the Union 2.1 SuperNova (SN) data (Suzuki et al., 2012). Halo catalogs were generated using the ROCKSTAR algorithm (Behroozi et al., 2013). Further details of the simulation data as well as the convergence tests done to ensure the reliability of these simulations are presented in DeRose et al. (2018).
We fit the simulation halo abundance data for the nuisance parameters and at each snapshot of each of the 40 simulations by computing the ratio . Next, we model the distribution of values as Gaussian, and fit for the mean values and covariance matrix describing the scatter of across all 320 snapshots. We find and , with a covariance matrix:
(11) 
The above matrix is the sum of the covariance matrix obtained from the fit of and , i.e. the statistical calibration uncertainty in the mean values of the parameters, and the covariance matrix describing the simulationtosimulation scatter in and . That is, the above covariance matrix fully accounts for any possible deviations of the Tinker et al. (2008) halo mass function to the test set of 40 Nbody simulations employed here. Note too the variance of corresponds to a 6% uncertainty in the amplitude of the halo mass function, consistent with the quoted precision in Tinker et al. (2008). The above covariance matrix and best fit values define the bivariate Gaussian priors used for the parameters and in our cosmological analysis. Future analyses should be benefit from the significantly higher precision that can be achieved using emulators (e.g. McClintock et al., 2018a).
The above analysis does not account for the impact of baryons on the halo mass function. Several recent works have estimated the impact of baryonic feedback on total masses of halos and, thereby, the mass function (e.g. Cui et al., 2014; Velliscig et al., 2014; Bocquet et al., 2016; Springel et al., 2017). These works all find that baryonic impact decreases with increasing radial aperture. For the specific mass definition we adopt, namely a 200 overdensity criterion relative to mean, the impact of AGNbased full physics is modest for massive () halos. In the IllustrisTNG simulation, full physics leads to a mean enhancement of while multiple methods analyzed by Cui et al. (2014) produce mainly decrements in of roughly similar magnitude. Due to the current uncertainty in modeling baryonic effects, we defer its inclusion into the error budget of the halo mass function to future work. For now, we simply note that a 3% systematic uncertainty in the halo mass function is already subdominant to the precision of the Tinker mass function, as characterized by our model parameters and .
3.3 Covariance Matrix
Having specified the expectation values for our observables (cf. Section 3.1) we need to define the covariance matrix of our data vector in order to fully specify the likelihood function. Here, we assume that the abundance and weak lensing data are uncorrelated. This assumption is well justified: the weak lensing error budget is strongly dominated by shape noise, and the dominant systematic is the overall multiplicative shear and photo bias of the source catalog. None of these errors affect the abundance data, and are therefore clearly uncorrelated.
Our Gaussian likelihood model takes the form
(12) 
where is the total covariance matrix detailed below, and and are respectively the data vectors (see Table 1) and the expectation values for the number counts and (Equation 2 and 5, respectively).
In reality, the likelihood for the abundance data is a convolution of a Poisson error on the counts and a Gaussian error due to density fluctuations within the survey area (e.g. Hu & Cohn, 2006; Takada & Spergel, 2014). Such a convolution does not have an analytic closed form. Here, we take care to ensure that all of our richness bins are well populated — our least populated richness bin contains over 200 galaxy clusters — so that the Poisson component can be adequately modeled with a Gaussian distribution. Consequently, our likelihood for the abundance data can be modeled as a Gaussian with a total covariance matrix having three distinct contributions:

A Poisson contribution due to the Poisson fluctuation in the number of halos at given mass in the survey volume.

A sample variance contribution due to the fluctuations of the density field in the survey volume.

A contribution due to uncertainty in the miscentering corrections detailed in section 2.2.
The Poisson and sample variance contributions to the covariance matrix are computed analytically at each step of the chain to properly account for their dependence on cosmology and model parameters. At high richness, the Poisson contribution dominates, with sample variance becoming increasingly important at low richness (Hu & Kravtsov, 2003). The analytical expression used to derive these two terms is provided in Appendix A, and it is validated by means of mock catalogs. In particular, we use the mock cluster catalogs generated by Costanzi et al. (2018) to calculate the variance in the abundance for each of the richness bins employed in our analysis. To do so, we split the full simulated survey area into patches, ranging from to , and use these patches to estimate the variance within the survey using a Jackknife method. Note that jackknife methods cannot adequately resolve modes that are larger than the jackknife region under consideration. Consequently, for the purposes of comparing our analytic estimates to the jackknife covariance matrix, we truncate our sample variance estimates by setting matter power spectrum to zero below the wavenumber , where is the radius of a spherical volume equal to the volume of the jackknife patch considered (see Eq. 23). As we can see from Figure 3, our analytic estimates are a good match to the observed variance in the simulations. We have also explicitly verified using our synthetic cluster catalog that the sample variance contribution to the mean cluster mass is negligible relative to the statistical and systematic uncertainties of the weak lensing analysis. We also find that the crosscorrelation between number counts and the mean weak lensing cluster mass obtained from our jackknife estimates is consistent with zero.
Turning to the uncertainty in the weak lensing mass estimates, we use the posteriors from the stacked weak lensing analysis described in Section 2.3 (see Table 1). Critically, these posteriors are found to be nearly Gaussian in the log, that is, the posterior is well described by a lognormal distribution. Consequently, we model the weak lensing likelihood as a Gaussian distribution in . The errors include not just the statistical uncertainty of the measurement, but also systematic uncertainties due to shear and photo biases (multiplicative shear bias), cluster projections, halo triaxiality, and miscentering effects. The overall shear and photo multiplicative bias is shared across all richness bins. Consequently, the variance associated with this error should be perfectly correlated across all bins.
To determine the overall mass error in each bin associated with the multiplicative bias we proceed as in Melchior et al. (2017, section 5.7) and McClintock et al. (2018a, section 5.6). First, we repeat our weak lensing analysis while fixing the systematic parameters. We then subtract the variance obtained from these fixedbias run to the total variance measured while floating the multiplicative bias parameters. The resulting variance is the variance attributable to the multiplicative bias. We do this for all five richness bins, and verify that, as expected, the variance due to multiplicative bias corrections to the profile, , is nearly richnessindependent, and roughly consistent with the expectation . The systematic uncertainty on the log of the masses is . This value is larger than the systematic quoted in Simet et al. (2017) because of the new contribution to the systematic error budget from blended sources, and the scatter correction. The error is perfectly correlated across all richness bins, so we set along the nondiagonal elements of the covariance matrix. This corresponds to the 6% total systematic error in Table 2. The diagonal entries of the covariance matrix are taken directly from the posteriors of our weak lensing analysis (see Simet et al. (2017) for details). Finally, as noted in Section 2.3, using numerical simulations we estimated that our fitting procedure introduces a bias that varies from 2% at low richness to 3% at high richness. We adopt a systematic uncertainty in this correction equal to half the magnitude of the correction, and add this (subdominant) contribution to arrive at the final covariance matrix employed for the weak lensing data.
In summary, our likelihood:

Is Gaussian in the abundances, and includes Poisson, sample variance, and miscentering uncertainties.

Is Lognormal in the weak lensing masses, as per the posteriors of Simet et al. (2017). It also accounts for the covariance due to shared multiplicative shear and photo biases, blended sources, and cluster triaxiality and projection effects. These systematics are assumed to be perfectly correlated across all richness bins.

Has no covariance between the abundance and weak lensing data. We expect this to be an excellent approximation.
The posteriors from our analysis are fully marginalized over all sources of systematic uncertainty described above.
4 Validation Tests
We validate our likelihood framework by placing cosmological constraints from a synthetic cluster catalog whose cosmology and richness–mass relation is known a priori. The mock data are generated starting from an body simulation run with dark matter particles in a box of comoving size . The code used is LGadget, a variant of Gadget (Springel, 2005). The simulation assumes a flatCDM model with , , , , and . Lightcone data, including a halo catalog down to , is constructed from the simulation on the fly. The halo finder is ROCKSTAR (Behroozi et al., 2013). For further details, see DeRose et al, in preparation; this is one realization of the L1 box described in that work.
We build a synthetic cluster catalog as follows. First, each halo is assigned a richness drawn from the distribution detailed in Section 3.1. We set our fiducial model parameters to: , , and . These fiducial values have been chosen by inverting the mass–richness relation of Simet et al. (2017) to arrive at parameters for the richness–mass relation (Eq. 9). The assigned richnesses are then modified to account for observational noise plus projection effects as described in Costanzi et al. (2018). Specifically:

Starting at the richest halo, we compute the projected richness of each halo via
(13) Here, the index runs over all clusters, is the fraction of the area of the th object which overlaps with the th cluster, and is the datacalibrated weight which accounts for the redshift distance between the cluster and .

A Gaussian noise is added to the projected richness to account for observational errors (including the effects of background subtraction).

Having calculated the observed projected richness of the richest galaxy cluster, we update the richnesses of all remaining clusters via
(14) That is, the richness of cluster is reduced by the amount it contributed to the richest cluster.

We now move on to the next richest cluster in the catalog, and iterate the above procedure to compute the observed projected richness of every halo in the simulation.

Finally, we assign a photometric redshift to every cluster based on their spectroscopic redshifts using the model in section 3.1.
Given our synthetic cluster catalog, we now compute our observable data vectors. The cluster counts are measured using the same binning scheme as for the real data (see Table 1). As for the synthetic weak lensing masses, these are obtained from the mean halo mass in each richness bin. We also add a lognormal noise to the mean weak lensing masses by drawing from the covariance matrix described in Section 3.3. To mimic the fact that the weak lensing masses are estimated assuming while the simulation use the fiducial value , we invert Equation 1 to arrive to our final mock data vector:
Figure 4 shows the constraints obtained by analyzing one realization of our synthetic data set with our pipeline. This analysis varies the same parameters and assumes the same priors as the real data analysis (see Section 6). We can see that our analysis successfully recovers the true values of the parameters of our synthetic data set (red lines in the triangle plots). Because this test relied on a single simulation, we could not use it to validate the width of our posteriors.
5 Blinding and Unblinding
In order to avoid confirmation biases our cosmological analysis was performed blinded. By “blinded”, we mean that the following protocols were followed:

The cosmological parameters in the MCMC were randomly displaced before being stored, and the random displacement was stored in binary (i.e. a nothumanreadable format).

All modeling choices — specifically which set of cosmological models and models for the scaling relations we would consider — were made before unblinding. In this work, we chose to focus exclusively on flat CDM cosmologies with massive neutrinos. Our choice to let neutrino mass vary follows the practice of the DES Year 1 combined probe analysis (DES Collaboration et al., 2017a).

The set of scaling relation models considered in our analysis was fully specified before unblinding. In addition to our fiducial model for the scaling relation, we considered three additional sets of models for the scaling relation (see Figure 8). "Randompointinjection" refers to the distribution calibrated by inserting galaxy clusters at random locations in the sky, i.e. neglecting the contribution from correlated structures to the observed richness. Since we know this model underestimates projection effects, we take half the difference between our fiducial model and the "Randompointinjection" model as our estimate of the systematic uncertainty in our cosmological posteriors. We further tested how allowing the intrinsic scatter of the richness–mass relation to scale with mass affected our cosmological posterior constraints. Finally, to compare to existing works, we considered both a powerlaw scaling relation with constant scatter and a scaling relation modeled after the work of Murata et al. (2017).

All priors were set before unblinding.

The set of external data sets we would compare against was decided before unblinding. Here, we distinguish between data sets with which we intend to combine our analysis and further test the flat CDM model (see next item), and other largescalestructure data sets that constrain the parameter. For the latter, we simply compare the recovered values, rather than test for consistency across the full parameter space. Specifically, we compare our constraints to those of Planck (Planck Collaboration et al., 2016a)^{9}^{9}9After we unblind our analysis new Planck results have been released (Planck Collaboration et al., 2018). To provide the most updated results we decided to consider the latest Planck constraints (hereafter Planck DR18) as CMB data set to compare against the SDSS clusters results., the DES Year 1 3x2 point analysis (DES Collaboration et al., 2017a), the KiDS450+2dFLenS analysis (Joudaki et al., 2018), and the KiDS450+GAMA analysis (van Uitert et al., 2018). In addition, we compare our constraint to several additional recent cluster abundance constraints, namely those from the PlanckSZ catalog (Planck Collaboration et al., 2016b), those from the South Pole Telescope (de Haan et al., 2016) and Atacama Cosmology Telescope (Hasselfield et al., 2013), and those from the Weighing the Giants team (Mantz et al., 2015). Finally, we compare our constraints to those derived from the SDSS maxBCG cluster sample (Rozo et al., 2010). For these purposes, we consider analyses and to be consistent if their central values of deviate by no more than , where .

The set of external data sets with which we intend to combine with the SDSS analysis (provided the two are deemed to be consistent, see metrics for consistency below) were fully specified. In addition to Planck (Planck Collaboration et al., 2016a), we decided to compare our constraints to those derived from Baryon Acoustic Oscillations (BAO) as measured in the 6dF galaxy survey (Beutler et al., 2011), the SDSS Data Release 7 Main Galaxy sample (Ross et al., 2015), and the BOSS Data Release 12 (Alam et al., 2017). As for the Planck constraints, we relied on the TT+lowP data set. With the expectation that our combined BAO and cluster abundance analysis (which included a prior on ) would result in a tight constraint on the Hubble constant (see e.g. DES Collaboration et al., 2017b), we also planned on comparing our posterior in to that of Riess et al. (2016).

The metrics for consistency with external data sets were selected before unblinding. Specifically, consistency between two data sets and was established by testing whether the hypothesis is acceptable (see method ‘3’ in Charnock et al., 2017). Here and are the model parameters of interest as constrained by data sets and respectively. The two data sets were deemed to be consistent if the point falls within the 99% confidence level of the multidimensional distribution of . If two data sets were found to be inconsistent with one another, we did not consider the combined analysis. We note that in order for this test to be the sharpest possible test, it is important to restrict oneself to parameter subspaces that are well constrained in both data sets. To that end, in all cases we restrict the parameter space for comparison to the set of parameters whose posterior is at most half as uncertain as the prior of each individual data set.

No comparison of our cosmological constraints to any other data sets were performed prior to unblinding.
The weak lensing analysis upon which our work relies was not performed blind (Simet et al., 2017), though our forthcoming analysis using data from the Dark Energy Survey data will have benefited from a blind weak lensing analysis. Importantly, all relevant weak lensing priors — specifically the multiplicative shear bias, photometric redshift correction, source dilution, etc — were finalized well before the advent of our particular analysis: no tuning of any input catalog was done in response to the weak lensing analysis of Simet et al. (2017) or our own cosmological analysis.
During the DES internal review process the question arose whether past experiences from people in our team — particularly those involved with the development of the maxBCG cluster catalog (Koester et al., 2007) — had in any way informed the development of redMaPPer, possibly leading to unintended biases in this way. The redMaPPer algorithm was not tuned based on previous experiences. The algorithm iteratively selftrains the colors of redsequence galaxies as a function of redshift without any cosmological assumptions. The parameters chosen for the richness estimation — particularly the aperture used — were selected by minimizing the scatter in Xray luminosity at fixed richness (Rykoff et al., 2012). We do adopt a fiducial cosmology when defining our richness measurements (this is necessary to turn physical apertures into angular apertures), but since the richness–mass relation is not known a priori, it carries no cosmological information.
Our unblinding protocol was defined by the set of requirements detailed below.

Our inference pipeline had to successfully recover the input cosmology in a synthetic data set. For details, see Section 4.

All SDSSonly chains (including alternative models) were run demanding the fulfillment of the GelmanRubin criteria (Gelman & Rubin, 1992) with being our convergence criteria.

We had to demonstrate systematics uncertainties in our model for did not appreciably bias our cosmological posteriors. Since the randompointinjection model for this distribution is an extreme model, we adopted half of the systematic shift in the values of the cosmological parameters between our fiducial model and this extreme model as our estimate of the associated systematic uncertainty. Note this definition implies that the extreme randompointinjection model is consistent with our fiducial model at despite being clearly extreme. We demanded that these systematic shifts be less than the corresponding statistical uncertainties.

We explicitly verified that the priors of all parameters which we expected would be wellconstrained a priori are not informative, that is the posteriors of such parameters did not run into the priors within the 95% confidence region. In case this condition was not fulfilled we planned to extend the relevant prior ranges until the requirement was met. Parameters that are prior dominated (i.e. their posterior runs into the prior) are , , , , , , , and . All of these were expected to be prior dominated a priori, and all prior ranges were purposely conservative. Of these, the two that might be most surprising to the reader are and , as these parameters help govern the richness–mass relation. However, notice that is the mass at which halos begin to host a single central galaxy; since our cluster sample is defined with the richness threshold , the mass regime of halos which host a single galaxy is not probed by our data set. Likewise, our data vector is comprised only of the mean mass of galaxy clusters in a given richness bin, a quantity that is essentially independent of the scatter in the richness–mass relation (see McClintock et al., 2018b).

The of the data for our bestfit model must be acceptable. To this end we considered the bestfit distribution recovered from 100 mock data realizations generated from the bestfit model of the data. We assumed these 100 trials were distributed according to a distribution, and fit for the effective number of degrees of freedom. The number of degrees of freedom is not obvious a priori due to the presence of priors in the analysis. We considered the of our data to be not acceptable if the probability to exceed the observed value was less than 1%, after marginalization over the posterior for the effective number of degrees of freedom. We emphasize that verifying an acceptable does not unblind the cosmological parameters. While our model did indeed have an acceptable (see section 6.1 for details), our plan was to revisit our model and covariance matrix estimation procedures in the case of an unacceptable value. This proved unnecessary.

Finally, this paper underwent internal review by the collaboration prior to unblinding. All members of the DES cluster working group, as well as all internal reviewers, agreed that our analysis was ready to unblind before we proceeded.
After all of our unblinding requirements were satisfied, we proceeded to unblind our analysis. No changes to the analyses were made post unblinding.
Parameter  Description  Prior  Posterior 

Mean matter density  
Amplitude of the primordial curvature perturbations  
Amplitude of the matter power spectrum  
Cluster normalization condition  
Minimum halo mass to form a central galaxy  
Characteristic halo mass to acquire one satellite galaxy  
Powerlaw index of the richness–mass relation  
Intrinsic scatter of the richness–mass relation  
Slope correction to the halo mass function  
Amplitude correction to the halo mass function  
Hubble rate  
Baryon density  
Energy density in massive neutrinos  
Spectral index 
6 Results
6.1 SDSS Cluster Abundances and Weak Lensing Data
We model the abundance of galaxy clusters and their weak lensing masses assuming a flat CDM cosmological model, allowing for massive neutrinos. The full set of cosmological parameters we consider is: , , , , , and . Neutrinos are included assuming three degenerate neutrino species. We adopt the same priors as in the DES Year 1 analysis of galaxy clustering and weak lensing (DES Collaboration et al., 2017a), with the exception of , where we adopt the slightly more restrictive prior . There are also two parameters ( and ) associated with systematic uncertainties in the halo mass function, and 4 parameters governing the richness–mass relation : , , , and . The priors for all parameters are summarized in Table 3. We have explicitly verified that increasing the range of the priors adopted for the richness–mass relation parameters does not adversely impact the recovered cosmological constraints.
The result of our MCMC fitting procedure is shown in Figure 5, while the marginalized posterior values are reported in Table 3. Parameters not shown in Figure 5 and without a quoted posterior in Table 3 are those whose posterior is equal to their prior. Also shown in the table is the posterior for the socalled cluster normalization condition parameter . In practice, the – degeneracy in our cosmology analysis corresponds to
(15) 
Nevertheless, unless otherwise specified in the text, from this point on we will focus on the cluster normalization condition as it has become a standard parameter in the literature.
A comparison of our bestfit model with the data is shown in Figure 1. The of our bestfit model is . To assess the goodness of the fit we generated mock data vectors from our bestfit model of the data and, for each of them, we recovered the bestfit value. Assuming a distribution for the trials we fit for the effective number of degree of freedom, finding (see Figure 6). This corresponds to a probability to exceed of . Thus, we conclude that our model provides a good fit to the data.
We wish to determine whether the error budget in the cosmological parameter is dominated by uncertainties in the abundance data or the weak lensing data. To determine this, we first compute the predicted abundance and weak lensing data using our bestfit model. We then run two additional chains using the predicted expectation values as a synthetic data vector. The key difference between the two chains is that for one we reduce the abundance covariance matrix by a factor of 100, while for the other chain we reduce the covariance matrix of the weak lensing data by a factor of 100. By comparing the cosmological posteriors for these two chains we can determine if there is one observable which dominates our error budget. In both cases the corresponding posterior on have an error bar (see left panels of Figure 7). Thus we conclude that our error budget on cosmological parameters is not dominated by either of the uncertainties associated with the two observables. Rather, both observables contribute in comparable ways to the total error budget of .
The balance between weak lensing errors and abundance uncertainties is surprising in light of the fact that all analyses to date have been dominated by uncertainties in the calibration of the mean of the observable–mass relation. Nevertheless, this feature of our results can be easily understood. The left panel of Figure 7 demonstrates that there is a strong degeneracy between and . Unlike previous analysis, which have incorporated wellmotivated simulationbased priors on the scatter of the observable–mass relation, our analysis adopts a very broad prior on the scatter of the richness–mass relation. This broad prior reflects the difficulty inherent to predicting properties of the richness–mass relation a priori. Since the scatter parameter impacts the detailed shape of the abundance function — larger scatter leads to flatter abundance functions — exquisitely precise measurements of the abundance function can break the degeneracy between scatter and , leading to significant improvements in the constraints. Conversely, even modest constraints of the scatter parameter can break the – degeneracy, leading to tighter constraints.
We demonstrate the impact that a modestlyprecise prior can have on our cosmological posteriors by redoing our analysis while imposing a flat prior . The corresponding posterior for is . If we now repeat our sensitivity analysis, and shrink the weak lensing mass errors, the width of the posterior decreases to , while decreasing the abundance errors while holding the weak lensing errors fixed has a negligible impact on the posterior. These trends are illustrated in the right panel of Figure 7. Evidently, external constraints on the scatter of the richness–mass relation of redMaPPer clusters are extremely valuable from a cosmological perspective.
6.2 Robustness to Assumptions About the Richness–Mass Relation
Figure 5 shows that there is strong covariance between cosmological parameters and parameters governing the richness–mass relation. Consequently, it is of critical importance to simultaneously fit for the richness–mass relation as part of the cosmological analysis. Likewise, one may ask to what extent are our cosmological constraints sensitive to the detailed assumptions we have made about the richness–mass relation. To address this question, we have repeated our analysis with a range of richness–mass relation models as summarized in Figure 8. The models considered are:

A model that allows for the intrinsic scatter to vary with mass via
(16) 
A model that neglects the perturbations on the observed richness due to correlated structures. To this end we set equal to the probability distributions recovered from injecting synthetic clusters at random positions in the survey mask. This calibration provides a very strict lower limit on the scatter of due to projection effects: clusters do contain correlated large scale structure. The difference in the posteriors of the cosmological parameters between our fiducial run and the randompointinjection model places a strict upper limit on the systematic associated with our modeling of projection effects.

A model in which the richness–mass relation is a simple power law — — and is a lognormal distribution in which the total scatter it the sum of a Poissonlike term and an intrinsic scatter term — .

The richness–mass relation model of Murata et al. (2017). This model assumes is lognormal, the mean richness–mass relation is a powerlaw, and the intrinsic scatter is massdependent, and given by . According to Murata et al. (2017) all integrals are truncated at . Reassuringly, when we mirror the analysis of Murata et al. (2017) and fix our cosmological parameters to Planck Collaboration et al. (2016a) we reproduce their results despite significant methodological differences in how the weak lensing data is incorporated into the likelihood.
As can been seen in Figure 8 our cosmological posteriors are all consistent with one another. It is clear that the more restrictive parameterizations (e.g. lognormal + powerlaw) result in somewhat tighter constraints. Notably, our standard result — which we believe is the most appropriate model — results in the most conservative posteriors. In particular, we see that opening up the freedom of a mass–dependent intrinsic scatter does not negatively impact the posterior on . Interestingly, we noticed that our randompoint injection model, which grossly underestimates the impact of projection effects, had a small impact on the posterior of the intrinsic scatter, modifying instead the best fit value for the slope of the richness–mass relation, a degeneracy that should be easily broken via multiwavelength analyses of the redMaPPer clusters. Finally, we note that the Muratalike parameterization has the largest impact on our posteriors, with the systematic shift in being comparable to the width of the posterior. We will address the origin of this shift in the next section.
6.3 The Observable–Mass Relation for redMaPPer Clusters
The left panel of Figure 9 shows the observed richness–mass relation for our fiducial model at the mean sample redshift . The error bars reflect the posterior of the mean relation at each mass. These are computed as follows: for each point in our chains we evaluate along a grid of masses. The mean and variance of across the chain are then recorded. In the left panel of Figure 9 we use these quantities to plot the 68% confidence interval for the posterior of the mean of the richness–mass relation. The central panel of Figure 9 is computed in a similar way, only now we show the posterior for the scatter .^{10}^{10}10A reader might find useful to have simple powerlaw fits to the data shown in Figure 9. We provide such fits below: The fits correspond to the best fitrelations. No errors are provided since these are meant to be “quicklook” references. Detailed quantitative analyses should rely on the full posterior of our model, which will be made available at http://risa.stanford.edu/redmapper/ when the paper is published.
For comparison we include in the two panels the richness–mass relation and scatter derived in Murata et al. (2017), who analyzed this sample of clusters using the same weak lensing shear catalog we employed. There are significant methodological differences between the two analyses. Specifically, Murata et al. (2017):

Uses emulators instead of an analytic model for the weak lensing profile of clusters, effectively holding the concentration of the galaxy clusters fixed. They also place no priors on the miscentering parameters.

Adopts a lognormal model for .

Adopts a powerlaw relation for both the mean and variance of at fixed mass.
The yellow bands in Figure 9 are obtained by estimating the 68% confidence region of their posterior from 10,000 Gaussian random draws of the richness–mass relation parameters quoted in Murata et al. (2017). Despite methodological differences, the two results are remarkably similar over the mass range probed by the survey (, see below). It is especially impressive how well the scatter found in Murata et al. (2017) agrees with our finding, confirming their argument that the masstrend in the scatter of the richness–mass relation they recovered is due to contamination from projection effects.
There is, however, one notable difference between our results and those of Murata et al. (2017). Figure 7 of Murata et al. (2017) shows that the mass distribution for clusters in the richness bin extends to masses as low as , the mass cut imposed in that analysis. Figure 10 shows the posterior of the mass distribution for clusters in our analysis, as labeled. Unlike Murata et al. (2017), we do not see any evidence for a population of lowmass () clusters. We believe the large number of lowmass halos in the Murata et al. (2017) analysis is driven by the combination of a lognormal model whose scatter increases with decreasing mass, and a model that ignores the central/satellite dichotomy at low masses. Our model avoids this problem by: 1) enforcing the appropriate Poisson limit in the limit of low ; and 2) developing a simulationbased model for projection effects that adequately characterizes nonGaussian tails in the distribution .
As an independent check of our conclusions, we estimate the halo masses of individual redMaPPer clusters using the stellar content of the redMaPPer central galaxy. Specifically, we fit the photometric SDSS data using a stellar population synthesis model to derive the stellar mass of each of the assigned redMaPPer central galaxies (Moustakas et al., in preparation). We then use the UniverseMachine algorithm (Behroozi et al., 2018) to determine the relation between the stellar mass of the central galaxy of a halo and the mass of the parent halo. Using the relation between halo mass and the stellar mass of the central galaxy, we can readily assign a mass estimate to each redMaPPer cluster. We find that the stellar mass estimates of 95% of redMaPPer central galaxies correspond to a halo mass of or higher. While the stellarmass to halomass relation of central galaxies is relatively uncertain, this result disfavors the existence of a significant population of lowmass () redMaPPer clusters, in agreement with our results.
We suspect that the tail of lowmass halos recovered in the Murata et al. (2017) model is responsible for the shift in seen in Figure 8 when adopting a powerlaw lognormal model for the richness–mass relation (i.e. the Muratalike analysis): the artificial boost in the abundance of low richness clusters is compensated by a decrease in the predicted halo abundance, which is in turn achieved by lowering the cluster normalization condition.
Finally, the right panel of Figure 9 shows the massselection function of the galaxy clusters selected in our experiment. That is, it shows the probability that a halo of mass at redshift is included in the SDSS redMaPPer sample. This probability is given by
(17) 
The probability is evaluated at a grid of masses for each point in the chain, and the corresponding mean and uncertainty is calculated. Using linear interpolation over this grid, we find that the mass value for which the detection probability is is .
Using our bestfit cosmological model we can combine the halo mass function with our recovered richness–mass relation to arrive at our bestfit mass–richness relation. In particular, for each point in the chain we can readily compute along a grid of richnesses, and calculate mean and variance of these quantities at each value as we sample our posterior. We estimated the mean mass at a grid of values, and computed the corresponding covariance matrix, and then fit the data with a powerlaw to arrive at the corresponding mass–richness relation. The resulting mass–richness relation is shown in Figure 11 along with the relation derived in Simet et al. (2017). The posteriors for the mass–richness relation in our analysis is
(18) 
where we used only data points (black dots in the figure) to fit the powerlaw relation. This is to be compared to the Simet et al. (2017) relation, . Note that the exponent has been obtained correcting the bestfit value derived in Simet et al. (2017) assuming via Equation 1 using the mean value derived in this analysis. Moreover, we expect a small difference in amplitude due to the updates to our multiplicative shear bias model and the redMaPPer centering fraction. All together, these corrections should boost the amplitude of the mass–richness relation by , in good agreement with our results.
7 Comparison to External Data Sets
Data sets  

Planck15  
BAO  
SDSS 