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

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


Introduction
The thermal Sunyaev-Zeldovich (tSZ) effect (Sunyaev & Zeldovich 1972) is the spectral distortion of the cosmic microwave background (CMB) through the inverse Compton scattering of CMB photons by hot thermal electrons in galaxy clusters. The tSZ angular power spectrum is sensitive to the amplitude of matter density fluctuations characterized by σ8 and Ωm (Komatsu & Kitayama 1999;Komatsu & Seljak 2002). These parameters are, however, strongly degenerate with the so-called "mass bias" parameter B, which is defined as the ratio of the mass proxy used by the Planck team and the true mass of galaxy clusters, where M500c is the mass enclosed by the radius r500c within which the average mass density is 500 times of the critical density of the Universe. The mass bias B is related to the more commonly used parameter b as B = (1 − b) −1 (Planck Collaboration et al. 2016b). Specifically, the tSZ power spectrum depends primarily on F ≡ σ8(Ωm/B) 0.40 h −0.21 (Bolliet et al. 2018), where Ωm is the matter density parameter, σ8 the r.m.s. matter density fluctuation smoothed over a 8 h −1 Mpc sphere, and h the dimensionless Hubble constant defined by H0 = 100 h km s −1 Mpc −1 .
The mass bias was introduced to account for the cluster mass uncertainty in the Planck analysis. The galaxy cluster masses used by the Planck team were calibrated against a local cluster sample observed by XMM-Newton, assuming the hydrostatic equilibrium (HSE) with thermal pressure (Arnaud et al. 2010). The Planck team reported that 20-40% mass bias (i.e., 1.25 < B < 1.67) is required to reconcile the power spectra of tSZ clusters with the joint result of Planck primordial CMB, CMB lensing and Baryonic Acoustic Oscillations (BAO) (Planck Collaboration et al. 2016b). Several authors performed follow-up analyses of the Planck tSZ power spectrum and found similar results (Horowitz & Seljak 2017;Hurier & Lacasa 2017;Salvati et al. 2018;Bolliet et al. 2018;Makiya et al. 2018;Salvati et al. 2019).
On the other hand, state-of-the-art cosmological hydrodynamic simulations show that the HSE mass underestimates the true mass due to non-thermal pressure support (e.g., Kay et al. 2004;Rasia et al. 2006Rasia et al. , 2012Nagai et al. 2007;Piffaretti & Valdarnini 2008;Lau et al. 2009;Meneghetti et al. 2010;Lau et al. 2013;Nelson et al. 2014bNelson et al. , 2014a. The dominant contribution to the non-thermal pressure support seen in the simulations is the mass accretion of galaxy clusters via structure formation (Shi & Komatsu 2014;Shi et al. 2015), which yields B = 1.28 ± 0.20 (68% C.L.) over a wide range of dynamical states of galaxy clusters (see Table 1 of Shi et al. 2016, for the mass-limited sample and the fitting range of (0.1,1.5)r true 500 ). The predicted HSE mass bias is therefore lower than, or at least on the lower side of, the Planck inferred value. This indicates that other sources of non-thermal pressure (such as cosmic rays and magnetic fields) are larger than expected; there might be other systematic effects such as the calibration error in X-ray observations; and/or that new physics beyond the standard Λ Cold Dark Matter (CDM) model, e.g., dark energy different from the cosmological constant (Bolliet et al. 2018), and/or modified gravity, is required to resolve a tension between cosmological parameters inferred from the tSZ clusters and those from primordial CMB based on ΛCDM (Planck Collaboration et al. 2016b). Massive neutrinos, which also modify the evolution of matter density fluctuations, have been shown not to help resolve this tension (Salvati et al. 2018;Sakr et al. 2018;Bolliet et al. 2019).
To gain more insights into this potential tension, we need an estimate of the mass bias that is independent of the primordial CMB. To this end, in this paper we perform a joint analysis of power spectra of the late-time Universe probes: tSZ and the cosmic weak lensing shear, to obtain a new constraint on the mass bias that is independent of the primordial CMB. Cosmic shear is a unique probe of the growth of total matter distribution including dark matter. The angular power spectrum of the cosmic shear is sensitive to S8 ≡ σ8(Ωm/0.3) α with α ∼ 0.5. This information is useful to break the degeneracy between the cosmological parameters and B.
Here we use the cosmic shear measurements from the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP, hereafter HSC) (Aihara et al. 2018b(Aihara et al. , 2018a. HSC is a wide-field imaging survey using a 1.5 deg diameter field-of-view camera mounted on the prime focus of the 8.2m Subaru telescope . The unique property of HSC is a combination of the depth (i ∼ 26) and the excellent image quality (typical iband seeing is ∼ 0.58"), which enables us to measure cosmic shear signals with unprecedentedly high precision. Hikage et al. (2019) measured tomographic lensing power spectra using the first-year shear catalog (Mandelbaum et al. 2018;Oguri et al. 2018) where the total sky coverage is 137 deg 2 and the effective number density is 17 arcmin −2 in the range of photometric redshifts from 0.3 to 1.5. The first-year shear catalog has been made publicly available as a part of the second public data release of HSC (Aihara et al. 2019).
Constraints on the mass bias from the tSZ and cosmic shear power spectra can be estimated as follows. From the definition of F , the scaling of mass bias with other parameters can be written as Combining the scaling relation of σ8 from the cosmic shear, σ8 = S8(0.3/Ωm) α , we find Similarly we can replace Ωm in Eq.
(2) with σ8. In this case we obtain These quick estimates are closer to B = 1.28 ± 0.20 expected from non-thermal pressure due to mass accretion (Shi et al. 2016), which motivates our calculating B from a joint analysis of the tSZ and cosmic shear power spectra in this paper. This paper is organized as follows. In Section 2 we describe the datasets of tSZ and cosmic shear that we use. In Section 3 we outline the model of the tSZ power spectrum and the details of our likelihood analysis. We discuss our results in Section 4 and conclude in Section 5. Throughout the paper, we adopt a flat ΛCDM cosmology with the minimal total neutrino mass of mν = 0.06 eV.

tSZ
We use the tSZ power spectrum data before marginalizing over the foreground components. Specifically, we take the data from Table 3 of Bolliet et al. (2018), which are based on Planck Collaboration et al. (2016b). The tSZ power spectrum was calculated by cross-correlating the first-half data of the Needlet Internal Linear Combination (NILC) map and last-half of the Modified Internal Linear Combination Algorithm (MILCA) map, where NILC and MILCA are two different methods for reconstructing the tSZ map. 1 We also take into account contributions from residual contaminating sources such as the cosmic infrared background (CIB), IR and radio point sources, and the correlated noise. The amplitudes of CIB, IR and radio point source power spectra are treated as free parameters, while their shapes, which are also taken from Table 3 of Bolliet et al. (2018), are fixed. The amplitude of correlated noise is fixed to reproduce the tSZ spectrum at = 2742, following Bolliet et al. (2018).

Cosmic shear
We use a sample of the posterior distributions of the cosmological parameters of a flat ΛCDM model with a fixed minimal total neutrino mass mν = 0.06 eV obtained from the HSC lensing power spectra alone (Hikage et al. 2019) 2 . The parameters include the five basic cosmological parameters (Ω b h 2 , Ωch 2 , h, As and ns) and nine nuisance parameters regarding modelling errors of point spread functions (PSF), shear biases, photo-z errors, and intrinsic alignments. The reionization optical depth τ is not used in the cosmic-shear-alone analysis. The range of flat priors of cosmological parameters adopted in the HSC cosmic shear analysis are listed in Table 1.

Analysis
We model and analyze the tSZ power spectrum using the halomodel-based approach that we have established in Bolliet et al. (2018) and Makiya et al. (2018), which are based on Komatsu & Seljak (2002). There are two important changes from the methodology established in these papers. First, we properly take into account the effect of massive neutrinos in the modelling of the power spectrum and mass function of collapsed objects (Section 3.1). Second, we improve the method to calculate the likelihood when the covariance matrix includes a non-Gaussian term (Section 3.2).

Massive neutrinos
We include the effect of massive neutrinos by following the socalled "CDM prescription" (Ichiki & Takada 2012;Costanzi et al. 2013;Villaescusa-Navarro et al. 2014;Castorina et al. 2014). The basic idea is to remove the contribution of neutrinos to the mass of collapsed objects (halos) when computing statistics of halos, as neutrino stream out of them. In the calculation of the dark matter halo mass function, we modify the matter density ρm and the matter power spectrum P (k, z) as and where As,ns and k * = 0.05 Mpc −1 are the amplitude, the spectral index and the pivot scale of the primordial power spectrum, Ωc and Ω b are the density parameters of CDM and baryons, ρcrit is the critical density of the Universe, and Tc and T b are the transfer functions of CDM and baryons respectively. We refer the reader to Bolliet et al. (2019) for more details of the effects of massive neutrios on the tSZ power spectrum.

Likelihood analysis
To quickly perform a joint likelihood analysis we use the importance sampling technique (Lewis & Bridle 2002) as follows.
We leave a full likelihood analysis as a future work. First we read a set of cosmological parameters from the likelihood chains of the HSC cosmic shear analysis line by line and calculate the tSZ power spectrum and its likelihood LtSZ for the same set of cosmological parameters. Then we reweigh the chain by multiplying a weight by exp(− ln LtSZ). Applying this procedure to the entire sample of HSC chains, we obtain this paper. Since their mass functions are for the overdensity with respect to the mean mass density (rather than the critical density), we use the spline interpolation of the parameters at various overdensities to obtain the mass function for M500c. As Tinker et al. (2008) only provide the mass function parameters up to the mean mass overdensity of ∆m = 3200, we linearly extrapolate the parameters beyond this bound. We have checked the extrapolation against the fitting function provided by Tinker et al. (2008) and found that the extrapolation method does not have significant effects on our results. the tSZ-shear joint posterior distributions.
Our model has four nuisance parameters: the mass bias B and the amplitudes of the power spectra of CIB, IR and radio point sources, ACIB, AIR and A Rad , respectively. The nuisance parameters are randomly picked from a flat prior with the range listed in Table 1. Following Bolliet et al. (2018), the total power of the contaminating sources are restricted not to exceed the residual of the total tSZ power and the sum of the contributions from resolved clusters. We iterate the above procedure by changing the random seed for nuisance parameters until the chains are converged, and then combine multiple chains to obtain the posterior distributions. We judge that the chains converge when the Gelman-Rubin estimator R − 1, where R is the ratio of the variance between chains and within chains (Gelman & Rubin 1992), is less than 0.05.
The tSZ likelihood is calculated as where ∆ is the difference vector between the observed and the model tSZ spectra, and Cov is the covariance matrix including the non-Gaussian term calculated from the model tSZ trispectrum. The Gaussian term of the covariance matrix is taken from Table 3 of Bolliet et al. (2018).
In Makiya et al. (2018) we calculated the non-Gaussian term at each step of parameter inference. However Carron (2013) pointed out that such a parameter-dependent covariance matrix adds extra artificial information and can bias the parameter constraints. Thus we adopt the new procedure in this paper. First, we perform a likelihood analysis without the non-Gaussian term and find the best-fitting parameters. Then we calculate the non-Gaussian term from the best-fitting parameters and repeat the likelihood analysis by including the fixed non-Gaussian term in the covariance matrix. From the new best-fitting parameters, we recalculate the covariance matrix and redetermine the bestfitting parameters. We iterate this procedure until the best-fitting parameters converge. Figure 1 shows the constraints on the parameters from the Planck tSZ alone, the HSC shear alone, and the joint analysis of them. The results from the joint likelihood analysis are summarized in Table 1. The nuisance parameters not shown in Figure 1 ( ACIB, AIR, A Rad ) are well determined within the prior ranges, as shown in Table 1. The constraints on S8 and F are consistent with those obtained from the HSC cosmic shear and the Planck tSZ power spectrum alone, S8 = 0.780 +0.030 −0.033 (Hikage et al. 2019) and F = 0.460 ± 0.012 (Bolliet et al. 2018), respectively. Since the cosmic shear and tSZ power spectra are not sensitive to individual parameters (i.e., all cosmological parameters and mass bias parameter) but only sensitive to the combination of them, S8 (for the cosmic shear) and F (for the tSZ), the best-fitting model from the joint analysis also gives a good fit to both data sets.

Results
As shown in Figure 1, the mass bias cannot be determined from the tSZ data alone because of the degeneracy with σ8, Ωm and h. The HSC cosmic shear lifts this degeneracy and helps to determine the mass bias. We find that the constraints from the joint analysis do not fully overlap with those from the tSZ alone or the shear alone analysis. This is due to a large difference between the best-fitting values of the tSZ alone or shear alone analysis and those from the joint analysis. The contour shows the likelihood distance from the best-fitting point. When the best-fitting values of two different distributions are largely different, direct comparison of two contours is not meaningful because those contours measure the likelihood distances from different points.
We find that the Planck tSZ and the HSC cosmic shear power spectra constrain the mass bias as B = 1.54 +0.20 −0.35 (mean and 68% C.L.). The black lines in Figure 2 are the predicted scaling relation shown in Eq. (4) and (5) with α = 0.5. We used h = 0.76 for B-Ωm and B-σ8 relations and Ωm = 0.23 for B-h relation, which are the mean value from the joint analysis, respectively. The measurements are consistent with the predicted relations. Figure 2 shows that the Hubble parameter h is not well constrained within the range of prior. While the limited range of the prior of h has an effect on the constraint on B, the effect is not significant since dependence of B on h is weak. We have also examined a Gaussian prior of h = 0.74 ± 0.014 taken from the distance ladder method (Riess et al. 2019) and found that the constraint on B does not improve significantly; B = 1.60 +0.21 −0.33 . This is again due to the weak h dependence of B.
Constraint on the mass bias from our analysis is mainly limited by the uncertainty of σ8. Currently, most (if not all) data sets indicate and are consistent with σ8 < 0.9 (Planck Collaboration et al. 2016a;Alam et al. 2017;Abbott et al. 2018;Burenin 2018). To incorporate this knowledge into our analysis, in Figure 3 we show marginalized joint constraints on B and σ8 with a prior of σ8 < 0.9. This σ8 cut roughly corresponds to Ωm > 0.2 (see the σ8-Ωm panel of Figure 1), which is also consistent with most (if not all) data sets. In this case the mass bias is constrained to be B = 1.37 +0.15 −0.23 or (1−b) = B −1 = 0.73 +0.08 −0.13 (mean and 68% C.L.). This is lower than that from, but is still consistent at 68% C.L. with, the joint analysis of the tSZ and CMB including CMB lensing, B = 1.64 ± 0.19. 4 On the other hand, our result is consistent with the value expected from cosmological hydrodynamical simulations, B = 1.28 ± 0.20 (Shi et al. 2016) and also with that estimated from weak lensing mass, B = 1.25 ± 0.13 (Salvati et al. 2018 and references therein; see also Miyatake et al. 2019;Stern et al. 2019;Dietrich et al. 2019 for recent attempts).

Conclusion
In this paper we have performed a joint likelihood analysis of the tSZ angular power spectrum obtained by Planck and the cosmic shear angular power spectra obtained by Subaru HSC. We have found that the mass bias of Planck clusters is constrained to be B = 1.37 +0.15 −0.23 or (1 − b) = B −1 = 0.73 +0.08 −0.13 (mean and 68% C.L.) for σ8 < 0.9.
Our result is consistent with the HSE mass bias estimated from the cosmological hydrodynamical simulations, which predict B = 1.28 ± 0.20 over a wide range of dynamical states of galaxy clusters (Shi et al. 2016). The origin of this bias is nonthermal motion arising from mass accretion of galaxy clusters via structure formation (Shi & Komatsu 2014;Shi et al. 2015). Therefore, as long as we adopt the cosmological parameters inferred from the HSC cosmic shear with a weak prior of σ8 < 0.9, the origin of the mass bias can be mostly understood. In other words, the cosmological parameters inferred from two different probes of a late-time Universe, tSZ and cosmic shear, are consistent with each other when the mass bias agrees with the expectations from cosmological hydrodynamical simulations. It has been known that the cosmological parameters inferred from the primordial CMB and those from the late-time probes are in a mild tension (e.g., Riess et al. 2018a, Joudaki et al. 2017Köhlinger et al. 2017;Troxel et al. 2018;Burenin 2018;Hikage et al. 2019). Our results suggest that the tSZ power spectrum may also be in tension with the primordial CMB; a higher value of B reported by Planck may come from the tension of σ8 (or S8) between Planck and the late-time Universe probes of the cosmic shear and SZ clusters.
More accurate measurements of cosmic shear is required to obtain a robust conclusion on this issue. In this paper we have used the HSC year 1 data, which are based on only 11% of the planned HSC survey data. Full HSC survey will put tighter constraints on S8 and also improve a constraint on the mass bias. Combining other probes such as galaxy-galaxy lensing and galaxy clustering will improve the constraints.  Mass bias B tSZ+Shear (σ 8 < 0.90) tSZ+CMB Fig. 3. Marginalized constraints on σ8 and the mass bias B from the joint analysis of the Planck tSZ and HSC cosmic shear power spectra (blue), and those from the Planck tSZ, primordial CMB and CMB lensing power spectra (red). The tSZ+shear contours are the same as those in the lower middle panel of Figure 1 except a prior of σ8 < 0.9. The contours show the 68% and 95% confidence levels and are smoothed by a Gaussian of 0.3 times standard deviations in each parameter. The horizontal dashed lines show the 68% confidence region of the mass bias estimated from the hydrodynamical simulations, B = 1.28 ± 0.20 (Shi et al. 2016).
and Astronomy Data Center at National Astronomical Observatory of Japan.