© ESO, 2017

1. Introduction

Galaxy clusters are the largest gravitationally bound structures in the universe. They are consequently a tailored probe of the evolution of the universe, in particular the growth of structure with cosmic time. For example, using numerical simulations (Tinker et al. 2008; Watson et al. 2013), the number counts of galaxy clusters has been shown to scale tightly with cosmological parameters. Using galaxy clusters abundance or correlation functions is now a well-known activity to constrain cosmological parameters (Hasselfield et al. 2013; Planck Collaboration XX 2014; de Haan et al. 2016). Galaxy clusters can be observed through a large number of observational probes: over-density of galaxies (Wen et al. 2012; Rykoff et al. 2014), lensing of background galaxies (Heymans et al. 2012; Erben et al. 2013), X-ray emission from the hot intra-clusters electrons (Böhringer et al. 2001), and the thermal Sunyaev-Zel’dovich (tSZ) effect from the same electron populations (Sunyaev & Zeldovich 1972).

The main limitation to the cosmological use of galaxy clusters is the necessity to properly calibrate the relation between the cluster masses and the probe used to detect them (Planck Collaboration XX 2014). To solve this mass-observable problem, it is possible to combine different probes, so as to simultaneously constrain the cosmological parameters and the mass-observable relations.

In this context, the tSZ effect (Sunyaev & Zeldovich 1972) is a tailored mass proxy. This effect is a distortion of the cosmic microwave background (CMB) blackbody radiation through inverse Compton scattering. CMB photons receive an average energy boost when scattering off hot (a few keV) ionized electrons of the intra-cluster medium (see e.g., Birkinshaw 1999; Carlstrom et al. 2002, for reviews). The intensity of the tSZ effect in a given direction on the sky is measured by the thermal Compton parameter, y, which is related to the electron density along the line of sight by $\mathit{y}\mathrm{\left(}{n}\mathrm{\right)}\mathrm{=}\mathrm{\int }{\mathit{n}}_{\mathrm{e}}\frac{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{e}}}{{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}}\hspace{0.17em}{\mathit{\sigma }}_{\mathrm{T}}\mathrm{d}\mathit{s,}$(1)where ds is the distance along the line of sight, n_e and T_e are respectively the electron number density and temperature. In units of CMB temperature the contribution of the tSZ effect for a given observation frequency ν is $\frac{\mathrm{\Delta }{\mathit{T}}_{\mathrm{CMB}}}{{\mathit{T}}_{\mathrm{CMB}}}\mathrm{=}\mathit{g}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathit{y,}$(2)where, neglecting relativistic corrections, we have the frequency factor $\mathit{g}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\left[\mathit{x}\coth \left(\frac{\mathit{x}}{\mathrm{2}}\right)\mathrm{-}\mathrm{4}\right] \mathrm{with} \mathit{x}\mathrm{=}\frac{\mathit{h\nu }}{{\mathit{k}}_{\mathrm{B}}\hspace{0.17em}{\mathit{T}}_{\mathrm{CMB}}}\mathit{,}$(3)where T_CMB = 2.726 ± 0.001 K, the tSZ effect is negative below 217 GHz and positive for higher frequencies.

Tight cosmological constraints have been obtained from the tSZ signal using cluster number counts (Hasselfield et al. 2013; Reichardt et al. 2013; Planck Collaboration XX 2014; de Haan et al. 2016), the tSZ angular power spectrum (Reichardt et al. 2012; Sievers et al. 2013; George et al. 2015; Planck Collaboration XXII 2016), the tSZ skewness (Wilson et al. 2012), or more recently the tSZ bispectrum (Crawford et al. 2014; Planck Collaboration XXI 2014; Planck Collaboration XXII 2016). The present work aims at combining the constraining power from the cluster number counts, tSZ angular power spectrum and bispectrum, in order to set joint constraints on the cosmological parameters and the mass-observable relation.

The paper is organized as follows, Sect. 2 details the modeling of the tSZ effect angular power spectrum, then Sect. 3 present a new approach to measure the tSZ angular power spectrum with a low level of contamination and the derived cosmological constraints. Then, Sect. 4 presents the modeling of the tSZ bispectrum and Sect. 5 the derived cosmological constraints. Finally, in Sect. 6, we present a combined analysis of all probes: cluster number count, tSZ angular power spectrum and bispectrum, and present the resulting cosmological and mass-observable relation constraints.

2. Power spectra

2.1. General formalism

In recent decades, the statistics of the tSZ signal have been well studied, through analytical developements (see e.g., Cooray 2000; Zhang & Pen 2001; Zentner & Bullock 2003; Cohn & Kadota 2005; Rephaeli & Sadeh 2008; Shaw et al. 2009) and numerical simulations (see e.g., Shaw et al. 2010; Trac et al. 2011; Battaglia et al. 2012; Horowitz & Seljak 2017).

The tSZ angular power spectrum reads ${\mathit{C}}_{\mathit{\ell }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}\sum _{\mathit{m}}\mathrm{|}{\mathit{y}}_{\mathit{\ell m}}{\mathrm{|}}^{\mathrm{2}}\mathit{,}$(4)with y_ℓm the harmonic coefficients of the tSZ map. In the context of large-scale structure tracers, we model these correlations, assuming the following general expression ${\mathit{C}}_{\mathit{\ell }}\mathrm{=}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{1}\mathrm{h}}\mathrm{+}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{2}\mathrm{h}}\mathit{,}$(5)where ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{1}\mathrm{h}}$ is the one-halo contribution and ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{2}\mathrm{h}}$ is the two-halo contribution. These terms can be computed considering a halo model formalism. A key ingredient is the mass function, dn_h/ dM, which gives the abundance of dark matter halos depending on their mass and redshift. In this article we have used the fitting formula from Tinker et al. (2008).

The 1-halo or Poissonian term can be computed using the following ingredients: the Fourier transform of the normalized halo projected profiles in the tSZ map, the mass function, and the tSZ flux of the halo (see e.g., Komatsu & Seljak 2002, for a derivation of the tSZ auto-correlation angular power spectrum); ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{1}\mathrm{h}}\mathrm{=}\mathrm{4}\mathit{\pi }\mathrm{\int }\mathrm{d}\mathit{z}\frac{\mathrm{d}\mathit{V}}{\mathrm{d}\mathit{z}\mathrm{d\Omega }}\mathrm{\int }\mathrm{d}\mathit{M}\frac{\mathrm{d}{\mathit{n}}_{\mathrm{h}}}{\mathrm{d}\mathit{M}}\hspace{0.17em}{\mathit{Y}}_{\mathrm{500}}^{\mathrm{2}}\hspace{0.17em}{\mathit{y}}_{\mathit{\ell }}^{\mathrm{2}}\mathit{,}$(6)where Y₅₀₀ is the spherical tSZ halo flux in the R₅₀₀ radius. This flux depends on M₅₀₀ and z and can be obtained with scaling relations. dV/ dzdΩ is the comoving volume element. The line-of-sight projected Fourier transform of the 3-D profile reads ${\mathit{y}}_{\mathit{\ell }}\mathrm{\left(}{\mathit{M}}_{\mathrm{500}}\mathit{,z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }{\mathit{r}}_{\mathrm{s}}}{{\mathit{l}}_{\mathrm{s}}^{\mathrm{2}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{d}\mathit{x}\hspace{0.17em}{\mathit{x}}^{\mathrm{2}}P\mathrm{\left(}\mathit{x}\mathrm{\right)}\frac{\sin \mathrm{(}\mathit{\ell x}\mathit{/}{\mathit{\ell }}_{\mathrm{s}}\mathrm{)}}{\mathit{\ell x}\mathit{/}{\mathit{\ell }}_{\mathrm{s}}}\mathit{,}$(7)where P(x) is the tSZ halo 3-D profile, x = r/r_s, ℓ_s = D_ang(z) /r_s, and r_s is the scale radius of the halo.

The two-halo term is due to large scale fluctuations of the dark matter field, that induce correlations in the halo distribution over the sky. It can be computed as (see e.g., Komatsu & Kitayama 1999; Diego & Majumdar 2004; Taburet et al. 2011) $\begin{array}{ccc}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{2}\mathrm{h}}& \mathrm{=}& \mathrm{4}\mathit{\pi }\mathrm{\int }\mathrm{d}\mathit{z}\frac{\mathrm{d}\mathit{V}}{\mathrm{d}\mathit{z}\mathrm{d\Omega }}\left(\mathrm{\int }\mathrm{d}\mathit{M}\frac{\mathrm{d}{\mathit{n}}_{\mathrm{h}}}{\mathrm{d}\mathit{M}}{\mathit{Y}}_{\mathrm{500}}{\mathit{y}}_{\mathit{\ell }}\mathit{b}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\right)\mathit{,}\\ & & \mathrm{\times }\left(\mathrm{\int }\mathrm{d}\mathit{M}\frac{\mathrm{d}{\mathit{n}}_{\mathrm{h}}}{\mathrm{d}\mathit{M}}{\mathit{Y}}_{\mathrm{500}}{\mathit{y}}_{\mathit{\ell }}\mathit{b}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\right){\mathit{P}}_{\mathrm{m}}\mathrm{\left(}{\mathit{k}}_{\mathit{\ell }}\mathit{,z}\mathrm{\right)}\end{array}$(8)with k_ℓ = (ℓ + 1/2) /r(z), b(M,z) the linear bias, that relates the matter power spectrum, P_m(k,z), to the power spectrum of the cluster distribution. Following Mo & White (1996), Komatsu & Kitayama (1999) we adopt $\begin{array}{ccc}\mathit{b}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}\mathrm{\left(}{\mathit{\nu }}^{\mathrm{2}}\mathrm{\right(}\mathit{M,z}\mathrm{)}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{/}{\mathit{\delta }}_{\mathrm{c}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{,}& & \end{array}$with ν(M,z) = δ_c(z)/ [D_g(z)σ(M)], D_g(z) is the linear growth factor and δ_c(z) is the over-density threshold for spherical collapse.

Fig. 1 Power density, dlnX/ dlnz, where X stands for either Bℓℓℓ or Cℓ at ℓ = 500 as a function of the redshift for the tSZ power spectrum in dark blue and for the bispectrum in red.

Fig. 2 Power density, dlnX/ dlnz, where X stands for either Bℓℓℓ or Cℓ at ℓ = 500 as a function of the mass, M500, of dark matter halos for the tSZ power spectrum in dark blue and for the bispectrum in red.

In Fig. 1, we present the power distribution as a function of the redshift, showing that the tSZ power spectrum is dominated by objects at z ≤ 1. Figure 2 presents the same power distribution as a function of the galaxy cluster masses. The tSZ power spectrum is dominated by halos with M₅₀₀> 10¹⁴M_⊙, but also receives contribution from smaller halos down to M₅₀₀ = 10¹³M_⊙. The power distributions for the tSZ power spectrum are presented as solid blue lines, solid red lines show the same distributions for the tSZ bispectrum and will be discussed in Sect. 4.

2.2. The tSZ scaling relation

A key step in the modeling of the tSZ effect is to relate the mass, M₅₀₀, and the redshift, z, of a given cluster to its tSZ flux, Y₅₀₀. This relation has to be calibrated on a representative sample of galaxy clusters, with careful propagation of statistical and systematic uncertainties. In this work we have used the M₅₀₀−Y₅₀₀ scaling laws presented in Planck Collaboration XX (2014), ${\mathit{E}}^{\mathrm{-}{\mathit{\beta }}_{\mathrm{sz}}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\left[\frac{{\mathit{D}}_{\mathrm{ang}}^{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}{\mathit{Y}}_{\mathrm{500}}}{{\mathrm{10}}^{-4}\hspace{0.17em}{\mathrm{Mpc}}^{\mathrm{2}}}\right]\mathrm{=}{\mathit{Y}}_{\mathit{\star }}{\left[\frac{\mathit{h}}{\mathrm{0.7}}\right]}^{\mathrm{-}\mathrm{2}\mathrm{+}{\mathit{\alpha }}_{\mathrm{sz}}}{\left[\frac{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{b}\mathrm{)}{\mathit{M}}_{\mathrm{500}}}{\mathrm{6}\mathrm{\times }{\mathrm{10}}^{\mathrm{14}}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}}\right]}^{{\mathit{\alpha }}_{\mathrm{sz}}}\mathit{,}$(9)with E(z) = Ω_m(1 + z)³ + Ω_Λ. The coefficients Y_⋆, α_sz, and β_sz from Planck Collaboration XX (2014), are given in Table 1. We used b = 0.2 for the bias between the X-ray estimated mass and the effective mass of the clusters (hydrostatic mass bias).

2.3. Pressure and density profiles

The tSZ effect is directly proportional to the pressure integrated along the line of sight. In this work, we model the galaxy cluster pressure profile by a generalized Navarro Frenk and White (GNFW, Navarro et al. 1997; Nagai et al. 2007) profile of the form $P\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{0}}}{{\mathrm{(}{\mathit{c}}_{\mathrm{500}}{\mathit{x}}^{\mathrm{)}}}^{\mathit{\gamma }}{\left[\mathrm{1}\mathrm{+}\mathrm{(}{\mathit{c}}_{\mathrm{500}}\mathit{x}{\mathrm{)}}^{\mathit{\alpha }}\right]}^{\mathrm{(}\mathit{\beta }\mathrm{-}\mathit{\gamma }\mathrm{)}\mathit{/}\mathit{\alpha }}}\mathrm{·}$(10)For the parameters c₅₀₀, α, β, and γ, we used the best-fitting values from Arnaud et al. (2010) presented in Table 1. The absolute normalization of the profile P₀ is set assuming the scaling laws Y₅₀₀−M₅₀₀ and presented in Sect. 2.2.

Fig. 3 tSZ power spectra reconstructed by cross-correlating Planck intensity map and MILCA y-map (black sample). From left to right and top to bottom: at 44, 70, 100, and 143 GHz. The best fitting model is presented by the solid red lines. We display the spectra in Compton parameter units, $\begin{array}{c}{}_{\mathrm{􏽢}}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{y\nu }}\end{array}\mathit{/}\mathit{g}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}$, for ease of comparison between cross-spectra.

3. Cosmological constraints from the tSZ power spectrum

In this section we revisit the cosmological constraints derived from the tSZ angular power spectrum by Planck Collaboration XXII (2016). The construction of tSZ maps (e.g., Remazeilles et al. 2011; Hurier et al. 2013; Bobin et al. 2008) is now a well-developed subject, however these maps still often suffer from significant contamination by different foregrounds, in particular by extra-galactic radio-sources and by the cosmic infrared background (CIB; Planck Collaboration XXII 2016). Previous works have handled this sources of contamination using a complete modeling and propagation of the CIB component from the frequency maps to the tSZ power-spectrum (Planck Collaboration XXII 2016; Planck Collaboration XXIII 2016). We present here a new approach to disentangle the tSZ signal from the contribution produced by other astrophysical components. Then, we use these measurements to put constraints on cosmological parameters σ₈ and Ω_b.

3.1. Measurement

We began by computing a tSZ y-map at 7 arcmin resolution using the MILCA procedure (Hurier et al. 2013) from the nine Planck frequencies (30, 44, 70, 100, 143, 217, 353, 545 and 957 GHz) (Planck Collaboration I 2014). By construction the instrumental noise contribution in MILCA maps is correlated with the instrumental noise in Planck intensity maps. In order to cancel this noise contribution, we used the so-called jacknife methodology: given two independent sub-datasets containing half of the Planck data, we computed ${\mathit{C}}_{\mathit{\ell }}^{\mathit{y\nu }}$ as $\begin{array}{ccc}{\mathit{C}}_{\mathit{\ell }}^{\mathit{y\nu }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{C}}_{\mathit{\ell }}^{{\mathit{y}}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}}\mathrm{+}{{\mathit{C}}_{\mathit{\ell }}^{{\mathit{y}}_{\mathrm{2}}{\mathit{\nu }}_{\mathrm{1}}}}^{\mathrm{)}}\mathit{,}& & \end{array}$(11)where y_i is the MILCA tSZ map coming from sub-dataset i, and ν_i is the corresponding Planck intensity map at frequency ν. We specified that the two MILCA maps, have been built using the same linear combination determined on the complete dataset and then applied to the two sub-datasets. We thus obtained nine angular cross power spectra, ${\mathit{C}}_{\mathit{\ell }}^{\mathit{y\nu }}$, each spectra being corrected for the beam and the mask effect following Tristram et al. (2005). In these cross spectra, the tSZ contribution follows the g(ν) spectral energy distribution. It allowed us to separate the tSZ contribution from other residual emissions in the MILCA map (mainly radio sources and CIB contributions).

The variance of ${\mathit{C}}_{\mathit{\ell }}^{\mathit{y\nu }}$ spectra reads, $\begin{array}{ccc}⟨{\mathrm{(}{{\mathit{C}}_{\mathit{\ell }}^{\mathit{y\nu }}}^{\mathrm{)}}}^{\mathrm{2}}⟩& \mathrm{=}& \frac{{\mathrm{(}{{\mathit{C}}_{\mathit{\ell }}^{{\mathit{y}}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}}}^{\mathrm{)}}}^{\mathrm{2}}\mathrm{+}{\mathrm{(}{{\mathit{C}}_{\mathit{\ell }}^{{\mathit{y}}_{\mathrm{2}}{\mathit{\nu }}_{\mathrm{1}}}}^{\mathrm{)}}}^{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathit{C}}_{\mathit{\ell }}^{{\mathit{y}}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}}{\mathit{C}}_{\mathit{\ell }}^{{\mathit{y}}_{\mathrm{2}}{\mathit{\nu }}_{\mathrm{1}}}}{\mathrm{4}\mathrm{(}\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{f}}_{\mathrm{sky}}}\\ & & \end{array}$(12)where, ${\mathit{C}}_{\mathit{\ell }}^{\mathit{\nu \nu }}$ is the auto correlation of the Planck intensity map at frequency ν, ${\mathit{C}}_{\mathit{\ell }}^{\mathit{yy}}$ is the auto-correlation of the tSZ Compton parameter map, and f_sky the covered sky fraction. In the present analysis we used the same mask as used in Planck Collaboration XXII (2016).

The noise is dominated by the second term in Eq. (12), and especially the contribution from astrophysical emissions. Thus, the cross-spectra are highly correlated from frequency to frequency. The covariance can be computed as, $\begin{array}{ccc}⟨{\mathit{C}}_{\mathit{\ell }}^{\mathit{y\nu }}{\mathit{C}}_{\mathit{\ell }}^{\mathit{y}{\mathit{\nu }}^{\mathrm{\prime }}}⟩& \mathrm{=}& \frac{{\mathit{C}}_{\mathit{\ell }}^{{\mathit{y}}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}}{\mathit{C}}_{\mathit{\ell }}^{{\mathit{y}}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}^{\mathrm{\prime }}}\mathrm{+}{\mathit{C}}_{\mathit{\ell }}^{{\mathit{y}}_{\mathrm{2}}{\mathit{\nu }}_{\mathrm{1}}}{\mathit{C}}_{\mathit{\ell }}^{{\mathit{y}}_{\mathrm{2}}{\mathit{\nu }}_{\mathrm{1}}^{\mathrm{\prime }}}\mathrm{+}{\mathit{C}}_{\mathit{\ell }}^{{\mathit{y}}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}^{\mathrm{\prime }}}{\mathit{C}}_{\mathit{\ell }}^{{\mathit{y}}_{\mathrm{2}}{\mathit{\nu }}_{\mathrm{1}}}\mathrm{+}{\mathit{C}}_{\mathit{\ell }}^{{\mathit{y}}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}}{\mathit{C}}_{\mathit{\ell }}^{{\mathit{y}}_{\mathrm{2}}{\mathit{\nu }}_{\mathrm{1}}^{\mathrm{\prime }}}}{\mathrm{4}\mathrm{(}\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{f}}_{\mathrm{sky}}}\\ & & \end{array}$(13)we then performed two additional steps in order to clean the cross-spectra from non-tSZ signal. First, we used the cross-spectra at 217 GHz, ${\mathit{C}}_{\mathit{\ell }}^{\mathit{y,}\mathrm{217}}$, to clean for CMB induced variance¹ in the power spectra estimates, $\begin{array}{ccc}\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{y\nu }}\end{array}\mathrm{=}\frac{\mathit{g}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}}{\mathit{g}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{-}\mathit{g}\mathrm{\left(}\mathrm{217}\mathrm{\right)}}\mathrm{(}{\mathit{C}}_{\mathit{\ell }}^{\mathit{y\nu }}\mathrm{-}{{\mathit{C}}_{\mathit{\ell }}^{\mathit{y,}\mathrm{217}}}^{\mathrm{)}}\mathit{.}& & \end{array}$(14)Additionally, we used the 857 GHz cross-spectra, ${\mathit{C}}_{\mathit{\ell }}^{\mathit{y,}\mathrm{857}}$, where the tSZ power spectrum is expected to be negligible to clean for CIB contribution at first order, $\begin{array}{ccc}\begin{array}{c}\mathrm{􏽢}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{y\nu }}\end{array}\mathrm{=}\frac{\mathit{g}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}}{\mathit{g}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{-}{\mathit{\alpha }}_{\mathit{\nu ,}\mathrm{857}}\mathit{g}\mathrm{\left(}\mathrm{857}\mathrm{\right)}}\mathrm{(}\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{y\nu }}\end{array}\mathrm{-}{\mathit{\alpha }}_{\mathit{\nu ,}\mathrm{857}}\stackrel{}{{{\mathit{C}}_{\mathit{\ell }}^{\mathit{y,}\mathrm{857}}}_{\mathrm{˜}}}\mathrm{)}\mathit{,}& & \end{array}$(15)where, $\begin{array}{ccc}{\mathit{\alpha }}_{\mathit{\nu ,}\mathrm{857}}\mathrm{=}\frac{{\mathrm{⟨}\mathit{\ell }\stackrel{}{{{\mathit{C}}_{\mathit{\ell }}^{\mathit{\nu ,}\mathrm{857}}}_{\mathrm{˜}}}\mathrm{⟩}}_{\mathit{\ell }\mathrm{\in }\mathrm{[}\mathrm{1000}\mathrm{-}\mathrm{2000}\mathrm{]}}}{{\mathrm{⟨}\mathit{\ell }\stackrel{}{{{\mathit{C}}_{\mathit{\ell }}^{\mathrm{857}\mathit{,}\mathrm{857}}}_{\mathrm{˜}}}\mathrm{⟩}}_{\mathit{\ell }\mathrm{\in }\mathrm{[}\mathrm{1000}\mathrm{-}\mathrm{2000}\mathrm{]}}}\mathit{,}& & \end{array}$(16)is an estimation of the infrared contamination SED with $\stackrel{}{{{\mathit{C}}_{\mathit{\ell }}^{\mathit{\nu ,}\mathrm{857}}}_{\mathrm{˜}}}\mathrm{=}{\mathit{C}}_{\mathit{\ell }}^{\mathit{\nu ,}\mathrm{857}}\mathrm{-}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{217}\mathit{,}\mathrm{857}}$.

Such cleaning procedure presents a limitation as described in Hurier et al. (2014), however at low-frequency (<217 GHz), the CIB intensity is faint, and thus will not bias significantly the amplitude of $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{y\nu }}\end{array}$, except at 353 and 545 GHz. These two frequencies have thus been excluded from the analysis. We also excluded the 30 GHz cross-spectra due to its poor angular resolution (≃ 30′) and the contamination from galactic and extra-galactic radio sources.

The two cleaning steps presented above are completely linear, thus we can propagate the ${\mathit{C}}_{\mathit{\ell }}^{\mathit{y\nu }}$ covariance matrix through the linear cleaning operation to determine the $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{y\nu }}\end{array}$ covariance matrix.

In Fig. 3, we present the derived cross-power spectra at 44, 70, 100, and 143 GHz. The global best fitting model on the four spectra is presented as a solid red line. We observe that the spectra from 44 to 143 GHz are consistent with a tSZ spectral behavior, and do not present evidence of a significant contamination from other astrophysical emissions. Such contamination would indeed appear as a frequency dependent bias on the spectra. We stress that, due to the cleaning process and astrophysical sources of noise, the four cross-spectra present highly correlated uncertainties. Thus, the main advantage of having access to this four power spectra is to control and access contamination by other astrophysical sources.

3.2. Cosmological constraints

We fit for cosmological parameters using each of the cross spectra from 44 to 143 GHz. As shown in a previous analysis (Hurier et al. 2014) the shape of the angular power spectrum is essentially sensitive to the Y−M mapping and does not significantly depend on cosmological parameters for the angular scales observed with Planck (ℓ< 2000). Cosmological parameters only affect the overall normalization of the tSZ angular power spectrum. Consequently, there is a degeneracy between parameters σ₈ and Ω_m, and we can only fit for the amplitude S₈ = σ₈ ∗ (Ω_m/ 0.28)^3.2/8.1. Figure 4, shows the likelihood function for the four cross-spectra we used, for two possible priors on the bias of the hydrostatic mass: (i) a Gaussian prior (1−b) = 0.8 ± 0.05, and (ii) a flat prior 0.7 < (1−b) < 0.9. For illustration, these constraints are compared with the cosmological constraints derived from the CMB power spectrum by Planck Collaboration XVI (2014). We find that the uncertainties on S₈ are completely dominated by the uncertainties on the Y−M mapping.

We summarize in Table 2 the constraints derived from each spectra. Combining all spectra we obtain S₈ = 0.77 ± 0.02, consistently with previous tSZ analysis using the same assumptions on the Y−M mapping. In that case the CMB best fitting parameters can be recovered for an hydrostatic mass bias of (1−b) ≃ 0.6 (see, Planck Collaboration XX 2014, for a detailed discussion).

Fig. 4 Cosmological parameters likelihood function derived from MILCA-Planck maps cross-correlation, at 44 (dark blue), 70 (green), 100 (yellow), and 143 GHz (red). Solid lines presents the constraints assuming a Gaussian prior on (1−b) = 0.8 ± 0.05, the dashed lines assume a flat prior 0.7 < (1−b) < 0.9. The solid black line presents the constraints provided by the CMB angular power spectrum.

4. Bispectrum modeling

The angular tSZ equilateral bispectrum is a projection of the corresponding 3D equilateral bispectrum, with Limber’s approximation giving ${\mathit{b}}_{\mathit{\ell \ell \ell }}^{\mathrm{tSZ}}\mathrm{=}\mathrm{\int }\mathrm{d}\mathit{z}\frac{\mathrm{d}\mathit{V}}{\mathrm{d}\mathit{z}\mathrm{d\Omega }}{\mathit{B}}_{\mathrm{tSZ}}\mathrm{\left(}{\mathit{k}}_{\mathit{\ell }}\mathit{,z}\mathrm{\right)}\mathit{,}$(17)with $\mathrm{d}\mathit{V}\mathrm{=}{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\hspace{0.17em}\frac{\mathrm{d}\mathit{r}}{\mathrm{d}\mathit{z}}$ the comoving volume per unit solid angle. The 3D bispectrum is composed of 1-halo, 2-halo and 3-halo terms, with (Lacasa 2014) $\begin{array}{ccc}{\mathit{B}}_{\mathrm{1}\mathrm{h}}\mathrm{\left(}{\mathit{k}}_{\mathit{\ell }}\mathit{,z}\mathrm{\right)}& \mathrm{=}& \\ {\mathit{B}}_{\mathrm{2}\mathrm{h}}\mathrm{\left(}{\mathit{k}}_{\mathit{\ell }}\mathit{,z}\mathrm{\right)}& \mathrm{=}& \mathrm{3}\mathrm{\int }\mathrm{d}{\mathit{M}}_{\mathrm{ab}}{\frac{\mathrm{d}{\mathit{n}}_{\mathrm{h}}}{\mathrm{d}\mathit{M}}|}_{{\mathit{M}}_a}{\frac{\mathrm{d}{\mathit{n}}_{\mathrm{h}}}{\mathrm{d}\mathit{M}}|}_{{\mathit{M}}_b}\hspace{0.17em}({\mathit{Y}}_{\mathrm{500}}\mathrm{\left(}{\mathit{M}}_a\mathit{,z}\mathrm{\right)}{\mathit{y}}_{\mathit{\ell }}\mathrm{\left(}{\mathit{M}}_a\mathit{,z}\mathrm{\right)}{)}^{\mathrm{2}}\\ & & \\ {\mathit{B}}_{\mathrm{3}\mathrm{h}}\mathrm{\left(}{\mathit{k}}_{\mathit{\ell }}\mathit{,z}\mathrm{\right)}& \mathrm{=}& \int \mathrm{d}{\mathit{M}}_{\mathrm{123}}\hspace{0.17em}\left(\underset{\mathit{i}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}{\mathrm{􏽙}}{\frac{\mathrm{d}{\mathit{n}}_{\mathrm{h}}}{\mathrm{d}\mathit{M}}|}_{{\mathit{M}}_{\mathit{i}}}{\mathit{Y}}_{\mathrm{500}}\mathrm{\left(}{\mathit{M}}_{\mathit{i}}\mathit{,z}\mathrm{\right)}{\mathit{y}}_{\mathit{\ell }}\mathrm{\left(}{\mathit{M}}_{\mathit{i}}\mathit{,z}\mathrm{\right)}\right)\\ & & \end{array}$We take the halo power spectrum and bispectrum at the lowest order in bias and perturbation theory (tree-level): $\begin{array}{ccc}{\mathit{P}}_{\mathrm{halo}}\mathrm{\left(}\mathit{k}\mathrm{\right|}{\mathit{M}}_a\mathit{,}{\mathit{M}}_b\mathit{,z}\mathrm{)}& \mathrm{=}& {\mathit{b}}_{\mathrm{1}}\mathrm{\left(}{\mathit{M}}_a\mathit{,z}\mathrm{\right)}\hspace{0.17em}{\mathit{b}}_{\mathrm{1}}\mathrm{\left(}{\mathit{M}}_b\mathit{,z}\mathrm{\right)}{\mathit{P}}_{\mathrm{m}}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}\mathit{,}\\ {\mathit{B}}_{\mathrm{halo}}\mathrm{\left(}\mathit{k}\mathrm{\right|}{\mathit{M}}_{\mathrm{123}}\mathit{,z}\mathrm{)}& \mathrm{=}& {\mathit{B}}_{\mathrm{halo}}^{\mathrm{2}\mathrm{PT}}\mathrm{+}{\mathit{B}}_{\mathrm{halo}}^{\mathrm{b}\mathrm{2}}\mathit{,}\\ {\mathit{B}}_{\mathrm{halo}}^{\mathrm{2}\mathrm{PT}}\mathrm{\left(}\mathit{k}\mathrm{\right|}{\mathit{M}}_{\mathrm{123}}\mathit{,z}\mathrm{)}& \mathrm{=}& \mathrm{6}{\mathit{b}}_{\mathrm{1}}\mathrm{\left(}{\mathit{M}}_{\mathrm{1}}\mathit{,z}\mathrm{\right)}\hspace{0.17em}{\mathit{b}}_{\mathrm{1}}\mathrm{\left(}{\mathit{M}}_{\mathrm{2}}\mathit{,z}\mathrm{\right)}\hspace{0.17em}{\mathit{b}}_{\mathrm{1}}\mathrm{\left(}{\mathit{M}}_{\mathrm{3}}\mathit{,z}\mathrm{\right)}\\ & & \mathrm{\times }{\mathit{F}}_{\mathrm{equi}}{\mathit{P}}_{\mathrm{m}}\mathrm{(}\mathit{k,z}{\mathrm{)}}^{\mathrm{2}} \mathrm{with} {\mathit{F}}_{\mathrm{equi}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{7}}\mathit{,}\\ {\mathit{B}}_{\mathrm{halo}}^{\mathrm{b}\mathrm{2}}\mathrm{\left(}\mathit{k}\mathrm{\right|}{\mathit{M}}_{\mathrm{123}}\mathit{,z}\mathrm{)}& \mathrm{=}& ({\mathit{b}}_{\mathrm{1}}\mathrm{\left(}{\mathit{M}}_{\mathrm{1}}\mathit{,z}\mathrm{\right)}\hspace{0.17em}{\mathit{b}}_{\mathrm{1}}\mathrm{\left(}{\mathit{M}}_{\mathrm{2}}\mathit{,z}\mathrm{\right)}\hspace{0.17em}{\mathit{b}}_{\mathrm{2}}\mathrm{\left(}{\mathit{M}}_{\mathrm{3}}\mathit{,z}\mathrm{\right)}\\ & & \mathrm{+}\mathrm{2}\mathrm{perm}\mathit{.})\mathrm{\times }{\mathit{P}}_{\mathrm{m}}\mathrm{(}\mathit{k,z}{\mathrm{)}}^{\mathrm{2}}\mathit{.}\end{array}$The 1-halo term of the bispectrum, given by the combination of Eqs. (17) and (18), recovers the result from Bhattacharya et al. (2012). The 2-halo and 3-halo terms, derived by Cooray (2000) and Lacasa (2014), are additional contributions due to spatial modulations of the cluster distribution by the dark matter large scale structure. We however found these terms to be subdominant, as the 2-halo power spectrum is subdominant compared to the 1-halo power spectrum, except on the lowest multipoles.

The red lines in Figs. 1 and 2 show the power density of the tSZ bispectrum for equilateral triangle at ℓ = 500 as a function of redshift and mass respectively. We observe that the tSZ bispectrum power density as a function of redshift is similar to the tSZ power spectrum power density, however favoring slightly lower redshift objects. The power density as a function of the mass shows that the tSZ bispectrum is dominated by objects with M₅₀₀ ≃ 10¹⁵M_⊙, and favors more massive halos compared to the tSZ angular power spectrum. These two findings are consistent with our expectations. Indeed the bispectrum, being a higher order quantity, is sensitive to more luminous objects than the power spectrum. Consequently, the bispectrum is dominated by a smaller number of objects and presents a higher sensitivity to the cosmic variance.

5. Cosmological constraints from the tSZ bispectrum

5.1. Measurement

To measure the bispectrum of a y(n) map, we use the following estimator $\begin{array}{ccc}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}\mathit{,}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\mathrm{=}\mathrm{\int }\frac{{\mathrm{d}}^{\mathrm{2}}{n}}{\mathrm{4}\mathit{\pi }}{\mathit{T}}_{{\mathit{\ell }}_{\mathrm{1}}}\mathrm{\left(}{n}\mathrm{\right)}\hspace{0.17em}{\mathit{T}}_{{\mathit{\ell }}_{\mathrm{2}}}\mathrm{\left(}{n}\mathrm{\right)}\hspace{0.17em}{\mathit{T}}_{{\mathit{\ell }}_{\mathrm{3}}}\mathrm{\left(}{n}\mathrm{\right)}\mathit{,}& & \end{array}$(25)where T_ℓ is the so-called scale maps that only contain harmonic coefficients of order ℓ, that is, T_ℓ(n) = ∑ _my_ℓmY_ℓm(n) with y_ℓm the harmonic coefficients of the Compton parameter map.

We refer to Planck Collaboration XXII (2016) for a more detailed description of the bispectrum estimation.

Fig. 5 tSZ bispectrum measured on the MILCA y-map (black sample). The best fitting model is shown as a red solid line.

The uncertainties on the bispectrum are usually estimated under the weak non-Gaussian limit and can be expressed as $\begin{array}{ccc}\mathrm{⟨}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\mathit{,}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}^{\mathrm{\prime }}{\mathit{\ell }}_{\mathrm{2}}^{\mathrm{\prime }}{\mathit{\ell }}_{\mathrm{3}}^{\mathrm{\prime }}}\mathrm{⟩}& \mathrm{=}& \frac{{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{3}}}}{{\mathit{N}}_{{\mathit{\ell }}_{\mathrm{1}}\mathit{,}{\mathit{\ell }}_{\mathrm{2}}\mathit{,}{\mathit{\ell }}_{\mathrm{3}}}}({\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{1}}^{\mathrm{\prime }}}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{2}}^{\mathrm{\prime }}}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{\ell }}_{\mathrm{3}}^{\mathrm{\prime }}}\\ & & \mathrm{+}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{1}}^{\mathrm{\prime }}}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}^{\mathrm{\prime }}}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{\ell }}_{\mathrm{2}}^{\mathrm{\prime }}}\mathrm{+}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}^{\mathrm{\prime }}}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{1}}^{\mathrm{\prime }}}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{\ell }}_{\mathrm{3}}^{\mathrm{\prime }}}\\ & & \mathrm{+}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}^{\mathrm{\prime }}}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}^{\mathrm{\prime }}}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{\ell }}_{\mathrm{1}}^{\mathrm{\prime }}}\mathrm{+}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{3}}^{\mathrm{\prime }}}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{2}}^{\mathrm{\prime }}}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{\ell }}_{\mathrm{1}}^{\mathrm{\prime }}}\\ & & \mathrm{+}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{3}}^{\mathrm{\prime }}}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{1}}^{\mathrm{\prime }}}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{\ell }}_{\mathrm{2}}^{\mathrm{\prime }}})\mathit{,}\end{array}$(26)with N_ℓ1ℓ2ℓ3, being the number of modes for the (ℓ₁,ℓ₂,ℓ₃) triangle.

To estimate the total uncertainty level in the bispectrum, we produced 100 tSZ simulated maps, using a Poissonian sampling of the cluster mass function and putting the corresponding halos randomly in the sky with tSZ fluxes following a log-normal distribution consistent with the tSZ Y−M scaling relation. The resulting total covariance is presented in the upper right panel of Fig. 7.

We note that the MILCA tSZ map is contaminated by non-Gaussian astrophysical component that could in principle bias the measured tSZ bispectrum. First, a contamination by radio sources would appear as a negative contribution in the measured bispectrum at high-ℓ. We avoided this contamination by applying aggressive radio source masks, and indeed find no trace of it in the measured bispectrum. Second, contamination by the cosmic infrared background (CIB) is more delicate, as it cannot be masked and it produces a power spectrum similar to the tSZ one. Nevertheless, the CIB contribution in the y-map is essentially high-z CIB (Planck Collaboration XXIII 2016) and consequently is near-Gaussian, text with a low amplitude bispectrum. Furthermore, the (weak) bispectrum produced by the CIB is significantly less steep than the tSZ bispectrum (e.g., Lacasa et al. 2014), and would thus appear at high multipoles, which we do not see in our measurements. In Fig. 5, we present the measured tSZ bispectrum on the MILCA tSZ-map and the best fitting model for the equilateral configuration of the bispectrum. Uncertainties displayed in Fig. 5 only account for uncertainties induced by the instrumental noise and CIB-leakage. They do not include cosmic variance, that dominates the error budget at low-ℓ.

5.2. Cosmological constraints

Fig. 6 Cosmological parameters likelihood function derived from MILCA y-map bispectrum. The red line shows the likelihood function assuming a Gaussian prior (1−b) = 0.8 ± 0.05 and the blue line a flat prior 0.7 < (1−b) < 0.9. The black line shows the likelihood function derived from CMB angular power spectrum constraints.

Figure 6 presents the cosmological constraints derived from the tSZ bispectrum for two cases: (i) assuming a Gaussian prior of (1−b) = 0.8 ± 0.05 and (ii) a flat prior 0.7 < (1−b) < 0.9 for the bias on hydrostatic mass. We derive ${\mathit{\sigma }}_{\mathrm{8}}{{}^{\mathrm{[}}\left({\mathrm{\Omega }}_{\mathrm{m}}\mathit{/}\mathrm{0.28}\right)^{\mathrm{3.9}}{\left({\mathit{H}}_{\mathrm{0}}\mathit{/}\mathrm{67.1}\right)}^{-1.1}{}^{\mathrm{]}}}^{\mathrm{1}\mathit{/}\mathrm{12.9}}\mathrm{=}\mathrm{0.765}\mathrm{\pm }\mathrm{0.025}$. This constraint is consistent with previous work (Planck Collaboration XXII 2016) and with our derivation of cosmological parameters from the tSZ angular power spectrum in Sect. 3.

6. Combined analysis of number counts, power spectrum, and bispectrum

In this section, we combine our measurement of the tSZ power spectrum (from 44 to 143 GHz) and tSZ bispectrum with cluster number count analysis using the Planck cosmology sample presented in Planck Collaboration XX (2014). We refer to Planck Collaboration XX (2014) for a detailed description of the number count analysis.

6.1. Covariance matrix

Fig. 7 Correlation matrix between galaxy cluster number count as a function of redshift Ncl, tSZ angular power spectrum, Cℓ, and tSZ bispectrum, Bℓ, for the cosmic variance contribution to the uncertainties.

The covariance matrix of galaxy cluster number count, tSZ power spectrum, and tSZ bispectrum is particularly challenging to estimate. A complete analytic derivation would involve the computation of one to six points correlation functions. Consequently, we estimated the covariance between probes by performing simulation of the galaxy cluster mass-function. The tSZ effect is sensitive to very high-mass galaxy cluster, thus we assumed that the mass-function covariance matrix is diagonal with respect to the galaxy cluster masses.

In Fig. 7, we present the correlation matrix derived from the 100 simulations used already in Sect. 5. Similarly to Planck Collaboration XX (2014), we consider ten redshift bins from z = 0 to 1 for galaxy cluster number count, N_cl. For the tSZ power spectrum and bispectrum we consider multipoles from ℓ = 0 to 1000. We verified that we derive similar results when using an analytic computation of the covariance matrix, described in Appendix A. The cosmic variance contribution to the total uncertainties is dominated by non-Gaussian terms that induces significant non-diagonal terms in the tSZ power spectrum and bispectrum covariance matrices. We observe that the tSZ power spectrum is highly correlated with the tSZ bispectrum, with a correlation factor higher than ≃ 0.8 for the presented ℓ range. Consequently, the combination of the tSZ power spectrum and bispectrum will not reduce significantly the cosmic variance contribution to the error on cosmological parameters.

Fig. 8 Likelihood function on Y⋆, Ωm, and σ8 derived from the combined analysis of galaxy cluster number count, tSZ power spectrum, and tSZ bispectrum. Blue contours indicates the 1, 2, and 3σ confidence levels. The red solid lines-sample shows the constraints from CMB power spectrum, the dashed red lines shows the 1σ confidence level. The solid green line shows (1−b) = 0.8 and the dashed green lines show (1−b) = 0.7 and (1−b) = 0.9.

We also observe that the galaxy cluster number count is not significantly correlated with the tSZ power spectrum and bispectrum. The highest level of correlation, ~ 50%, is obtained for low redshift galaxy cluster number count bins and the tSZ power spectrum at high ℓ. The variance in tSZ power spectrum and bispectrum is essentially produced by very massive galaxy clusters (M₅₀₀> 10¹⁵M_⊙), whereas the variance in galaxy cluster number counts is dominated by object with a lower mass, at the edge of the sample selection function. The small overlap between this two populations of galaxy clusters explains the overall small correlation of galaxy cluster number counts with other probes. A more detailed description of the amount of correlation between probes is discussed in Appendix A. At higher ℓ the tSZ power spectrum and bispectrum receives significant contribution from undetected galaxy clusters at higher redshift, inducing an overall small correlation level between number count and angular spectra.

We stress that this correlation matrix only represents the contribution from the cosmic variance to the total uncertainties. The total uncertainties also receive contributions from the instrumental noise and CIB residuals. To account for these additional uncertainties, we added random realizations of them to our 100 tSZ sky simulations.

We estimated the instrumental noise properties using Planck half-dataset difference for each frequency. We assumed that the noise in Planck frequency map is uncorrelated. Then, we propagated the noise through the MILCA linear combination. Using this half-dataset noise map, we estimated the noise inhomogeneities by computing the local standard deviation in a four-degree FWHM Gaussian beam. For the noise in the MILCA map, we also estimated the noise angular power-spectrum that can not be considered spatially uncorrelated due to the component separation process. We performed 100 realistic simulations of the instrumental instrumental noise.

For the CIB component, we also performed 100 homogeneous correlated Gaussian realizations using power spectra from the best fitting model from Planck Collaboration XXIII (2016) and propagate them through the MILCA linear combination.We note that, given that these two sources of uncertainties are Gaussian, they do not add correlations between the power spectrum and bispectrum measurements.

6.2. Cosmological constraints

In Fig. 8, we present the constraints we obtained on Y_⋆, σ₈, and Ω_m. In this case we have not added priors on the Y−M normalisation. And we interpreted the Y−M calibration modification as an adjustement of (1−b). We observe that the combined analysis favors a lower calibration for the Y−M relation, leading to a best fitting value of (1−b) = 0.71 ± 0.07. We also obtain σ₈ = 0.79 ± 0.02 and Ω_m = 0.29 ± 0.02.

We checked that individual results derived from each probe are consistent together as well as consistent with the total combination. In order to investigate where the information is coming from, we examined the constraints from the three possible combination of two probes. We find that the combination of galaxy cluster number counts and power spectra allows us to break degeneracies between the three considered parameters, giving optimal results on σ₈, σ₈ = 0.79 ± 0.02, but information is still missing for the other parameters, Ω_m = 0.29 ± 0.04 and (1−b) = 0.74 ± 0.14. The combination of the tSZ angular power-spectrum and bispectrum is efficient to break the degeneracy between Ω_m and ((1−b),σ₈), though still having a degeneracy between (1−b) and σ₈); it gives Ω_m = 0.28 ± 0.03 and σ₈((1−b)/0.7)^-0.42 = 0.81 ± 0.02. Finally the combination of number counts and the tSZ bispectrum gives Ω_m = 0.29 ± 0.03, σ₈ = 0.79 ± 0.03, and (1−b) = 0.70 ± 0.07. This last combination is thus particularly efficient to determine (1−b) and drives our constraint on the hydrostatic mass bias when combining the three probes.

7. Conclusion and discussion

We have revisited the cosmological constraints derived from the tSZ angular power spectrum by performing a combined analysis of a tSZ y-map and Planck intensity map per frequency. This approach provided us with a measurement of the tSZ angular power spectrum robust with respect to the CIB contamination in the Planck tSZ y-maps, which was previously an important issue. From this analysis, we derived robust cosmological constraints, which come out consistent with previous works from the Planck collaboration.

We presented a halo-model description of the tSZ bispectrum and compared it with a measurement of the equilateral bispectrum of a Planck-derived tSZ y-map. Using the measurement we were able to set cosmological constraints, also robust with respect to contamination by near-Gaussian signals such as the CIB. These constraints are found to be both competitive and consistent with those from the tSZ angular power spectrum.

By computing their joint covariance, we demonstrated that the tSZ power spectrum and bispectrum present a high degree of correlation for the cosmic variance contribution to the uncertainties. However, the number counts of galaxy clusters, as performed in the analysis by Planck Collaboration XX (2014), are not significantly correlated with the tSZ power spectrum nor bispectrum. Combining the number count analysis with our measurement of the tSZ angular power spectrum and bispectrum, we have been able to set tight constraints on the hydrostatic mass bias and cosmological parameter simultaneously.

These results favor a value for the hydrostatic mass bias (1−b) = 0.71 ± 0.07, consistent with the prior used in (Planck Collaboration XX 2014), (1−b) ∈ [0.7,1.0]. It is particularly interesting to note that our combined analysis enables us to break the degeneracy between cosmological parameters and the normalization of the scaling relation. It is the combination of the number counts and the bispectrum that drives the constraint on (1−b), with the power spectrum helping to further reduce the cosmological errors. Finally, comparing these results with cosmological parameters derived from the Planck CMB analysis, we obtain an agreement between the two probes at 1σ level.

Acknowledgments

We acknowledge the support of the Agence Nationale de la Recherche through grant ANR-11-BS56-015.

References

Appendix A: Correlation matrix between galaxy cluster number count, tSZ power spectrum, and tSZ bispectrum

We compared the correlation matrix presented in Fig. 7, that has been derived though simulations, with an analytical derivation. We compute the correlation, ρ, matrix as follows: $\begin{array}{ccc}\mathit{\rho }\mathrm{\left(}\mathit{X,Y}\mathrm{\right)}& \mathrm{=}& \frac{\mathit{I}\mathrm{\left(}\mathit{X,Y}\mathrm{\right)}}{\sqrt{\mathit{I}\mathrm{\left(}\mathit{X,X}\mathrm{\right)}\mathit{I}\mathrm{\left(}\mathit{Y,Y}\mathrm{\right)}}}\\ \mathit{I}\mathrm{\left(}\mathit{X,Y}\mathrm{\right)}& \mathrm{=}& \end{array}$(A.1)where X and Y stand either for number counts redshift bins, power spectrum multipole bins, or bispectrum multipole bins, W^X and W^Y are the window functions for the X and Y quantities, and $\frac{\mathrm{d}{\mathit{n}}_{\mathrm{eff}}}{\mathrm{d}\mathit{M}}\mathrm{=}\frac{\mathrm{d}{\mathit{n}}_{\mathrm{h}}}{\mathrm{d}\mathit{M}}$ if the number count is not involved and $\frac{\mathrm{d}{\mathit{n}}_{\mathrm{eff}}}{\mathrm{d}\mathit{M}}\mathrm{=}\frac{\mathrm{d}{\mathit{n}}_{\mathrm{h}}}{\mathrm{d}\mathit{M}}{\mathit{F}}_{\mathrm{sel}}$ if the number counts is involved, with F_sel is the cosmological sample selection function. We computed these window functions as, $\begin{array}{ccc}{\mathit{W}}^{\mathrm{PS}}& \mathrm{=}& {\mathit{Y}}_{\mathrm{500}}^{\mathrm{2}}{\mathit{y}}_{\mathit{\ell }}^{\mathrm{2}}\\ {\mathit{W}}^{\mathrm{BS}}& \mathrm{=}& {\mathit{Y}}_{\mathrm{500}}^{\mathrm{3}}{\mathit{y}}_{\mathit{\ell }}^{\mathrm{3}}\\ {\mathit{W}}^{\mathrm{NC}}& \mathrm{=}& \end{array}$(A.2)with the superscripts PS, BS, and NC corresponding to the power spectrum, bispectrum, and number counts respectively, and w(z_m) is equal to 1 for z_m−0.05 <z<z_m + 0.05, and 0 otherwise. We note that that this derivation is an approximation for the number count analysis, as it does not account for M−Y relation intrinsic scatter and redshift uncertainties. We present the derived covariance matrix in Fig. A.1, finding an excellent agreement with the covariance matrix derived through simulations and shown in Fig. 7.

To illustrate the origin of this low level of correlation between galaxy cluster number counts and tSZ power spectrum and bispectrum, we present in Figs. A.2−A.5 a few examples of the covariance density in the M₅₀₀−z plane, d²ln(Cov)/dlnMdlnz. We observe that the variance in the power spectrum and bispectrum is produced by very massive halos (M₅₀₀> 10¹⁵M_⊙) at low redshift. The variance in number count is produced by less massive halos (M₅₀₀< 10¹⁵M_⊙). This effect can be explained through the weighting applied to the halo mass function (${\mathit{Y}}_{\mathrm{500}}^{\mathrm{2}}$ for the power spectrum, ${\mathit{Y}}_{\mathrm{500}}^{\mathrm{3}}$ for the bispectrum, F_sel for the number counts). This results illustrates that the tSZ-based sample already contains enough low-mass galaxy clusters to be in a regime where the correlation between number count and spectra is negligible.

Fig. A.1 As Fig. 7, but derived through analytical considerations. See Eqs. (A.1) and (A.2) for more details.

Fig. A.2 Covariance density, $\frac{{\mathrm{d}}^{\mathrm{2}}\ln \mathrm{\left(}\mathrm{Cov}\mathrm{\right)}}{\mathrm{d}\ln \mathrm{\left(}{\mathit{M}}_{\mathrm{500}}\mathrm{\right)}\mathrm{d}\mathit{ln}\mathrm{\left(}\mathit{z}\mathrm{\right)}}$, as a function of halo mass, M500, and redshift, z. Number count is considered for z< 0.1, power spectrum and bispectrum are at at ℓ = 250. The plots are organized as for Fisher forecasts. On the diagonal are covariance density for the auto-probe variance, with from left to right: contribution to the variance of number counts, power spectrum and bispectrum. Off-diagonal plots represent cross-probes covariance density, with from left to right and top to bottom, covariance between: number counts and spectrum, number counts and bispectrum, and between power spectrum and bispectrum. The color scale is consistent between panels and represents the total amount of correlation between probes.

Fig. A.3 As Fig. A.2, with the power spectrum and bispectrum at ℓ = 1000.

Fig. A.4 As Fig. A.2, with the number counts at 0.1 <z< 0.2.

Fig. A.5 As Fig. A.2, with the number counts at 0.1 <z< 0.2 and the power spectrum and bispectrum at ℓ = 1000.

All Tables

All Figures

	Fig. 1 Power density, dlnX/ dlnz, where X stands for either Bℓℓℓ or Cℓ at ℓ = 500 as a function of the redshift for the tSZ power spectrum in dark blue and for the bispectrum in red.
In the text

	Fig. 2 Power density, dlnX/ dlnz, where X stands for either Bℓℓℓ or Cℓ at ℓ = 500 as a function of the mass, M500, of dark matter halos for the tSZ power spectrum in dark blue and for the bispectrum in red.
In the text

Fig. 3 tSZ power spectra reconstructed by cross-correlating Planck intensity map and MILCA y-map (black sample). From left to right and top to bottom: at 44, 70, 100, and 143 GHz. The best fitting model is presented by the solid red lines. We display the spectra in Compton parameter units, $\begin{array}{c}{}_{\mathrm{􏽢}}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{y\nu }}\end{array}\mathit{/}\mathit{g}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}$, for ease of comparison between cross-spectra.

In the text

	Fig. 4 Cosmological parameters likelihood function derived from MILCA-Planck maps cross-correlation, at 44 (dark blue), 70 (green), 100 (yellow), and 143 GHz (red). Solid lines presents the constraints assuming a Gaussian prior on (1−b) = 0.8 ± 0.05, the dashed lines assume a flat prior 0.7 < (1−b) < 0.9. The solid black line presents the constraints provided by the CMB angular power spectrum.
In the text

	Fig. 5 tSZ bispectrum measured on the MILCA y-map (black sample). The best fitting model is shown as a red solid line.
In the text

	Fig. 6 Cosmological parameters likelihood function derived from MILCA y-map bispectrum. The red line shows the likelihood function assuming a Gaussian prior (1−b) = 0.8 ± 0.05 and the blue line a flat prior 0.7 < (1−b) < 0.9. The black line shows the likelihood function derived from CMB angular power spectrum constraints.
In the text

	Fig. 7 Correlation matrix between galaxy cluster number count as a function of redshift Ncl, tSZ angular power spectrum, Cℓ, and tSZ bispectrum, Bℓ, for the cosmic variance contribution to the uncertainties.
In the text

Fig. 8 Likelihood function on Y⋆, Ωm, and σ8 derived from the combined analysis of galaxy cluster number count, tSZ power spectrum, and tSZ bispectrum. Blue contours indicates the 1, 2, and 3σ confidence levels. The red solid lines-sample shows the constraints from CMB power spectrum, the dashed red lines shows the 1σ confidence level. The solid green line shows (1−b) = 0.8 and the dashed green lines show (1−b) = 0.7 and (1−b) = 0.9.

In the text

	Fig. A.1 As Fig. 7, but derived through analytical considerations. See Eqs. (A.1) and (A.2) for more details.
In the text

Fig. A.2 Covariance density, $\frac{{\mathrm{d}}^{\mathrm{2}}\ln \mathrm{\left(}\mathrm{Cov}\mathrm{\right)}}{\mathrm{d}\ln \mathrm{\left(}{\mathit{M}}_{\mathrm{500}}\mathrm{\right)}\mathrm{d}\mathit{ln}\mathrm{\left(}\mathit{z}\mathrm{\right)}}$, as a function of halo mass, M500, and redshift, z. Number count is considered for z< 0.1, power spectrum and bispectrum are at at ℓ = 250. The plots are organized as for Fisher forecasts. On the diagonal are covariance density for the auto-probe variance, with from left to right: contribution to the variance of number counts, power spectrum and bispectrum. Off-diagonal plots represent cross-probes covariance density, with from left to right and top to bottom, covariance between: number counts and spectrum, number counts and bispectrum, and between power spectrum and bispectrum. The color scale is consistent between panels and represents the total amount of correlation between probes.

In the text

	Fig. A.3 As Fig. A.2, with the power spectrum and bispectrum at ℓ = 1000.
In the text

	Fig. A.4 As Fig. A.2, with the number counts at 0.1 <z< 0.2.
In the text

	Fig. A.5 As Fig. A.2, with the number counts at 0.1 <z< 0.2 and the power spectrum and bispectrum at ℓ = 1000.
In the text