New constraints on the mass bias of galaxy clusters from the power spectra of the thermal Sunyaev-Zeldovich effect and cosmic shear

TL;DR

This study combines thermal Sunyaev-Zeldovich and cosmic shear power spectra to constrain the galaxy cluster mass bias, providing insights into cluster physics and cosmological consistency without relying on primordial CMB data.

Contribution

It presents a joint analysis of tSZ and cosmic shear data to constrain the mass bias parameter independently of CMB priors, aligning observations with hydrodynamical simulations.

Findings

Measured mass bias B = 1.37^{+0.15}_{-0.23} for _8 < 0.9

Mass bias consistent with hydrodynamical simulations

Results support non-thermal pressure support in galaxy clusters

Abstract

Thermal Sunyaev-Zeldovich (tSZ) power spectrum is a powerful probe of the present-day amplitude of matter density fluctuations, and has been measured up to $\ell\approx 10^3$ from the {\it Planck} data. The largest systematic uncertainty in the interpretation of this data is the so-called "mass bias" parameter $B$, which relates the true halo mass to the mass proxy used by the {\it Planck} team as $M_{\rm 500c}^{Planck}=M_{\rm 500c}^{\rm true}/B$. Since the power spectrum of the cosmic weak lensing shear is also sensitive to the amplitude of matter density fluctuations via $S_8\equiv \sigma_8 \Omega_m^{\alpha}$ with $\alpha\sim 0.5$, we can break the degeneracy between the mass bias and the cosmological parameters by combining the tSZ and cosmic shear power spectra. In this paper, we perform a joint likelihood analysis of the tSZ power spectrum from {\it Planck} and the cosmic shear power spectrum from Subaru Hyper Suprime-Cam. Our analysis does not use the primordial cosmic microwave background (CMB) information. We obtain a new constraint on the mass bias as $B = 1.37 ^{+0.15}_{-0.23}$ or $(1-b) = B^{-1}=0.73^{+0.08}_{-0.13}$ (68\%~C.L.), for $\sigma_8 < 0.9$. This value of $B$ is lower than that needed to reconcile the tSZ data with the primordial CMB and CMB lensing data, i.e., $B = 1.64 \pm 0.19$, but is consistent with the mass bias expected from hydrodynamical simulations, $B = 1.28 \pm 0.20$. Our results thus indicate that the mass bias is consistent with the non-thermal pressure support from mass accretion of galaxy clusters via the cosmic structure formation, and that the cosmologies inferred from the tSZ and the cosmic shear are consistent with each other.

Equations

$$B\equiv M_{500c}^{\mathrm{true}}/M_{500c}^{Planck},$$

$$B=F^{-2.5}{\sigma }_8^{2.5}{\Omega }_m{h}^{-0.53}.$$

$$B=0.3{\left(\frac{S_8}{F}\right)}^{2.5}{\left(\frac{{\Omega }_m}{0.3}\right)}^{1-2.5\alpha }{h}^{-0.53}.$$

$$B=(1.36\pm 0.17){\left(\frac{0.3}{{\Omega }_m}\right)}^{0.25}{\left(\frac{0.7}{h}\right)}^{0.53}.$$

$$B=(1.37\pm 0.15){\left(\frac{{\sigma }_8}{0.8}\right)}^{0.5}{\left(\frac{0.7}{h}\right)}^{0.53}.$$

$${\rho }_m=({\Omega }_c+{\Omega }_b){\rho }_{\mathrm{crit}},$$

$$P(k,z)=\frac{2{\pi }^3}{k^3}{A}_s{\left(\frac{k}{k_{\ast }}\right)}^{n_s-1}{\left(\frac{{\Omega }_c{\mathcal{T}}_c(k,z)+{\Omega }_b{\mathcal{T}}_b(k,z)}{{\Omega }_c+{\Omega }_b}\right)}^2,$$

$$-2\ln {\mathcal{L}}_{\mathrm{tSZ}}(\mathbf{d}|\mathbf{\vartheta })={\Delta }^{\mathrm{T}}{\mathrm{Cov}}^{-1}\Delta ,$$

Full text

1 Introduction

The thermal Sunyaev-Zeldovich (tSZ) effect (Sunyaev & Zeldovich 1972) is the spectral distortion of the cosmic microwave background (CMB) through the inverse Compton scattering of CMB photons by hot thermal electrons in galaxy clusters. The tSZ angular power spectrum is sensitive to the amplitude of matter density fluctuations characterized by $\sigma_{8}$ and $\Omega_{m}$ (Komatsu & Kitayama 1999; Komatsu & Seljak 2002). These parameters are, however, strongly degenerate with the so-called “mass bias” parameter $B$, which is defined as the ratio of the mass proxy used by the Planck team and the true mass of galaxy clusters,

[TABLE]

where $M_{\rm 500c}$ is the mass enclosed by the radius $r_{\rm 500c}$ within which the average mass density is 500 times of the critical density of the Universe. The mass bias $B$ is related to the more commonly used parameter $b$ as $B=(1-b)^{-1}$ (Planck Collaboration et al. 2016b). Specifically, the tSZ power spectrum depends primarily on $F\equiv\sigma_{8}(\Omega_{m}/B)^{0.40}h^{-0.21}$ (Bolliet et al., 2018), where $\Omega_{m}$ is the matter density parameter, $\sigma_{8}$ the r.m.s. matter density fluctuation smoothed over a $8~{}h^{-1}~{}{\rm Mpc}$ sphere, and $h$ the dimensionless Hubble constant defined by $H_{0}=100~{}h~{}{\rm km~{}s^{-1}~{}Mpc^{-1}}$.

The mass bias was introduced to account for the cluster mass uncertainty in the Planck analysis. The galaxy cluster masses used by the Planck team were calibrated against a local cluster sample observed by XMM-Newton, assuming the hydrostatic equilibrium (HSE) with thermal pressure (Arnaud et al. 2010). The Planck team reported that 20–40% mass bias (i.e., $1.25<B<1.67$) is required to reconcile the power spectra of tSZ clusters with the joint result of Planck primordial CMB, CMB lensing and Baryonic Acoustic Oscillations (BAO) (Planck Collaboration et al. 2016b). Several authors performed follow-up analyses of the Planck tSZ power spectrum and found similar results (Horowitz & Seljak 2017; Hurier & Lacasa 2017; Salvati et al. 2018; Bolliet et al. 2018; Makiya et al. 2018; Salvati et al. 2019).

On the other hand, state-of-the-art cosmological hydrodynamic simulations show that the HSE mass underestimates the true mass due to non-thermal pressure support (e.g., Kay et al. 2004; Rasia et al. 2006, 2012; Nagai et al. 2007; Piffaretti & Valdarnini 2008; Lau et al. 2009; Meneghetti et al. 2010; Lau et al. 2013; Nelson et al. 2014b, a). The dominant contribution to the non-thermal pressure support seen in the simulations is the mass accretion of galaxy clusters via structure formation (Shi & Komatsu, 2014; Shi et al., 2015), which yields $B=1.28\pm 0.20$ (68% C.L.) over a wide range of dynamical states of galaxy clusters (see Table 1 of Shi et al., 2016, for the mass-limited sample and the fitting range of $(0.1,1.5)r^{\rm true}_{500}$). The predicted HSE mass bias is therefore lower than, or at least on the lower side of, the Planck inferred value. This indicates that other sources of non-thermal pressure (such as cosmic rays and magnetic fields) are larger than expected; there might be other systematic effects such as the calibration error in X-ray observations; and/or that new physics beyond the standard $\Lambda$ Cold Dark Matter (CDM) model, e.g., dark energy different from the cosmological constant (Bolliet et al., 2018), and/or modified gravity, is required to resolve a tension between cosmological parameters inferred from the tSZ clusters and those from primordial CMB based on $\Lambda$CDM (Planck Collaboration et al., 2016b). Massive neutrinos, which also modify the evolution of matter density fluctuations, have been shown not to help resolve this tension (Salvati et al., 2018; Sakr et al., 2018; Bolliet et al., 2019).

To gain more insights into this potential tension, we need an estimate of the mass bias that is independent of the primordial CMB. To this end, in this paper we perform a joint analysis of power spectra of the late-time Universe probes: tSZ and the cosmic weak lensing shear, to obtain a new constraint on the mass bias that is independent of the primordial CMB. Cosmic shear is a unique probe of the growth of total matter distribution including dark matter. The angular power spectrum of the cosmic shear is sensitive to $S_{8}\equiv\sigma_{8}(\Omega_{m}/0.3)^{\alpha}$ with $\alpha\sim 0.5$. This information is useful to break the degeneracy between the cosmological parameters and $B$.

Here we use the cosmic shear measurements from the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP, hereafter HSC) (Aihara et al. 2018b, a). HSC is a wide-field imaging survey using a 1.5 deg diameter field-of-view camera mounted on the prime focus of the 8.2m Subaru telescope (Miyazaki et al., 2018). The unique property of HSC is a combination of the depth ($i\sim 26$) and the excellent image quality (typical $i$-band seeing is $\sim 0.58"$), which enables us to measure cosmic shear signals with unprecedentedly high precision. Hikage et al. (2019) measured tomographic lensing power spectra using the first-year shear catalog (Mandelbaum et al., 2018; Oguri et al., 2018) where the total sky coverage is 137 deg2 and the effective number density is 17 arcmin*-2* in the range of photometric redshifts from 0.3 to 1.5. The first-year shear catalog has been made publicly available as a part of the second public data release of HSC (Aihara et al., 2019).

Constraints on the mass bias from the tSZ and cosmic shear power spectra can be estimated as follows. From the definition of $F$, the scaling of mass bias with other parameters can be written as

[TABLE]

Combining the scaling relation of $\sigma_{8}$ from the cosmic shear, $\sigma_{8}=S_{8}(0.3/\Omega_{m})^{\alpha}$, we find

[TABLE]

If we put $F=0.460\pm 0.012$ from Planck (Bolliet et al., 2018) and $S_{8}=0.780^{+0.030}_{-0.033}$ with $\alpha=0.5$ from HSC (Hikage et al., 2019), we obtain

[TABLE]

Similarly we can replace $\Omega_{m}$ in Eq. (2) with $\sigma_{8}$. In this case we obtain

[TABLE]

These quick estimates are closer to $B=1.28\pm 0.20$ expected from non-thermal pressure due to mass accretion (Shi et al., 2016), which motivates our calculating $B$ from a joint analysis of the tSZ and cosmic shear power spectra in this paper.

This paper is organized as follows. In Section 2 we describe the datasets of tSZ and cosmic shear that we use. In Section 3 we outline the model of the tSZ power spectrum and the details of our likelihood analysis. We discuss our results in Section 4 and conclude in Section 5. Throughout the paper, we adopt a flat $\Lambda$CDM cosmology with the minimal total neutrino mass of $\sum m_{\nu}=0.06\;{\rm eV}$.

2 Data

2.1 tSZ

We use the tSZ power spectrum data before marginalizing over the foreground components. Specifically, we take the data from Table 3 of Bolliet et al. (2018), which are based on Planck Collaboration et al. (2016b).

The tSZ power spectrum was calculated by cross-correlating the first-half data of the Needlet Internal Linear Combination (NILC) map and last-half of the Modified Internal Linear Combination Algorithm (MILCA) map, where NILC and MILCA are two different methods for reconstructing the tSZ map.111In Makiya et al. (2018) we used the cross-power spectrum of NILC’s fist- and last-half maps as the Compton-Y auto spectrum, in order to be consistent with the galaxy–Compton-Y cross spectrum analysis. In this paper we decided to follow Planck Collaboration et al. (2016b) and Bolliet et al. (2018) for simplicity. We also take into account contributions from residual contaminating sources such as the cosmic infrared background (CIB), IR and radio point sources, and the correlated noise. The amplitudes of CIB, IR and radio point source power spectra are treated as free parameters, while their shapes, which are also taken from Table 3 of Bolliet et al. (2018), are fixed. The amplitude of correlated noise is fixed to reproduce the tSZ spectrum at $\ell=2742$, following Bolliet et al. (2018).

2.2 Cosmic shear

We use a sample of the posterior distributions of the cosmological parameters of a flat $\Lambda$CDM model with a fixed minimal total neutrino mass $\sum m_{\nu}=0.06\;{\rm eV}$ obtained from the HSC lensing power spectra alone (Hikage et al., 2019)222The likelihood chain can be found at http://gfarm.ipmu.jp/~surhud/#hikage. The parameters include the five basic cosmological parameters ($\Omega_{b}h^{2}$, $\Omega_{c}h^{2}$, $h$, $A_{s}$ and $n_{s}$) and nine nuisance parameters regarding modelling errors of point spread functions (PSF), shear biases, photo-$z$ errors, and intrinsic alignments. The reionization optical depth $\tau$ is not used in the cosmic-shear-alone analysis. The range of flat priors of cosmological parameters adopted in the HSC cosmic shear analysis are listed in Table 3.1.

3 Analysis

We model and analyze the tSZ power spectrum using the halo-model-based approach that we have established in Bolliet et al. (2018) and Makiya et al. (2018), which are based on Komatsu & Seljak (2002).333While we used the mass function of Bocquet et al. (2016) in Makiya et al. (2018), we follow Bolliet et al. (2018) to use that of Tinker et al. (2008) in this paper. Since their mass functions are for the overdensity with respect to the mean mass density (rather than the critical density), we use the spline interpolation of the parameters at various overdensities to obtain the mass function for $M_{500c}$. As Tinker et al. (2008) only provide the mass function parameters up to the mean mass overdensity of $\Delta_{m}$ = 3200, we linearly extrapolate the parameters beyond this bound. We have checked the extrapolation against the fitting function provided by Tinker et al. (2008) and found that the extrapolation method does not have significant effects on our results. See Section 3.1.1 of Makiya et al. (2018) for details of the implementation.

There are two important changes from the methodology established in these papers. First, we properly take into account the effect of massive neutrinos in the modelling of the power spectrum and mass function of collapsed objects (Section 3.1). Second, we improve the method to calculate the likelihood when the covariance matrix includes a non-Gaussian term (Section 3.2).

3.1 Massive neutrinos

We include the effect of massive neutrinos by following the so-called “CDM prescription” (Ichiki & Takada, 2012; Costanzi et al., 2013; Villaescusa-Navarro et al., 2014; Castorina et al., 2014). The basic idea is to remove the contribution of neutrinos to the mass of collapsed objects (halos) when computing statistics of halos, as neutrino stream out of them. In the calculation of the dark matter halo mass function, we modify the matter density $\rho_{m}$ and the matter power spectrum $P(k,z)$ as

[TABLE]

and

[TABLE]

where $A_{s},n_{s}$ and $k_{*}=0.05~{}{\rm Mpc}^{-1}$ are the amplitude, the spectral index and the pivot scale of the primordial power spectrum, $\Omega_{c}$ and $\Omega_{b}$ are the density parameters of CDM and baryons, $\rho_{\rm crit}$ is the critical density of the Universe, and $\mathcal{T}_{c}$ and $\mathcal{T}_{b}$ are the transfer functions of CDM and baryons respectively. We refer the reader to Bolliet et al. (2019) for more details of the effects of massive neutrios on the tSZ power spectrum.

3.2 Likelihood analysis

To quickly perform a joint likelihood analysis we use the importance sampling technique (Lewis & Bridle 2002) as follows. We leave a full likelihood analysis as a future work.

First we read a set of cosmological parameters from the likelihood chains of the HSC cosmic shear analysis line by line and calculate the tSZ power spectrum and its likelihood $\mathcal{L}_{\rm tSZ}$ for the same set of cosmological parameters. Then we reweigh the chain by multiplying a weight by $\exp(-\ln\mathcal{L}_{\rm tSZ})$. Applying this procedure to the entire sample of HSC chains, we obtain the tSZ-shear joint posterior distributions.

Our model has four nuisance parameters: the mass bias $B$ and the amplitudes of the power spectra of CIB, IR and radio point sources, $A_{\rm CIB},A_{\rm IR}$ and $A_{\rm Rad}$, respectively. The nuisance parameters are randomly picked from a flat prior with the range listed in Table 3.1. Following Bolliet et al. (2018), the total power of the contaminating sources are restricted not to exceed the residual of the total tSZ power and the sum of the contributions from resolved clusters. We iterate the above procedure by changing the random seed for nuisance parameters until the chains are converged, and then combine multiple chains to obtain the posterior distributions. We judge that the chains converge when the Gelman-Rubin estimator $R-1$, where $R$ is the ratio of the variance between chains and within chains (Gelman & Rubin 1992), is less than 0.05.

The tSZ likelihood is calculated as

[TABLE]

where $\Delta$ is the difference vector between the observed and the model tSZ spectra, and ${\rm Cov}$ is the covariance matrix including the non-Gaussian term calculated from the model tSZ trispectrum. The Gaussian term of the covariance matrix is taken from Table 3 of Bolliet et al. (2018).

In Makiya et al. (2018) we calculated the non-Gaussian term at each step of parameter inference. However Carron (2013) pointed out that such a parameter-dependent covariance matrix adds extra artificial information and can bias the parameter constraints. Thus we adopt the new procedure in this paper. First, we perform a likelihood analysis without the non-Gaussian term and find the best-fitting parameters. Then we calculate the non-Gaussian term from the best-fitting parameters and repeat the likelihood analysis by including the fixed non-Gaussian term in the covariance matrix. From the new best-fitting parameters, we recalculate the covariance matrix and redetermine the best-fitting parameters. We iterate this procedure until the best-fitting parameters converge.

4 Results

Figure 1 shows the constraints on the parameters from the Planck tSZ alone, the HSC shear alone, and the joint analysis of them. The results from the joint likelihood analysis are summarized in Table 3.1. The nuisance parameters not shown in Figure 1 ( $A_{\rm CIB}$, $A_{\rm IR}$, $A_{\rm Rad}$) are well determined within the prior ranges, as shown in Table 3.1. The constraints on $S_{8}$ and $F$ are consistent with those obtained from the HSC cosmic shear and the Planck tSZ power spectrum alone, $S_{8}=0.780^{+0.030}_{-0.033}$ (Hikage et al., 2019) and $F=0.460\pm 0.012$ (Bolliet et al., 2018), respectively. Since the cosmic shear and tSZ power spectra are not sensitive to individual parameters (i.e., all cosmological parameters and mass bias parameter) but only sensitive to the combination of them, $S_{8}$ (for the cosmic shear) and $F$ (for the tSZ), the best-fitting model from the joint analysis also gives a good fit to both data sets.

As shown in Figure 1, the mass bias cannot be determined from the tSZ data alone because of the degeneracy with $\sigma_{8}$, $\Omega_{m}$ and $h$. The HSC cosmic shear lifts this degeneracy and helps to determine the mass bias. We find that the constraints from the joint analysis do not fully overlap with those from the tSZ alone or the shear alone analysis. This is due to a large difference between the best-fitting values of the tSZ alone or shear alone analysis and those from the joint analysis. The contour shows the likelihood distance from the best-fitting point. When the best-fitting values of two different distributions are largely different, direct comparison of two contours is not meaningful because those contours measure the likelihood distances from different points.

We find that the Planck tSZ and the HSC cosmic shear power spectra constrain the mass bias as $B=1.54^{+0.20}_{-0.35}$ (mean and 68% C.L.). The black lines in Figure 2 are the predicted scaling relation shown in Eq. (4) and (5) with $\alpha=0.5$. We used $h=0.76$ for $B$–$\Omega_{m}$ and $B$–$\sigma_{8}$ relations and $\Omega_{m}=0.23$ for $B$–$h$ relation, which are the mean value from the joint analysis, respectively. The measurements are consistent with the predicted relations.

Figure 2 shows that the Hubble parameter $h$ is not well constrained within the range of prior. While the limited range of the prior of $h$ has an effect on the constraint on $B$, the effect is not significant since dependence of $B$ on $h$ is weak. We have also examined a Gaussian prior of $h=0.74\pm 0.014$ taken from the distance ladder method (Riess et al., 2019) and found that the constraint on $B$ does not improve significantly; $B=1.60^{+0.21}_{-0.33}$. This is again due to the weak $h$ dependence of $B$.

Constraint on the mass bias from our analysis is mainly limited by the uncertainty of $\sigma_{8}$. Currently, most (if not all) data sets indicate and are consistent with $\sigma_{8}<0.9$ (Planck Collaboration et al., 2016a; Alam et al., 2017; Abbott et al., 2018; Burenin, 2018). To incorporate this knowledge into our analysis, in Figure 3 we show marginalized joint constraints on $B$ and $\sigma_{8}$ with a prior of $\sigma_{8}<0.9$. This $\sigma_{8}$ cut roughly corresponds to $\Omega_{m}>0.2$ (see the $\sigma_{8}$–$\Omega_{m}$ panel of Figure 1), which is also consistent with most (if not all) data sets. In this case the mass bias is constrained to be $B=1.37^{+0.15}_{-0.23}$ or $(1-b)=B^{-1}=0.73^{+0.08}_{-0.13}$ (mean and 68% C.L.). This is lower than that from, but is still consistent at 68% C.L. with, the joint analysis of the tSZ and CMB including CMB lensing, $B=1.64\pm 0.19$.444This constraint is different from that shown in Makiya et al. (2018) due to the difference of dark matter mass function, the tSZ auto-power spectrum data and the treatment of the non-Gaussian term of the covariance matrix, as noted in Section 2 and 3. On the other hand, our result is consistent with the value expected from cosmological hydrodynamical simulations, $B=1.28\pm 0.20$ (Shi et al. 2016) and also with that estimated from weak lensing mass, $B=1.25\pm 0.13$ (Salvati et al. 2018 and references therein; see also Miyatake et al. 2019; Stern et al. 2019; Dietrich et al. 2019 for recent attempts).

5 Conclusion

In this paper we have performed a joint likelihood analysis of the tSZ angular power spectrum obtained by Planck and the cosmic shear angular power spectra obtained by Subaru HSC. We have found that the mass bias of Planck clusters is constrained to be $B=1.37^{+0.15}_{-0.23}$ or $(1-b)=B^{-1}=0.73^{+0.08}_{-0.13}$ (mean and $68\%$ C.L.) for $\sigma_{8}<0.9$.

Our result is consistent with the HSE mass bias estimated from the cosmological hydrodynamical simulations, which predict $B=1.28\pm 0.20$ over a wide range of dynamical states of galaxy clusters (Shi et al., 2016). The origin of this bias is non-thermal motion arising from mass accretion of galaxy clusters via structure formation (Shi & Komatsu, 2014; Shi et al., 2015). Therefore, as long as we adopt the cosmological parameters inferred from the HSC cosmic shear with a weak prior of $\sigma_{8}<0.9$, the origin of the mass bias can be mostly understood. In other words, the cosmological parameters inferred from two different probes of a late-time Universe, tSZ and cosmic shear, are consistent with each other when the mass bias agrees with the expectations from cosmological hydrodynamical simulations. It has been known that the cosmological parameters inferred from the primordial CMB and those from the late-time probes are in a mild tension (e.g., Riess et al. 2018a, b, Joudaki et al. 2017; Köhlinger et al. 2017; Troxel et al. 2018; Burenin 2018; Hikage et al. 2019). Our results suggest that the tSZ power spectrum may also be in tension with the primordial CMB; a higher value of $B$ reported by Planck may come from the tension of $\sigma_{8}$ (or $S_{8}$) between Planck and the late-time Universe probes of the cosmic shear and SZ clusters.

More accurate measurements of cosmic shear is required to obtain a robust conclusion on this issue. In this paper we have used the HSC year 1 data, which are based on only 11$\%$ of the planned HSC survey data. Full HSC survey will put tighter constraints on $S_{8}$ and also improve a constraint on the mass bias. Combining other probes such as galaxy-galaxy lensing and galaxy clustering will improve the constraints.

{ack}

R. M. thanks Boris Bolliet for useful discussion. Kavli IPMU is supported by World Premier International Research Center Initiative (WPI), MEXT, Japan. This work was supported in part by JSPS KAKENHI Grant Number JP15H05896, JP15K21733, JP16K17684 and JP18H04348.

The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University.

This paper makes use of software developed for the Large Synoptic Survey Telescope. We thank the LSST Project for making their code available as free software at http://dm.lsst.org

The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under Grant No. AST-1238877, the University of Maryland, and Eotvos Lorand University (ELTE) and the Los Alamos National Laboratory.

Based [in part] on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by Subaru Telescope and Astronomy Data Center at National Astronomical Observatory of Japan.