The mass-richness relation of optically-selected clusters from weak gravitational lensing and abundance with Subaru HSC first-year data

TL;DR

This study measures the relation between cluster richness and halo mass using weak lensing and abundance data from Subaru HSC, providing insights into the mass-richness relation across redshifts and its implications for cosmology.

Contribution

It presents the first detailed measurement of the richness-mass relation for optically-selected clusters over a wide redshift range using Subaru HSC data, including modeling of redshift-dependent scatter.

Findings

Richness-mass relation is consistent with models when including redshift-dependent terms.

Scatter in mass-richness relation varies with richness and differs between Planck and WMAP cosmologies.

Constraints on the mean relation are consistent across cosmological models.

Abstract

Constraining the relation between the richness $N$ and the halo mass $M$ over a wide redshift range for optically-selected clusters is a key ingredient for cluster-related science in optical surveys, including the Subaru Hyper Suprime-Cam (HSC) survey. We measure stacked weak lensing profiles around 1747 HSC CAMIRA clusters over a redshift range of $0.1\leq z_{\rm cl}\leq 1.0$ with $N\geq 15$ using the HSC first-year shear catalog covering $\sim$$140$ ${\rm deg^2}$. The exquisite depth and image quality of the HSC survey allow us to measure lensing signals around the high-redshift clusters at $0.7\leq z_{\rm cl}\leq 1.0$ with a signal-to-noise ratio of 19 in the comoving radius range $0.5\lesssim R\lesssim 15 h^{-1}{\rm Mpc}$. We constrain richness-mass relations $P(\ln N|M,z)$ of the HSC CAMIRA clusters assuming a log-normal distribution without informative priors on model parameters, by jointly fitting to the lensing profiles and abundance measurements under both Planck and WMAP cosmological models. We show that our model gives acceptable $p$-values when we add redshift dependent terms which are proportional to $\ln (1+z)$ and $[\ln (1+z)]^{2}$ into the mean and scatter relations of $P(\ln N|M,z)$. Such terms presumably originate from the variation of photometric redshift errors as a function of the redshift. We show that the constraints on the mean relation $\langle M|N \rangle$ are consistent between the Planck and WMAP models, whereas the scatter values $\sigma_{\ln M|N}$ for the Planck model are systematically larger than those for the WMAP model. We also show that the scatter values for the Planck model increase toward lower richness values, whereas those for the WMAP model are consistent with constant values as a function of richness. This result highlights the importance of the scatter in the mass-richness relation for cluster cosmology.

Tables

Table 1. Table 1: Binning scheme for the CAMIRA clusters and characteristics of each bin for the abundance and lensing measurements. ∗

Abundance bin	$\alpha$	$\beta$	$N_{\min }$	$N_{\max }$		$z_{\mathrm{cl},\min }$	$z_{\mathrm{cl},\max }$		$N_{\alpha ,\beta }^{\mathrm{uncorr}}$	${\hat{N}}_{\alpha ,\beta }$
1	1	1	15	20	17.1	0.1	0.4	0.27	255	258.3
2	2	1	20	30	24.1	0.1	0.4	0.27	208	210.7
3	3	1	30	60	38.2	0.1	0.4	0.26	92	93.3
4	4	1	60	200	77.4	0.1	0.4	0.26	9	9.1
5	1	2	15	20	17.2	0.4	0.7	0.56	301	301.7
6	2	2	20	30	24.0	0.4	0.7	0.53	210	210.6
7	3	2	30	60	38.9	0.4	0.7	0.54	79	79.2
8	4	2	60	200	73.1	0.4	0.7	0.50	7	7.0
9	1	3	15	20	17.0	0.7	1.0	0.84	339	339.2
10	2	3	20	30	23.8	0.7	1.0	0.83	181	181.1
11	3	3	30	60	36.3	0.7	1.0	0.84	65	65.0
12	4	3	60	200	64.7	0.7	1.0	1.0	1	1.0
Lensing bin	$\alpha$	$\beta$	$N_{\min }$	$N_{\max }$		$z_{\mathrm{cl},\min }$	$z_{\mathrm{cl},\max }$		$N_{\alpha ,\beta }^{\mathrm{uncorr}}$	${\hat{N}}_{\alpha ,\beta }$
1	1	1	15	20	17.1	0.1	0.4	0.27	255	258.3
2	2	1	20	30	24.1	0.1	0.4	0.27	208	210.7
3	3	1	30	200	41.7	0.1	0.4	0.26	101	102.4
4	1	2	15	20	17.2	0.4	0.7	0.56	301	301.7
5	2	2	20	30	24.0	0.4	0.7	0.53	210	210.6
6	3	2	30	200	41.6	0.4	0.7	0.53	86	86.2
7	1	3	15	20	17.0	0.7	1.0	0.84	339	339.2
8	2	3	20	30	23.8	0.7	1.0	0.83	181	181.1
9	3	3	30	200	36.7	0.7	1.0	0.84	66	66.0

Table 2. Table 2: The model parameters, priors, and parameter estimations from our joint analysis of lensing and abundance measurements. ∗

Parameter	Description	Prior	Median and Error	Median and Error
			Planck	WMAP
$A$	$⟨\ln N⟩$ at pivot mass scale and pivot redshift	$\mathrm{flat}[2,5]$	${3.15}_{-0.08}^{+0.07}$	${3.36}_{-0.06}^{+0.05}$
$B$	Coefficient of halo mass dependence in $⟨\ln N⟩$	$\mathrm{flat}[0,2]$	${0.86}_{-0.05}^{+0.05}$	${0.83}_{-0.03}^{+0.03}$
$B_z$	Coefficient of linear redshift dependence in $⟨\ln N⟩$	$\mathrm{flat}[-50,50]$	$-{0.21}_{-0.42}^{+0.35}$	$-{0.20}_{-0.34}^{+0.26}$
$C_z$	Coefficient of square redshift dependence in $⟨\ln N⟩$	$\mathrm{flat}[-50,50]$	${3.61}_{-2.23}^{+1.96}$	${3.51}_{-1.59}^{+1.32}$
${\sigma }_0$	${\sigma }_{\ln N|M,z}$ at pivot mass scale and pivot redshift	$\mathrm{flat}[0,1.5]$	${0.32}_{-0.06}^{+0.06}$	${0.19}_{-0.07}^{+0.07}$
$q$	Coefficient of halo mass dependence in ${\sigma }_{\ln N|M,z}$	$\mathrm{flat}[-1.5,1.5]$	$-{0.08}_{-0.04}^{+0.05}$	$-{0.02}_{-0.03}^{+0.03}$
$q_z$	Coefficient of linear redshift dependence in ${\sigma }_{\ln N|M,z}$	$\mathrm{flat}[-50,50]$	${0.03}_{-0.30}^{+0.31}$	${0.23}_{-0.35}^{+0.37}$
$p_z$	Coefficient of square redshift dependence in ${\sigma }_{\ln N|M,z}$	$\mathrm{flat}[-50,50]$	${0.70}_{-1.60}^{+1.71}$	${1.26}_{-1.49}^{+1.77}$
$f_0$	Centering fraction at pivot richness and redshift		${0.68}_{-0.06}^{+0.05}$	${0.68}_{-0.06}^{+0.05}$
$f_N$	Coefficient of richness dependence in $f_{\mathrm{cen}}^{\alpha ,\beta }$		${0.33}_{-0.09}^{+0.10}$	${0.33}_{-0.09}^{+0.10}$
$f_z$	Coefficient of redshift dependence in $f_{\mathrm{cen}}^{\alpha ,\beta }$		$-{0.14}_{-0.34}^{+0.35}$	$-{0.19}_{-0.34}^{+0.34}$
$R_{1,\mathrm{off}}$	Off-centering radius for the first redshift bin	$\mathrm{flat}[{10}^{-3},0.64]$	${0.39}_{-0.09}^{+0.10}$	${0.38}_{-0.09}^{+0.10}$
$R_{2,\mathrm{off}}$	Off-centering radius for the second redshift bin	$\mathrm{flat}[{10}^{-3},0.78]$	${0.55}_{-0.11}^{+0.12}$	${0.52}_{-0.10}^{+0.12}$
$R_{3,\mathrm{off}}$	Off-centering radius for the third redshift bin	$\mathrm{flat}[{10}^{-3},0.92]$	${0.59}_{-0.17}^{+0.17}$	${0.58}_{-0.15}^{+0.17}$
$m_{\mathrm{lens}}$	Marginalization parameter of the lensing amplitudes	$\mathrm{Gauss}(0,0.01)$	${0.00}_{-0.01}^{+0.01}$	${0.00}_{-0.01}^{+0.01}$
${\chi }_{\min }^2/\mathrm{dof}$			$107.0/97$	$106.6/97$

Table 3. Table 3: The median and 68 % percent \% percentile uncertainties of the model parameters with the Planck cosmological parameters when we use the measurements only from one redshift bin. ∗

Parameter	Low-$z$	Middle-$z$	High-$z$
$A$	${3.34}_{-0.20}^{+0.25}$	${3.19}_{-0.15}^{+0.20}$	${3.31}_{-0.26}^{+0.15}$
$B$	${0.85}_{-0.07}^{+0.08}$	${0.94}_{-0.07}^{+0.09}$	${0.88}_{-0.05}^{+0.08}$
${\sigma }_0$	${0.36}_{-0.21}^{+0.07}$	${0.33}_{-0.21}^{+0.06}$	${0.27}_{-0.20}^{+0.14}$
$q$	$-{0.06}_{-0.11}^{+0.09}$	$-{0.08}_{-0.09}^{+0.09}$	$-{0.03}_{-0.11}^{+0.04}$
$f_{\mathrm{cen}}^{1,\beta }$	${0.52}_{-0.27}^{+0.18}$	${0.43}_{-0.21}^{+0.12}$	${0.64}_{-0.17}^{+0.21}$
$f_{\mathrm{cen}}^{2,\beta }$	${0.75}_{-0.35}^{+0.13}$	${0.62}_{-0.16}^{+0.13}$	${0.39}_{-0.21}^{+0.20}$
$f_{\mathrm{cen}}^{3,\beta }$	${0.73}_{-0.45}^{+0.13}$	${0.86}_{-0.12}^{+0.09}$	${0.62}_{-0.24}^{+0.19}$
$R_{1,\mathrm{off}}$	${0.33}_{-0.24}^{+0.20}$	—	—
$R_{2,\mathrm{off}}$	—	${0.47}_{-0.20}^{+0.19}$	—
$R_{3,\mathrm{off}}$	—	—	${0.50}_{-0.25}^{+0.22}$
$m_{\mathrm{lens}}$	${0.00}_{-0.01}^{+0.01}$	${0.00}_{-0.01}^{+0.01}$	${0.00}_{-0.01}^{+0.01}$
${\chi }_{\min }^2/\mathrm{dof}$	$30.6/29$	$25.2/29$	$34.4/29$

Equations

$$w_{\mathrm{rand}}(z_{\mathrm{cl}})=\frac{n_{\mathrm{sample}}(z_{\mathrm{cl}})}{n_{\mathrm{keep}}(z_{\mathrm{cl}})},$$

$${\Omega }_{\mathrm{eff}}(z_{\mathrm{cl}})=\frac{{\Omega }_{\mathrm{tot}}}{w_{\mathrm{rand}}(z_{\mathrm{cl}})},$$

$${\hat{N}}_{\alpha ,\beta }=\sum _{l;N_l\in N_{\alpha },z_l\in z_{\beta }}\frac{{\Omega }_{\mathrm{tot}}}{{\Omega }_{\mathrm{eff}}(z_l)},$$

$${\widehat{\Delta \Sigma }}_{l,\alpha ,\beta }(R)=\frac{1}{1+{\hat{m}}_{l,\alpha ,\beta }(R)}[\frac{1}{2{\hat{\mathcal{R}}}_{l,\alpha ,\beta }(R)}$$

$$\times \frac{1}{N_{ls}^{\alpha ,\beta }(R)}{\sum _{l,s;N_l\in N_{\alpha },z_l\in z_{\beta }}{w}_{ls}{⟨{\Sigma }_{\mathrm{cr}}^{-1}⟩}_{ls}^{-1}{e}_{+}({\mathbf{\theta }}_s)|}_{R={\chi }_l|{\mathbf{\theta }}_l-{\mathbf{\theta }}_s|}$$

$${-\frac{1}{N_{ls}^{\alpha ,\beta }(R)}\sum _{l,s;N_l\in N_{\alpha },z_l\in z_{\beta }}{w}_{ls}{⟨{\Sigma }_{\mathrm{cr}}^{-1}⟩}_{ls}^{-1}{c}_{+}({\mathbf{\theta }}_s)|}_{R={\chi }_l|{\mathbf{\theta }}_l-{\mathbf{\theta }}_s|}],$$

$${\Sigma }_{\mathrm{cr}}^{-1}(z_l,z_s)=4\pi G(1+z_l){\chi }_l\left(1-\frac{{\chi }_l}{{\chi }_s}\right)(z_l\le z_s)$$

$${⟨{\Sigma }_{\mathrm{cr}}^{-1}⟩}_{ls}={\int }_{z_l}^{\infty }\mathrm{d}{z}_s{\Sigma }_{\mathrm{cr}}^{-1}(z_l,z_s)P(z_s).$$

$$w_{ls}=\frac{{\Omega }_{\mathrm{tot}}}{{\Omega }_{\mathrm{eff}}(z_l)}{⟨{\Sigma }_{\mathrm{cr}}^{-1}⟩}_{ls}^2{w}_s,$$

$$w_s=\frac{1}{{\sigma }_e^2+e_{\mathrm{rms}}^2},$$

$$N_{ls}^{\alpha ,\beta }(R)={\sum _{l,s;N_l\in N_{\alpha },z_l\in z_{\beta }}{w}_{ls}|}_{R={\chi }_l|{\mathbf{\theta }}_l-{\mathbf{\theta }}_s|}.$$

$${\hat{m}}_{l,\alpha ,\beta }(R)=\frac{{\sum _{l,s;N_l\in N_{\alpha },z_l\in z_{\beta }}{m}_s{w}_{ls}|}_{R={\chi }_l|{\mathbf{\theta }}_l-{\mathbf{\theta }}_s|}}{{\sum _{l,s;N_l\in N_{\alpha },z_l\in z_{\beta }}{w}_{ls}|}_{R={\chi }_l|{\mathbf{\theta }}_l-{\mathbf{\theta }}_s|}},$$

$${\hat{\mathcal{R}}}_{l,\alpha ,\beta }(R)=1-\frac{{\sum _{l,s;N_l\in N_{\alpha },z_l\in z_{\beta }}{e}_{\mathrm{rms}}^2{w}_{ls}|}_{R={\chi }_l|{\mathbf{\theta }}_l-{\mathbf{\theta }}_s|}}{{\sum _{l,s;N_l\in N_{\alpha },z_l\in z_{\beta }}{w}_{ls}|}_{R={\chi }_l|{\mathbf{\theta }}_l-{\mathbf{\theta }}_s|}},$$

$${\hat{m}}_{\mathrm{sel},l,\alpha ,\beta }(R)=A_{\mathrm{sel}}{p}_{l,\alpha ,\beta }(R_2=0.3;R),$$

$${\widehat{\Delta \Sigma }}_{\alpha ,\beta }(R)={\widehat{\Delta \Sigma }}_{l,\alpha ,\beta }(R)-{\widehat{\Delta \Sigma }}_{r,\alpha ,\beta }(R),$$

$$P(\ln N|M,z)=\frac{1}{\sqrt{2\pi }{\sigma }_{\ln N|M,z}}\exp \left(-\frac{x^2(N,M,z)}{2{\sigma }_{\ln N|M,z}^2}\right),$$

$$x(N,M,z)$$

$$⟨\ln N⟩(M,z)$$

$${\sigma }_{\ln N|M,z}$$

$$N_{\alpha ,\beta }$$

$$S(M,z|N_{\alpha ,\min },N_{\alpha ,\max })$$

$$\equiv {\int }_{\ln N_{\alpha ,\min }}^{\ln N_{\alpha ,\max }}\mathrm{d}\ln NP(\ln N|M,z)$$

$$=\frac{1}{2}\left[\mathrm{erf}\left(\frac{x(N_{\alpha ,\max },M,z)}{\sqrt{2}{\sigma }_{\ln N|M,z}}\right)-\mathrm{erf}\left(\frac{x(N_{\alpha ,\min },M,z)}{\sqrt{2}{\sigma }_{\ln N|M,z}}\right)\right],$$

$${\rho }_{\mathrm{hm}}(r;M,z_l)={\bar{\rho }}_{m0}\left[1+{\xi }_{\mathrm{hm}}(r;M,z_l)\right].$$

$${\xi }_{\mathrm{hm}}(r;M,z_l)={\int }_0^{\infty }\frac{k^2\mathrm{d}k}{2{\pi }^2}{P}_{\mathrm{hm}}(k;M,z_l)j_0(kr),$$

$$\Sigma (R;M,z_l)$$

$$\Delta \Sigma (R;M,z_l)$$

$$\Delta {\Sigma }_{\alpha ,\beta }(R)$$

$$w_l^{\alpha ,\beta }(z)=\sum _{l,s;N_l\in N_{\alpha },z_l\in z_{\beta },z_l\in z}{w}_{ls}.$$

$$N_{\Delta \Sigma }^{\alpha ,\beta }$$

$$P_{\mathrm{hm}}(k;M,z_l)$$

$$f_{\mathrm{cen}}^{\alpha ,\beta }=f_0+f_N\ln \left(\frac{⟨N{⟩}_{\alpha ,\beta }}{N_{\mathrm{pivot}}}\right)+f_z\ln \left(\frac{1+⟨z_{\mathrm{cl}}{⟩}_{\alpha ,\beta }}{1+z_{\mathrm{pivot}}}\right),$$

Full text

1 Introduction

Clusters of galaxies are dominated by dark matter and therefore are useful sites for cosmological studies since $N$-body simulations can predict cluster observables reasonably well. The abundance and clustering of massive clusters and their time evolution are known to be sensitive to cosmological parameters such as the matter density ($\Omega_{\rm m}$), the normalization of the matter power spectrum ($\sigma_{8}$), and dark energy (see e.g., White et al., 1993; Eke et al., 1996; Kitayama & Suto, 1997; Haiman et al., 2001; Voit, 2005; Vikhlinin et al., 2009; Mantz et al., 2010; Rozo et al., 2010; Allen et al., 2011; Oguri & Takada, 2011; Weinberg et al., 2013; Planck Collaboration et al., 2016b). Clusters of galaxies also play an important role in understanding galaxy formation physics via their possible large environmental effect (e.g., Renzini et al., 2006; Kravtsov & Borgani, 2012).

Clusters of galaxies can be identified in optical, X-ray, and radio/mm/submm wavelengths. The recent development of wide-field optical imaging surveys makes an optical selection of clusters particularly powerful. This is because optical surveys take wide-field images with multiple photometric passbands from which we can select clusters of galaxies efficiently via the enhancement of galaxy number counts and derive photometric redshifts of clusters (e.g., Gladders & Yee, 2000). Ongoing and upcoming wide-field optical or infrared galaxy surveys allow us to study the cosmology and galaxy formation physics in great detail if systematic errors are under control. Such surveys include the Kilo-Degree Survey (KiDS; de Jong et al., 2013; Kuijken et al., 2015), the Dark Energy Survey (DES; Flaugher, 2005; Dark Energy Survey Collaboration et al., 2016), and the Subaru Hyper Suprime-Cam (HSC) survey (Aihara et al., 2018a, b) for ongoing surveys, and the Large Synoptic Survey Telescope (LSST; Ivezić et al., 2008), Euclid (Laureijs et al., 2011), and Wide-Field Infrared Survey Telescope (WFIRST; Spergel et al., 2015) for upcoming surveys.

Since theoretical predictions of cluster observables are primarily determined with respect to the halo mass for a given cosmological model, we need a valid statistical model to connect the halo mass and observed mass proxy in order to make full use of cluster samples. In optical surveys, a commonly-used mass proxy is an optical richness, which roughly corresponds to the number of red-sequence member galaxies above some luminosity threshold in each cluster (e.g., Rozo et al., 2009). Well-calibrated and unbiased mass-richness relations allow us to infer cluster masses from observed richness values.

Weak gravitational lensing provides a powerful means to constrain mass-observable relations of clusters. It is the deflection of light due to the intervening matter density field along the line-of-sight direction to produce a coherent distortion pattern in background galaxy shapes (for reviews, see e.g., Bartelmann & Schneider, 2001; Kilbinger, 2015; Mandelbaum, 2018b). Stacked weak lensing measurements statistically probe the projected average mass distribution of clusters with equal sensitivity to the dark and baryonic matter. Stacking shapes of background galaxies for a sample of clusters enhances the signal-to-noise ratio of the measurements. Previous studies have utilized the weak gravitational lensing effect to constrain mass-observable relations, including mass-richness relations (e.g., Johnston et al., 2007; Leauthaud et al., 2010; Okabe et al., 2013; von der Linden et al., 2014; Hoekstra et al., 2015; Battaglia et al., 2016; Simet et al., 2017; Melchior et al., 2017; Murata et al., 2018; Medezinski et al., 2018a; Miyatake et al., 2019; McClintock et al., 2019).

The Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) is a wide-field optical imaging survey with a 1.77 ${\rm deg}^{2}$ field-of-view camera on the 8.2-meter Subaru telescope (Miyazaki et al., 2012, 2015, 2018a; Komiyama et al., 2018; Furusawa et al., 2018; Kawanomoto et al., 2018). The HSC survey is unique in its combination of depth and high-resolution image quality, which allows us to detect clusters of galaxies over a wide redshift range up to $z_{\rm cl}\sim 1$ and to measure stacked lensing profiles around such high-redshift clusters with a lower shape noise than that of other ongoing surveys such as KiDS and DES.

In this paper, we present constraints on the relation between the optical richness and halo mass of optically-selected HSC clusters in the redshift range $0.1\leq z_{\rm cl}\leq 1.0$ (Oguri et al., 2018a) detected by the CAMIRA cluster-finding algorithm (Oguri, 2014). For this purpose, we conduct a joint analysis of the abundance and stacked lensing profiles from the first-year data catalogs of the Subaru HSC-SSP survey. For this work, we adapt and apply a pipeline developed in Murata et al. (2018), which was used for clusters at $0.1\leq z_{\rm cl}\leq 0.33$ in the Sloan Digital Sky Survey (SDSS) selected by the redMaPPer cluster-finding algorithm (Rykoff et al., 2014; Rozo & Rykoff, 2014; Rozo et al., 2015a, b; Rykoff et al., 2016). We model the probability distribution of the richness for a given halo mass and redshift $P(\ln N|M,z)$ without informative prior distributions for richness-mass parameters. We then use Bayes theorem to calculate the mass-richness relation $P(\ln M|N)$ in each redshift bin. In order to accurately model the abundance and stacked lensing profiles, we use Dark Emulator (Nishimichi et al., 2018), which is constructed from a suite of high-resolution $N$-body simulations. We employ an analytic model for the covariance matrix describing statistical errors for the abundance and stacked lensing profiles. We validate the analytic covariance matrix against realistic mock shear and halo catalogs (Shirasaki et al., 2019).

The structure of this paper is as follows. In Section 2, we briefly describe the Subaru HSC data and catalogs used in our richness-mass relation analysis. In Section 3, we describe measurements of the cluster abundance and stacked cluster lensing profiles. In Section 4, we summarize model ingredients for our richness-mass relation analysis. In Section 5, we show the resulting constraints on the mass-richness relation. We discuss the robustness of our results in Section 6. We conclude in Section 7.

Throughout this paper we use natural units where the speed of light is set equal to unity, $c=1$. We use $M\equiv M_{\rm 200m}=4\pi(R_{200{\rm m}})^{3}\bar{\rho}_{\rm m0}\times 200/3$ for the halo mass definition, where $R_{\rm 200m}$ is the spherical halo boundary comoving radius within which the mean mass density is $200$ times the present-day mean mass density. We note that in this paper we use a radius and density in comoving coordinates rather than in physical coordinates. We adopt the standard flat $\Lambda$-dominated cold dark matter model as the fiducial cosmological model with the parameters from the Planck15 (hereafter the Planck) result (Planck Collaboration et al., 2016a): $\Omega_{\rm b0}h^{2}=0.02225$ and $\Omega_{\rm c0}h^{2}=0.1198$ for the density parameters of baryon and cold dark matter, respectively, $\Omega_{\Lambda}=0.6844$ for the cosmological constant, $\sigma_{8}=0.831$ for the normalization of the matter fluctuation, and $n_{s}=0.9645$ for the spectral index. We also use cosmological parameters consistent with those from WMAP9 (hereafter the WMAP) results (Hinshaw et al., 2013) to compare with the results when assuming the Planck cosmological parameters: $\Omega_{\rm b0}h^{2}=0.02254$, $\Omega_{\rm c0}h^{2}=0.1142$, $\Omega_{\Lambda}=0.721$, $\sigma_{8}=0.82$, and $n_{s}=0.97$.

2 HSC first-year dataset

2.1 HSC-SSP survey

HSC is a wide-field prime focus camera with a $1.5$ ${\rm deg}$ diameter field-of-view mounted on the 8.2-meter Subaru telescope (Miyazaki et al., 2012, 2015, 2018a; Komiyama et al., 2018; Furusawa et al., 2018; Kawanomoto et al., 2018). With its unique combination of a wide field-of-view, a large aperture of the primary mirror, and excellent image quality, HSC enables us to measure lensing signals out to a relatively high redshift. Under the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP; Aihara et al., 2018a), the HSC is conducting a multi-band wide-field imaging survey over six years with 300 nights of Subaru time. The HSC-SSP survey consists of three layers: Wide, Deep, and UltraDeep. The Wide layer is designed for weak lensing science and aims at covering 1400 deg2 of the sky with five broad bands, ${\it grizy}$, with a $5\sigma$ point-source depth of $r\sim 26$. The $i$-band images are taken when the seeing is better since $i$-band images are used for galaxy shape measurements for weak lensing analysis, resulting in a median PSF FWHM of \sim$$0_{.}^{{}^{\prime\prime}}58 for $i$-band images for the HSC first-year shear catalog described in Section 2.3. The software pipelines that reduce the data are described in Bosch et al. (2018).

While the HSC-SSP Data Release 1 (Aihara et al., 2018b) is based on data taken on 61.5 nights between March 2014 and November 2015, in this paper we use HSC cluster, lensing shear, and photometric redshift catalogs based on the S16A internal release data of the HSC-SSP survey that was taken during March 2014 through April 2016, about 90 nights in total.

2.2 HSC CAMIRA cluster catalog

We use the CAMIRA (Cluster finding Algorithm based on Multi-band Identification of Red-sequence gAlaxies) cluster catalog from the S16A internal release data of the HSC-SSP Wide dataset presented in Oguri et al. (2018a), which was constructed using the CAMIRA algorithm (Oguri, 2014). The CAMIRA algorithm is a red-sequence cluster finder based on a stellar population synthesis model (Bruzual & Charlot, 2003) to predict colors of red-sequence galaxies at a given redshift and to compute likelihoods of being red-sequence galaxies as a function of redshift. In addition, Oguri et al. (2018a) calibrated the stellar population synthesis model with spectroscopic galaxies to improve its accuracy. In the CAMIRA algorithm, the richness corresponds to the number of red member galaxies with stellar mass $M_{*}\gtrsim 10^{10.2}M_{\odot}$ (roughly corresponding to a luminosity range of $L\gtrsim 0.2L_{*}$) within a circular aperture with radius $R\lesssim 1\ h^{-1}{\rm Mpc}$ in physical coordinates. The CAMIRA algorithm does not include a richness-dependent scale radius to define the richness, unlike the redMaPPer algorithm (Rykoff et al., 2012). The HSC images are deep enough to detect cluster member galaxies down to $M_{*}\sim 10^{10.2}M_{\odot}$ even at a cluster redshift $z_{\rm cl}\sim 1$, which allows a reliable cluster detection at such high redshifts without a richness incompleteness correction. The algorithm employs a spatially-compensated filter such that the background level is automatically subtracted in deriving the richness. The CAMIRA algorithm identifies a brightest cluster galaxy (BCG) for each cluster candidate that is defined by a peak in the three-dimensional richness map in RA, dec, and redshift space (Oguri, 2014; Oguri et al., 2018a). The cluster centers are defined as the locations of the identified BCGs. The mask-corrected richness $N$ and cluster photometric redshift $z_{\rm cl}$ are refined iteratively during the BCG identification process. The offset of BCG positions from matched X-ray cluster centers is investigated in Oguri et al. (2018a). The bias and scatter in photometric cluster redshifts of $\Delta z_{\rm cl}/(1+z_{\rm cl})$ are shown to be better than $0.005$ and $0.01$ respectively with 4$\sigma$ clipping over most of the redshift range by using available spectroscopic redshifts of BCGs. We use the cluster catalog without applying a bright-star mask (Oguri et al., 2018a).

The catalog contains $1921$ clusters of $0.1\leq z_{\rm cl}\leq 1.1$ and richness $N\geq 15$ with almost uniform completeness and purity over the sky region. We use 1747 CAMIRA clusters with $15\leq N\leq 200$ and $0.1\leq z_{\rm cl}\leq 1.0$. The total area for the CAMIRA clusters is estimated to be $\Omega_{\rm tot}=232.8$ ${\rm deg}^{2}$ (Oguri et al., 2018a), which is larger than the area for the lensing shear catalog presented in Section 2.3. The CAMIRA algorithm calculates the mask area fraction $f_{\rm mask}$ to correct for the richness and adopts the minimum mask fraction value to reject detections (Oguri, 2014). Thus areas for clusters with lower redshifts are slightly smaller since the richness is defined within a circular aperture with a radius $R\lesssim 1\ h^{-1}{\rm Mpc}$ in physical coordinates and clusters at lower redshift have a higher rejection rate from the mask cut. We use a random catalog to estimate this effect by injecting clusters at the catalog level into the footprint to calculate the rejection rate from the masking cut as a function of cluster redshift. We then define a weighting function to quantify the masking effect as

[TABLE]

where $n_{\rm sample}(z_{\rm cl})$ and $n_{\rm keep}(z_{\rm cl})$ are the number of injected clusters and the number of clusters not rejected by the masking cut, respectively (similarly defined in Murata et al., 2018). In addition, we define the effective area of the CAMIRA cluster catalog at a given redshift as

[TABLE]

to account for the detection efficiency as a function of redshift in measurements below, although this effect is not very large for the CAMIRA clusters with $w_{\rm rand}(z_{\rm cl}=0.1)=1.02$ at most. Here we note that $\Omega_{\rm eff}(z_{\rm cl})\leq\Omega_{\rm tot}$. We use the same random catalog to measure lensing signals for the subtraction of systematics in Section 3.2.

2.3 HSC weak lensing shear catalog

We employ the HSC first-year shear catalog (Mandelbaum et al., 2018a, c) based on the S16A internal release data for weak lensing measurements around HSC CAMIRA clusters described in Section 2.2. The galaxy shapes are measured on coadded $i$-band images with the re-Gaussianization moment-based method (Hirata & Seljak, 2003) and fully described in Mandelbaum et al. (2018c)111In this method, galaxy shapes are defined in terms of distortion: $(e_{1},e_{2})=(e\cos 2\phi,e\sin 2\phi)$ with $e=(a^{2}-b^{2})/(a^{2}+b^{2})$ where $a$ and $b$ are the major and minor axes of galaxy shape respectively, and $\phi$ indicates the position angle with respect to the RA/dec coordinate system (Bernstein & Jarvis, 2002). . This method has been applied extensively to SDSS data, and thus the systematics of the method are well understood (Mandelbaum et al., 2005; Reyes et al., 2012; Mandelbaum et al., 2013). Both shape uncertainties and biases are estimated per galaxy with simulations created using an open source software package GalSim (Rowe et al., 2015) with galaxy samples from the Hubble Space Telescope COSMOS survey (see Mandelbaum et al., 2018c, for more details). More specifically, Mandelbaum et al. (2018c) estimate multiplicative bias $m$, additive bias $(c_{1},c_{2})$, intrinsic root-mean-square ellipticity $e_{\rm rms}$, and shape measurement error $\sigma_{e}$ for the galaxy ensemble in the simulation, and define interpolation functions to produce an estimate of that quantity for each galaxy in the real data. The systematic uncertainty in the overall shear calibration is estimated to be $0.01$ in Mandelbaum et al. (2018c).

Mandelbaum et al. (2018a) applied a number of cuts to satisfy the requirements for carrying out first-year weak lensing cosmology analyses. For example, the catalog is constructed using regions of sky with approximately full depth in all five bands to ensure the homogeneity of the sample. Mandelbaum et al. (2018a) also limited the cmodel magnitude (see Bosch et al., 2018, for the definition of cmodel magnitude) with ${\it i}_{\rm cmodel}<24.5$, which is conservative compared to the i-band magnitude limit of \sim$$26.4 ($5\sigma$ for point sources; Aihara et al., 2018a). As a result, the first-year shear catalog covers $\Omega_{\rm lens}=136.9~{}{\rm deg}^{2}$ in total with six distinct fields (XMM, GAMA09H, GAMA15H, HECTOMAP, VVDS, and WIDE12H). We note that the area coverage of the shear catalog is smaller than that of the CAMIRA clusters ($\Omega_{\rm tot}=232.8$ ${\rm deg}^{2}$) due to the conservative cuts for the shape measurements. Mandelbaum et al. (2018a), Oguri et al. (2018b), and Hikage et al. (2019) performed extensive null tests of the shear catalog to show that the shear catalog satisfies the requirements for HSC first-year weak lensing analyses using cosmic shear and galaxy-galaxy lensing. Even after relatively conservative cuts, the HSC first-year shear catalog includes galaxy shapes with a high source number density, $24.6$ (raw) and $21.8$ (effective) ${\rm arcmin}^{-2}$ (Mandelbaum et al., 2018a), which enables us to measure lensing signals around high redshift HSC CAMIRA clusters ($z_{\rm cl}\leq 1.0$).

2.4 HSC photometric redshift catalog

We use a photometric redshift (photo-$z$) catalog (Tanaka et al., 2018) for the source galaxies in Section 2.3 estimated from the S16A internal release data of the HSC five broadband photometry. In the HSC survey, several different codes are employed to estimate the photometric redshifts: a machine-learning code based on a self-organizing map (MLZ), a classical template-fitting code (Mizuki), an empirical polynomial fitting code (DEmP; Hsieh & Yee, 2014), an extended (re)weighing method to find the nearest neighbors in color/magnitude space from a reference spectroscopic redshift sample (NNPZ), a neural network code using PSF-matched aperture (afterburner) photometry (Ephor AB), and a hybrid code combining machine learning with template fitting (FRANKEN-Z; Speagle et al., 2019, Speagle et al. in prep.) The codes are trained, validated, and tested with spectroscopic and grism redshifts as well as COSMOS 30-band data with high accuracy photo-$z$ (Ilbert et al., 2009; Laigle et al., 2016) in Tanaka et al. (2018). The photo-$z$ estimation is most accurate in the range $0.2\lesssim z\lesssim 1.5$, where the HSC filter set straddles the 4000${\rm\AA}$ break (Tanaka et al., 2018).

Among these catalogs from different codes, we choose MLZ as the fiducial photometric redshift catalog for the lensing measurement in Section 3.2, while we also use the other photometric redshift catalogs to check the robustness of our results to our choice of MLZ (see Section 6.1.2). In this paper, we use the redshift probability distribution functions (PDFs), $P(z)$, and randomly (i.e., Monte Carlo) sampled point estimates drawn from the full PDFs, $z_{\rm mc}$ from Tanaka et al. (2018) for the lensing measurements and the source galaxy selection. We also use the best point estimates $z_{\rm best}$ (see section 4.2 of Tanaka et al., 2018) only for lensing covariance estimation as described in Appendix A.1. We correct for the effect of photometric redshift bias on the lensing measurements (More et al. in prep) using COSMOS 30-band photo-$z$ data (see Section 3.2).

3 Measurement

We describe the measurement method for the cluster abundance in Section 3.1 and the stacked cluster lensing profile in Section 3.2.

3.1 Cluster abundance

We use the abundance of CAMIRA clusters in given richness and redshift bins as cluster observables to constrain the richness-mass relation of the clusters. We divide the CAMIRA clusters into 12 bins for the abundance measurements with four richness and three redshift bins as shown in Table 1. We use the point estimate of richness and redshift of the clusters to calculate the abundance (i.e., we ignore any errors in the richness estimation and the cluster photometric redshift).

We measure the abundance of the clusters in each bin corrected for the detection efficiency (Murata et al., 2018) as

[TABLE]

where $N_{\alpha}$ and $z_{\beta}$ denote the $\alpha$-th richness bin and $\beta$-th redshift bin in Table 1, respectively. The summation runs over all clusters in each richness and redshift bin. The factor $\Omega_{\rm tot}/\Omega_{\rm eff}$ corrects for the detection efficiency as discussed in Section 2.2. Equation (3) gives an estimate of the abundance of the clusters we could observe for the survey area of $\Omega_{\rm tot}$ without the mask effect. We then do not need to include the mask effect in the model prediction given in Section 4.2.

The numbers of clusters in the $\alpha$-th richness and $\beta$-th redshift bin before and after the detection efficiency correction are given in Table 1 by $N_{\alpha,\beta}^{\rm uncorr}$ and $\widehat{N}_{\alpha,\beta}$, respectively. The correction is less than \sim$$1\% for all bins.

3.2 Stacked cluster lensing profile

We cross-correlate the positions of CAMIRA clusters with the shapes of background galaxies to measure the average excess mass density profile around the clusters (hereafter the stacked lensing profile). We follow the procedure in Mandelbaum et al. (2018a, c) to estimate the stacked lensing profile for a sample of CAMIRA clusters for $\alpha$-th richness and $\beta$-th redshift bin in Table 1 as

[TABLE]

where the subscripts $l$ and $s$ stand for lens (cluster) and source, respectively, and $e_{+}$ and $c_{+}$ are the tangential component of the source galaxy ellipticity and the additive bias with respect to the cluster center, respectively. The summation runs over all pairs of clusters and source galaxies after a source selection cut described below in a given comoving transverse separation bin of $R=\chi_{l}|\bm{\theta}_{l}-\bm{\theta}_{s}|$, where $\chi_{l}$ is the comoving distance to each cluster, and $\bm{\theta}_{l}$ and $\bm{\theta}_{s}$ are angular positions of the lenses and the sources, respectively. We use 11 radial bins that are equally spaced logarithmically from $0.42$ $h^{-1}{\rm Mpc}$ to $14.0$ $h^{-1}{\rm Mpc}$ in comoving coordinates. We use area-weighted mean values of comoving radii for the representative radial values. The value for the first inner bin is $0.51~{}h^{-1}{\rm Mpc}$. We do not use lensing profiles at $<0.42h^{-1}$Mpc to avoid a possible dilution effect (Medezinski et al., 2018b) and an increased blending effect (Murata et al. in prep.) at such small radii.

The critical surface mass density is defined for a system of lens and source for a flat universe as

[TABLE]

and [math] when $z_{l}>z_{s}$, where $G$ is the gravitational constant. We average this over the source photometric redshift PDF for each lens-source pair as

[TABLE]

In addition, the lens-source pair weight is given as

[TABLE]

and $w_{s}$ is the source weight defined as

[TABLE]

where $\sigma_{e}$ is the measurement error of galaxy ellipticity and $e_{\rm rms}$ is the intrinsic per-component root-mean-square ellipticity estimated in Mandelbaum et al. (2018c) for each source galaxy. The lensing weight factor $\Omega_{\rm tot}/\Omega_{\rm eff}(z_{l})$ corrects for the effective area difference as a function of cluster redshift as in the abundance estimator in equation (3). We note that this correction to the lensing estimator is small (less than $1\%$) for all bins. The denominator in equation (4) is the weighted number of lens-source pairs in each separation bin computed as

[TABLE]

The multiplicative bias correction is estimated as

[TABLE]

where $m_{s}$ is the estimated multiplicative bias for each source galaxy (Mandelbaum et al., 2018c). The shear responsivity factor represents the statistically-averaged response of galaxy distortions to small shears (Kaiser et al., 1995; Bernstein & Jarvis, 2002), and is measured as

[TABLE]

which is found to be around $0.83$ for all the bins.

Here we describe the source selection cut. As shown in e.g., Medezinski et al. (2018b), it is important to use a secure background galaxy sample for cluster weak lensing in order to minimize the dilution effect by cluster member galaxies especially at inner radii. We use the Pcut method (Oguri, 2014) as our fiducial background source selection cut. In the Pcut method, we select source galaxies whose photometric redshift PDFs lie mostly beyond the cluster redshift plus some threshold $\Delta z$: $\int_{z_{l}+\Delta z}^{\infty}{\rm d}zP(z)>0.98$. We apply an additional cut on the randomly sampled redshift from PDFs, $z_{\rm mc}<2.5$, for each lens-source pair in equation (4). We use $\Delta z=0.1$ for the fiducial cut since Medezinski et al. (2018b) shows the possible dilution effect is negligible for the Pcut method with this threshold value at $R\gtrsim 0.5h^{-1}{\rm Mpc}$. In Section 6.1.2, we also check the robustness of our choice of this source selection cut by using the Pcut method with $\Delta z=0.2$ and the color-color cut in Medezinski et al. (2018b). The color-color cut uses the color-color space of $g-i$ vs $r-z$ for HSC, where cluster red-sequence member galaxies can be well isolated from background galaxies.

In addition, we correct for selection bias from the lower cut on the resolution factor at $R_{2}=0.3$ in the shape catalog (Mandelbaum et al., 2018c) as

[TABLE]

where $A_{\rm sel}=0.0087$ and $p_{l,\alpha,\beta}(R_{2}=0.3;R)$ is calculated by summation of lens-source weights, $w_{ls}$, in each of the radial, richness, and redshift bins. This correction is found to be $m_{\rm sel}\sim 0.01$ for all bins. We then correct for the redshift variation of the intrinsic shape noise (Mandelbaum et al., 2018c) by employing $m=0.03$ for $1.0\leq z_{\rm best}\leq 1.5$ and $m=-0.01$ for the other redshift ranges, and by averaging this over lens-source pairs in each bin. This correction is found to be a positive multiplicative bias, but less than $0.01$ for all bins.

We also correct for systematic bias effects of the photometric redshift on the lensing measurements by comparing the critical surface mass density from the photometric redshift estimates against that from the accurate photometric redshifts in the COSMOS 30-band catalog (Ilbert et al., 2009; Laigle et al., 2016) with the weak lensing weight $w_{ls}$ multiplied by a self-organizing map weight $w_{\rm SOM}$ which adjusts the COSMOS 30-band photo-$z$ sample to mimic the HSC source galaxy sample (More et al. in prep.). We assume that the COSMOS 30-band photometric redshift estimates are sufficiently accurate due to the larger numbers of bands. This method is based on Nakajima et al. (2012) in which the method was applied to SDSS data, and was also applied in Miyatake et al. (2019) (see their equation 11 for more details). The debias factor is found to be $m\sim 0.01$ for CAMIRA clusters with $0.1\leq z_{\rm cl}\leq 0.7$ and $m\sim 0.02$ for CAMIRA clusters with $0.7\leq z_{\rm cl}\leq 1.0$ for the fiducial photo-$z$ code MLZ and the fiducial Pcut with $\Delta z=0.1$. The debias factor is similar for the other photo-$z$ catalogs and photo-$z$ cuts. Our result indicates that the photo-$z$ bias correction is not very large when we average the critical surface mass density over the photo-$z$ PDF in equation (6) with the Pcut method or the color-color cut for the shape catalog. We also estimate the uncertainties of these correction factors for each richness and redshift bin using the jackknife resampling technique (Efron, 1982) with ten subsamples of the COSMOS 30-band catalog, where we recalculate the SOM weight for each jackknife resampled subsample. We estimate the photo-$z$ bias correction uncertainties as $\sigma_{\alpha,\beta,{\rm photoz}}\sim 0.001,0.002,0.004$ for $\beta=1,2,3$, respectively, for the fiducial photo-$z$ catalog and source selection cut. We note that we ignore the impact of photo-$z$ biases, outliers, and the limited field variance of galaxies due to the small area in the COSMOS 30-band catalog, thus the uncertainties might be underestimated. These values are used for marginalization together with the systematic uncertainty in the overall calibration of the shear $\sigma_{\rm shear}=0.01$ (Mandelbaum et al., 2018c) in Section 4.3 222 We note that we apply this photo-$z$ correction for all photo-$z$ catalogs except for FRANKEN-Z since we do not have the photo-$z$ estimates for the galaxies with HSC photometry in the COSMOS 30-band catalog from FRANKEN-Z. .

After the above corrections for $\widehat{\Delta\Sigma}_{l,\alpha,\beta}(R)$, we also subtract the lensing measurement around random points as

[TABLE]

where $\widehat{\Delta\Sigma}_{r,\alpha,\beta}(R)$ replaces the clusters with random points in the estimator for $\widehat{\Delta\Sigma}_{l,\alpha,\beta}(R)$. This subtraction allows us to measure the excess mass density profile with respect to the background density as stressed in Sheldon et al. (2004) and Mandelbaum et al. (2005) (also see Singh et al., 2017, for a recent detailed study). The random subtraction can also correct for an additive bias due to shear systematics including point-spread function ellipticity errors (Mandelbaum et al., 2005). We use a random catalog of the CAMIRA clusters presented in Section 2.2 with the same richness and redshift distributions as the data, and there are 100 times as many random points as real clusters in each redshift and richness bin.

4 Forward modeling of cluster observables

We adopt a forward modeling approach to model the abundance and stacked lensing profiles (e.g., Zu et al., 2014; Murata et al., 2018; Costanzi et al., 2018) for a fixed cosmological model. In this approach, we model the probability distribution of the richness for a given halo mass and redshift, $P(\ln N|M,z)$. An alternative approach models the probability distribution of halo mass for a given richness and redshift, $P(\ln M|N,z)$, as in some previous works (e.g., Baxter et al., 2016; Simet et al., 2017; Melchior et al., 2017; McClintock et al., 2019).

4.1 Richness-mass relation

Following Lima & Hu (2005), we assume that the probability distribution of the observed richness for halos with a fixed mass and redshift is given by a log-normal distribution as

[TABLE]

where $x(N,M,z)$ models the mean relation of $\ln N$ parametrized by four model parameters, $A$, $B$, $B_{z}$, and $C_{z}$ as

[TABLE]

Hence $x(N,M,z)=0$ gives the mean relation of $\ln N$ (also the median relation for $N$) as

[TABLE]

We adopt $M_{\rm pivot}=3\times 10^{14}h^{-1}M_{\odot}$ for the pivot mass scale and $z_{\rm pivot}=0.6$ for the pivot cluster redshift. We discuss the validity of including redshift evolution parameters $B_{z}$ and $C_{z}$ in Section 6.3. In this work, we ignore errors on $z_{\rm cl}$ for simplicity given the small bias and scatter of the cluster redshifts (Oguri et al., 2018a) as described in Section 2.2.

In addition, we assume that the scatter of the richness around the mean relation at a fixed halo mass and redshift can be parametrized by four parameters, $\sigma_{0}$, $q$, $q_{z}$, and $p_{z}$ as

[TABLE]

We also discuss the validity of including redshift evolution parameters, $q_{z}$ and $p_{z}$, in Section 6.3. We only consider the parameter regions of $\{\sigma_{0},q,q_{z},p_{z}\}$ that result in $\sigma_{\ln N|M,z}>0$ for all halo mass and redshift for the parameter estimation. In this treatment, $\sigma_{\ln N|M,z}$ should be effectively regarded as a total scatter, including contributions from the intrinsic scatters, the richness measurement errors, the halo orientation effect (e.g., Dietrich et al., 2014), the projection effect (e.g., Costanzi et al., 2019), and any other source of observational scatter that may be present.

In summary, we model the richness-mass relation in equation (14) with eight model parameters. We constrain these parameters and check whether this model can reproduce measurements of the cluster abundance and stacked lensing profiles simultaneously with an acceptable value of $\chi^{2}_{\rm min}/{\rm dof}$.

4.2 Abundance in richness and redshift bins

Once we fix the richness-mass relation parameters of $P(\ln N|M,z)$, we can predict the abundance of CAMIRA clusters for a given cosmology. For the $\alpha$-th richness bin ($N_{\alpha,{\rm min}}\leq N\leq N_{\alpha,{\rm max}}$) and the $\beta$-th redshift bin ($z_{\beta,{\rm min}}\leq z\leq z_{\beta,{\rm max}}$), the abundance of the clusters for the total survey area is given as

[TABLE]

where $\chi^{2}(z)/H(z)$ denotes the comoving volume per unit redshift interval and unit steradian, and $\mathrm{d}n/\mathrm{d}M$ is the halo mass function at redshift $z$ for a fixed cosmological model. The selection function of halo mass at a fixed redshift in the richness bin is calculated by integrating the log-normal distribution of $P(\ln N|M,z)$ over the richness range as

[TABLE]

where erf$(x)$ is the error function.

4.3 Stacked cluster lensing profile in richness and redshift bins

The stacked lensing profile of halos with mass $M$ at redshift $z_{l}$ probes the average radial profile of the matter distribution around halos, $\rho_{\rm hm}(r;M,z_{l})$. Assuming statistical isotropy in the cluster detections, the average matter distribution around the halos is one-dimensional as a function of separation from the halo center $r$, where $r$ is in comoving coordinates. We express the average matter density profile with the cross-correlation function between the halo distribution and the matter density fluctuation field, $\xi_{\rm hm}(r;M,z_{l})$, as

[TABLE]

We note that we use the present-day mean matter density, $\bar{\rho}_{\rm m0}$, since we use the comoving density. The cross-correlation function is connected to the cross-power spectrum, $P_{\rm hm}(k;M,z_{l})$, via the Fourier transform as

[TABLE]

where $j_{0}(x)$ is the zeroth-order spherical Bessel function. The surface mass density profile is obtained from a projection of the three-dimensional matter profile along the line-of-sight direction as

[TABLE]

where $J_{0}(x)$ is the zeroth-order Bessel function and $R$ is the projected separation from the halo center in comoving coordinates. The excess surface mass density profile around halos, which is the direct observable from the stacked cluster lensing measurement, is given as

[TABLE]

where $\left\langle\Sigma(R;M,z_{l})\right\rangle_{<R}$ is the average of $\Sigma(R;M,z_{l})$ within a circular aperture of radius $R$, and $J_{2}(x)$ is a second-order Bessel function.

We can compute the model prediction for the stacked lensing profile accounting for the distribution of halo masses and redshifts for CAMIRA clusters in the $\alpha$-th richness and $\beta$-th redshift bin as

[TABLE]

The term in square brackets accounts for the non-linear contribution of reduced shear, which might not be negligible at very small radii (e.g., Johnston et al., 2007), where $\langle 1/\Sigma_{\mathrm{cr}}\rangle_{\alpha,\beta}(R)$ is measured from pairs of CAMIRA clusters and source galaxies in each radial bin for the $\alpha$-th richness and $\beta$-th redshift cluster bin. The lens redshift weight of $w_{l}^{\alpha,\beta}(z)$ is introduced to account for the weight distribution of lens redshift in the lensing measurement, and is calculated as

[TABLE]

More specifically, we compute $w_{l}^{\alpha,\beta}(z)$ as follows. First we divide the lens-source pairs into 12 lens redshift bins for each $\alpha$-th richness and $\beta$-th redshift bin, which are linearly spaced in $z_{l}\in[z_{\beta,{\rm min}},z_{\beta,{\rm max}}]$. We then estimate the weight in equation (25) over all sources in each lens redshift bin. We interpolate the weight values linearly as a function of lens redshift. The normalization factor in the denominator of equation (24) is similar to the abundance prediction in equation (18), but is defined accounting for the lens redshift weight as

[TABLE]

This lens redshift weight changes the model prediction by less than $1$–$2\%$ in the lensing profile amplitude compared to that without accounting for the lens redshift weight.

The identified BCGs as cluster centers can be off-centered from the true halo centers (e.g., Lin et al., 2004; Rozo & Rykoff, 2014; Oguri, 2014; Oguri et al., 2018a). We marginalize over the effect of off-centered clusters on the lensing profiles in equation (24) by modifying the halo-matter cross-power spectrum in equations (22) and (23) for CAMIRA clusters in the $\alpha$-th richness and $\beta$-th redshift bin as

[TABLE]

following Oguri & Takada (2011) and Hikage et al. (2012, 2013). Here $f_{{\rm cen}}^{\alpha,\beta}$ is a parameter to model the fraction of centered clusters in the $\alpha$-th richness and $\beta$-th redshift bin, and $(1-f_{{\rm cen}}^{\alpha,\beta})$ is the fraction of off-centered clusters. While Murata et al. (2018) assigned a model parameter describing the centering fraction independently for each lensing bin, we employ an empirical parametrization $f_{\rm cen}^{\alpha,\beta}$ that depends on the average richness and redshift in each bin (see Table 1) to reduce the number of model parameters as

[TABLE]

with three model parameters ($f_{0}$, $f_{N}$, and $f_{z}$), where we use $N_{\rm pivot}=25$ and $z_{\rm pivot}=0.6$ and we restrict the parameter regions of $f_{0}$, $f_{N}$, and $f_{z}$ such that $0<f_{{\rm cen}}^{\alpha,\beta}<1$ for all richness and redshift bins. The function $p_{\rm off}(r;R_{\beta,\rm off})$ is the normalized one-dimensional radial profile of detected centers by the CAMIRA algorithm with respect to the true halo center, for which we assume a Gaussian distribution (e.g., Johnston et al., 2007; Oguri & Takada, 2011) given as $p_{\rm off}(r;R_{\beta,\rm off})\propto\exp(-r^{2}/2R^{2}_{\beta,\rm off})$, where $R_{\beta,\rm off}$ is a parameter to model the typical off-centering radius in the $\beta$-th lens redshift bins. The Fourier transform of this function is denoted as $\tilde{p}_{\rm off}(k;R_{\beta,\rm off})=\exp(-k^{2}R_{\beta,\rm off}^{2}/2)$. Since the aperture radius of CAMIRA clusters is independent of the richness, with aperture size of \sim$$1\ h^{-1}{\rm Mpc} in physical coordinates (Oguri et al., 2018a), we use one model parameter ($R_{\beta,\rm off}$) for each $\beta$-th redshift bin, which is common for all richness bins, with a flat prior from $10^{-3}$ to $0.5\times(1+\langle z_{\rm cl}\rangle_{\beta})h^{-1}{\rm Mpc}$. We note that Oguri et al. (2018a) investigated the offset distribution of centers detected using the CAMIRA algorithm from centers of matched X-ray clusters, and showed that in most cases the offset is less than $0.5h^{-1}{\rm Mpc}$ in physical transverse distances.

We also marginalize over the shape calibration and photometric redshift bias uncertainties by introducing a nuisance parameter $m_{\rm lens}$ following a Gaussian prior distribution with a mean of zero and a standard deviation of $0.01$ to change the lensing model prediction as

[TABLE]

where $m_{\alpha,\beta}=m_{\rm lens}\times\sqrt{\sigma_{\rm shear}^{2}+\sigma_{\alpha,\beta,\rm photoz}^{2}}/0.01$. Mandelbaum et al. (2018c) show that the systematic uncertainty in the overall calibration of the shear is $\sigma_{\rm shear}=0.01$, and we use the photometric redshift bias uncertainty of $\sigma_{\alpha,\beta,\rm photoz}$ derived in Section 3.2. We expect a large correlation in the residual shape calibration and photometric redshift bias among different richness and redshift bins, and thus we conservatively use just one parameter $m_{\rm lens}$ for all bins.

Once the halo mass function ${\rm d}n(M,z)/{\rm d}M$ and the three-dimensional halo-matter cross-correlation $\xi_{\rm hm}(r;M,z)$ or $P_{\rm hm}(k;M,z)$ are provided for a given cosmological model (see Section 4.4), we can calculate the model prediction of the stacked lensing profile in each richness and redshift bin. The model is specified by 15 parameters in total for the entire cluster sample of $15\leq N\leq 200$ and $0.1\leq z_{\rm cl}\leq 1.0$: eight parameters $\{A,B,B_{z},C_{z},\sigma_{0},q,q_{z},p_{z}\}$ for the richness-mass relation $P(\ln N|M,z)$, six parameters $\{f_{0},f_{N},f_{z},R_{1,\rm off},R_{2,\rm off},R_{3,\rm off}\}$ for the off-centering effect, and one nuisance parameter ($m_{\rm lens}$) for the uncertainty in lensing amplitudes. We use the FFTLog algorithm (Hamilton, 2000) for Fourier transforms, which allows a fast, but sufficiently accurate and precise computation of the model prediction.

4.4 $N$-body simulation based halo emulator for the mass function and the lensing profile

The model predictions must be accurate in order to estimate the model parameters in an unbiased way. For this purpose, cosmological $N$-body simulations are one of the methods commonly used in the literature. Here, we use the database generated by the Dark Quest campaign (Nishimichi et al., 2018) to predict the halo mass function and the halo-matter cross-correlation function.

Nishimichi et al. (2018) develop a scheme called Dark Emulator for predicting statistical quantities of halos, including the mass function, the halo-matter cross-correlation function, the halo-halo auto-correlation function as a function of halo mass, redshift, separation length, and cosmological parameters, based on a series of high-resolution, cosmological $N$-body simulations. The simulation suite is composed of cosmological N-body simulations for 101 cosmological models within a flat $w$CDM framework, which are sampled around the Planck cosmological parameters. The simulations trace the nonlinear evolution of $2048^{3}$ particles in a box size of $1~{}$or$~{}2h^{-1}{\rm Gpc}$ on a side with mass resolution of \sim$$10^{10}h^{-1}M_{\odot} or $\sim 8\times 10^{10}h^{-1}M_{\odot}$, respectively. There are 21 redshift bins for each simulation realization across the range $0\leq z\leq 1.47$. To identify dark matter halos, Rockstar (Behroozi et al., 2013) is employed. The halo mass defined for the simulations includes all particles within a distance $R_{\rm 200m}$ from the halo center. The minimum halo mass of the emulator is $10^{12}h^{-1}M_{\odot}$. Throughout this paper, we set $M_{\rm min}=10^{12}h^{-1}M_{\odot}$ and $M_{\rm max}=2\times 10^{15}h^{-1}M_{\odot}$ for the minimum and maximum halo masses, respectively, to evaluate the halo mass integration in the model predictions of abundances and stacked lensing profiles.

For the Planck cosmological model, instead of the prediction of the Dark Emulator code, we use a simpler interpolation scheme for a fixed cosmological model. Specifically, we employ exactly the same cosmological parameters as their fiducial cosmological model, for which 24 independent realizations of high-resolution runs are available. The relevant statistics are interpolated as a function of mass and redshift (see Murata et al., 2018). On the other hand, we use Dark Emulator presented in Nishimichi et al. (2018) to interpolate over the cosmological parameter space and compute the predictions for the WMAP cosmological model. This relies on data compression based on Principal Component (PC) Analysis followed by Gaussian Process Regression for each PC coefficient. This is done for a subset of 80 cosmological models. Note that the realizations for the Planck cosmology are not used in the regression, but rather are used as part of the validation set to assess the performance of the emulator.

We estimate impacts of the uncertainties in the emulator on constraints of the mass-richness parameters as follows. For the Planck cosmology, we shift the halo mass function and the halo-matter cross-correlation function in each bin of halo mass, redshift, and separation length by one standard deviation uncertainty from the 24 realizations to opposite sides, and we calculate $\chi^{2}$ at the best-fit mass-richness relation parameters for the fiducial analysis below. We find that $\chi^{2}=106.9$ while we have $\chi^{2}=107.0$ for the fiducial emulator as shown in Table 2, suggesting that the emulator precision for the Planck cosmology is high enough for our analysis. For the WMAP cosmology, we estimate the impacts using the outputs of Dark Emulator for the Planck cosmology by comparing it with the fiducial emulator for the Planck cosmology. Since the realizations for the Planck cosmology are not used in the regression, differences between the 24 realizations and the Dark Emulator reveal the typical impact of the uncertainties in the emulator from the Dark Emulator for the WMAP cosmology. The errors for the model parameters are consistent between the two emulators. On the other hand, a shift for the median value of $A$ is $\sim 0.35$ compared to the error width, while shifts for the other parameters are smaller than $\sim 0.15$ compared to the error widths. This suggests that systematic uncertainties from the emulator for the WMAP cosmology increase the error for $A$ by $\sim 6\%$ and the errors for other parameters by $\sim 1\%$. Since these impacts of the emulator precision do not change the constraints on the mass-richness relation parameters very significantly, we ignore these uncertainties in this paper for simplicity.

4.5 Covariance

We must estimate the covariance describing the statistical uncertainties of the stacked lensing profiles and the abundance. The covariance consists of shape noise covariance from the finite number of lens-source pairs, Poisson noise for the abundance from the finite number of the clusters, and sample covariance from an imperfect sampling of the fluctuations in large-scale structure within a finite survey volume (e.g., Hu & Kravtsov, 2003; Takada & Bridle, 2007; Oguri & Takada, 2011; Takada & Hu, 2013; Hikage & Oguri, 2016; Takahashi et al., 2018; Shirasaki & Takada, 2018). We estimate the shape noise covariance for the stacked lensing profiles directly from the data (e.g., Murata et al., 2018) by repeating lensing measurements (equation 4) for source galaxies with their orientations randomized $15,000$ times. The shape noise covariance estimated in this manner accounts for the survey geometry and the inhomogeneous distribution of source galaxies. As discussed in Appendix A.1, we adopt an analytic halo model (Cooray & Sheth, 2002) to compute the sample covariance and the covariance for the abundance assuming that the distribution of clusters and lensing fields obey a Gaussian distribution. We note that we compute the shape noise covariance and the sample covariance for each setup of cosmological parameters, photo-$z$ catalog, and source selection cut as described in Appendix A.1.

In Appendix A.2, we validate the analytic covariance model for the sample covariance by comparing it with the covariance estimation from 2268 mock catalogs of the source galaxies and the CAMIRA clusters for the HSC footprint, which are generated from full-sky ray-tracing simulations with halo catalogs (Shirasaki et al., 2019) based on methods described in Shirasaki et al. (2017) (see also Shirasaki & Yoshida, 2014). The mock catalogs are based on full-sky, light-cone cosmological simulations constructed from sets of $N$-body simulations (Takahashi et al., 2017) with WMAP cosmological parameters (Hinshaw et al., 2013). The lensing effects at a given angular position are computed by a ray-tracing simulation through the foreground matter distribution from the multiple lens-plane algorithm (Hamana & Mellier, 2001; Shirasaki et al., 2015). Each source plane is given in HEALPix format (Górski et al., 2005) with an angular resolution of about $0.43$ arcmin. The mock galaxy shape catalog accounts for various effects as in the real data, including survey geometry, the inhomogeneous angular distribution of source galaxies, statistical uncertainties in the photo-$z$ estimation of each galaxy, and variations in the lensing weight from observational conditions and galaxy properties, since the mock catalog is constructed based on the real shape catalog (see Shirasaki et al., 2019, for more details). We also construct mock catalogs of the CAMIRA clusters by using the catalog of halos in each light-cone simulation realization in Takahashi et al. (2017), which are identified with Rockstar (Behroozi et al., 2013). We refer the readers to Appendix A.2 for more details.

Figure 1 shows the diagonal components of the covariance matrix for the fiducial analysis with Planck cosmological parameters, the fiducial photo-$z$ catalog, and the fiducial source selection cut. The left panel compares the full covariance with the shape noise covariance for the stacked lensing profiles. The sample covariance starts to become the dominant contribution even from $R=2h^{-1}{\rm Mpc}$ for the lowest redshift bin ($0.1\leq z_{\rm cl}\leq 0.4$) due to the high source density of the HSC shear catalog (see also Miyatake et al., 2019). The right panel shows diagonal components of the covariance for the abundance compared to the Poisson term. The sample variance contributions in the abundance for the lower redshift bins are higher than those for the higher redshift bins.

Figure 2 shows the correlation coefficient matrix defined as

[TABLE]

where $D_{i}$ is the $i$-th component of the data vector D. There are large correlations among neighboring bins especially for large radii and lower redshift, since the same large-scale structure causes spatially-correlated fluctuations in the stacked lensing profile and the abundance measurements especially among neighboring radii or different richness bins for each redshift bin (Takada & Bridle, 2007; Oguri & Takada, 2011).

We then calculate the signal-to-noise ratio with the measurements and the covariance for the fiducial analysis with the Planck cosmological parameters. The total signal-to-noise ratio from the stacked lensing profiles and abundance measurements is $61.0$, and the signal-to-noise ratios for the stacked lensing profiles and the abundance are $53.0$ and $30.3$, respectively. The lensing signal-to-noise ratios in each redshift bin are $36.6$, $34.2$, and $19.3$ from the lowest to highest redshift bins.

5 Results

In this section, we show the posterior distribution of the parameters, the joint probability $P_{\beta}(\ln M,\ln N)$ in each redshift bin, the mass-richness relation $P_{\beta}(\ln M|N)$ in each redshift bin, and the richness-mass relation $P(\ln N|M,z)$, from the joint analysis of the model parameters based on the abundance and lensing profiles under the fiducial setups. In this work, we fix the cosmological parameters to the Planck cosmology (Planck Collaboration et al., 2016a) or the WMAP cosmology (Hinshaw et al., 2013) as described in Section 1 to investigate how and to what degree the difference of the cosmological parameters affects the constraints on the mass-richness relation without informative priors on the parameters, and to gain insight into cluster cosmology by comparing the results.

5.1 Posterior distribution of parameters

We constrain the model parameters by comparing the model predictions with the measurements of the abundance and the lensing profiles. We perform Bayesian parameter estimation assuming a Gaussian form for the likelihood, $\mathcal{L}\propto|{\bf C}|^{-1/2}\exp(-\chi^{2}/2)$, with

[TABLE]

where ${\bf D}$ is a data vector that consists of the lensing profiles and the abundance in different radial, richness, and redshift bins, ${\bf D}^{\rm model}$ is the model prediction, and ${\bf C}^{-1}$ is the inverse of the covariance matrix (see Section 4.5). We use 11 radial bins ${\bf R}$ in each richness and redshift bin for the stacked lensing profile. The indices $i$ and $j$ run from $1$ to the total number of data points ($111$ for the fiducial analysis). We perform the parameter estimation with the affine-invariant Markov Chain Monte Carlo (hereafter MCMC) sampler of Goodman & Weare (2010) as implemented in the python package emcee (Foreman-Mackey et al., 2013).

In Table 2, we summarize the results of the parameter estimation, including a description of each parameter, priors, the median and 68$\%$ credible level interval after removing the burn-in chains and marginalizing over the other parameters, and $\chi^{2}_{\rm min}/{\rm dof}$ to show goodness-of-fit under ${\it Planck}$ or ${\it WMAP}$ cosmology. In fitting we employ uninformative flat priors for all of the richness-mass relation parameters and the off-centering parameters. We also show the 68$\%$ and 95$\%$ credible level contours in each two-parameter subspace, and one-dimensional posterior distributions in Appendix B. Each richness-mass relation parameter is well constrained by the joint analysis compared to its prior. From Table 2 and Appendix B, we find that constraints on $A$, $\sigma_{0}$ and $q$ are systematically different between Planck and WMAP cosmologies mainly due to the differences in their halo mass functions. On the other hand, constraints on the other parameters are similar to each other. The result for the WMAP cosmology prefers a higher mean normalization $A$ and lower scatter normalization $\sigma_{0}$ than for the Planck cosmology. In addition, the result for Planck cosmology prefers negative $q$ values (i.e., larger scatter at the lower halo mass) more significantly than the WMAP cosmology, although it is still a $<2\sigma$ preference. For both cosmological models, off-centering parameters are constrained well compared to their priors. The centering fraction at the pivot richness and redshift is constrained as $f_{0}=0.68^{+0.05}_{-0.06}$, which is consistent with $f_{\rm cen}=0.68\pm 0.09$ in the analysis of Oguri et al. (2018a), who estimated the centering fraction without redshift or richness dependences by comparing CAMIRA cluster centers with X-ray centroids. The richness dependence of the off-centering fraction $f_{N}$ prefers positive values with a high significance, indicating that the higher richness clusters are centered better than the lower richness clusters. This might be partly because CAMIRA obtained lower richness values for off-centered clusters given that the richness is computed around the identified BCG using a circular aperture. On the other hand, the redshift dependence parameter $f_{z}$ is consistent with zero. The off-centering radius parameter for each redshift bin $R_{\beta,{\rm off}}$ is constrained compared to its prior, but only marginally.

We check the validity of our model by monitoring the $\chi^{2}$ of the best-fitting models. Since the posterior distribution of the nuisance parameter $m_{\rm lens}$ for lensing amplitudes is strongly determined by its prior, we do not include it as a parameter when calculating the number of degrees-of-freedom, thus ${\rm dof}=111-14=97$ where $111$ is the total number of data points for the fitting and $14$ is the total number of richness-mass relation parameters and off-centering parameters. We find that $\chi^{2}_{\rm min}/{\rm dof}=107.0/97$ ($p$-${\rm value}=0.23$) for the Planck cosmology, and $\chi^{2}_{\rm min}/{\rm dof}=106.6/97$ ($p$-${\rm value}=0.24$) for the WMAP cosmology, both of which are acceptable $p$-values. This indicates that we cannot distinguish between Planck and WMAP cosmologies from the abundance and lensing measurements, partly because of our adoption of a flexible richness-mass relation.

We also show the comparison of the model predictions from the MCMC chains with the measurements of the lensing profiles and abundance in Figures 3 and 4 for Planck and WMAP cosmologies, respectively. The figures show that the model predictions reproduce the lensing profiles and abundance simultaneously for both sets of cosmological parameters with the fiducial richness-mass relation model.

5.2 Joint probability $P_{\beta}(\ln M,\ln N)$ and mass-richness relation $P_{\beta}(\ln M|N)$

We then calculate the joint probability distribution of halo mass and richness in each redshift bin after averaging over the redshift range with volume weight ${\rm d}^{2}V/{\rm d}z{\rm d}\Omega=\chi^{2}(z)/H(z)$ as

[TABLE]

where $P(\ln M|z)$ is the probability distribution of the halo mass for a given redshift, and thus is proportional to the halo mass function ${\rm d}n(M,z)/{\rm d}\ln M$. The normalization factor is determined in the range of $10^{12}\leq M/[h^{-1}M_{\odot}]\leq 2\times 10^{15}$ and $15\leq N\leq 200$, and we restrict the domain of the joint probability to this range. We use the best-fit model parameters of the fiducial analysis in Table 2 for $P(\ln N|M,z)$. The contours in Figure 5 show the joint probability distribution (equation 32) in each redshift bin for the Planck or WMAP cosmological parameters. We obtain the distribution of $\ln N$ and $\ln M$ in the range of richness ($N_{\rm min}=15$, $N_{\rm max}=200$) for each redshift bin by integrating the joint probability distribution along the halo mass and the richness directions respectively as

[TABLE]

and

[TABLE]

The richness distributions $P_{\beta}(\ln N|M_{\rm min}\leq M\leq M_{\rm max})$ shown in Figure 5 indicate that the model reproduces the observed richness function at much finer bins than those used for the analysis. The joint probability contours in Figure 5 show that scatter widths of the mass at a given fixed richness for the WMAP cosmology are smaller than those for the Planck cosmology for all redshift bins, and that the widths for the middle redshift bin ($0.4\leq z_{\rm cl}\leq 0.7$) are the smallest among the redshift bins for both cosmologies. We discuss the origin of these results in Section 6.3.

Since we constrain the richness-mass relation $P(\ln N|M,z)$ from a joint analysis, we can compute the mass-richness relation $P(\ln M|N,z)$ using Bayes theorem as

[TABLE]

where we use the halo mass function for $P(\ln M|z)$. We then average the mass-richness relation over the redshift range with volume weight for each redshift bin as

[TABLE]

We obtain the mean mass for a given richness in each redshift bin as

[TABLE]

Figure 6 shows the median, mean, and 16$\%$ and 84$\%$ percentile region of the mass distribution at fixed richness values in each redshift bin for Planck and WMAP cosmologies, using the best-fit richness-mass relation parameters of the fiducial analysis.

We also show the constraint on the mean relation computed in equation (37) and the scatter relation $\sigma_{\ln M|N}$ in Figure 7 from the MCMC chains after marginalizing over the model parameters (i.e., not only at the best-fit parameters as shown in Figure 6). Since the model generally predicts a skewed distribution of halo mass for a fixed richness value in $\ln M$ space, we define $\sigma_{\ln M|N}$ as the half width of the 68$\%$ percentile region of $P_{\beta}(\ln M|N)$ for each redshift bin as $\sigma_{\ln M|N}=(\ln M_{84}-\ln M_{16})/2$, where $M_{84}$ and $M_{16}$ are masses corresponding to the 84th and 16th percentiles of $P_{\beta}(\ln M|N)$ at a fixed richness, respectively. The mean relations for Planck and WMAP cosmologies are consistent with each other given the error bars for all the redshift bins, and the mean relation of $0.4\leq z_{\rm cl}\leq 0.7$ has a larger amplitude than in the other redshift bins for both cosmologies with relatively high significance given the error bars. We constrain the mean relations at $N=25$ with $\sim 4\%$ precision for $0.1\leq z_{\rm cl}\leq 0.4$ and $0.4\leq z_{\rm cl}\leq 0.7$, and $\sim 8\%$ precision for $0.7\leq z_{\rm cl}\leq 1.0$. The scatter relation for the Planck cosmology increases toward lower richness values for all redshift bins, whereas the scatter for the WMAP cosmology is consistent with a constant value as a function of richness for all redshift bins. The scatter values for the Planck cosmology are systematically larger than those for the WMAP cosmology. This result is qualitatively consistent with the one obtained from cosmological analysis of SDSS redMaPPer clusters in Costanzi et al. (2018), which shows larger scatter values for larger $S_{8}=\sigma_{8}(\Omega_{\rm m}/0.3)^{0.5}$ values ($S_{8}=0.85$ for the Planck and $S_{8}=0.79$ for the WMAP cosmology in this work), although the scatter modeling method is different from ours. The scatter values in the middle redshift bin ($0.4\leq z_{\rm cl}\leq 0.7$) are also lower than those in the other redshift bins. We discuss the origin of this result in Section 6.3.

5.3 Richness-mass relation $P(\ln N|M,z)$

As complementary results to Figure 7 on $P_{\beta}(\ln M|N)$, we show the mean and scatter relations in equations (LABEL:eq:mean_relation) and (17) of the richness-mass relation $P(\ln N|M,z)$ in Figure 8 for typical masses as a function of redshift. We find that, for both Planck and WMAP cosmologies, the mean relation has a minimum around $z\sim 0.5$ for all typical masses. This is also true for the scatter relation especially for WMAP cosmology. We also discuss the origin of these constraints in Section 6.3.

6 Discussion

We discuss the robustness of the fiducial results in Section 6.1, and the complementarity of lensing profile and abundance measurements to constrain the richness-mass relation parameters in Section 6.2. In Section 6.3, we discuss redshift evolution in the richness-mass relation and the difference between the middle redshift bin versus the lower and higher redshift bins shown in Sections 5.2 and 5.3.

6.1 Robustness of our results

6.1.1 Shape and photometric redshift measurement uncertainties

In the fiducial analysis, we marginalize over the shape and photometric redshift uncertainties on the lensing measurements by including the nuisance parameter $m_{\rm lens}$ in equation (29). To check how these uncertainties affect the model parameter constraints, we repeat the MCMC analysis ignoring these errors (i.e., setting $m_{\rm lens}=0$) for the Planck cosmology with the same measurements and covariance as in the fiducial analysis. We find no significant shift in the best-fit parameters with $\chi^{2}_{\rm min}=106.9$ compared to $\chi^{2}_{\rm min}=107.0$ for the fiducial analysis with $m_{\rm lens}$. We also find the 68$\%$ percentile error widths do not change significantly from the fiducial analysis. Specifically, the difference between the error widths in this modified versus the fiducial analysis is smaller than 5$\%$ of the error widths in the fiducial analysis for all parameters. This result shows that the impact of these shape and photometric redshift measurement errors on the joint analysis of current lensing and abundance measurements is negligible when constraining the richness-mass relation.

6.1.2 Different photometric redshift catalogs or different photometric redshift cuts

We use the photometric catalog MLZ and the photometric redshift cuts of Pcut with $\Delta z=0.1$ for the fiducial analysis in Section 5. We check the robustness of our results by using different photometric redshift catalogs (see Section 2.4) or different photometric redshift cuts (see Section 3.2) for the lensing measurements. Here we assume the Planck cosmological parameters for this test, using the same model parameters and priors as shown in Table 2. We repeat the same procedure for the lensing measurements and lensing covariance calculation for different photometric redshift catalogs or cuts to constrain the model parameters by jointly fitting to the lensing and abundance measurements.

In Figure 9, we show the median and the 16th and 84th percentiles of the mean relation $\langle M|N\rangle$ for each redshift bin with respect to the median from the fiducial result of MLZ and the Pcut method with $\Delta z=0.1$. The top panels in Figure 9 show the results with different photometric redshift catalogs, but with the same photometric redshift cuts of the Pcut method with $\Delta z=0.1$. We show that the mean relation $\langle M|N\rangle$ is consistent between the difference photometric redshift catalogs for all redshift bins. The lower panels in Figure 9 show the results with different photometric redshift cuts, but with the same photometric redshift catalog MLZ, showing that the results are consistent with each other given the error bars for all of the redshift bins. This is partly due to our conservative choice of the radial range from $0.5h^{-1}{\rm Mpc}$ in comoving coordinates for the lensing measurements to avoid possible dilution effects on the lensing measurement based on Medezinski et al. (2018b). These results show the robustness of our result to photo-$z$ differences.

6.2 Complementarity of abundance and stacked lensing profile measurements

We constrained the model parameters by jointly fitting to the lensing profiles and abundance measurements in Section 5. Here we study how lensing profile or abundance measurements alone constrain the model parameters, which helps explain how the joint analysis lifts the model parameter degeneracies. We note that we use the same measurements and covariance matrix as in the fiducial analysis.

Figure 10 shows the 68$\%$ and 95$\%$ credible level contours in each two-parameter subspace of the richness-mass relation parameters. The figure shows that the two observables are complementary to each other, which is why their combination can efficiently lift model parameter degeneracies. The lensing measurements constrain the mean normalization parameter $A$ better than the abundance measurements, mainly because the lensing measurements are more sensitive to the mass scale in each redshift and richness bin. The lensing measurements also constrain the scatter normalization parameter $\sigma_{0}$ better than the abundance. On the other hand, other parameters such as the mass dependence parameter $B$ and the redshift evolution parameters $B_{z}$, $C_{z}$ in the mean relation are relatively better constrained by the abundance measurements. In summary, the two observables combine effectively to break complex degeneracies in the model parameters.

6.3 Redshift evolution in the richness-mass relation

In this paper, we use the richness-mass relation model with linear and square redshift evolution parameters ($B_{z}$, $C_{z}$, $q_{z}$, and $p_{z}$) in the fiducial analysis. When we use the richness-mass relation model without any redshift evolution parameters (i.e., without $B_{z}$, $C_{z}$, $q_{z}$, and $p_{z}$) we obtain $\chi^{2}_{\rm min}/{\rm dof}=151.3/101$ ($p$-${\rm value}=8\times 10^{-4}$), and when we adopt the model without square redshift evolution parameters (i.e., without $C_{z}$ and $p_{z}$) we obtain $\chi^{2}_{\rm min}/{\rm dof}=138.6/99$ ($p$-${\rm value}=5\times 10^{-3}$), both of which are unacceptable $p$-values. Here we use the Planck cosmological parameters with the same covariance as the fiducial analysis for these analyses. Given the acceptable $p$-value of $\chi^{2}_{\rm min}/{\rm dof}=107.0/97$ in Table 2, we use the richness-mass model with both linear and square redshift evolution parameters in the fiducial analysis.

With such linear and square redshift evolution parameters, the model prediction can include a non-monotonic dependence on redshift for a fixed mass in $P(\ln N|M,z)$ or a fixed richness in $P_{\beta}(\ln M|N)$. These non-monotonic behaviors in $P(\ln N|M,z)$ and $P_{\beta}(\ln M|N)$ are preferred given the significant improvements of $\chi^{2}_{\rm min}$ by adding the redshift evolution parameters as shown in Figures 5, 6, 7 and 8. More specifically, the mean relation $\langle M|N\rangle$ for the middle redshift bin ($0.4\leq z_{\rm cl}\leq 0.7$) has the \sim$$20\% higher amplitude than those for the other redshift bins for both the Planck and WMAP cosmologies as shown in Figure 7 with relatively high significance, whereas the scatter relation $\sigma_{\ln M|N}$ for the middle redshift bin is slightly smaller than those for the other redshift bins. In addition, Figure 8 shows the non-monotonic behaviors in the mean and scatter relations of $P(\ln N|M,z)$ as a function of redshift.

A possible explanation for the non-monotonic behavior as a function of redshift is different impacts of projection effects at different cluster redshifts. To illustrate this point, in Figure 11 we show a comparison between the photometric cluster redshifts measured from the CAMIRA algorithm and the available spectroscopic redshifts of BCGs without the $4\sigma$ clipping done in Oguri et al. (2018a). The standard deviations of the difference between the photometric and spectroscopic redshifts are smaller around $0.4\leq z_{\rm cl}\leq 0.55$ than for lower and higher cluster redshifts, which is also a non-monotonic behavior in terms of redshift. The larger errors at high redshifts can be understood by larger errors on galaxy magnitudes, whereas the larger errors at low redshifts are most likely due to the lack of $u$-band in the HSC survey as well as too-bright galaxy magnitudes (leading to the saturation in some cases) for such low redshifts.

Optical cluster-finding algorithms in imaging surveys essentially use photometric redshifts of individual galaxies to identify cluster member galaxies. The larger cluster photometric redshift errors imply that photometric redshift errors of individual galaxies are also larger, leading to larger contaminations of non-member galaxies along the line-of-sight direction in estimating the cluster richness. Since the photometric redshift errors in the middle redshift bin ($0.4\leq z_{\rm cl}\leq 0.7$) are smaller than in the other redshift bins, we expect that CAMIRA separates true member galaxies from non-member galaxies along the line-of-sight direction more effectively in the middle redshift bin. In this case, the mean richness values $\langle\ln N(M,z)\rangle$ should be smaller than those in the other redshift bins for a fixed mass even when non-monotonic behaviors in terms of redshift do not exist for intrinsic richness values (i.e., without non-member galaxies along the line-of-sight direction). This is consistent with the result in Figure 8 for the mean relation $\langle\ln N(M,z)\rangle$ for the typical masses. In addition, the non-monotonic behaviors in the order of the mean relation $P_{\beta}(\ln M|N)$ in Figure 7 might also be interpreted by the explanation above. Specifically, the smaller observed richness for a fixed mass leads to the higher mass for a fixed richness, since we expect a smaller contribution of non-member galaxies to observed richness values at the redshift with smaller photometric redshift errors.

The projection effect modifies not only the mean relation but also the scatter of the richness-mass relation. In particular, the larger projection effect implies larger scatter because the projection effect depends sensitively on the projection direction such that a large projection effect is expected when it is projected along the direction of the filamentary structure. This is also consistent with the results in Figures 7 and 8, which show slightly smaller scatter values $\sigma_{\ln M|N}$ in the middle redshift bin than the lower and higher redshift bins, and smaller scatter values $\sigma_{\ln N|M,z}$ around $z_{\rm cl}\sim 0.5$ for typical masses, respectively.

To check the robustness of the fiducial result to its parametrization for the redshift evolution, we repeat the MCMC analysis by using only one of the three redshift bins with a simpler richness-mass relation model without redshift evolution parameters (i.e., only $A$, $B$, $\sigma_{0}$, and $q$ for the richness mass relation), assuming Planck cosmology with the same covariances as the fiducial analysis. This model is similar to the one used in Murata et al. (2018) for SDSS redMaPPer clusters over $0.10\leq z_{\rm cl}\leq 0.33$. Table 3 shows the parameter constraint from this model for each redshift bin. We find that the $p$-values are acceptable for all redshift bins. Figure 12 shows the median and the 16th and 84th percentiles of $\langle M|N\rangle$ and $\sigma_{\ln M|N}$ from the parameter constraints shown in Table 3. The mean and scatter constraints shown in Figure 12 are consistent with the fiducial result shown in Figure 7 within the errors. It is worth noting that the mean relation in the middle redshift bin also favors higher values than for the lower and higher redshift bins, which is similar to the fiducial results. This result supports the non-monotonic redshift evolution of the richness-mass relation found in our fiducial analysis, although the significance is not very high given the larger errors.

7 Conclusion

In this paper, we present the results of the richness-mass relation analysis of 1747 HSC CAMIRA clusters in a wide redshift range ($0.1\leq z_{\rm cl}\leq 1.0$) with a richness range of $N\geq 15$ by jointly fitting to the stacked weak lensing profiles and abundance measurements from the HSC-SSP first-year data ($\sim$232 ${\rm deg^{2}}$ for the cluster catalog and $\sim$140 ${\rm deg^{2}}$ for the shear catalog). The exquisite depth and image quality of the HSC survey enables us to measure stacked weak lensing signals even for high-redshift clusters at $0.7\leq z_{\rm cl}\leq 1.0$ with a total signal-to-noise ratio of $19$.

We constrain the richness-mass relation defined in equations (LABEL:eq:mean_relation) and (17) assuming a log-normal distribution for the relation $P(\ln N|M,z)$ for both the Planck and WMAP cosmological parameters, based on a forward modeling method. We constrain the richness-mass relation parameters without informative priors when marginalizing over off-centering effects on the stacked lensing profiles. We employ an $N$-body simulation-based halo emulator for theoretical predictions of the halo mass function and the lensing profiles. We also use an analytic model for the sample covariance matrix, which is validated against the HSC mock shear and halo catalogs. We find that our model simultaneously fits the stacked lensing profiles and abundance measurements quite well for both the Planck and WMAP cosmological parameters with $\chi^{2}_{\rm min}/{\rm dof}=107.0/97$ and $\chi^{2}_{\rm min}/{\rm dof}=106.6/97$, respectively. We check the robustness of the results against the choice of different photo-$z$ catalogs and source selection cuts. We also show how the stacked lensing and abundance measurements individually constrain the model parameters, and show that the joint analysis efficiently breaks the richness-mass parameter degeneracies.

We then derive the mass-richness relation $P_{\beta}(\ln M|N)$ in each redshift bin, using Bayes theorem from the constraint on $P(\ln M|N,z)$. We show that the mean relations $\langle M|N\rangle$ in each redshift bin are consistent between the Planck and WMAP cosmological parameters within the errors, but the scatter relation values $\sigma_{\ln M|N}$ for the Planck cosmological model are larger than those for the WMAP model. In addition, scatter values for the Planck model increase toward lower richness values, whereas those for the WMAP model are consistent with constant values as a function of richness.

We also show that we need to include the linear and square redshift-dependent parameters in terms of $\ln(1+z)$ for the mean and scatter relations in $P(\ln N|M,z)$ to have acceptable $p$-values. The models without such redshift-dependent parameters resulted in much worse $p$-values. By including the square redshift-dependent parameters, we show that the mean relation $\langle M|N\rangle$ in the middle redshift bin has \sim$$20\% higher amplitude than in the lower and higher redshift bins, whereas the scatter relation $\sigma_{\ln M|N}$ in the middle redshift bin is slightly smaller than in the other bins. We ascribe this non-monotonic redshift dependence to the non-monotonic behavior of the projection effect as a function of redshift, which is supported by the redshift dependence of cluster photometric redshift errors. Redshift evolution in the mean relation $\langle M|N\rangle$ should be properly accounted for when one uses the stacked weak lensing signals around the HSC CAMIRA clusters to validate shear and photo-$z$ catalogs, or to define source selection cuts (e.g., Medezinski et al., 2018b), by matching the cluster weight distributions in terms of cluster redshift and richness values. We also check the consistency of our fiducial results based on the richness-mass relation including the redshift-dependent parameters with those from the analysis of individual redshift bins without redshift-dependent parameters.

Our results indicate that we cannot distinguish between Planck and WMAP cosmological models from the current abundance and lensing profile measurements. This is partly because of our choice of a flexible richness-mass relation model without any informative prior constraints on the model parameters. However, we find that the predicted scatter values are clearly different between Planck and WMAP cosmologies, suggesting that any additional constraints on the scatter of the mass-richness relation may break the degeneracy between these two cosmological models. For instance, one could add other independent probes (e.g., X-ray temperature, the Sunyaev–Zel’dovich effect, or lensing magnification effect) for constraining the scatter values in order to distinguish cosmological models from cluster observables. The analysis of lensing magnification effect might provide complementary information, and will constrain the richness-mass relation of the HSC CAMIRA clusters (Chiu et al. in prep.).

Another possible observable to break the degeneracy is spatial clustering of clusters. Here we briefly investigate the difference of the real-space three-dimensional halo-halo correlation function $\xi_{\rm hh}(r)$ in each redshift and richness bin between the Planck and WMAP models by using the best-fit richness-mass relation parameters from the fiducial analyses in this paper, to roughly assess its power to break the degeneracy. We find that the predicted amplitudes of the halo-halo correlation function for the WMAP model are \sim$$20–$30\%$ larger than those for the Planck model for all richness and redshift bins at $3\lesssim r\lesssim 50h^{-1}{\rm Mpc}$. This implies that clustering of clusters adds useful information that is complementary to abundance and lensing, although careful investigations of cluster photo-$z$ accuracy, redshift-space distortions, and projection effects should be conducted in order to combine our results with the cluster clustering analysis to obtain tight cosmological constraints. We leave this exploration for future work.

In addition, our result should be useful for other cluster-related observational studies (e.g., Lin et al., 2017; Jian et al., 2018; Nishizawa et al., 2018; Miyaoka et al., 2018; Ota et al., 2018; Hashimoto et al., 2019), including galaxy formation in cluster regions and the mass scale estimation for clusters detected via strong-lensing, X-ray, and Sunyaev–Zel’dovich effect with richness values. We can also use our constraint on $P(\ln N|M,z)$ to construct mock CAMIRA cluster catalogs with richness values by using halo mass and redshift in $N$-body simulations, which may be useful for testing the performance of cluster-finding algorithms with simulations (e.g., Dietrich et al., 2014; Oguri et al., 2018a; Costanzi et al., 2019).

Our analysis involves several assumptions. Most critically, we have assumed that the CAMIRA clusters are randomly oriented with respect to the line-of-sight direction in the forward modeling method to compute the cluster observables from the mass function and the spherically averaged halo-matter cross-correlation function from $N$-body simulation outputs. However, this assumption is inaccurate if the CAMIRA clusters are affected by projections effects such as misidentification of non-member galaxies along the line-of-sight direction as member galaxies in the richness estimation (e.g., Cohn et al., 2007; Zu et al., 2017; Busch & White, 2017; Costanzi et al., 2019; Sunayama & More, 2019). Since our cluster selection is based on richness values, CAMIRA could preferentially detect clusters with filamentary structure along the line-of-sight direction. Given the correlation between the halo orientation and surrounding large-scale structure, this effect can change the lensing profile from the spherically-symmetric case (Osato et al., 2018). This investigation is beyond the scope of this paper. Further careful investigations of projection effects for the mass-richness relation are warranted. In order to properly address the projection effects, we need to construct realistic mock catalogs of the CAMIRA clusters with intrinsic richness values and to derive cluster observables accounting for the projection effects. Our results might be useful to check the consistency of such simulation setups with the observations by comparing our measurements and constraints on the mean and scatter relations with ones from such mock catalogs with projection effects.

We have presented a richness-mass relation analysis from HSC first-year data. We will have $\sim$1000 ${\rm deg^{2}}$ area for cluster and shear catalogs when the HSC-SSP survey is complete in 2020. The final HSC cluster analysis has the potential to provide stronger constraints on the richness-mass relation and its better applications for cosmological, galaxy formation, and cluster-related studies particularly when combined with other probes. There is also room for improving the measurement methods e.g., improving the shear measurement technique to include more galaxy shapes, and increasing the sample of galaxies with spectroscopic redshifts independent from the COSMOS 30-band catalog to improve and understand the accuracy of cluster photometric redshifts, and also improving the model framework by further accounting for possible systematic effects such as projection effects.

{ack}

We thank the anonymous referee for helpful comments that improved the quality of this work.

RyM acknowledges financial support from the University of Tokyo-Princeton strategic partnership grant, Advanced Leading Graduate Course for Photon Science (ALPS), Research Fellowships of the Japan Society for the Promotion of Science for Young Scientists (JSPS), and JSPS Overseas Challenge Program for Young Researchers. This work was supported by JSPS KAKENHI Grant Numbers JP15H05892, JP17J00658, JP17K14273, JP18K03693, and JP18H04358, and by Japan Science and Technology Agency CREST JPMHCR1414. This work was supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. Numerical computations were in part carried out on Cray XC30 and XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.

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 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, National Astronomical Observatory of Japan.

Appendix A Covariance

We use analytic calculations of the sampling variance contribution to the covariances, assuming that the distribution of clusters and lensing fields obeys the Gaussian statistics. We describe the analytic model and the detailed estimation procedure in Appendix A.1, and validate it against realistic shear and cluster mock catalogs in Appendix A.2. We also use an analytic model for Poisson shot noise in the abundance covariance. On the other hand, we do not use an analytic model for the shape noise covariance in the lensing profiles, but rather estimate it directly from the data catalogs as described below.

A.1 Analytic model of the covariance matrix

We employ an analytic covariance model for cluster abundances (Hu & Kravtsov, 2003; Takada & Bridle, 2007; Oguri & Takada, 2011) as

[TABLE]

where $\delta_{\alpha\alpha^{\prime}}^{K}$ denotes the Kronecker delta function. The first term denotes the Poisson shot noise from the finite number of available clusters and the second term gives the sample covariance as

[TABLE]

where $\tilde{W}_{s}(\ell\Theta_{s})$ is the Fourier transform of the survey window function, for which we assume a circular survey geometry with survey area $\Omega_{\rm tot}=\pi\Theta_{s}^{2}$ for simplicity: $\tilde{W}_{s}(\ell\Theta_{s})=2J_{1}(\ell\Theta_{s})/(\ell\Theta_{s})$. We use CAMB (Lewis et al., 2000) for computing the linear matter power spectrum $P_{\rm mm}^{\rm L}(k;z)$. The halo weight function is defined as

[TABLE]

Here $b_{\rm h}(M;z)$ is the bias parameter for halos with mass $M$ at redshift $z$, for which we employ a halo bias function presented in Tinker et al. (2010) calculated using the colossus package (Diemer, 2018).

We calculate the covariance model for the stacked lensing profiles among different redshift, richness, and radial bins by decomposing it into the shape noise covariance ${\bf C}^{\rm SN}$ and the sample covariance ${\bf C}^{\rm SV}$ as

[TABLE]

where we compute the shape noise covariance by randomly rotating background galaxies (e.g., Murata et al., 2018). More specifically, we measure the lensing estimators around the clusters in the data catalog with all the multiplicative biases after randomly rotating background shapes in the HSC data catalog, repeating the process 15,000 times. We then calculate the covariance among different richness, redshift, and radial bins based on these measurements.

For the sample covariance of the lensing profiles, we use a Gaussian covariance (Oguri & Takada, 2011; Shirasaki & Takada, 2018) as

[TABLE]

where the power spectrum of convergence fields with $\langle\Sigma_{\rm cr}\rangle_{ls}$ terms are defined as

[TABLE]

and $\langle\chi_{l,\beta}\rangle$ is the average of $\chi(z)$ in $\beta$-th cluster redshift bin from the data with $N\geq 15$. Here $z_{s,\beta}$ is the weighted mean of $z_{s,{\rm best}}$ ($z_{\rm best}$ for a source galaxy, $s$) in $\beta$-th cluster redshift bin over all the radial bins as

[TABLE]

We find $z_{s,1}=1.09$, $z_{s,2}=1.30$, and $z_{s,3}=1.57$ for the fiducial photo-$z$ catalog and source selection cut with the Planck cosmological parameters. The weighted mean critical surface mass density is calculated as

[TABLE]

We use Halofit (Smith et al., 2003) for the nonlinear matter power spectrum based on Takahashi et al. (2012). The terms $C_{\rm hh,\alpha\alpha^{\prime}\beta}^{\rm obs}$ and $C_{\Delta\Sigma,\alpha\beta}$ are defined as

[TABLE]

with

[TABLE]

and

[TABLE]

The second-order Bessel function after averaging within radial bins is given as

[TABLE]

We do not account for the window function effect of the cluster and shear catalogs in this analytic model. In Appendix A.2, we validate this analytical covariance by using realistic mock shear and cluster catalogs. We also ignore the cross-covariance between the stacked lensing profiles and abundance measurements since this cross-covariance does not have a large impact in the parameter estimation, which we confirm by using the cross-covariance estimated from the mock catalogs. Specifically, we repeat the MCMC analysis based on the fiducial covariance with the cross-covariance from the mock catalogs derived in Appendix A.2 to find that the $68\%$ percentile widths are consistent with the fiducial ones and the shift of the $\chi^{2}_{\rm min}$ value from the fiducial value is $\sim 0.1$.

In the parameter estimation, we fix the richness-mass relation parameters for the analytic covariance model to reduce the model calculation time (especially of the lensing covariance). For each setup of the photo-$z$ catalog, source selection cut, and the cosmological parameters, we estimate the analytic covariance model as follows. First, we perform the MCMC analysis with a simpler covariance model which does not include the richness-mass relation dependent terms of $C_{\Delta\Sigma,\alpha\beta}(k)C_{\Delta\Sigma,\alpha^{\prime}\beta^{\prime}}(k)$ in equation (42) and $C_{\rm hh,\alpha\alpha^{\prime}\beta}(k)$ in equation (46). We do not fix the richness-mass relation parameters for other terms in the abundance and lensing profiles. We obtain the best-fit parameters as $\{A,B,B_{z},C_{z},\sigma_{0},q,q_{z},p_{z}\}=\{3.16,0.92,-0.13,4.17,0.29,-0.12,-0.02,0.52\}$ for the Planck model, and $\{A,B,B_{z},C_{z},\sigma_{0},q,q_{z},p_{z}\}=\{3.37,0.84,-0.14,4.47,0.17,-0.02,0.19,0.50\}$ for the WMAP model with the fiducial photo-$z$ catalog (MLZ) and source selection cut (Pcut with $\Delta z=0.1$). Second, we calculate the covariance with the analytic model based on these parameters from the simpler covariance to derive the parameter constraints. For the case of the fiducial photo-$z$ catalog and source selection cut, the richness-mass relation parameters for the covariance calculation above are consistent with our final results shown in Table 2, and the $\chi^{2}_{\rm min}$ values are not very different from the final values shown in Table 2. We show our covariance matrix for the fiducial setup with the Planck cosmological parameters in Figures 1 and 2.

A.2 Validation against realistic mock shear and cluster catalogs

We validate our model of the covariance matrix presented in Appendix A.1 against the realistic HSC shear and halo catalogs (Shirasaki et al., 2019, see Section 4.5 for more details). We use 2268 realizations of the mock catalogs that share the same footprints of the shear and halo catalogs as the real data catalogs. The cosmological parameters for the mock catalogs are the same as those for WMAP used in this paper. We assign richness values for halos with $M\geq 10^{12}h^{-1}M_{\odot}$ to create the mock catalogs of the CAMIRA clusters with the richness values. The richness values are assigned according to the richness-mass relation parameters consistent with our results shown in Table 2 for the WMAP model as $\{A,B,B_{z},C_{z},\sigma_{0},q,q_{z},p_{z}\}=\{3.37,0.85,-0.14,4.47,0.18,-0.05,0.19,0.50\}$. We repeat the measurements of the abundance and lensing profiles for each realization to calculate the covariance matrix from the 2268 realizations. Here we measure the lensing profiles from the shear values without shape noise as we use the shape noise covariance estimated in Appendix A.1. We also calculate the covariance contribution ${\bf C}^{\rm R}$ from the random subtraction (Singh et al., 2017). As shown in Section 2.3 of Singh et al. (2017), we subtract ${\bf C}^{\rm R}$ from the covariance above to account for the random subtraction. Murata et al. (2018) found that this term is negligible (${\bf C}^{\rm R}/{\bf C}\sim 0.01$ for the diagonal terms) for the SDSS redMaPPer clusters. Similarly, we find that these values for the HSC CAMIRA clusters are similar to those for the SDSS redMaPPer clusters and thus are negligible for the HSC CAMIRA clusters.

In Figure 13, we show the comparison between the covariance estimated from mock catalogs and the covariance with the analytic model for the WMAP model using the same richness-mass relation parameters. Since the resolution of the lensing shear in the mock catalogs is limited to $0.43$ arcmin, we compare the lensing covariance only above an effective resolution limit for each cluster redshift bin. Here we set the resolution limit by comparing the mean of the lensing profiles from the mock catalogs with the model prediction (see Section 4.3). The figure shows that the diagonal components of the covariance with the analytic model agree well with those from the mock catalogs at better than the \sim$$10\% level for both the abundance and lensing profiles measurements. We use the covariance with the analytic model for our parameter estimation because the covariance matrix from the mock catalogs is affected by the resolution effect as mentioned above.

Appendix B Model parameter constraint contours

We show the model parameter constraint contours in Figure 14 from the fiducial analysis to show the marginalized one-dimensional posterior distributions for each parameter and the 68$\%$ and 95$\%$ credible levels contours for each two-parameter subspace from the MCMC chains.