COSMOLOGICAL CONSTRAINTS FROM GALAXY CLUSTERS IN THE 2500 SQUARE-DEGREE SPT-SZ SURVEY

The SPT is a 10 m telescope located at the geographic South Pole. With a 1 degree field of view and ∼1 arcmin resolution, it was designed to rapidly map large areas of sky while being well matched to the angular size of high-redshift clusters. The SPT-SZ camera operated from 2007 through 2011 and consisted of a 960-element, photon-noise-limited bolometer array observing in three frequency bands centered at 95, 150, and 220 GHz, though this work uses only the first two bands. The observation strategy and analysis for the 2500 square-degree SPT-SZ survey are described in many previous SPT papers (e.g., Schaffer et al. 2011), and the analysis specific to obtaining the cluster sample is described in B15. Briefly, the majority of the 2500 square-degree survey was performed using the following scan strategy. The telescope is moved using a right-going scan in azimuth, followed by a left-going scan, followed by a step in elevation. This pattern is repeated until a ≳100 square degree patch of sky is observed. One such iteration is termed an observation, whereas the patch of sky is termed a field. In the only exception to this observing strategy, for two thirds of the data taken on one field (ra21hdec-50), the scan pattern instead consisted of scans in elevation at a series of fixed azimuth positions. Each field was observed until a depth of ≲$18\,\mu {\rm{K}}$-arcmin at 150 GHz was reached. Table 1 lists the 19 fields that comprise the full survey. For more information on the properties of the fields, see Table 1 of B15.

Table 1.  The 19 SPT-SZ Fields with the Simulation-derived Inputs to the Cosmological Analysis of the 2500 Square-degree Survey

Name	${\alpha }_{\mathrm{field}}$	${\beta }_{\mathrm{field}}$	${\gamma }_{\mathrm{field}}$
ra5h30dec-55	16.79	4.60	1.33
ra23h30dec-55	17.58	4.03	1.39
ra21hdec-60	25.64	4.07	1.29
ra3h30dec-60	20.53	4.70	1.25
ra21hdec-50	25.28	4.14	1.11
ra4h10dec-50	16.75	5.48	1.27
ra0h50dec-50	20.76	5.11	1.14
ra2h30dec-50	14.98	4.78	1.19
ra1hdec-60	17.25	5.38	1.18
ra5h30dec-45	15.91	4.81	1.08
ra6h30dec-55	17.77	4.58	1.16
ra3h30dec-42.5	16.85	4.31	1.20
ra23hdec-62.5	14.90	4.92	1.18
ra21hdec-42.5	17.11	4.49	1.15
ra1hdec-42.5	18.41	5.55	1.19
ra22h30dec-55	16.45	5.23	1.13
ra23hdec-45	17.00	5.20	1.19
ra6h30dec-45	14.78	4.23	1.16
ra6hdec-62.5	16.53	4.70	1.18

Note. The parameter ${\alpha }_{\mathrm{field}}$ describes the number of false detections expected above $\xi =5$ scaled to 2500 square degrees, while ${\beta }_{\mathrm{field}}$ describes the scaling of the number of false detections with ξ as defined in Equation (3). ${\gamma }_{\mathrm{field}}$ describes the renormalization of A_SZ in the SZ–mass scaling relation of Equation (2).

Download table as:  ASCIITypeset image

Two-dimensional maps of the sky are made by binning the time-ordered data (TOD), to which mild time-domain filtering has been applied, into 0.25 arcmin pixels. This produces estimates of the 95 and 150 GHz sky, to which we apply a simultaneous spatial-spectral filter, yielding a filtered estimate of the tSZ sky, optimized for extracting cluster candidates. The candidates are identified using a peak-finding algorithm.

Due to the exact choice of field extent, as well as the finite footprint of the SPT-SZ array on the sky, the fields overlap slightly. The overlap regions only comprise $\sim 2 \%$ of the total survey area, and for simplicity we choose to treat the fields as fully independent. As a result of this treatment, we would expect to double-count roughly 2% of cluster candidates, or seven to eight of our total number of candidates. In fact we find nine cases of candidates with a very high probability of having been detected in two different fields (i.e., 18 total candidates that correspond to nine physical clusters). There is no bias to our cosmological constraints from this treatment. Our uncertainties will be very slightly underestimated, but this has a negligible effect on our final cosmological constraints.

2.1.1. SZ–Mass Scaling Relation Parameterization

We use the same functional form for the SZ–mass scaling relation as previous SPT cluster cosmology analyses (V10; Benson et al. 2013; hereafter B13; Reichardt et al. 2013; hereafter R13; Bocquet et al. 2015).

Briefly, we introduce two SZ parameters related to the cluster detection process. The first is the detection significance ξ. After generating filtered synthesized SZ maps at a series of filter scales, each map pixel value is divided by the rms of the map in a strip that spans 90 arcmin in declination. Then the maximum peak height over all filter scales is identified and defined as ξ, the SZ observable used in this paper. Because this observable allows for a very well understood selection function, and cosmological constraints are dominated by the unknown normalization of the observable-mass scaling relation, we find that this observable is preferred over going to a different SZ observable such as Y_SZ.

The second parameter is the unbiased significance ζ. It is defined as the value of ξ that would be found in the absence of instrumental noise and astrophysical contaminants (including the SZ background). Due to the fact that ξ is determined by maximizing the significance after searching in two-dimensional position space and source template size, the average ξ found across many noise realizations is enhanced by those three degrees of freedom, resulting in the approximate relation

$$\begin{eqnarray}&&{\langle \xi \rangle }^{2}={\zeta }^{2}+3.\end{eqnarray} \tag{ 1 }$$

Since this relation cannot hold down to very low ζ, we only model this maximization bias for $\zeta \gt 2$. We find that changing the location of this cutoff to $\zeta \gt 1.5$ or $\zeta \gt 2.5$ has negligible impact on the results presented in this work. Motivated by the definition of ξ, and the fact that the astrophysical contaminants and instrument noise are Gaussian to a high degree, we model ξ as related to $\langle \xi \rangle$ by unit-width Gaussian scatter. We then parameterize the ζ–M scaling relation as

$$\begin{eqnarray}&&\zeta ={A}_{\mathrm{SZ}}{\left(\displaystyle \frac{M}{3\times {10}^{14}{M}_{\odot }{h}^{-1}}\right)}^{{B}_{\mathrm{SZ}}}{\left(\displaystyle \frac{E(z)}{E(0.6)}\right)}^{{C}_{\mathrm{SZ}}},\end{eqnarray} \tag{ 2 }$$

with an additional log-normal intrinsic scatter parameter ${\sigma }_{\mathrm{ln}\zeta }$. The dimensionless Hubble parameter is denoted by E(z). We follow B13 and R13 and apply Gaussian priors to these scaling relation parameters. However, the mean values and widths of the priors are updated to reflect results from the more recent cosmo-OWLS hydrodynamic simulations (Le Brun et al. 2014), which will be discussed in Section 2.2. Specifically, we use the Gaussian priors ${A}_{\mathrm{SZ}}=5.38\pm 1.61$, ${B}_{\mathrm{SZ}}=1.340\pm 0.268$, ${C}_{\mathrm{SZ}}=0.49\pm 0.49$, and ${\sigma }_{\mathrm{ln}\zeta }=0.13\pm 0.13$, while requiring ${\sigma }_{\mathrm{ln}\zeta }\gt 0.05$.

Each of the 19 fields that comprises the 2500 square degree survey has slightly different properties, which we account for using two similar types of simulations. The first is to account for the fact that the map noise is slightly different in each of the fields. This causes the detection significance ξ (the selection observable) of a cluster to relate to true underlying mass differently depending on the field in which the cluster was found. We model this by simulating the relation between true cluster mass and unbiased significance ζ separately for each field. The simulations contain several components.

• 1.  
  A Gaussian random field with a power spectrum equal to the CMB power spectrum calculated using the best-fit ΛCDM parameters from Keisler et al. (2011).
• 2.  
  A Gaussian random field meant to approximate the background of emissive point sources after the brightest sources are masked. The power spectrum is generated using three components. The first is a Poisson component modeling the radio point source population with an amplitude at ${\ell }=3000$ of ${D}_{3000}^{r}=1.28\mu {{{\rm{K}}}_{\mathrm{CMB}}}^{2}$ at 150 GHz and spectral index ${\alpha }_{r}=-0.6$ (defined by flux $\propto {\nu }^{\alpha }$). The second and third components model the Poisson and clustered dusty star-forming galaxy (DSFG) populations, respectively. The assumed spectral index is ${\alpha }_{\mathrm{DSFG}}=3.6$ and the amplitudes are the best-fit values from Shirokoff et al. (2011).
• 3.  
  Atmospheric and instrumental noise. These are simulated using the actual TOD. The coherent nature of the sky signal, modulated by the scan strategy, and the incoherent nature of the instrumental and atmospheric noise contributions allow us to estimate the map noise by simply subtracting the right-going scan TOD from the left-going scan TOD (up- and down-going in the case of most of the ra21hdec-50 field) in the map-making process. This removes any coherent signal present on the sky. We furthermore apply a sign change to half the observations before they are coadded into the final map. By changing which signs are assigned, we generate different realizations of realistic noise maps.
• 4.  
  Maps of the Compton y parameter from the cosmo-OWLS simulations in Le Brun et al. (2014). We use maps from the AGN8.0 model and convert the Comptonization maps to CMB temperature units by integrating the measured frequency response of the instrument against a relativistic 5 keV tSZ spectrum.

A modified version of the standard cluster-finding process is run on the sum of these four components, and the cluster candidates are recorded. In addition, the same spatial-spectral filter is applied to the maps only containing SZ signal. The latter is used to identify cluster candidates and define their SZ center. The amplitude in these filtered, SZ-only maps, divided by the standard deviation in a 90 arcmin strip in the filtered maps that include noise, is precisely the unbiased significance ζ introduced in Section 2.1.1. These candidates are then compared to the underlying halo catalog used to generate the mock SZ maps, which provide the cluster redshift. For each SZ-identified cluster candidate, we choose the nearest dark matter halo with ${M}_{200}\gt 5\times {10}^{13}{M}_{\odot }$.

With mass, unbiased significance and redshift in hand, we use the method of least absolute deviations to fit the ζ–M relation from Equation (2) over the range $z\gt 0.25$, ${M}_{500}\gt 1\times {10}^{14}{M}_{\odot }$. We find nearly identical results when using a linear least-squares method. The normalization parameter A_SZ differs from field to field by up to 30%, while the other scaling relation parameters only vary at the few percent level. We therefore model field variation with one parameter ${\gamma }_{\mathrm{field}}$ that renormalizes the overall scaling: ${A}_{\mathrm{SZ}}\to {\gamma }_{\mathrm{field}}{A}_{\mathrm{SZ}}$. For consistency with previous publications, we normalize ${\gamma }_{\mathrm{field}}$ so that the weighted average over the three fields which share 100% of the raw data used in R13 matches the weighted average in that work. We note that any other choice of normalization would simply alter the definition of A_SZ, leaving all other results unaffected. The values of ${\gamma }_{\mathrm{field}}$ are shown in Table 1.

The other field-specific simulation is used to compute the expected false detection rate. At a detection threshold of $\xi \gt 5$, approximately 5% of cluster candidates are expected to be false. In order to simulate the rate of false detections, we perform the cluster-finding process on the simulated maps described above, omitting the SZ component. The resulting number of cluster candidates as a function of ξ is recorded for 100 such simulations. To reduce shot noise from the finite number of simulated realizations, we model the false detection rate with the empirically chosen fitting function

$$\begin{eqnarray}&&{N}_{\mathrm{false}}(\gt \xi )={\alpha }_{\mathrm{field}}\exp (-{\beta }_{\mathrm{field}}(\xi -5)).\end{eqnarray} \tag{ 3 }$$

The values of ${\alpha }_{\mathrm{field}}$ and ${\beta }_{\mathrm{field}}$ are shown in Table 1. The total number of expected false detections per 2500 square degrees with $\xi \gt 5$ from these simulations is 18 ± 4, which is consistent with the number of optically unconfirmed cluster candidates, 21 (B15).

Following V09, B13, R13, and Bocquet et al. (2015), we use Y_X as a proxy for the total cluster mass. We write the scaling relation as

$$\begin{eqnarray}\displaystyle \frac{{M}_{500}}{{10}^{14}{M}_{\odot }/h} & = & \left({A}_{{\rm{X}}}{h}^{3/2}{\left(\displaystyle \frac{h}{0.72}\right)}^{\tfrac{5}{2}{B}_{{\rm{X}}}-\displaystyle \frac{3}{2}}\right)\\ & & \times {\left(\displaystyle \frac{{Y}_{{\rm{X}}}}{3\times {10}^{14}{M}_{\odot }\mathrm{keV}}\right)}^{{B}_{{\rm{X}}}}E{(z)}^{{C}_{{\rm{X}}}},\end{eqnarray} \tag{ 4 }$$

where the parameters A_X, B_X, C_X describe the normalization, mass dependence and redshift dependence of the relation, respectively.⁵⁰ As we did for the ζ–M relation, we introduce a parameter ${\sigma }_{\mathrm{lnYx}}$ which models a log-normal intrinsic scatter in the ${Y}_{{\rm{X}}}\mbox{--}M$ relation. We allow the intrinsic scatter in the ${Y}_{{\rm{X}}}\mbox{--}M$ relation and the ζ–M relation to be correlated, parameterizing this with the correlation coefficient ${\rho }_{\zeta ,\mathrm{Yx}}$. We follow B13 and R13 and apply Gaussian priors of ${B}_{{\rm{X}}}=0.57\pm 0.03$, ${C}_{{\rm{X}}}=-0.40\pm 0.20$, and ${\sigma }_{\mathrm{lnYx}}=0.12\pm 0.08$, respectively, as well as a uniform prior between −0.98 and 0.98 on ${\rho }_{\zeta ,\mathrm{Yx}}$.

In B13 and R13, we used a prior on the normalization of the ${Y}_{{\rm{X}}}\mbox{--}M$ relation, A_X, motivated from V09. In V09, the normalization was cross-checked against WL mass estimates from Hoekstra (2007) for the 10 clusters at $z\lt 0.3$ that, at the time, also had sufficient Chandra observations to measure Y_X. In B13, this normalization was remeasured using more recent WL mass estimates from Hoekstra et al. (2012), and found to be consistent with the assumed calibration from V09. In this work, we revisit the ${Y}_{{\rm{X}}}\mbox{--}M$ normalization, using the most recent WL mass estimates from H15 to constrain any mean offset in the normalization of the original hydrostatic mass calibration of the Y_X-mass scaling relation from V09. In addition, we also consider any systematic offset in the H15 WL mass estimates, by comparing to alternative WL mass estimates from WtG.

In Figure 1, we compare WL-based mass estimates from H15 to Y_X-based mass estimates and to weak-lensing masses from WtG. We have remeasured Y_X from archival Chandra X-ray data for the 14 clusters from H15 with sufficiently deep X-ray measurements, which we use to estimate a cluster mass (${M}_{\mathrm{Yx}}$) using the Y_X-mass scaling relation from V09. For these clusters, we also remeasure the deprojected weak-lensing aperture mass from the H15 data set using the X-ray implied ${r}_{500}$. Using the X-ray implied ${r}_{500}$ approximates what would happen if the WL mass estimates were included in a joint Y_X–WL scaling relation fit because of the relatively low scatter of Y_X with cluster mass and, as we will show later, the agreement in the resulting ${Y}_{{\rm{X}}}\mbox{--}M$ relation. To compare H15 and WtG mass estimates at the same radius, we extract a WtG mass for each of the 18 clusters in common at the value of ${r}_{500}$ implied by the H15 analysis. We specifically use the mass estimates computed by fitting an NFW model to the observed shear profiles from the WtG "color-cut" method (Applegate et al. 2014).

Zoom In Zoom Out Reset image size

In Table 2, we give three different measures of the relative cluster mass estimates between H15, WtG, and V09. First, we fit a scaling relation of the form $\mathrm{ln}({M}_{\mathrm{WL}})=A+\mathrm{ln}({M}_{\mathrm{Yx}})$ with a free log-normal intrinsic scatter, using the Bayesian linear regression fitting code from Kelly (2007). For comparison, we also estimate the bootstrap mean and median of the log of the mass ratio, which makes no assumption about the underlying scatter.

First considering the mass ratio between H15 and V09, we find all three estimates imply a mass ratio near unity with the scaling relation, bootstrap mean, and bootstrap median, giving ratios of 1.01 ± 0.07, 1.03 ± 0.06, and 1.12 ± 0.05, respectively. This implies that the Y_X-based masses from V09 are in relatively good agreement with H15, with the H15 masses ∼1%–12% larger, depending on how the mass ratio is estimated. Next, we consider the mass ratio between WtG and H15. We also find that all three estimates imply a mass ratio near unity with the scaling relation, bootstrap mean, and bootstrap median giving ratios of 1.06 ± 0.07, 1.07 ± 0.07, and 1.11 ± 0.05, respectively. Therefore the WL-based mass estimates from H15 and WtG are also in relatively good agreement, with the WtG masses ∼6%–11% larger, comparable to the 7% estimated systematic uncertainty from H15. The overlapping cluster sample between our V09-Y_X sample and WtG is smaller (eight clusters), and therefore has larger uncertainty on the mass ratio. We find that the three different mass ratio estimates between WtG and V09 range from 1.15–1.24. This range is somewhat larger than, but comparable with, the expectation when considering the mass ratios of the larger samples used to calculate the WtG/H15 and H15/V09 mass ratios, which would predict a mass ratio between WtG and V09 of ∼1.1–1.2.

In summary, comparing to the original ${Y}_{{\rm{X}}}\mbox{--}M$ calibration from V09, we find that the WL measurements imply a normalization that increases the Y_X-based cluster masses by a factor between 1.0 and 1.2, depending on the WL analysis and data set (H15, WtG), and the fit assumptions. We therefore choose a Gaussian prior on the normalization of our ${Y}_{{\rm{X}}}\mbox{--}M$ relation of ${A}_{{\rm{X}}}=6.35\pm 0.61$, which increases the normalization from V09 by a factor of 1.1 with an overall uncertainty of 10%. We expect that this uncertainty includes the sum of our statistical and systematic uncertainty, because it brackets the range of values in the above comparison with the WL measurements and is comparable to our naive expectation given the statistical uncertainty in the fit above and the systematic uncertainty estimated in H15. This effectively gives us a purely WL-derived mass scale, independent of the hydrostatic mass estimates. We expect our constraints on the mass scale to improve in future work using WL observations of SPT-SZ survey clusters, which we will then be able to directly include in our cosmological analysis.

Finally, we have also estimated the dependence of the weak-lensing and X-ray-derived mass estimates on cosmological parameters. Over the range of cosmological parameters explored, the ratio of the two varies negligibly compared to the uncertainty of 10%.

For a given set of cosmological parameters, we use the Code for Anisotropies in the Microwave Background (CAMB, Lewis et al. 2000) to generate the matter power spectrum as a function of redshift. We use this power spectrum as input to the Tinker et al. (2008) halo mass function, with which we calculate the number of dark matter halos as a function of spherical overdensity mass M₅₀₀ and redshift.

More recent results (e.g., Bhattacharya et al. 2011; Skillman et al. 2014) confirm that, over the range of cluster masses and redshifts considered in this work, the Tinker mass function is accurate to a level that is more than sufficient for the cosmological parameter constraints presented. For instance, even a shift as large as 10% in the halo mass function at ${M}_{500}\sim 3\times {10}^{14}\,{M}_{\odot }$ would only affect the cosmological constraints presented in Section 6.2 at the $\lesssim 0.2\sigma$ level.

Given the Tinker mass function, we relate the cluster mass to the observed quantities ξ and Y_X using the scaling relations described in Sections 2.1.1 and 4.1.

In our default cosmological analysis, we use the standard six-parameter ΛCDM model, parameterized by the physical densities of baryons and cold dark matter at z = 0, ${{\rm{\Omega }}}_{b}{h}^{2}$ and ${{\rm{\Omega }}}_{c}{h}^{2}$, the angular scale of the sound horizon at last scattering ${\theta }_{s}$, the tilt of the scalar power spectrum n_s, the amplitude of the scalar power spectrum A_s, and the optical depth to reionization τ. We will often refer to values of parameters derived from combinations of these original six, such as ${\sigma }_{8}$, the amplitude of linear matter fluctuations on 8 Mpc/h scales at z = 0, ${{\rm{\Omega }}}_{m}$, the total matter density at z = 0, and ${H}_{0}=100h\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$, the value of the Hubble parameter at z = 0. We explore three extensions to ΛCDM. We first explore a model in which the sum of the neutrino masses ${\rm{\Sigma }}{m}_{\nu }$ is a free parameter. In this work we use the prescription from Costanzi et al. (2013) to predict the effect of massive neutrinos on the halo mass function. We also consider a cosmological model in which ${\rm{\Sigma }}{m}_{\nu }$ and ${N}_{\mathrm{eff}}$, the effective number of relativistic particle species, are both free parameters, and a model in which w, the parameter describing the equation of state of dark energy, is allowed to vary.

To compare theoretical predictions for cluster abundance to the observations in this work, we closely follow the likelihood derivation of V10, B13, and R13. We differ from those analyses in this work by presenting an efficient numerical technique that scales linearly rather than exponentially with the number of mass proxies.

As is appropriate for a cluster abundance measurement (Hu & Kravtsov 2003; Holder 2006), we start with the binned Poisson statistic in cluster observable space, which consists of redshift z and an arbirary number of mass proxies. In this work, we use the SZ mass proxy ξ which was described in Section 2.1.1, and the X-ray mass proxy Y_X as described in Section 4.1. For a given observable-space bin x_i, the probability of observing n events, with expectation value $y({x}_{i})$ is

$$\begin{eqnarray}&&{P}_{i}=\displaystyle \frac{{e}^{-y({x}_{i})}y{({x}_{i})}^{n}}{n!}.\end{eqnarray} \tag{ 5 }$$

We choose to take the limit of small bins where y becomes arbitrarily small at all x_i, and n is zero except at the locations of x_i where clusters have been observed (Mantz et al. 2008, 2010, V09). If we let x_j denote the bin that contains the jth cluster, we obtain

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=\displaystyle \sum _{i}\mathrm{ln}{P}_{i}=-\displaystyle \sum _{i}y({x}_{i})+\displaystyle \sum _{j}\mathrm{ln}y({x}_{j}),\end{eqnarray} \tag{ 6 }$$

where i runs over all observable space bins, and j runs over detected clusters. In this work, we are dealing specifically with the three-dimensional cluster observable space $\{z,\xi ,{Y}_{{\rm{X}}}\}$, so we can write the model expectation value $y({x}_{i})=N({z}_{i},{\xi }_{i},{Y}_{{\rm{X}},i})$ and express the likelihood as

$$\begin{eqnarray}\mathrm{ln}{ \mathcal L } & = & -\displaystyle \sum _{{i}_{1},{i}_{2},{i}_{3}}N({z}_{{i}_{1}},{\xi }_{{i}_{2}},{Y}_{{\rm{X}},{i}_{3}})\\ & & +\displaystyle \sum _{j}\mathrm{ln}N({z}_{j},{\xi }_{j},{Y}_{{\rm{X}},j}),\end{eqnarray} \tag{ 7 }$$

where the ${i}_{1},{i}_{2},{i}_{3}$ sums again run over all possible values and j only runs over the bins where clusters were detected. Going to the continuous limit would result in a divergent likelihood expression. This can be understood by the fact that the model is increasingly unlikely to produce our particular realization of the cluster catalog as we go to finer binning. The divergence can be removed by simply adding $-\mathrm{ln}{\rm{\Delta }}z{\rm{\Delta }}\xi {\rm{\Delta }}{Y}_{{\rm{X}}}$ to Equation (7). This quantity (the logarithm of the observable-space bin volume) depends only on the choice of bin size, so that ${\rm{\Delta }}\mathrm{ln}{ \mathcal L }$ for different values of cosmological or scaling relation parameters remains meaningful. This procedure corresponds to substituting the number density ${dN}(z,\xi ,{Y}_{{\rm{X}}})/{dzd}\xi {{dY}}_{{\rm{X}}}$ for the number contained in a bin $N({z}_{j},{\xi }_{j},{Y}_{{\rm{X}},j})$. We then obtain the expression

$$\begin{eqnarray}\mathrm{ln}{ \mathcal L } & = & -\,\displaystyle \int {dzd}\xi {{dY}}_{{\rm{X}}}\displaystyle \frac{{dN}(z,\xi ,{Y}_{{\rm{X}}})}{{dzd}\xi {{dY}}_{{\rm{X}}}}\\ & & +\displaystyle \sum _{j}\mathrm{ln}\displaystyle \frac{{dN}({z}_{j},{\xi }_{j},{Y}_{{\rm{X}},j})}{{dzd}\xi {{dY}}_{{\rm{X}}}}.\end{eqnarray} \tag{ 8 }$$

V10, B13, and R13 evaluated this expression on either a two-dimensional or three-dimensional uniformly spaced grid. Each dimension was gridded into several hundred points, so that the computational cost rises exponentially with the number of mass proxies. The computational cost was trivial for V10, challenging for B13 and R13, and likely intractable in future work if another mass proxy is added. The techniques used in Mantz et al. (2010) and Vikhlinin et al. (2009b) also scale exponentially with the number of mass proxies.

Instead of computing Equation (8) on a uniformly spaced grid, we use the following numerical techniques. First, we note that the follow-up mass proxy ${Y}_{{\rm{X}}}$ immediately integrates out of the first term. This applies to any mass proxy that does not appear in the survey selection function. Therefore, we only need to perform a two-dimensional integral

$$\begin{eqnarray}&&-\int {dzd}\xi {{dY}}_{{\rm{X}}}\displaystyle \frac{{dN}(z,\xi ,{Y}_{{\rm{X}}})}{{dzd}\xi {{dY}}_{{\rm{X}}}}=-{\int }_{{z}_{\mathrm{cut}}}^{\infty }{dz}{\int }_{{\xi }_{\mathrm{cut}}}^{\infty }d\xi \displaystyle \frac{{dN}}{{dzd}\xi }.\end{eqnarray} \tag{ 9 }$$

We evaluate this expression with

$$\begin{eqnarray}&&{\displaystyle \int }_{{z}_{\mathrm{cut}}}^{\infty }{dz}{\displaystyle \int }_{{\xi }_{\mathrm{cut}}}^{\infty }d\xi \displaystyle \frac{{dN}}{{dzd}\xi }\hspace{4cm}\\ &&\quad ={\displaystyle \int }_{{z}_{\mathrm{cut}}}^{\infty }{dz}{\displaystyle \int }_{0}^{\infty }{dM}\displaystyle \frac{{dN}}{{dzdM}}P(\xi \gt {\xi }_{\mathrm{cut}}| M),\end{eqnarray} \tag{ 10 }$$

where $P(\xi \gt {\xi }_{\mathrm{cut}}| M)$ is simply the significance cut modeled through the scaling relation as described in Section 2.1.1.

The second term of Equation (8) is more challenging to evaluate. It can be written as an integral over the mass function,

$$\begin{eqnarray}\displaystyle \frac{{dN}({z}_{j},{\xi }_{j},{Y}_{{\rm{X}},j})}{{dzd}\xi {{dY}}_{{\rm{X}}}} & = & \displaystyle \int {dz}^{\prime} d\langle \xi \rangle {{dY}}_{{{\rm{X}}}^{\prime }}d\mathrm{ln}M\\ & & \times P({z}_{j}| z^{\prime} )P({\xi }_{j}| \langle \xi \rangle )P({Y}_{{\rm{X}},j}| {Y}_{{{\rm{X}}}^{\prime }})\\ & & \times P(z^{\prime} ,\langle \xi \rangle ,{Y}_{{{\rm{X}}}^{\prime }}| z,\mathrm{ln}M)\displaystyle \frac{{dN}}{{dzd}\mathrm{ln}M}.\end{eqnarray} \tag{ 11 }$$

The probability density functions $P({z}_{j}| z^{\prime} )$, $P({\xi }_{j}| \langle \xi \rangle )$, and $P({Y}_{{\rm{X}},j}| {Y}_{{{\rm{X}}}^{\prime }})$ describe the (independent) measurement error of each of the observables. $P(z^{\prime} ,\langle \xi \rangle ,{Y}_{{{\rm{X}}}^{\prime }}| z,\mathrm{ln}M)$ describes the joint scaling relations, implemented with multi-dimensional log-normal intrinsic scatter. In our case of two mass proxies (ξ, ${Y}_{{\rm{X}}}$), this is implemented with three parameters: two describing the marginal variances ${\sigma }_{\mathrm{ln}\zeta }$ and ${\sigma }_{\mathrm{Yx}}$, as well as one correlation coefficient ${\rho }_{\zeta ,\mathrm{Yx}}$, each of which is marginalized over in this work. Due to the treatment of maximization bias in Section 2.1.1, the scaling relation is defined in terms of ζ rather than $\langle \xi \rangle$. We perform the change of variables using

$$\begin{eqnarray}&&P(z^{\prime} ,\langle \xi \rangle ,Y{{\rm{X}}}^{\prime }| z,\mathrm{ln}M)=\displaystyle \int d\zeta ^{\prime} \\ &&\quad \displaystyle \frac{1}{\zeta ^{\prime} }\exp \left(-\displaystyle \frac{1}{2}{(\langle \xi \rangle -\zeta ^{\prime} )}^{2}\right)P(z^{\prime} ,\zeta ^{\prime} ,Y{{\rm{X}}}^{\prime }| z,\mathrm{ln}M).\end{eqnarray} \tag{ 12 }$$

The computational bottleneck for evaluating the cluster likelihood lies in evaluating Equation (11). We use Monte Carlo integration, randomly drawing from the probability distribution $P({z}_{j}| z^{\prime} )$ $P({\xi }_{j}| \langle \xi \rangle )$ $P({Y}_{{\rm{X}},j}| Y{{\rm{X}}}^{\prime })$ $P(z^{\prime} ,\langle \xi \rangle ,Y{{\rm{X}}}^{\prime }| z,\mathrm{ln}M)$, weighted by the mass function $\tfrac{{dN}}{d\mathrm{ln}M}$. Let $\mathrm{ln}{m}_{k}^{\prime }$ denote the kth integration variable corresponding to its flat-prior mass estimate (i.e., $\zeta ^{\prime}$ and $Y{{\rm{X}}}^{\prime }$ substituted into their respective observable-mass scaling relations). We can then write

$$\begin{eqnarray}&&P(\mathrm{ln}{m}_{k}^{\prime }| \mathrm{ln}M^{\prime} )=\displaystyle \frac{1}{\sqrt{{(2\pi )}^{{N}_{\mathrm{obs}}}| {\rm{\Psi }}| }}\\ &&\quad \exp \left(-\displaystyle \frac{1}{2}\displaystyle \sum _{l}(\mathrm{ln}{m}_{k}^{\prime }-\mathrm{ln}M^{\prime} ){({{\rm{\Psi }}}^{-1})}_{{kl}}(\mathrm{ln}{m}_{l}^{\prime }-\mathrm{ln}M^{\prime} )\right),\end{eqnarray} \tag{ 13 }$$

where ${\rm{\Psi }}$ is the mass proxy covariance matrix containing the intrinsic scatter and correlation coefficient parameters introduced in Sections 2.1.1 and 4.1 so that ${{\rm{\Psi }}}_{{kl}}=\langle (\mathrm{ln}{m}_{k}^{\prime }\mbox{--}\mathrm{ln}M^{\prime} )(\mathrm{ln}{m}_{l}^{\prime }\mbox{--}\mathrm{ln}M^{\prime} )\rangle$.

Now, we wish to obtain samples of the integration variable $\mathrm{ln}M^{\prime}$, given the location at which we are attempting to evaluate ${dN}({x}_{j})/{dx}$. In order to efficiently draw samples we first, for each k, draw samples from the (inverse) measurement error $P({x}_{k,j}| {x}_{k}^{\prime })$, which is assumed to be independent for each k. We then substitute these values into the scaling relations to obtain an ensemble of $\mathrm{ln}{m}_{k}^{\prime }$ values. The remaining task is then to draw random deviates $\mathrm{ln}M^{\prime}$ that follow the probability distribution explicitly shown in Equation (13), given each value of $\mathrm{ln}{m}_{k}^{\prime }$. To do so, we note that

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{kl}}(\mathrm{ln}{m}_{k}^{\prime }-\mathrm{ln}M^{\prime} )(\mathrm{ln}{m}_{l}^{\prime }-\mathrm{ln}M^{\prime} ){({{\rm{\Psi }}}^{-1})}_{{kl}}\\ &&\quad ={{\rm{\Sigma }}}_{{kl}}\mathrm{ln}{m}_{k}^{\prime }\mathrm{ln}{m}_{l}^{\prime }{({{\rm{\Psi }}}^{-1})}_{{kl}}\\ &&\qquad -2\mathrm{ln}M^{\prime} {{\rm{\Sigma }}}_{k}\mathrm{ln}{m}_{k}^{\prime }{{\rm{\Sigma }}}_{l}{({{\rm{\Psi }}}^{-1})}_{{kl}}+\mathrm{ln}M{^{\prime} }^{2}{{\rm{\Sigma }}}_{{kl}}{({{\rm{\Psi }}}^{-1})}_{{kl}}\\ &&\quad ={T}_{0}-2{T}_{1}\mathrm{ln}M^{\prime} +{T}_{2}\mathrm{ln}M{^{\prime} }^{2}\\ &&\quad ={\left(\tfrac{\mathrm{ln}M^{\prime} -{T}_{1}/{T}_{2}}{1/\sqrt{{T}_{2}}}\right)}^{2}+{T}_{0}-\tfrac{{T}_{1}^{2}}{{T}_{2}},\end{eqnarray} \tag{ 14 }$$

where ${T}_{0}={{\rm{\Sigma }}}_{{kl}}\mathrm{ln}{m}_{k}^{\prime }\mathrm{ln}{m}_{l}^{\prime }{({{\rm{\Psi }}}^{-1})}_{{kl}}$, ${T}_{1}={{\rm{\Sigma }}}_{k}\mathrm{ln}{m}_{k}^{\prime }{{\rm{\Sigma }}}_{l}{({{\rm{\Psi }}}^{-1})}_{{kl}}$, and ${T}_{2}={{\rm{\Sigma }}}_{{kl}}{({{\rm{\Psi }}}^{-1})}_{{kl}}$, which is quadratic in $\mathrm{ln}M^{\prime}$, so that Equation (13) is a log-normal distribution in $M^{\prime}$ with a known mean, width and normalization. We compute T₀, T₁, and T₂ explicitly and sample from the resulting distribution.

Having obtained samples of $\mathrm{ln}M^{\prime}$, which we denote as $\mathrm{ln}\tilde{M^{\prime} }$, we average the mass function

$$\begin{eqnarray}&&\displaystyle \frac{1}{{N}_{\mathrm{samples}}}\displaystyle \sum _{i}{\left.\displaystyle \frac{{dN}}{d\mathrm{ln}M^{\prime} }\right|}_{\mathrm{ln}{\tilde{M^{\prime} }}_{i}}\sqrt{\displaystyle \frac{{T}_{2}}{2\pi }}\displaystyle \frac{\exp -\tfrac{1}{2}\left({T}_{0}-\tfrac{{T}_{1}^{2}}{{T}_{2}}\right)}{\sqrt{{(2\pi )}^{{N}_{\mathrm{obs}}}| {\rm{\Psi }}| }}\end{eqnarray} \tag{ 15 }$$

and combine with Equation (14) to obtain an estimate for Equation (11). Our implementation of this estimator has been demonstrated to be unbiased through extensive simulations, with well-behaved residuals. The error on the mean is found to decrease as the inverse square root of the number of deviates drawn. In practice, for two mass proxies and 377 cluster candidates, we draw 10⁴ deviates per cluster candidate, resulting in rms noise on the likelihood surface of ${\rm{\Delta }}(-2\mathrm{ln}{ \mathcal L })\lesssim 0.1$ near the maximum likelihood. For the cluster catalog used in this work, the execution time is approximately one second on a single CPU thread. This likelihood module has been checked against an independent implementation, based on the estimator presented in Bocquet et al. (2015). The comparison shows agreement to well within the 1σ uncertainties on the cosmological parameter constraints found in this work.

In Section 6, we will discuss the cosmological constraints obtained using the analysis laid out so far. We will evaluate the compatibility of the cluster data with other data sets, as well as show the improvements in parameter constraints when the cluster data set is combined with other data sets. The primary external data we use are CMB power spectrum measurements, measurements of BAO from galaxy surveys, and distances to Type Ia supernovae (SNe).

The canonical CMB power spectrum data we use is the temperature–temperature power spectrum from the Planck 2013 release, combined with low-ℓ polarization information from WMAP (Planck Collaboration et al. 2014a, hereafter Planck+WP). The qualitative results from this work will be very similar to those that would have been obtained by considering the Planck 2015 data set, since the constraints on ${\sigma }_{8}$ and ${{\rm{\Omega }}}_{m}$ from the temperature and low-ℓ polarization power spectrum are very similar between the two Planck releases. In Section 6.2, we will briefly review the effect of substituting the best CMB power spectrum data from before the Planck 2013 release, using instead the combination of WMAP9 (Hinshaw et al. 2013) and SPT (Story et al. 2013, hereafter S13) CMB power spectrum data.

We sometimes use a prior on the angular scale of the sound horizon ${\theta }_{s}$. This is a powerful piece of cosmological information, and relies only very weakly on the details of the CMB analysis, since it affects the peak positions. We use the Planck+WP measurement, though we conservatively increase the uncertainty by a factor of five to $100{\theta }_{\mathrm{MC}}=1.0413\pm 0.0034$.⁵¹ Since ${\theta }_{s}$ is so well-measured, none of the results presented in this work that use a prior on ${\theta }_{s}$ are sensitive to the specific value of the assumed uncertainty.

When BAO data are used in this work, we use the SDSS-III BOSS results from data release 11. Specifically, we use the measurements of the parameter combination ${D}_{V}/{r}_{s}$ at z = 0.32 (LOWZ) and z = 0.57 (CMASS) from Anderson et al. (2014).

When considering cosmologies with a free parameter describing the dark energy equation of state, we compare to and contrast with the SN results from the joint likelihood analysis of the SDSS-II and SNLS SN samples from Betoule et al. (2014).

Cluster abundances do weakly constrain the baryon density, but for our main results we choose to adopt a prior based on big-bang nucleosynthesis calculations and deuterium abundance measurements of ${{\rm{\Omega }}}_{b}{h}^{2}=0.02202\pm 0.00045$ (Cooke et al. 2014). Rather than fixing this parameter, we marginalize over the uncertainty, though that analysis choice has no significant effect on the results presented in this work.

Finally, where stated, we adopt a prior of ${H}_{0}=73.8\pm 2.4\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$ from the direct H₀ measurements of Riess et al. (2011). We find that changing this prior by $1\sigma$ affects the resulting ${\sigma }_{8}$ constraint by $\lesssim 0.2\sigma$, while the ${{\rm{\Omega }}}_{m}$ constraint is affected by $\sim 0.4\sigma$. These two effects are anti-correlated so that the commonly used parameter combination ${\sigma }_{8}{({{\rm{\Omega }}}_{m}/0.27)}^{0.3}$ is highly insensitive to the assumed value of H₀.

A measurement of cluster abundance as a function of redshift provides constraints on cosmological parameters through several mechanisms. First, the total number of clusters found strongly depends on the matter density and the amplitude of the matter power spectrum (e.g., Bahcall & Cen 1992; White et al. 1993). As shown in, e.g., White et al. (1993), a simple spherical-collapse model for the halo mass function (Press & Schechter 1974) indicates that the total number of clusters should most strongly constrain the parameter combination ${\sigma }_{8}{{\rm{\Omega }}}_{m}^{\alpha }$, where α is related to the local slope of the matter power spectrum at the mean mass of the cluster catalog. This prediction has been empirically borne out in many cluster abundance studies, including this one.

The redshift dependence of the cluster abundance contains information on the growth function, as well as a dependence on the cosmic volume surveyed. This combination of sensitivity to growth and volume provides a unique constraint on parameters that affect the expansion history, notably w, the equation of state of dark energy (e.g., Haiman et al. 2001). In the case of w, the effect of changing this parameter on the cluster abundance actually changes sign at $z\sim 1$ (e.g., Figure 1 of Mohr 2005). The 377 cluster candidates presented in this work provide a large enough sample to constrain cosmological parameters by measuring the evolution of cluster abundance. However, this constraint is also potentially limited by knowledge of the observable-mass scaling relations, in this case their evolution with redshift. This is further discussed in the context of dark energy in Section 6.5.

In addition to their effect on the true abundance of clusters as a function of mass and redshift, cosmological parameters also affect the cluster observables through their influence on the SZ and X-ray scaling relations. Specifically, while the mass estimate inferred from the ζ–M relation depends weakly on redshift, the mass estimate inferred from the ${Y}_{{\rm{X}}}\mbox{--}M$ scaling relation depends strongly on the angular diameter distance to the cluster.

In this section, we present constraints on the parameters of the ΛCDM model. Because not all parameters of this model are well constrained by cluster counts alone, we adopt priors on some of the six cosmological parameters (${{\rm{\Omega }}}_{b}{h}^{2}$, ${{\rm{\Omega }}}_{c}{h}^{2}$, ${\theta }_{s}$, n_s, A_s, and τ, defined in Section 5.1), in addition to the priors on the scaling relation parameters (A_SZ, B_SZ, C_SZ, and ${\sigma }_{\mathrm{ln}\zeta }$, described in Section 2.1.1, as well as A_X, B_X, C_X, ${\sigma }_{\mathrm{lnYx}}$, and ${\rho }_{\zeta ,\mathrm{Yx}}$, described in Section 4.1). The cluster likelihood is insensitive to the optical depth to reionization τ as well as to the primordial scalar spectrum power-law index n_s, once an appropriate pivot point is chosen. Therefore, when discussing cluster constraints without the inclusion of CMB temperature power spectrum data, we fix these parameters to the best-fit values from Planck Collaboration et al. (2014a).

These two prior constraints, in addition to the H₀ and BBN priors discussed in 5.3, leave us two remaining degrees of freedom in the ΛCDM model, which we choose to express as ${{\rm{\Omega }}}_{m}$ and ${\sigma }_{8}$ (though in exploring the likelihood surface we actually vary the base parameters ${{\rm{\Omega }}}_{b}{h}^{2}$, ${{\rm{\Omega }}}_{c}{h}^{2}$, $\mathrm{ln}{10}^{10}{A}_{s}$ and H₀). We fix the species-summed neutrino mass ${\rm{\Sigma }}{m}_{\nu }$ to 0.06 eV, the minimum allowed value from terrestrial measurements of squared neutrino mass differences (see e.g., Gonzalez-Garcia et al. 2012 for a review). Figure 2 shows the constraints on ${\sigma }_{8}$ and ${{\rm{\Omega }}}_{m}$ from B13, R13, and this work. The large increase in the number of clusters primarily reduces the uncertainty on ${{\rm{\Omega }}}_{m}$ by constraining the shape of the halo mass function. The uncertainty on the parameter combination ${\sigma }_{8}{({{\rm{\Omega }}}_{m}/0.27)}^{0.3}$ is set by the knowledge of the overall mass scale, rather than the number of clusters, since Poisson errors on this number are subdominant. The improvement in the ${\sigma }_{8}{({{\rm{\Omega }}}_{m}/0.27)}^{0.3}$ constraint results from the increase in the number of clusters with Chandra follow-up data. These two effects result in tighter joint constraints on the two parameters than we obtained in previous cluster analyses, if we use the same prior on the observable-mass relations. This can be seen in Figure 2: the innermost, purple contours use the current cluster sample and the same observable-mass priors as B13 and R13. However, in contrast to B13 and R13 we choose to use a WL-based prior on the overall mass scale for our baseline results, as discussed in 4.2. This results in slightly degraded constraints, in particular on ${\sigma }_{8}$, as can be seen from the light-blue contours in Figure 2. With the current cluster sample, the stated priors on cosmological parameters, and the updated priors on the observable-mass relations, we obtain

$$\begin{eqnarray}&&{\sigma }_{8}=0.784\pm 0.039,\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{m}=0.289\pm 0.042,\end{eqnarray} \tag{ 17 }$$

and

$$\begin{eqnarray}&&{\sigma }_{8}{({{\rm{\Omega }}}_{m}/0.27)}^{0.3}=0.797\pm 0.031.\end{eqnarray} \tag{ 18 }$$

The ΛCDM parameter constraints, including scaling relation parameters, are shown in Table 3.

Zoom In Zoom Out Reset image size

Table 3.  Constraints on Cosmological and Scaling Relation Parameters Assuming a ΛCDM Cosmology

Parameter	Prior	SPT${}_{{\rm{CL}}}$+H0+BBN	SPT${}_{{\rm{CL}}}$+BAO+BBN	SPT${}_{{\rm{CL}}}$+Planck+WP+BAO
ASZ	5.38 ± 1.61	4.842 ± 0.913	4.936 ± 0.955	3.531 ± 0.273
BSZ	1.340 ± 0.268	1.668 ± 0.083	1.660 ± 0.064	1.661 ± 0.060
CSZ	0.49 ± 0.49	0.550 ± 0.315	0.864 ± 0.159	0.733 ± 0.123
${\sigma }_{\mathrm{ln}\zeta }$	0.13 ± 0.13	0.199 ± 0.069	0.201 ± 0.070	0.203 ± 0.066
AX	6.38 ± 0.61	6.235 ± 0.514	6.316 ± 0.505	7.030 ± 0.341
BX	0.57 ± 0.03	0.491 ± 0.023	0.493 ± 0.023	0.498 ± 0.021
CX	−0.40 ± 0.20	−0.251 ± 0.127	−0.280 ± 0.122	−0.174 ± 0.102
${\sigma }_{\mathrm{lnYx}}$	0.12 ± 0.08	0.162 ± 0.070	0.160 ± 0.069	0.154 ± 0.064
${\rho }_{\zeta ,\mathrm{Yx}}$	$[-0.98,0.98]$	−0.147 ± 0.458	−0.136 ± 0.465	−0.204 ± 0.443
${{\rm{\Omega }}}_{m}$	⋯	0.289 ± 0.042	0.306 ± 0.010	0.304 ± 0.007
${\sigma }_{8}$	⋯	0.784 ± 0.039	0.768 ± 0.030	0.820 ± 0.009
${\sigma }_{8}{({{\rm{\Omega }}}_{m}/0.27)}^{0.3}$	⋯	0.797 ± 0.031	0.797 ± 0.030	0.850 ± 0.013

Download table as:  ASCIITypeset image

In the left- and right-hand panels of Figure 3 we show the cluster abundance as a function of redshift and detection significance, respectively. Both show one-dimensional representations of the observable-space mass function. The data points, with approximate ($\sqrt{N}$) error bars shown, are independent of cosmological and scaling relation parameters. The points are independent between bins for ${dN}/d\xi$, and nearly independent for dN/dz, where clusters without spectroscopic redshift information can contribute to multiple bins.

Zoom In Zoom Out Reset image size

Shifts in the parameters ${\sigma }_{8}$ and A_SZ result in roughly global shifts in the amplitude of both curves, simultaneously. The scaling relation parameter B_SZ induces a roughly power-law tilt in ${dN}/d\xi$, and parameters such as C_SZ and ${{\rm{\Omega }}}_{m}$ induce tilts in dN/dz. This visualization shows two important ways in which the model, marginalized over a large number of scaling relation and cosmological parameters, is tested for agreement with the data.

Figure 4 shows the cluster constraints on the ΛCDM model in the ${\sigma }_{8}$-${{\rm{\Omega }}}_{m}$ plane, when combined with either the H₀+BBN prior or the BAO$+{\theta }_{s}$+BBN prior. We also show the constraints from CMB power spectrum measurements from Planck+WP and WMAP9+S13 data. The 68% confidence regions from the cluster constraints and the CMB power spectrum constraints overlap. We proceed to adopt the Planck+WP data set as the baseline CMB data set for the remainder of this work. We also combine the SPT${}_{{\rm{CL}}}$+Planck+WP+BAO data sets to obtain joint ΛCDM parameter constraints. Finally, we note that the ${\sigma }_{8}$ constraint obtained using the SPT${}_{{\rm{CL}}}$ data in Table 3 is within $1\sigma$ of the value recently reported by the Planck collaboration for the full-mission data (their 2015 TT+lowP+lensing data set) of ${\sigma }_{8}=0.815\pm 0.009$.

Zoom In Zoom Out Reset image size

6.2.1. Constraints on Scaling Relation Parameters

While the main focus of this work is on the cosmological constraints, the nuisance parameters are of interest themselves, both in terms of the cluster scaling relation parameter constraints and their degeneracy with the cosmological parameters. In Figure 5, we show the marginalized posterior in the multi-dimensional parameter space.

Zoom In Zoom Out Reset image size

The dominant systematic uncertainty limiting the cosmological constraints from the SPT-SZ cluster sample is the uncertainty of the overall cluster mass scale. This can be seen as a strong degeneracy between A_X and ${\sigma }_{8}{({{\rm{\Omega }}}_{m}/0.27)}^{0.3}$, which are 85% correlated given the SPT${}_{{\rm{CL}}}$+H₀+BBN data set. The second most important source of systematic uncertainty in the cosmological constraints is from the parameter characterizing the redshift evolution of the scaling relation, C_SZ. In a ΛCDM cosmology, it is highly degenerate with ${{\rm{\Omega }}}_{m}$, correlated at 87% when considering the SPT${}_{{\rm{CL}}}$+H₀+BBN data set.

The parameters that shift most significantly away from their priors are the B_SZ and B_X parameters which encode the power-law slopes of the scaling relations. We find that the preference for high B_SZ persists when the X-ray data are not used. In addition, this preference is not localized to any particular region of the data; when considering half the cluster sample at a time, either by redshift, ξ, or ${\gamma }_{\mathrm{field}}$, the high B_SZ persists, albeit at a lower significance.

This high B_SZ implies that the measured observable-space mass function, ${dN}/d\xi$, is shallower than expected given the scaling relation found from the simulations described in Section 2.2. We can approximately quantify this by assuming ${dN}/d\xi$ follows a power law. The data prefer a power-law index of ∼4.0, compared to the simulation prediction of ∼5.0, which is disfavored by the measured ${dN}/d\xi$ at $\sim 4.5\sigma$.

The preference for a higher B_SZ, i.e., a steeper ζ–M relation, is statistically weak, approximately $1\sigma$, due to the relatively weak 20% width assumed on the B_SZ prior. By contrast, the assumed width of the prior on the X-ray scaling relation slope B_X is 5%, leading to a tension between the likelihood and the prior of approximately $4\sigma$. There is some evidence for a similarly steep ${Y}_{{\rm{X}}}\mbox{--}M$ scaling relation from the comparison of the V09 and H15 WL mass estimates in Figure 1, which prefers ${B}_{X}\sim 0.42$. Such a slope would be in $\sim 4\sigma$ tension with the expected self-similar cluster slope of ${B}_{{\rm{X}}}=0.6$.

In this work, we will assume that our observational priors on the ${Y}_{{\rm{X}}}\mbox{--}M$ slope and theoretical priors on the mass function slope are well motivated, but note that the shape of the observable-space mass function will need to be studied in more detail for future work.

We now consider a cosmological model in which the species-summed neutrino mass ${\rm{\Sigma }}{m}_{\nu }$ is a free parameter. Constraints on this model from the CMB power spectrum show a strong degeneracy between ${\sigma }_{8}$ and ${\rm{\Sigma }}{m}_{\nu }$ (e.g., Komatsu et al. 2009; Planck Collaboration et al. 2014a; Abazajian et al. 2015). This allows even modest measurements of ${\sigma }_{8}$ to improve on neutrino mass constraints from the CMB power spectrum data alone. Figure 6 shows the improvement in constraints when adding cluster abundance information to CMB power spectrum data alone (Planck+WP), and also when including BAO data. In both cases, the addition of cluster information tightens the constraint and causes the posterior to peak at positive values of ${\rm{\Sigma }}{m}_{\nu }$, though the 95% upper limit on ${\rm{\Sigma }}{m}_{\nu }$ remains largely unchanged. Allowing the ${\rm{\Sigma }}{m}_{\nu }$ to vary in the range of 0–2 eV, and using the combination of the SPT${}_{{\rm{CL}}}$+Planck+WP+BAO data sets, we find

$$\begin{eqnarray}&&{\rm{\Sigma }}{m}_{\nu }=0.14\pm 0.08\mathrm{eV}.\end{eqnarray} \tag{ 19 }$$

Zoom In Zoom Out Reset image size

We note that the preference for positive ${\rm{\Sigma }}{m}_{\nu }$, when we combine our cluster abundance measurements with Planck data, is driven by the small residual tension between the preferred values of ${\sigma }_{8}$ in the two data sets. Such a preference has been pointed out by several authors (e.g., Battye & Moss 2014; Planck Collaboration et al. 2014b; Wyman et al. 2014), but is in contrast to the preference for positive ${\rm{\Sigma }}{m}_{\nu }$ shown, for example, by the combination of WMAP+SPT CMB power spectrum data and SPT cluster data in Hou et al. (2014), in which the evidence for positive ${\rm{\Sigma }}{m}_{\nu }$ is not driven by the cluster data. In this work, we find relatively good agreement between the preferred ${\sigma }_{8}$ using the CMB and SPT${}_{{\rm{CL}}}$ data sets, so the preference for positive ${\rm{\Sigma }}{m}_{\nu }$ is weak and consistent at 1σ with the minimum expected value of ∼0.06 eV from neutrino oscillation experiments (Lesgourgues & Pastor 2006). Relative to the previous cluster-based constraints cited above, the updated WL-based calibration described in Section 4.2, has shifted the normalization and increased the uncertainty of the observable-mass relation in a way that relieves tension with the Planck CMB data (see Figure 2). We also note that our constraint on ${\rm{\Sigma }}{m}_{\nu }$ is largely independent of the change in the optical depth to reionization, τ, between the Planck 2013 and 2015 data release; the SPT${}_{{\rm{CL}}}$ constraints are independent of τ, and the Planck CMB constraints on ${\sigma }_{8}$ and ${{\rm{\Omega }}}_{m}$ negligibly changed between the two data releases.

The effective number of relativistic species, ${N}_{\mathrm{eff}}$, affects the CMB power spectrum by altering the time of matter–radiation equality, changing the apparent sound horizon at recombination (e.g., Hou et al. 2013). This mechanism results in strong degeneracies between ${N}_{\mathrm{eff}}$ and other cosmological parameters, notably H₀ and ${\sigma }_{8}$, when considering CMB data alone (e.g., Bashinsky & Seljak 2004). Thus, the addition of constraints on ${\sigma }_{8}$—such as from the cluster data in this work—and H₀ can improve upon CMB-only constraints on ${N}_{\mathrm{eff}}$.

Here, we consider simultaneously varying the species-summed neutrino mass and the effective number of relativistic species. In this cosmological model, the Planck+WP data alone constrain the ${\sigma }_{8}$-${{\rm{\Omega }}}_{m}$-H₀ parameter volume relatively poorly. Adding the cluster information improves on all three of those parameters by roughly a factor of two. Through parameter degeneracies, this improves the simultaneous constraints on ${N}_{\mathrm{eff}}$ and ${\rm{\Sigma }}{m}_{\nu }$ as shown in Figure 7. The simultaneous constraints are

$$\begin{eqnarray}&&{N}_{\mathrm{eff}}=3.25\pm 0.23\end{eqnarray} \tag{ 20 }$$

and

$$\begin{eqnarray}&&{\rm{\Sigma }}{m}_{\nu }=0.39\pm 0.20\,\mathrm{eV},\end{eqnarray} \tag{ 21 }$$

which represent factors of 1.3 and 1.6 respective improvement over the Planck+WP data alone.

Zoom In Zoom Out Reset image size

Adding BAO data reduces the remaining allowed parameter space significantly to ${N}_{\mathrm{eff}}=3.28\pm 0.20$ and ${\rm{\Sigma }}{m}_{\nu }=0.18\pm 0.09$, and results in a degeneracy between ${N}_{\mathrm{eff}}$ and ${\rm{\Sigma }}{m}_{\nu }$, allowing for larger values of ${N}_{\mathrm{eff}}$ for increasing ${\rm{\Sigma }}{m}_{\nu }$. These constraints can be further tightened with the addition of local H₀ measurements, with the caveat that the best-fit value of H₀ from the SPT${}_{{\rm{CL}}}$+Planck+WP+BAO data set, ${H}_{0}=68.6\pm 1.2\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$, is in mild tension with direct local measurements from Riess et al. (2011). Proceeding to add those local measurements, so that we consider a SPT${}_{{\rm{CL}}}$+Planck+WP+H₀+BAO data set, we find a preference for larger ${N}_{\mathrm{eff}}$, resulting in the marginalized constraints of

$$\begin{eqnarray}&&{N}_{\mathrm{eff}}=3.43\pm 0.16\end{eqnarray} \tag{ 22 }$$

and

$$\begin{eqnarray}&&{\rm{\Sigma }}{m}_{\nu }=0.16\pm 0.08\,\mathrm{eV}.\end{eqnarray} \tag{ 23 }$$

The combined data set has a 2.3σ preference for ${N}_{\mathrm{eff}}\gt 3.046$, the standard model prediction. This is partially driven by the weak tension between local H₀ measurements and the Planck+BAO data set, as has been noted by other authors (e.g., Battye & Moss 2014; Hou et al. 2014; Wyman et al. 2014). However, the sensitivity to the H₀ prior is relatively weak. The preference for ${N}_{\mathrm{eff}}$ exceeding the standard model prediction would still be $2.0\sigma$ if the central value of the H₀ prior was reduced by one standard deviation.

This mild preference for massive neutrinos and extra ${N}_{\mathrm{eff}}$ has also been noted by other authors. For example, by combining CMB power spectrum data, BAO, SZ cluster counts, cosmic shear data, redshift space distortions from BOSS, and CMB lensing from Planck, Battye et al. (2015) found an even stronger preference for ${N}_{\mathrm{eff}}$ being larger than the ΛCDM model value. In MacCrann et al. (2015), the tension between CFHTLenS cosmic shear data and Planck+WP in the ${{\rm{\Omega }}}_{m}-{\sigma }_{8}$ plane is shown, and the effective number of extra neutrino species (${\rm{\Delta }}{N}_{\mathrm{eff}}$) is favored to be non-zero in the joint fit, alleviating such tension to some extent. More concretely, Ruiz & Huterer (2015) found redshift-space distortion data in $\sim 3\sigma$ tension with ΛCDM but that this tension was alleviated when the species-summed neutrino mass was made a free parameter. Another similar study has been done by Bernal et al. (2016) where adding neutrino mass reduces the tension on dark energy parameters inferred from the combination of geometry and the growth of the large-scale structure.

With the increased number of clusters in this work we are able to place constraints not only on the local cluster abundance but also on the evolution of cluster abundance with redshift. In particular, we examine the constraints on the wCDM cosmology, where the equation of state of dark energy w is a free parameter. We assume that w is a constant (i.e., its value does not evolve with redshift). This additional parameter affects the cluster abundance and observables through its influence on the geometry of the universe and, more importantly, the growth of structure. The geometrical effects include the change in the survey volume element and the angular diameter distance that modifies the implied X-ray mass information. However, in contrast to other probes of dark energy, the cluster abundance measurement is very sensitive to the effect of w on the growth of structure, primarily ${\sigma }_{8}(z)$ (e.g., Wang & Steinhardt 1998; Haiman et al. 2001).

In Figure 8, we show constraints on the dark energy equation of state parameter w and the energy density of dark energy today ${{\rm{\Omega }}}_{\mathrm{DE}}$ for different cosmological probes. With the SPT${}_{{\rm{CL}}}$+H₀+BBN data set, we obtain

$$\begin{eqnarray}&&w=-1.28\pm 0.31\end{eqnarray} \tag{ 24 }$$

and

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{DE}}=0.738\pm 0.046.\end{eqnarray} \tag{ 25 }$$

This is in good agreement, and of comparable precision, with the constraints when considering other cosmological probes, including Planck+WP, BAO, and SNe, as can be seen in Figure 8. Since all these probes except the cluster abundance measurement are primarily geometrical tests in this plane, the consistency between the cluster-implied parameter constraints, and those from other data sets, offers an important systematic test of dark energy. This measurement is limited primarily by our knowledge of the redshift evolution of the ζ–M scaling relation, C_SZ. Specifically, both ${{\rm{\Omega }}}_{m}$ and w are correlated with C_SZ at the $\sim 60 \%$ level, whereas the correlation with A_X is only $\sim 30 \%$.

Zoom In Zoom Out Reset image size

As shown in Figure 9, the Planck+WP measurements of the primary CMB show a strong degeneracy between w and ${\sigma }_{8}$. The addition of cluster data breaks the degeneracy and results in the marginalized constraints

$$\begin{eqnarray}&&w=-1.04\pm 0.17\end{eqnarray} \tag{ 26 }$$

and

$$\begin{eqnarray}&&{\sigma }_{8}=0.803\pm 0.045.\end{eqnarray} \tag{ 27 }$$

This level of w-uncertainty is ∼2.5 times larger compared to the constraints when adding either of the BAO or SNe data sets to the Planck+WP measurements (Aubourg et al. 2015; Betoule et al. 2014).

Zoom In Zoom Out Reset image size

For the data set combination H₀+BAO+SNe, which does not include primary CMB data, ${\sigma }_{8}$ is not determined. Adding the cluster data improves the w constraint by 37% through a direct measurement of the dark energy parameters w and ${{\rm{\Omega }}}_{\mathrm{DE}}$, rather than by breaking the w-${\sigma }_{8}$ degeneracy. For this combination, we find

$$\begin{eqnarray}&&w=-1.08\pm 0.07.\end{eqnarray} \tag{ 28 }$$

Finally, when considering the Planck+WP+H₀+BAO+SNe data sets, w is constrained to −1.062 ± 0.048. The allowed parameter space shows a significant w-${\sigma }_{8}$ degeneracy, which allows the addition of the cluster data to improve this constraint by 14%, to

$$\begin{eqnarray}&&w=-1.023\pm 0.042,\end{eqnarray} \tag{ 29 }$$

consistent with ΛCDM where $w=-1$.