# Breaking degeneracies with the Sunyaev-Zeldovich full bispectrum

Andrea Ravenni, Matteo Rizzato, Slađana Radinović, Michele Liguori, Fabien Lacasa, Elena Sellentin

► To cite this version:

Andrea Ravenni, Matteo Rizzato, Slađana Radinović, Michele Liguori, Fabien Lacasa, et al. Breaking degeneracies with the Sunyaev-Zeldovich full bispectrum. Journal of Cosmology and Astroparticle Physics, 2021, 06, pp.026. 10.1088/1475-7516/2021/06/026. [hal-02946238](https://hal.science/hal-02946238v1)

HAL Id: hal-02946238

https://hal.science/hal-02946238v1

Submitted on 1 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Breaking degeneracies with the Sunyaev-Zeldovich full bispectrum

To cite this article: Andrea Ravenni et al JCAP06(2021)026

View the article online for updates and enhancements.

## You may also like

- Assessing the Importance of Noise from Thermal Sunyaev-Zeldovich Signals for CMB Cluster Surveys and Cluster Cosmology
Srinivasan Raghunathan

- DETECTION OF QUASAR FEEDBACK FROM THE THERMAL SUNYAEV-ZEL'DOVICH EFFECT IN PLANCK
John J. Ruan, Matthew McQuinn and Scott F. Anderson

- IMPROVED CONSTRAINTS ON COSMIC MICROWAVE BACKGROUND SECONDARY ANISOTROPIES FROM THE COMPLETE 2008 SOUTH POLE TELESCOPE DATA
E. Shirokoff, C. L. Reichardt, L. Shaw et al.
# Breaking degeneracies with the Sunyaev-Zeldovich full bispectrum

Andrea Ravenni,$^{a,*}$ Matteo Rizzato,$^{b,c,d}$ Slađana Radinović,$^e$ Michele Liguori,$^{f,g,h}$ Fabien Lacasa$^i$ and Elena Sellentin$^d$

Received September 2, 2020

Revised March 21, 2021

Accepted May 7, 2021

Published June 14, 2021

**Abstract.** Non-Gaussian (NG) statistics of the thermal Sunyaev-Zeldovich (tSZ) effect carry significant information which is not contained in the power spectrum. Here, we perform a joint Fisher analysis of the tSZ power spectrum and bispectrum to verify how much the full bispectrum can contribute to improve parameter constraints. We go beyond similar studies of this kind in several respects: first of all, we include the complete power spectrum and bispectrum (auto- and cross-) covariance in the analysis, computing all NG contributions; furthermore we consider a multi-component foreground scenario and model the effects of component separation in the forecasts; finally, we consider an extended set of both cosmological and intra-cluster medium parameters. We show that the tSZ bispectrum is very efficient

<footnote>$^a$Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, The University of Manchester, Manchester M13 9PL, U.K.</footnote>

<footnote>$^b$Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France</footnote>

<footnote>$^c$Sorbonne Universités, Institut Lagrange de Paris (ILP), 98 bis bd Arago, 75014 Paris, France</footnote>

<footnote>$^d$Leiden Observatory, Leiden University, Huygens Laboratory, Niels Bohrweg 2, NL-2333 CA Leiden, The Netherlands</footnote>

<footnote>$^e$Institute of Theoretical Astrophysics, University of Oslo, 0315 Oslo, Norway</footnote>

<footnote>$^f$Dipartimento di Fisica e Astronomia “G. Galilei”, Università degli Studi di Padova, via Marzolo 8, I-35131 Padova, Italy</footnote>

<footnote>$^g$INFN, Sezione di Padova, via Marzolo 8, I-35131 Padova, Italy</footnote>

<footnote>$^h$INAF-Osservatorio Astronomico di Padova, vicolo dell'Osservatorio 5, I-35122 Padova, Italy</footnote>

<footnote>$^i$Université Paris-Saclay, CNRS, Institut d'astrophysique spatiale, 91405 Orsay, France</footnote>

<footnote>E-mail: andrea.ravenni@manchester.ac.uk, rizzato@strw.leidenuniv.nl, sladana.radinovic@astro.uio.no, liguori@pd.infn.it, fabien.lacasa@universite-paris-saclay.fr, sellentin@strw.leidenuniv.nl</footnote>

<footnote>*Corresponding author.</footnote>
at breaking parameter degeneracies, making it able to produce even stronger cosmological
constraints than the tSZ power spectrum: e.g. the standard deviation on σ<sub>8</sub> shrinks from
σ<sup>PS</sup>(σ<sub>8</sub>) = 0.35 to σ<sup>BS</sup>(σ<sub>8</sub>) = 0.065 when we consider a multi-parameter analysis. We find
that this is mostly due to the different response of separate triangle types (e.g. equilateral and
squeezed) to changes in model parameters. While weak, this shape dependence is clearly non-
negligible for cosmological parameters, and it is even stronger, as expected, for intra-cluster
medium parameters.

**Keywords:** Sunyaev-Zeldovich effect, cosmological parameters from LSS

ArXiv ePrint: 2008.12947
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Theoretical tSZ correlation functions</td><td>2</td></tr><tr><td>2.1</td><td>Halo mass function</td><td>4</td></tr><tr><td>2.2</td><td>Electron pressure profile</td><td>5</td></tr><tr><td>3</td><td>Covariance matrix for the (binned) observables</td><td>5</td></tr><tr><td>3.1</td><td>Structure of the covariance matrix</td><td>6</td></tr><tr><td>3.2</td><td>Power spectrum covariance matrix</td><td>7</td></tr><tr><td>3.3</td><td>Bispectrum covariance matrix</td><td>7</td></tr><tr><td>4</td><td>Power spectrum and bispectrum dependence on the parameters</td><td>9</td></tr><tr><td>5</td><td>Gaussian experimental noise and foregrounds</td><td>11</td></tr><tr><td>5.1</td><td>Spectral components</td><td>11</td></tr><tr><td>5.2</td><td>Component separation</td><td>13</td></tr><tr><td>6</td><td>Forecast results</td><td>13</td></tr><tr><td>6.1</td><td>Noiseless survey with no foreground contamination</td><td>14</td></tr><tr><td>6.1.1</td><td>Analysis of a subset of configurations</td><td>19</td></tr><tr><td>6.2</td><td>Validation against Planck noise and foreground contamination</td><td>19</td></tr><tr><td>6.3</td><td>Forecast for future realistic surveys</td><td>20</td></tr><tr><td>7</td><td>Discussion and conclusions</td><td>22</td></tr><tr><td>A</td><td>Projection of the 3D tSZ field</td><td>24</td></tr><tr><td>A.1</td><td>Correlation functions for the projected Compton-y field</td><td>24</td></tr><tr><td>A.2</td><td>Data vector</td><td>26</td></tr><tr><td>A.3</td><td>Covariance matrix</td><td>27</td></tr><tr><td>B</td><td>Three-halo term</td><td>27</td></tr><tr><td>C</td><td>Noise and foreground spectral shape</td><td>28</td></tr><tr><td>D</td><td>Full triangle plots</td><td>29</td></tr><tr><td>E</td><td>Binning</td><td>30</td></tr><tr><td>F</td><td>Bispectrum derivatives and breaking the degeneracies</td><td>34</td></tr></table>
# 1 Introduction

The thermal Sunyaev-Zeldovich (tSZ) effect [1–3] [4, for a recent review] is a spectral distortion of the Cosmic Microwave Background (CMB), mostly generated in galaxy clusters by inverse Compton scattering of CMB photons off hot electrons. The tSZ effect is a powerful cosmological observable, mainly applied to the study of individual clusters, to build cluster catalogues and to extract number count statistics. A complementary possibility consists in the study of the tSZ angular power spectrum [5]. This approach has been adopted as a powerful probe of the low-redshift Universe, to test both the standard $\Lambda$CDM scenario [6, 7] and some extended models, which encompass primordial non-Gaussianity, massive neutrinos and dark energy [8–12]. One of the advantages of the tSZ power spectrum analysis is that it allows including also small, unresolved clusters, and it does not require direct measurements of cluster masses.

As it has been argued long before the tSZ was routinely measured across the sky [13], it is important to notice that the tSZ map is highly non-Gaussian, therefore only a part of the available tSZ information is actually captured by the power spectrum. A natural question which arises is therefore how much additional information can be extracted from higher order statistics, starting with the bispectrum (i.e. the 3-point multipole correlation function).

An initial theoretical study of the tSZ bispectrum was performed in [14]. There, it was pointed out that, besides having a different amplitude scaling with $\sigma_8$, compared to the power spectrum, the bispectrum signal also takes its main contributions from massive clusters at low redshift. This makes the impact of astrophysical uncertainties smaller in the bispectrum than in the power spectrum, since the latter has a sizeable contribution from less-well-understood, low-mass, high-redshift clusters. At the same time the tSZ skewness was detected using Atacama Cosmology Telescope (ACT) data [15, 16], while both the tSZ skewness and bispectrum were measured by the South Pole Telescope (SPT) [17]. An all-sky Compton-$y$ map, and subsequent measurements of both the skewness and the bispectrum were later obtained by Planck [6, 7]. In all cases it was shown that the bispectrum could be used to obtain significant constraints on $\sigma_8$. A later study [18] combined cluster counts with power spectrum and equilateral (all three $\ell$'s equal) bispectrum measurements, showing that the bispectrum — even just by considering the small equilateral triangle subset — can play a significant role in breaking (cosmological and astrophysical) parameter degeneracies. If this is the case, as we will show, the tSZ bispectrum would be an extremely valuable source of information, especially in light of how difficult it is proving to understand the astrophysics of halos and reconcile various measurement of gas parameters, particularly the hydrostatic mass bias $b_{\text{HSM}}$ [19–24].

These and other similar results clearly encourage further investigation. The aim of this work is to explore how much extra-information is contained in the 3-point correlation function, via a detailed joint Fisher analysis of the tSZ power spectrum and bispectrum. Our main goal is to make our forecasts as realistic and accurate as possible. For this purpose:

* We consider the entire bispectrum domain (not just specific triangles) and compute for the first time the full tSZ bispectrum covariance, beyond the Gaussian approximation. We therefore include all contributions in the bispectrum covariance, up to the connected 6-point correlator, and we will also account for correlations between the power spectrum and the bispectrum.

* We extend our parameter space with respect to previous forecasts and analyses, in order to model in greater detail the impact of uncertainties in the electron pressure profile.
* We account for foreground contamination and model (in a simple way) the effects of component separation, using Internal Linear Combination (ILC).

As we show in the following, our analysis reinforces the conclusion that the tSZ bispectrum is a powerful observable to constrain cosmology, especially due to its efficiency in breaking degeneracies between astrophysical and cosmological parameters, otherwise present in a power spectrum-only analysis.

As a note to our previous statements, let us also stress here that the bispectrum does not account for all the cosmological information which can be extracted from the NG component of the tSZ map, since of course relevant contributions also come from higher order correlators. As long as most of the information is captured by the amplitude — and not by the shape — of different tSZ correlation functions, it was shown that the tSZ 1-point probability-density-function is a near-optimal statistic to constrain cosmology, since it encodes information from all n-point amplitudes [25, 26]. This method has been used in [7, 25] to constrain cosmological parameters, outperforming the skewness analysis. This fact shows again the value of the NG information, in this case applied to even higher order correlators. However, the implementation of the 1-point statistic in parameter inference relies on simplifying assumptions that might not be valid for future datasets. Its generalization is work in progress at the moment, making bispectrum estimation still a viable and worth pursuing alternative approach in practice. Moreover, we will show in the remainder that, while indeed small, the tSZ bispectrum shape dependence on cosmology is not negligible, and it actually plays an important role in breaking parameter degeneracies; this argument becomes even stronger when we account for the shape dependence on IntraCluster Medium (ICM) parameters.

The paper is structured as follows: in section 2 we review the halo model approach applied to the theoretical calculation of tSZ n-point correlation function, we show our choice of halo mass function and bias and we describe the analytical model of electron pressure profile which we use in the analysis; in section 3 we illustrate the calculation of the full covariance matrix for our observables; in section 4 we analyze the dependence of the tSZ power spectrum and bispectrum on different cosmological and ICM parameters; in section 5 we discuss foreground contamination issues; in section 6 we describe in detail and comment the results of our analysis. Finally, we summarize and draw our conclusions in section 7.

## 2 Theoretical tSZ correlation functions

The temperature fluctuations associated with the tSZ effect appear in sky maps as

$$ \frac{\Delta T}{T}(\nu, \hat{n}) = g(\nu) y(\hat{n}), \quad g(\nu) \equiv x \coth(x/2) - 4, \quad y(\hat{n}) \equiv \int \sigma_T \frac{k_B n_e T_e}{m_e c^2} dr. \quad (2.1) $$

Here $g(\nu)$ encapsulates the spectral dependence of the tSZ effect in units of CMB temperature with $x \equiv h\nu/(k_B T)$, while $y$, called the Compton-$y$ parameter, is the tSZ intensity along the line of sight $\hat{n}$. $\sigma_T$ is the Thomson cross section and $m_e$, $n_e$ and $T_e$ are respectively the electron rest mass, number density and temperature.

The statistical properties of the $y$ field can be extracted from its n-point angular correlation functions that we predict using the halo model approach [27, for a review]. Under some simplifying assumptions, discussed later, the expressions for all n-point correlators are well known [5, 7, 14, 18, 28, 29]. They can be derived via the formalism detailed in appendix A, where we express the $y$-polyspectra as the projection onto the past light-cone of the three-dimensional correlators of the corresponding three-dimensional field. This formulation is very
useful in the mathematical derivation and numerical implementation of the general y field n-point correlator.

The tSZ power spectrum is given by the sum of the Poisson one-halo term, which account for intra-halo correlations, and the two-halo term which models inter-halo correlations [8, 30]. They read

$$C_{\ell}^{1h} = \int dz \frac{d^2V}{dz d\Omega} \int dM \frac{dn}{dM}(z, M) |\tilde{y}_{\ell}(z, \tilde{M})|^2, \quad (2.2)$$

$$C_{\ell}^{2h} = \int dz \frac{d^2V}{dz d\Omega} D_{+}^2(z) P_m(k_{\ell}) \left[ \int dM \frac{dn}{dM}(z, M) b(z, M) \tilde{y}_{\ell}(z, \tilde{M}) \right]^2. \quad (2.3)$$

Here $P_m(k_\ell)$ is the linear matter power spectrum calculated at $k_\ell = (\ell + 1/2)/\chi(z)$, $D_+(z)$ is the growth factor, and $d^2V/dz d\Omega$ is the comoving volume element per steradian, which can be calculated as $d^2V/dz d\Omega = c\chi^2(z)/H(z)$, where $\chi(z)$ is the radial comoving distance. The Halo Mass Function (HMF) $\frac{dn}{dM}$ and halo bias $b$ will be described in section 2.1. $\tilde{y}_\ell(z, M)$ is the 2D Fourier transform of the projected y parameter image of the halo defined in detail in section 2.2. Notice that in eq. (2.2), (2.3), and throughout the paper, we also consider a hydrostatic mass bias rescaling the true mass of the halo, i.e., we use $\tilde{M} = (1 - b_{\text{HSM}})M$ in the expression for the projected y.

Higher order correlators are a straightforward extension of the power spectrum. Throughout this work, we will adopt the flat-sky approximation. Therefore, we will have to include the reduced bispectrum in our NG analysis. The bispectrum is described by one-, two- and three-halo terms. Including the three halo term would require the second order halo bias, which is not provided in [31]. To overcome this shortcoming, one could in principle resort to the peak-background-split formalism [32–35]; however, since the three-halo contribution is negligible, as shown in appendix B, we just omit it here. The one- and two-halo terms read

$$b_{\ell_1 \ell_2 \ell_3}^{1h} = \int dz \frac{d^2V}{dz d\Omega} \int dM \frac{dn}{dM}(z, M) \tilde{y}_{\ell_1}(z, \tilde{M}) \tilde{y}_{\ell_2}(z, \tilde{M}) \tilde{y}_{\ell_3}(z, \tilde{M}), \quad (2.4)$$

$$b_{\ell_1 \ell_2 \ell_3}^{2h} = \left[ \int dz \frac{d^2V}{dz d\Omega} D_+^2(z) P_m(k_{\ell_1}) \int dM \frac{dn}{dM}(z, M) b(z, M) \tilde{y}_{\ell_1}(z, \tilde{M}) \times \int dM \frac{dn}{dM}(z, M) b(z, M) \tilde{y}_{\ell_2}(z, \tilde{M}) \tilde{y}_{\ell_3}(z, \tilde{M}) \right] + 2 \text{ permutations.} \quad (2.5)$$

Finally, for the higher order correlators, we will always employ only the one-halo term, which in general reads

$$P_{(n)}^{1h}(\ell_1, \dots, \ell_n) = \int dz \frac{d^2V}{dz d\Omega} \int dM \frac{dn}{dM}(z, M) \prod_i \tilde{y}_{\ell_i}(z, \tilde{M}). \quad (2.6)$$

For later use, we point out that the general one-halo term does depend solely on the magnitudes of the multipoles involved.

Unless specified otherwise, in all the spectra we integrate over redshift between $z_{\text{min}} = 10^{-6}$ and $z_{\text{max}} = 4.5$, and over masses between $M_{\text{min}} = 10^{10} M_{\odot} h^{-1}$ and $M_{\text{max}} = 10^{16} M_{\odot} h^{-1}$. These limits, that stretch beyond the values commonly used in the literature, allow us to ensure that even higher order correlators are integrated correctly. We also point out that with the used parametrization the integrals are stable at low redshift. We integrate directly
over the overdensity mass $M_{500,c}$, as in, e.g., [7, 10], so that a conversion to the virial masses is never required, and also assume that $d\ln M_{\Delta,c}/d\ln M_{\text{vir}} \approx 1$ [10]. We use $b_{\text{HSM}} = 0.2$ [7] as the fiducial value throughout the analysis. We choose the values from [36] for the cosmological parameters $h = 0.6711$, $n_s = 0.9624$, $w_0 = -1$, but we pick $\Omega_m = 0.28$, $\sigma_8 = 0.8$ to be in agreement with [7].<fnref n="1" />

The fitting functions we employ, that we are about to describe, are the most commonly used in the literature. A possible alternative, proposed in [38], would be to jointly fit both the pressure profile and the HMF from simulations.

Before moving on, we ought to make some remarks about the modelling of the tSZ field. Since this effect is generated by hot electrons that mostly sit in galaxy clusters, generally speaking the halo model prescription is a valid approximation to calculate the theoretical signal. However, this means contributions from filaments are neglected. Along the same lines, the scatter around the mean pressure profile, cluster substructure and asphericity are not captured by the analytical model. It is not trivial to assess at what level the n-point correlation functions are biased and how much their errors increase because of these effects. In [39] it has been shown using simulations that they amount to a 10–20% power spectrum increase, substructures being the main contribution. For the bispectrum we can therefore expect at least the same amount of deviation, most likely with the configurations coupling different scales being more affected than the equilaterals.

Since, as we will show, signal modulations along the harmonic configurations play an important role, it will be important to study the bispectrum of hydrodynamical n-body simulations, in order to quantify, and possibly account for, this kind of contributions.

## 2.1 Halo mass function

We use the HMF and bias from [31], converting the mass definition $M_{500,c}$ into $M_{\Delta,m}$ to fit the parameters of their table 4 to the appropriate value of $\Delta = 500\rho_c/\rho_m$.

$$ \frac{dn}{dM} = \nu f(\nu) \frac{\bar{\rho}_m}{M} \frac{d\ln(\sigma^{-1})}{dM}, \qquad f(\nu) = \alpha \left[1 + (\beta\nu)^{-2\phi}\right] \nu^{2\eta} e^{-\gamma\nu^2/2}. \qquad (2.7) $$

Here the HMF is parametrized in terms of the peak height $\nu = \delta_c/\sigma(M)$, $\delta_c$ is the critical density for collapse and $\sigma^2$ is the linear matter variance. We assume that the redshift scaling provided for $\Delta = 200$ applies to any $\Delta$:

$$ \beta = \beta_0(1+z)^{0.2}, \quad \phi = \phi_0(1+z)^{-0.08}, \quad \eta = \eta_0(1+z)^{0.27}, \quad \gamma = \gamma_0(1+z)^{-0.01}, \qquad (2.8) $$

where the parameters at $z = 0$ are taken from table 4 of [31], using a linear interpolation along $\Delta$. As recommended, for each $z > 3$ we use the value calculated at $z = 3$, and at each $z$ we calculate $\alpha(z)$ imposing $\int d\nu f(\nu, z) b(\nu, z) = 1$. At $z = 0$ we correctly recover the value of $\alpha$ interpolated from the aforementioned table.

We use a halo bias of the form

$$ b(\nu) = 1 - A \frac{\nu^a}{\nu^a + \delta_c^a} + B\nu^b + C\nu^c, \qquad (2.9) $$

with the parameters from table 2 of [31].

<footnote n="1">¹For consistency with [7], we do not use the updated parameter from [37].</footnote>
## 2.2 Electron pressure profile

In our analysis, we are interested in the projected two-dimensional Compton-y field, which is obtained via line-of-sight integration of the rescaled electron density profile $P_e$, in any given direction of the sky through

$$ \tilde{y}_{\ell}(z, M) = \frac{\sigma_{\text{T}}}{m_e c^2} \frac{4\pi r_s}{\ell_s^2} \int dx x^2 j_0\left(\frac{\ell + 1/2}{\ell_s \chi(x)} x\right) P_e(x, z, M). \quad (2.10) $$

Here $r_s = r_s(z, M)$ and $\ell_s = a(z)\chi(z)/r_s(z, M)$ are, respectively, the typical scale radius of the y-image of the halo, and the multipole moment associated with it. For the parametrization of the electron profile that we employ, given in [40], $r_s = r_{500}$ and, for a single cluster, $P_e$ relies on the generalised Navarro-Frenk-White profile

$$ P_e(x, z, M_{500,c}) = \frac{C(x, z, M_{500,c}) \times P_0}{(c_{500}x)^{\gamma_G} [1 + (c_{500}x)^{\alpha_G}]^{\frac{\gamma_G - \beta_G}{\alpha_G}}}, \quad x \equiv \frac{r}{r_{500}}, \quad (2.11) $$

for which we use the best fit parameters obtained in [40]:<fnref n="2" />

$$ \{P_0, c_{500}, \gamma_G, \alpha_G, \beta_G\} = \{8.403 h_{70}^{-3/2}, 1.177, 0.3081, 1.0510, 5.4905\}. \quad (2.12) $$

As we rely on the halo model for the description of matter clustering, $r_{500}$ is the distance from the center of the halo, delimiting a sphere containing a density of dark matter which is 500 times the critical density of the Universe $\rho_c$. The function $C(x, z, M_{500,c})$ has the expression

$$ C(x, z, M_{500,c}) = 1.65 h_{70}^2 E(z)^{\frac{8}{3}} \left( \frac{h_{70} M_{500,c}}{3 \times 10^{14} M_{\odot}} \right)^{\frac{3}{2} + \alpha_P + \alpha'(x)} \text{eV cm}^{-3}, \quad (2.13) $$

where $\alpha_P = 0.12$, and the exponent $\alpha'_P(x)$, also fitted in [40], has the following dependence

$$ \alpha'_{\text{P}}(x) = 0.10 - (\alpha_{\text{P}} + 0.10) \frac{\left(\frac{x}{0.5}\right)^3}{1 + \left(\frac{x}{0.5}\right)^3}. \quad (2.14) $$

These exponents parametrise deviation from the standard self-similar case. In eqs. (2.12)–(2.13) we also introduced the widely used parameter $h_{70} \equiv h/0.7$.

## 3 Covariance matrix for the (binned) observables

In this work, we consider flat-sky binned estimators for the observables of interest as they are smooth and slowly varying functions in multipole space. The bins are defined in harmonic space: $\ell' \in \ell^b$ if $\ell - \Delta\ell^b/2 \le \ell' \le \ell + \Delta\ell^b/2$, $\ell$ being the central value of the bin $\ell^b$ of width $\Delta\ell^b$. As we are working in flat-sky, $\ell \in \mathbb{R}^2$. We are also considering a survey observing a sky fraction corresponding to a solid angle of $\Omega_{\text{sky}}$ steradians. Similarly to the three-dimensional matter field case, the binned angular power spectrum is defined as [41–43]

$$ \hat{C}_{\ell^b} \equiv \sum_{\ell \in \ell^b} \frac{\delta(\ell + \ell') \tilde{y}_{\ell} \tilde{y}_{\ell'}}{\Omega_{\text{sky}} N(\ell^b)}, \quad (3.1) $$

<footnote n="2">²We use the non-self-similar profile derived in the main text (cf. their eq. (12)).</footnote>
where the sum runs over discrete modes which are integer multiples of the fundamental frequency of our survey $\ell_f \equiv 2\pi/\Theta_{\text{sky}}$, $\Theta_{\text{sky}}$ being the survey linear angular size. We recall that

$$ \Omega_{\text{sky}} = 2\pi (1 - \cos \Theta_{\text{sky}}). \qquad (3.2) $$

The field $\tilde{y}_\ell$ is the two-dimensional Fourier transform of the $y$ field (2.1). The quantity $N(\ell^b) \approx 2\ell\Delta\ell^b f_{\text{sky}}$ (in the limit $\ell \gg \ell_f$) gives the number of vector pairs $\hat{\ell}, -\hat{\ell}$ whose magnitude $\hat{\ell}$ is within the bin $\ell^b$, each pair being discriminated by a deviation in the module of the vectors of a unit of the survey fundamental mode $\ell_f$ [44–46]. In a similar fashion to the three-dimensional matter field case, we define the binned angular bispectrum as [46–50]

$$ \hat{b}_{\ell_1^b \ell_2^b \ell_3^b} \equiv \sum_{\ell_1 \in \ell_1^b} \sum_{\ell_2 \in \ell_2^b} \sum_{\ell_3 \in \ell_3^b} \frac{\delta^{(2)}(\ell_1 + \ell_2 + \ell_3) \tilde{y}_{\ell_1} \tilde{y}_{\ell_2} \tilde{y}_{\ell_3}}{\Omega_{\text{sky}} N_{\text{tri.}}(\ell_1^b, \ell_2^b, \ell_3^b)}, \qquad (3.3) $$

where $N_{\text{tri.}}(\ell_1^b, \ell_2^b, \ell_3^b)$ is the number of valid triplets (i.e. respecting the triangular delta) in the proposed bin triplets [44–46]

$$ N_{\text{tri.}}(\ell_1^b, \ell_2^b, \ell_3^b) \approx \frac{\Omega_{\text{sky}}^2 \ell_1 \ell_2 \ell_3 \Delta \ell_1^b \Delta \ell_2^b \Delta \ell_3^b}{2\pi^3 \sqrt{2\ell_1^2 \ell_2^2 + 2\ell_1^2 \ell_3^2 + 2\ell_2^2 \ell_3^2 - \ell_1^4 - \ell_2^4 - \ell_3^2}}, \qquad (3.4) $$

where the approximation assumed is $\Delta \ell_i^b \gg \ell_f$. To speed up the calculations of the Fisher matrix, we assume that the power spectrum and the bispectrum are constant in each bin $\ell^b$, and can therefore be described by a single suitably picked representative $\bar{\ell}^b$.<fnref n="3" /> In practice, it is possible to prove that the 2 estimators above are unbiased and $\langle \hat{C}_{\ell^b} \rangle = C_{\bar{\ell}^b}$ and $\langle \hat{b}_{\ell_1^b \ell_2^b \ell_3^b} \rangle = b_{\bar{\ell}_1^b \bar{\ell}_2^b \bar{\ell}_3^b}$. In the following we drop the bars to distinguish the bin from its representative, as they should be distinguishable from the context. The representative considered in this work will be the central magnitude value $\bar{\ell}^b \equiv \ell$ of the bin $\ell^b$.

## 3.1 Structure of the covariance matrix

For our computation, we chose the binned estimators $\hat{C}_{\ell^b}$ (3.1) and $\hat{b}_{\ell_1^b \ell_2^b \ell_3^b}$ (3.3) for which we present the covariance matrix here below. We split the joint covariance as [51–54]

$$ \text{Cov}[\hat{C}_{\ell^b}, \hat{C}_{\ell'^b}] = \text{Cov}[\dots]_{\text{Gauss}} + \text{Cov}[\dots]_{\text{NG}}, \qquad (3.5) $$

$$ \text{Cov}[\hat{b}_{\ell_1^b \ell_2^b \ell_3^b}, \hat{b}_{\ell_1'^b \ell_2'^b \ell_3'^b}] = \text{Cov}[\dots]_{\text{Gauss}} + \text{Cov}[\dots]_{\text{NG}}, \qquad (3.6) $$

$$ \text{Cov}[\hat{C}_{\ell^b}, \hat{b}_{\ell_1^b \ell_2^b \ell_3^b}] = \text{Cov}[\dots]_{\text{NG}}. \qquad (3.7) $$

In eqs. (3.5) and (3.6), the subscript Gauss labels the covariance terms containing only 2-point statistics, which are non-vanishing only for correlations within the same $\ell$-bin. The other covariance terms arise due to the non-Gaussian statistics of the $y$ field (2.1) and correlate modes in different $\ell$-bins and the power spectrum and bispectrum. The derivation of the flat-sky joint power spectrum-bispectrum covariance has already been outlined for other projected scalar fields, such as the weak lensing convergence [46, 53]. As the derivation of the covariance for the $y$ field follows the same mathematical and conceptual steps, we refer to the above cited works for more insights on its rigorous derivation. In the following, we report the final expressions employed in the present work, along with the most important concepts associated to their structure.

<footnote n="3">³We discuss and validate our choice of binning and bin representative in appendix E.</footnote>
## 3.2 Power spectrum covariance matrix

Starting from the power spectrum binned estimator (3.1), we can obtain its full covariance by applying the following standard definition

$$ \mathrm{Cov}[\hat{C}_{\ell^b}, \hat{C}_{\ell'^b}] = \langle \hat{C}_{\ell^b} \hat{C}_{\ell'^b} \rangle - \langle \hat{C}_{\ell^b} \rangle \langle \hat{C}_{\ell'^b} \rangle. \quad (3.8) $$

A detailed calculation leads to a non-connected correlator of 4 instances of the field $\tilde{y}_\ell$ which can be split via Wick theorem into a sum of products of 2-point correlators and one 4-point connected correlator. The former can be further simplified, leading to the following Gaussian term

$$ \mathrm{Cov}[\hat{C}_{\ell^b}, \hat{C}_{\ell'^b}]_{\mathrm{Gauss}} = \frac{2\delta_{\ell\ell'}^K}{N(\ell^b)} (C_{\ell}^{\mathrm{n.}})^2, \quad C_{\ell}^{\mathrm{n.}} \equiv C_{\ell} + N_{\ell}. \quad (3.9) $$

In eq. (3.9) we potentially account for sources of Gaussian noise via the associated power spectrum $N_\ell$. The 4-point connected component leads instead to the NG term in eq. (3.5)

$$ \mathrm{Cov}[\hat{C}_{\ell^b}, \hat{C}_{\ell'^b}]_{\mathrm{NG},T} \approx \frac{1}{\Omega_{\mathrm{sky}}} T(\ell, -\ell, \ell', -\ell'), \quad (3.10) $$

where $T(\ell, -\ell, \ell', -\ell')$ is the trispectrum (A.13) of the tSZ field. In eq. (3.10) the exact covariance evaluation would require an average of the trispectrum over the two bins $\ell^b, \ell'^b$ involving both angular and magnitude integrations. However, we work under the common approximation of slowly varying polyspectra within the bins chosen for our work [46, 53]. We therefore assumed the averaged trispectrum being the same as the trispectrum computed at the central values of the bins.

## 3.3 Bispectrum covariance matrix

Similarly to the procedure outlined in the previous section, we can compute the covariance matrix for the flat-sky bispectrum binned estimator (3.3), shown in figure 1, left panel [46, 55, 56]. The Wick theorem allows us to derive the Gaussian part of the bispectrum covariance

$$ \mathrm{Cov}[\hat{b}_{\ell_1^b \ell_2^b \ell_3^b}, \hat{b}_{\ell_1'^b \ell_2'^b \ell_3'^b}]_{\mathrm{Gauss}} = \frac{\Omega_{\mathrm{sky}}}{N_{\mathrm{tri.}} (\ell_1^b, \ell_2^b, \ell_3^b)} C_{\ell_1}^{\mathrm{n.}} C_{\ell_2}^{\mathrm{n.}} C_{\ell_3}^{\mathrm{n.}} \left[ \delta_{\ell_1 \ell_1'}^{\mathrm{K}} (\delta_{\ell_2 \ell_2'}^{\mathrm{K}} \delta_{\ell_3 \ell_3'}^{\mathrm{K}} + \delta_{\ell_2 \ell_3'}^{\mathrm{K}} \delta_{\ell_3 \ell_2'}^{\mathrm{K}}) + \right. \\ \left. + 2 \text{ permutations of } (\ell_1' \leftrightarrow \ell_2') + 2 \text{ permutations of } (\ell_1' \leftrightarrow \ell_3') \right]. \quad (3.11) $$

Meanwhile, the non-Gaussian (NG) components read

$$ \mathrm{Cov}[\hat{b}_{\ell_1^b \ell_2^b \ell_3^b}, \hat{b}_{\ell_1'^b \ell_2'^b \ell_3'^b}]_{\mathrm{NG,BB}} = \frac{2\pi}{\Omega_{\mathrm{sky}}} b_{\ell_1, \ell_2, \ell_3} b_{\ell_1, \ell_2, \ell_3'} \left[ \frac{1}{\ell_1 \Delta \ell_1^b} (\delta_{\ell_1 \ell_1'}^{\mathrm{K}} + \delta_{\ell_1 \ell_2'}^{\mathrm{K}} + \delta_{\ell_1 \ell_3'}^{\mathrm{K}}) + \right. \\ \left. + 3 \text{ permutations of } (\ell_1 \leftrightarrow \ell_2) + 3 \text{ permutations of } (\ell_1 \leftrightarrow \ell_3) \right], \quad (3.12) $$

$$ \mathrm{Cov}[\hat{b}_{\ell_1^b \ell_2^b \ell_3^b}, \hat{b}_{\ell_1'^b \ell_2'^b \ell_3'^b}]_{\mathrm{NG,PT}} = \frac{2\pi}{\Omega_{\mathrm{sky}}} \frac{C_{\ell_1}^{\mathrm{n.}}}{\ell_1 \Delta \ell_1^b} \left[ T(\ell_2, \ell_3, \ell_2', \ell_3') \delta_{\ell_1 \ell_1'}^{\mathrm{K}} + T(\ell_2, \ell_3, \ell_1', \ell_3') \delta_{\ell_1 \ell_2'}^{\mathrm{K}} + \right. \\ \left. + T(\ell_2, \ell_3, \ell_1', \ell_2') \delta_{\ell_1 \ell_3'}^{\mathrm{K}} \right] + 3 \text{ permutations of } (\ell_1 \leftrightarrow \ell_2) + 3 \text{ permutations of } (\ell_1 \leftrightarrow \ell_3). \quad (3.13) $$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Left) Full bispectrum covariance matrix. Right) Fractional amplitude of the bispectrum covariance terms different from $\text{Cov}[\hat{b}_{\ell_1^{\ell_2} \ell_3^{\ell_3}}, \hat{b}_{\ell_1^{\ell_2} \ell_3^{\ell_3}}]_{\text{NG},6\text{P}}$ compared to the full covariance. The white areas correspond to configurations for which the covariance is fully determined by $\text{Cov}_{\text{BB}}^{(6)}$ for which the plotted quantity is not defined. In both panels, the ticks mark equilateral configurations with sides of the labelled length.</figcaption></figure>

By employing the estimator definitions, we can also obtain an expression for the cross-covariance between the binned flat-sky observables used in our analyses

$$ \text{Cov}[C_{\ell^b}, \hat{b}_{\ell_1^b \ell_2^b \ell_3^b}]_{\text{NG,PB}} = \frac{4\pi}{\Omega_{\text{sky}}} \left[ \frac{C_{\ell}^{\text{n.}}(\ell) b_{\ell, \ell_2, \ell_3}}{\ell_1 \Delta \ell_1^b} \delta_{\ell \ell_1}^{\text{K}} + 1 \text{ permutation of } (\ell_1 \leftrightarrow \ell_2) + \right. \\ \left. + 1 \text{ permutation of } (\ell_1 \leftrightarrow \ell_3) \right]. \quad (3.14) $$

So far, all the equations listed in the present section are exact. The two terms still missing are those associated to the 6- and 5-point correlator arising from the bispectrum covariance and the cross-covariance calculation. In these cases, the spectra should be averaged over the bins involved. However, as it was the case for the NG component of the power spectrum covariance (3.10), we assume them to be almost constant over the bins thus avoiding this computationally expensive operation. Therefore, the final expressions are

$$ \text{Cov}[\hat{b}_{\ell_1^{\ell_2} \ell_3^{\ell_3}}, \hat{b}_{\ell_1^{\ell_2} \ell_3^{\ell_3}}]_{\text{NG},6\text{P}} \approx \frac{1}{\Omega_{\text{sky}}} P_{(6)}(\ell_1, \ell_2, \ell_3, \ell_1', \ell_2', \ell_3'), \quad (3.15) $$

$$ \text{Cov}[\hat{C}_{\ell^b}, \hat{b}_{\ell_1^b \ell_2^b \ell_3^b}]_{\text{NG},5\text{P}} \approx \frac{1}{\Omega_{\text{sky}}} P_{(5)}(\ell, -\ell, \ell_1, \ell_2, \ell_3). \quad (3.16) $$

In the right panel of figure 1 we show the fraction of the total covariance due to the 6-point connected term we just described, which is the dominant one: as shown in figure 2 and 3, it is possible to identify a hierarchy among the different components listed above. In these figures we use 27 logarithmically spaced bins spanning from $\ell = 10$ to $\ell = 5000$. The $\text{Cov}[\dots]_{\text{NG,T}}$ and the $\text{Cov}[\dots]_{\text{NG},6\text{P}}$ component dominating the power spectrum and the bispectrum covariance, respectively, clearly emphasises the level of non-Gaussianities in the tSZ field.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Different components of the diagonal part for the tSZ power spectrum covariance matrix. For comparison we show both the variance provided directly by *Planck*, and the one we adapted in our analysis for which we employ slightly finer bins.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Selected components of the tSZ bispectrum covariance matrix. Whenever of relevance, the impact of the Gaussian Planck-like noise is displayed through markers. Left) equilateral configurations. Right) squeezed configuration of type {~22, ℓ, ℓ}.</figcaption></figure>

# 4 Power spectrum and bispectrum dependence on the parameters

It is well known that varying cosmological and gas parameters results in an amplitude shift of the power spectrum [e.g. 10, 57] and a slight tilt at smaller scales, where the power spectrum has a stronger dependence on the shape of the halo profile. As a consequence, there are important degeneracies among the various parameters and, both with current and future data, one needs to consider additional observables to constrain various parameter jointly, either exploiting other datasets [e.g. 57] or other statistics extracted from the same tSZ maps [18].

Overall, the bispectrum also displays a similar dependence on different parameters, leading again to an amplitude shift as the main effect. This comes from the fact that tSZ correlators are dominated by the one-halo term, making the overall tSZ statistic nearly Poisson. However, for the bispectrum we can also appreciate a non-negligible, configuration dependent modulation, which differs for varying parameters.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Left and central column) logarithmic derivative of the power spectrum (red solid line) and of the bispectrum (dark blue line) with respect to different parameters of interest. Whenever the bispectrum is considered, we show equilateral configurations (solid line) and squeezed configurations of kind {~22, *l*, *l*} (dotted line). The parameters with respect to which the derivatives are taken are: $\beta_G$ (top-left panel), *h* (top-central panel), $\sigma_8$ (bottom-left panel) and $\Omega_m$ (bottom central panel). Right column) fractional difference between the two derivatives on the left (per row). Color and line style code is the same as above.</figcaption></figure>

Let us start by considering the bispectrum derivatives for the equilateral ($\ell_1 = \ell_2 = \ell_3$) and “squeezed” configurations ($\ell_1 \approx 22$, $\ell_2 = \ell_3$), and compare them to their power spectrum counterparts in figure 4.<fnref n="4" /> Much like the power spectrum, the derivatives of the equilateral bispectrum are similar to each other, most of them being somewhat flat up to $\ell \approx 10^3$ where the halo inner structure starts to be resolved, with a subsequent reduction (in absolute value) on smaller scales. On the other hand, the derivatives of the squeezed bispectrum have a slightly more variegated phenomenology among themselves and, on top of that, they differ from the equilateral bispectrum ones. For this reason we can intuitively expect squeezed configurations to be less affected by parameter degeneracies. While in section 6.1.1 we will confirm this expectation, we will also find that the conditional errors will be much bigger to begin with than what we find using only equilateral configurations, further motivating the

<footnote n="4">⁴We do not add the further requirement $\ell_{2,3} \gg \ell_1$, so not all the configurations in the squeezed set are technically “squeezed”.</footnote>
analysis of the full bispectrum. Using the bispectrum parametrization introduced in [58],

$$ 
\begin{align} 
P &\equiv \ell_1 + \ell_2 + \ell_3, \nonumber \\ 
\tilde{\sigma}_2 &\equiv 12 \frac{(\ell_1\ell_2 + \ell_1\ell_3 + \ell_2\ell_3)}{(\ell_1 + \ell_2 + \ell_3)^2} - 3, \tag{4.1} \\ 
\tilde{\sigma}_3 &\equiv 27 \frac{\ell_1\ell_2\ell_3}{(\ell_1 + \ell_2 + \ell_3)^3}, \nonumber 
\end{align} 
$$

in figure 5 we show how each configuration responds to parameter changes in a slightly different manner. Even though each one has lower significance than the power spectrum, this means that it is easier to discern the effects of the various parameter changes. We will quantify this statement in section 6. The reason why different bispectrum configurations respond differently to changes in the various parameters lies in the slightly different ranges of halo masses and redshifts that contribute to each. This is investigated in appendix F.

# 5 Gaussian experimental noise and foregrounds

To forecast the performance of real instruments we need to account for instrumental noise and foregrounds.

To assess the impact of their Gaussian component on the forecast we employ the methodology of [8, 13]. We assume to be using multi-channel measurements to run a simple internal linear combination (ILC) to separate foregrounds from the signal. To do so, we will assume perfect knowledge of the foreground spectral energy densities, which might be unrealistic in a real life scenario [59–61]. We leave the assessment of the impact of more refined foreground modelling to future work.

In terms of non-Gaussianity, the dominant contributions are galactic foregrounds — thermal dust and synchrotron — which peak on large angular scales. Quantitatively, the Planck analysis [62] showed that these contributions are negligible after removing the multipoles $\ell \le 64$. For this reason, we will later investigate the Fisher matrix dependence on this type of constraint. Other contributions might be relevant on smaller angular scales. In [63] it has been shown how the Cosmic Infrared Background (CIB) and tSZ auto- and cross-bispectra are correlated. Therefore, if unaccounted for, CIB might introduce biases. Given that the amplitudes of CIB bispectra are smaller than the tSZ auto-bispectrum, and the fact that the spectral energy density of the two effects differ, we neglect this contribution for the time being, assuming those contributions can be effectively removed. However, prior to any analysis of actual data, non-Gaussianity of the foregrounds will require a thorough assessment to verify that residuals in the component separated $y$ maps do not introduce sizeable biases in the estimated parameters.

## 5.1 Spectral components

A comprehensive list of the spectral components that affect tSZ measurements is given by instrumental noise, CMB, free-free, thermal dust, synchrotron, radio point sources, infrared point sources. We will consider all of these contributions even if we use approximate descriptions that are, in any case, commonly used in the literature.

We assume we can factorize the inter-frequency correlation in three ingredients: the spectral energy density expressed with respect to the CMB black-body $\Theta^{\text{sc}}(\nu)$, a spectral coherence factor that parametrizes the correlation between different instrumental frequency channels $R(\nu_i, \nu_j, \xi^{\text{sc}})$, and an angular scaling $C_{\ell}^{\text{sc}}$:

$$ C^{\text{sc}}(\nu_i, \nu_j) = \Theta^{\text{sc}}(\nu_i) \Theta^{\text{sc}}(\nu_j) R(\nu_i, \nu_j, \xi^{\text{sc}}) C_{\ell}^{\text{sc}}. \qquad (5.1) $$
<figure><img src="image_6.png" /><figcaption>Figure 5. Logarithmic derivative of the bispectrum with respect to $\Omega_m$ (top), $\sigma_8$ (middle), and $\alpha_P$ (bottom). On the right we also provide a rescaled plot for the perimeter value of P = 8130 to illustrate the differences in derivatives for different triangles.</figcaption></figure>

In [64] the spectral coherence in two channels $\nu_i$ and $\nu_j$, for each component, was parametrized as

$$R(\nu_i, \nu_j, \xi^{\text{sc}}) = \exp \left[ -\frac{1}{2} \left( \frac{\log(\nu_i) - \log(\nu_j)}{\xi^{\text{sc}}} \right)^2 \right], \quad \xi^{\text{sc}} \approx \frac{1}{\sqrt{2\Delta\alpha^{\text{sc}}}}, \quad (5.2)$$

where, for each spectral component, $\Delta\alpha^{\text{sc}}$ is the variance of the spectral index over the sky. Depending on the spectral component, $\xi^{\text{sc}}$ varies in the range $[0, \infty)$. In particular, we assume the instrument frequency channels to be independent ($\xi^{\text{sc}} \to 0$), and notice that the CMB has perfect correlation across all the frequencies ($\xi^{\text{sc}} \to \infty$).
$\Theta^{\text{sc}}$, $\Delta\alpha^{\text{sc}}$, and $C_{\ell}^{\text{sc}}$ vary for each spectral component. We will use the values and functional forms provided in [8, 60, 64, 65], that we report in appendix C.

## 5.2 Component separation

The ILC aims at recovering, multipole by multipole, a signal of known spectral dependence computing a weighted average of data collected at different frequencies.<fnref n="5" /> The weights are calculated in such a way to have unitary response to the desired spectral shape, while minimizing spurious contributions from other spectral components.

The internal linear combination weights and the effective noise term after multifrequency subtraction are respectively [13]

$$ \vec{w} = \frac{[\sum_{\text{sc}} C^{\text{sc}}(\nu_i, \nu_j)]^{-1} \vec{g}}{\vec{g} \cdot [\sum_{\text{sc}} C^{\text{sc}}(\nu_i, \nu_j)]^{-1} \vec{g}}, \quad N = \vec{g} \cdot [\sum_{\text{sc}} C^{\text{sc}}(\nu_i, \nu_j)]^{-1} \vec{g}, \quad (5.3) $$

where $\vec{g} = (g(\nu_1), \dots, g(\nu_n))^T$ is a vector of the value of the tSZ spectral shape evaluated on the observed channels. This result, valid for a single multipolar coefficient, for each $\ell$ can be averaged over all $m$: we can therefore express the effective noise term $N_{\ell} = N/(2\ell + 1)$.<fnref n="6" />

## 6 Forecast results

We now have all the ingredients to understand what is the information content of the power spectrum, of the full bispectrum, and of their combination. We use a Fisher forecast to estimate how well a particular experiment can constrain a parameter through some observed quantity. The components of the Fisher information matrix are defined as [66]:

$$ F_{ij} = \left\langle -\frac{\partial^2(\ln \mathcal{L})}{\partial \theta_i \partial \theta_j} \right\rangle, \quad (6.1) $$

where $\mathcal{L}$ is the likelihood function and the $\theta_i$ are the parameters we want to constrain.

In case of Gaussian distributed data and rotational invariance of the observable, the Fisher matrix can also be written as [67]:

$$ F_{ij} = \frac{\partial \mathcal{O}^T}{\partial \theta_i} \mathrm{Cov}^{-1} \frac{\partial \mathcal{O}}{\partial \theta_j} + \frac{1}{2} \mathrm{Tr} \left[ \mathrm{Cov}^{-1} \frac{\partial \mathrm{Cov}}{\partial \theta_i} \mathrm{Cov}^{-1} \frac{\partial \mathrm{Cov}}{\partial \theta_j} \right], \quad (6.2) $$

where $\mathrm{Cov} = \mathrm{Cov}[\mathcal{O}, \mathcal{O}]$ is the covariance matrix, described in section 3. We drop the second term as all the relevant information is contained in the mean of the observable, while the covariance dependence on the parameters would introduce spurious information that cannot be extracted by the considered estimator [67]. In principle neither the power spectrum nor the bispectrum follow a Gaussian distribution, but we expect it to be a good approximation for the binned spectra in virtue of the central limit theorem. This ansatz should however be validated against simulations or more refined analysis [68, 69] that we leave for future work.

For each parameter the conditional error is:

$$ \sigma_{i, \text{conditional}} = \frac{1}{\sqrt{F_{ii}}}, \quad (6.3) $$

<footnote n="5">⁵All quantities in this section (except $N_{\ell}$) are calculated separately for each $(\ell, m)$.</footnote>

<footnote n="6">⁶As a sanity check we successfully reproduced figure 12 of [8]. To do so one has to use $N_{\ell}^{\text{sc}} = \vec{w} \cdot \mathbb{C}^{\text{sc}}\vec{w}$, where the $\vec{w}$ are still calculated according to eq. (5.3).</footnote>
which is the error one gets fitting only one parameter while fixing all the others. What we are more interested in, however, is the marginalized error. The diagonal entries of the inverse of a Fisher matrix give the error for each parameter marginalized over all the others, which one would obtain by a multivariate fit:

$$ \sigma_i = \sqrt{\mathrm{Cov}(\theta_i, \theta_i)} = \sqrt{(F^{-1})_{ii}}. \qquad (6.4) $$

The off diagonal elements then give us the covariance between two parameters, Cov($\theta_i, \theta_j$).

We start by considering an ideal, noiseless experiment in the absence of foregrounds, to understand the general properties of the spectra in a simplified case, then we show a forecast for the already available Planck data, and finally we consider two realistic future surveys.

For the reader convenience, we remind that the parameter fiducial values we assume are: $\Omega_m = 0.28$, $\sigma_8 = 0.8$, $h = 0.6711$, $n_s = 0.9624$, $w_0 = -1$, $b_{\mathrm{HSM}} = 0.2$, $P_0 = 8.403 h^{-3/2}$, $c_{500} = 1.177$, $\gamma_G = 0.3081$, $\alpha_G = 1.0510$, $\beta_G = 5.4905$.

## 6.1 Noiseless survey with no foreground contamination

In order to set an upper limit on the amount of information that can be in principle extracted from the tSZ power spectrum and bispectrum, we start by considering the case of a noiseless survey while also neglecting any issue of foreground contamination. This will also help us to build understanding of the main factors which affect the final constraints. For brevity we will refer to this configuration, slightly improperly, as Cosmic Variance Limited (CVL) survey. This initial oversimplified analysis will then be generalized in the following sections, including all realistic effects.

In tables 1 and 2 we show the forecasted error bars on all the parameters for a full sky survey, with $\ell_{\mathrm{max}} = 1000$ and $\ell_{\mathrm{max}} = 5000$ respectively. To assess the impact of a change in $\ell_{\mathrm{min}}$ instead, we repeated both analysis switching from $\ell_{\mathrm{min}} = 10$ to $\ell_{\mathrm{min}} = 70$ and found quantitatively negligible differences. In all these cases we have employed bins logarithmically spaced in $\ell$, 20 when $\ell \in [10, 1000]$ and [70, 1000], and 27 when $\ell \in [10, 5000]$ and [70, 5000]. The choice of binning is further discussed in appendix E.

In both cases, if one considers the conditional errors, the power spectrum outperforms the bispectrum analyzed on its own by a factor $2 \sim 10$; obviously when the two signals are analyzed jointly, the errors shrink marginally with respect to the power spectrum case. If all but one parameters are perfectly known, the power spectrum, comparatively bigger than its noise with respect to the bispectrum, leads to tighter constraints on the last parameter. This outcome is overturned by the marginalization needed for a joint fit of the parameters.

As anticipated in section 4, at power spectrum level the parameters have serious degeneracies among all of them, whereas the bispectrum has a more diverse response to parameters change. To quantify this statement, in figure 6 we compare for each couple of parameters the Pearson correlation coefficient

$$ \mathrm{Corr}(\theta_i, \theta_j) = \frac{\mathrm{Cov}(\theta_i, \theta_j)}{\sqrt{\mathrm{Cov}(\theta_i, \theta_i) \mathrm{Cov}(\theta_j, \theta_j)}}, \qquad (6.5) $$

that one obtains using the power spectrum, the bispectrum, and their combination, in the case of $\ell_{\mathrm{max}} = 5000$. Correlations are generally much higher for the power spectrum and for this reason after marginalization the power spectrum loses much of its constraining power. The same information is also displayed in figure 7, where we show the triangle plot of cosmological and gas parameters.
<table><caption>**Table 1.** Conditional and marginalized error forecasts for a CVL experiment with $\ell_{max} = 1000$ and $f_{sky} = 1$.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Conditional</th><th colspan="3">Marginalized</th></tr><tr><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.0021</td><td>0.0092</td><td>0.0015</td><td>8.1</td><td>0.12</td><td>0.092</td></tr><tr><td>h</td><td>0.013</td><td>0.033</td><td>0.010</td><td>16</td><td>0.36</td><td>0.34</td></tr><tr><td>σ<sub>8</sub></td><td>0.0020</td><td>0.0075</td><td>0.0015</td><td>3.4</td><td>0.047</td><td>0.047</td></tr><tr><td>n<sub>S</sub></td><td>0.060</td><td>0.060</td><td>0.051</td><td>8.0</td><td>0.43</td><td>0.30</td></tr><tr><td>w<sub>0</sub></td><td>0.096</td><td>0.31</td><td>0.071</td><td>105</td><td>0.59</td><td>0.54</td></tr><tr><td>b<sub>HSM</sub></td><td>0.0060</td><td>0.023</td><td>0.0043</td><td>407</td><td>0.93</td><td>0.81</td></tr><tr><td>P<sub>0</sub></td><td>0.097</td><td>0.36</td><td>0.071</td><td>8212</td><td>14</td><td>12</td></tr><tr><td>α<sub>P</sub></td><td>0.12</td><td>0.071</td><td>0.069</td><td>60</td><td>0.20</td><td>0.18</td></tr><tr><td>c<sub>500,c</sub></td><td>0.0059</td><td>0.024</td><td>0.0041</td><td>4.4</td><td>0.065</td><td>0.064</td></tr><tr><td>α<sub>G</sub></td><td>0.0041</td><td>0.015</td><td>0.0029</td><td>64</td><td>0.065</td><td>0.064</td></tr><tr><td>β<sub>G</sub></td><td>0.023</td><td>0.030</td><td>0.014</td><td>1.8</td><td>0.034</td><td>0.034</td></tr><tr><td>γ<sub>G</sub></td><td>0.011</td><td>0.034</td><td>0.0085</td><td>245</td><td>0.24</td><td>0.20</td></tr></tbody></table>

<table><caption>**Table 2.** Conditional and marginalized error forecasts for a CVL experiment with $\ell_{max} = 5000$ and $f_{sky} = 1$.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Conditional</th><th colspan="3">Marginalized</th></tr><tr><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.00044</td><td>0.0017</td><td>0.00031</td><td>0.86</td><td>0.025</td><td>0.018</td></tr><tr><td>h</td><td>0.0040</td><td>0.0079</td><td>0.0033</td><td>5.7</td><td>0.083</td><td>0.076</td></tr><tr><td>σ<sub>8</sub></td><td>0.00040</td><td>0.0013</td><td>0.00030</td><td>0.43</td><td>0.011</td><td>0.010</td></tr><tr><td>n<sub>S</sub></td><td>0.0050</td><td>0.026</td><td>0.0030</td><td>0.77</td><td>0.080</td><td>0.063</td></tr><tr><td>w<sub>0</sub></td><td>0.0052</td><td>0.014</td><td>0.0039</td><td>11</td><td>0.10</td><td>0.092</td></tr><tr><td>b<sub>HSM</sub></td><td>0.0013</td><td>0.0044</td><td>0.00091</td><td>78</td><td>0.16</td><td>0.16</td></tr><tr><td>P<sub>0</sub></td><td>0.018</td><td>0.060</td><td>0.014</td><td>1754</td><td>2.7</td><td>2.6</td></tr><tr><td>α<sub>P</sub></td><td>0.0050</td><td>0.022</td><td>0.0030</td><td>2.81</td><td>0.040</td><td>0.037</td></tr><tr><td>c<sub>500,c</sub></td><td>0.0014</td><td>0.0054</td><td>0.00097</td><td>0.90</td><td>0.016</td><td>0.015</td></tr><tr><td>α<sub>G</sub></td><td>0.00082</td><td>0.0026</td><td>0.00059</td><td>7.84</td><td>0.015</td><td>0.014</td></tr><tr><td>β<sub>G</sub></td><td>0.0063</td><td>0.0067</td><td>0.0036</td><td>0.31</td><td>0.0077</td><td>0.0077</td></tr><tr><td>γ<sub>G</sub></td><td>0.0017</td><td>0.0045</td><td>0.0013</td><td>50</td><td>0.039</td><td>0.036</td></tr></tbody></table>

To show the dramatic impact of degeneracies and the incremental reduction of constraining power when more and more parameters are jointly fitted, in figure 8 we directly compare the error bars recovered in various scenarios. For each parameter the baseline (tightest possible constraints in our analysis) is the conditional error with a joint power spectrum and bispectrum fit. All other errors are shown in figure 8 as ratio to this value. For each parameter the top error bar is calculated with the power spectrum analysis, the middle one with bispectrum analysis and the bottom one with the joint fit of the two. For each bar the most saturated and least saturated color show the conditional and marginalized error, respectively. In the left panel the mid-tone bar is obtained fixing the value of all gas parameters and marginalizing over cosmological ones; in the right panel we did the opposite.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Comparison of the Pearson correlation coefficients for all couples of model parameters. Broadly speaking the correlations among the parameters are clearly higher when one considers the power spectrum with respect to the bispectrum.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Triangle plots for cosmological (left) and gas parameters (right). Here we consider a Cosmic variance limited experiment with perfect foreground separation and $l_{max} = 5000$. Notice that the power spectrum $1\sigma$ ellipses have been rescaled by a factor to fit in the same graph. The grey bands are the power spectrum conditional errors.</figcaption></figure>

Having made the point that the bispectrum allows us to separately fit multiple parameters, we now inquire what are the prospects considering tSZ observations combined to prior knowledge of cosmological and gas parameters. Generally, the gas parameters are fixed to the best-fit values of simulations [39], external X-Ray [40], or stacked tSZ clusters [70] measurements, and not varied throughout the analysis of cosmological parameters. On the other hand, one can also think of using the *Planck* primary anisotropies measurements.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Impact of degeneracies on the parameters error bars. (Left) for each cosmological parameter we show the error marginalized over all others parameter, the error marginalized over cosmological parameters, fixing the gas ones, and the conditional error; all divided by the conditional error of power spectrum and bispectrum combined. This is repeated for (top to bottom in each triplet of bars) power spectrum, bispectrum, and combination of the two. The constraining power of the power spectrum is hindered by the marginalization over other parameters; on the other hand the bispectrum, in principle less constraining, is less affected. (Right) same, but switching cosmological and gas parameters.</figcaption></figure>

ments to set a prior on the cosmological parameters, and exploit tSZ anisotropies data to cross-check the gas parameters. We explore both strategies in table 3. In the left side we fixed ($P_0, \alpha_P, c_{500,c}, \alpha_G, \beta_G, \gamma_G$), removing the respective rows and columns from the Fisher matrix. In the right side, instead, we add a (Gaussian) prior to the cosmological parameters, obtained from the $1\sigma$ errors in the column TT,TE,EE+lowE+lensing+BAO in table 2 of [37], with the additional assumption of the errors being uncorrelated. Since our paper is meant to be a proof of concept, and a comparison between the power spectrum and the bispectrum, we refrain from comparing the results from table 3 with state of the art observations for both cosmological and gas parameters, as our forecast are bound to be an upper limit and most likely over-optimistic. However, they serve as a comparison between the power spectrum and the bispectrum efficiency and to motivate the analysis of both.

Our results heavily depend, through the covariance, on high order correlations functions. In particular we study the impact of massive halos at low redshift, since higher order correlators are incrementally sensitive to the signal from these clusters. In principle, to calculate the y-distortion along the line of sight one has to integrate from 0 to infinity, and again in principle the HMF has been calibrated in [31] integrating clusters up to a mass of $10^{16} M_\odot$. However, the most massive cluster in the Planck SZ cluster catalogue [62] has a mass of $1.61 \times 10^{15} M_\odot$, whereas the closest detected cluster is at $z = 0.011$. To account for these observational constraints we repeat our analysis enforcing these values as upper limits in the redshift and mass integrations. The results that follow are shown in table 4 which is the main result of the paper.
<table><caption>**Table 3.** (Left) marginalized error forecasts for cosmological parameters assuming the gas parameters are fixed. (Right) marginalized error forecast for gas parameters using Planck priors on the cosmological parameters. In this table we always considered a CVL experiment with $\ell_{max} = 5000$ and $f_{sky} = 1$.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Gas fixed</th><th rowspan="2">P<sub>0</sub></th><th colspan="3">Planck prior</th></tr><tr><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.078</td><td>0.023</td><td>0.016</td><td>P<sub>0</sub></td><td>1230</td><td>2.2</td><td>2.0</td></tr><tr><td>h</td><td>0.89</td><td>0.074</td><td>0.068</td><td>α<sub>P</sub></td><td>1.5</td><td>0.029</td><td>0.012</td></tr><tr><td>σ<sub>8</sub></td><td>0.26</td><td>0.010</td><td>0.010</td><td>c<sub>500,c</sub></td><td>0.39</td><td>0.015</td><td>0.015</td></tr><tr><td>n<sub>S</sub></td><td>0.45</td><td>0.069</td><td>0.047</td><td>α<sub>G</sub></td><td>6.9</td><td>0.014</td><td>0.014</td></tr><tr><td>w<sub>0</sub></td><td>1.1</td><td>0.050</td><td>0.035</td><td>β<sub>G</sub></td><td>0.15</td><td>0.0076</td><td>0.0075</td></tr><tr><td>b<sub>HSM</sub></td><td>1.3</td><td>0.086</td><td>0.069</td><td>γ<sub>G</sub></td><td>41</td><td>0.039</td><td>0.036</td></tr></tbody></table>

<table><caption>**Table 4.** Same as table 2 (CVL, $\ell_{max} = 5000$, $f_{sky} = 1$) but with the halo redshift and mass boundaries observationally set by Planck.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Conditional</th><th colspan="3">Marginalized</th></tr><tr><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.00041</td><td>0.0016</td><td>0.00030</td><td>0.59</td><td>0.11</td><td>0.037</td></tr><tr><td>h</td><td>0.0038</td><td>0.0072</td><td>0.0033</td><td>13</td><td>0.60</td><td>0.49</td></tr><tr><td>σ<sub>8</sub></td><td>0.00038</td><td>0.0013</td><td>0.00029</td><td>0.35</td><td>0.065</td><td>0.035</td></tr><tr><td>n<sub>S</sub></td><td>0.0046</td><td>0.022</td><td>0.0030</td><td>11</td><td>0.50</td><td>0.39</td></tr><tr><td>w<sub>0</sub></td><td>0.0048</td><td>0.015</td><td>0.0040</td><td>5.3</td><td>0.51</td><td>0.25</td></tr><tr><td>b<sub>HSM</sub></td><td>0.0012</td><td>0.0039</td><td>0.00089</td><td>5.6</td><td>0.89</td><td>0.77</td></tr><tr><td>P<sub>0</sub></td><td>0.017</td><td>0.056</td><td>0.013</td><td>147</td><td>13</td><td>9.1</td></tr><tr><td>α<sub>P</sub></td><td>0.0046</td><td>0.021</td><td>0.0031</td><td>2.5</td><td>0.15</td><td>0.13</td></tr><tr><td>c<sub>500,c</sub></td><td>0.0013</td><td>0.0049</td><td>0.00092</td><td>0.33</td><td>0.12</td><td>0.10</td></tr><tr><td>α<sub>G</sub></td><td>0.00078</td><td>0.0024</td><td>0.00057</td><td>2.5</td><td>0.086</td><td>0.078</td></tr><tr><td>β<sub>G</sub></td><td>0.0058</td><td>0.017</td><td>0.0038</td><td>2.2</td><td>0.46</td><td>0.37</td></tr><tr><td>γ<sub>G</sub></td><td>0.0016</td><td>0.0043</td><td>0.0013</td><td>4.0</td><td>0.15</td><td>0.11</td></tr></tbody></table>

Removing massive nearby clusters has a deeper impact on higher order correlators with respect to the power spectrums and the bispectrum, therefore the signal to noise increases for both observables. For the power spectrum, this translates to proportionally tighter constraints on the parameters, compared to the fiducial model. For the bispectrum the picture is more complicated: conditional error bars shrink as expected but correlations among different parameters become more severe, enlarging the marginalized error bars. In any case, the marginalized errors obtained with the bispectrum, even with this model, still improve the PS one by approximately one order of magnitude. Therefore, our conclusions are qualitatively unchanged. Even more than in the fiducial case, combining power spectrum and bispectrum improves the constraining power and reliability of the results. Since the correlations are particularly severe for a subset of parameters, fixing even a limited number of them would allow for recovering most of the constraining power we forecasted with the Tinker fiducial model, as seen in figure 11. Nevertheless we conservatively quote as marginalized errors the results obtained without fixing any parameter.
<table><caption>**Table 5.** Ratio of the equilateral (Eq.), squeezed (Sq.), joint equilateral-plus-squeezed (Eq.+Sq.), and full bispectrum error to the power spectrum one. After marginalization, the constraining power of the bispectrum in one single limit is degraded comparably to the power-spectrum. The same does not apply when we consider the combination of equilateral and squeezed limit. Obviously, the full bispectrum — that encompasses many other configurations, orthogonals and all the intermediate ones — is even more stable under the marginalization.</caption><thead><tr><th rowspan="2"></th><th colspan="4">Conditional</th><th colspan="4">Marginalized</th></tr><tr><th>σ<sup>Eq</sup>/σ<sup>PS</sup></th><th>σ<sup>Sq</sup>/σ<sup>PS</sup></th><th>σ<sup>Eq+Sq</sup>/σ<sup>PS</sup></th><th>σ<sup>Full</sup>/σ<sup>PS</sup></th><th>σ<sup>Eq</sup>/σ<sup>PS</sup></th><th>σ<sup>Sq</sup>/σ<sup>PS</sup></th><th>σ<sup>Eq+Sq</sup>/σ<sup>PS</sup></th><th>σ<sup>Full</sup>/σ<sup>PS</sup></th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>5.9</td><td>86</td><td>5.9</td><td>4.0</td><td>10</td><td>5.7</td><td>0.65</td><td>0.029</td></tr><tr><td>h</td><td>2.6</td><td>29</td><td>2.6</td><td>2.0</td><td>13</td><td>23</td><td>0.99</td><td>0.015</td></tr><tr><td>σ<sub>8</sub></td><td>4.7</td><td>74</td><td>4.7</td><td>3.3</td><td>11</td><td>25</td><td>3.2</td><td>0.025</td></tr><tr><td>n<sub>S</sub></td><td>6.1</td><td>41</td><td>6.0</td><td>5.1</td><td>43</td><td>217</td><td>7.4</td><td>0.10</td></tr><tr><td>w<sub>0</sub></td><td>4.5</td><td>184</td><td>4.5</td><td>2.8</td><td>3.3</td><td>14</td><td>0.67</td><td>0.0094</td></tr><tr><td>b<sub>HSM</sub></td><td>4.9</td><td>71</td><td>4.9</td><td>3.5</td><td>0.98</td><td>2.0</td><td>0.11</td><td>0.0021</td></tr><tr><td>P<sub>0</sub></td><td>4.7</td><td>75</td><td>4.7</td><td>3.3</td><td>1.2</td><td>1.5</td><td>0.12</td><td>0.0015</td></tr><tr><td>α<sub>P</sub></td><td>6.0</td><td>65</td><td>6.0</td><td>4.3</td><td>3.7</td><td>6.0</td><td>0.55</td><td>0.014</td></tr><tr><td>c<sub>500,c</sub></td><td>5.9</td><td>70</td><td>5.9</td><td>3.9</td><td>8.4</td><td>11</td><td>2.0</td><td>0.018</td></tr><tr><td>α<sub>G</sub></td><td>5.2</td><td>81</td><td>5.2</td><td>3.1</td><td>1.1</td><td>2.3</td><td>0.21</td><td>0.0019</td></tr><tr><td>β<sub>G</sub></td><td>6.7</td><td>57</td><td>6.7</td><td>1.1</td><td>3.4</td><td>37</td><td>1.4</td><td>0.025</td></tr><tr><td>γ<sub>G</sub></td><td>3.9</td><td>73</td><td>3.9</td><td>2.7</td><td>1.9</td><td>0.85</td><td>0.085</td><td>0.00079</td></tr></tbody></table>

### 6.1.1 Analysis of a subset of configurations

To validate our statement that the bispectrum is effective at breaking degeneracies because of the shape-dependent response to changes in parameters, combined with the large number of available triangles, we repeat our analysis limiting ourselves to only equilateral configurations, only squeezed configurations, and a combination of the two.

If we consider conditional errors (left part of table 5), the power spectrum constraining power exceeds the bispectrum. The equilateral configurations contribute the most for the bispectrum, whereas the squeezed ones have one order of magnitude less constraining power. Upon marginalization (right part of the table), the constraining power of the squeezed and equilateral bispectrum degrades, as it does for the power spectrum, with the squeezed one degrading noticeably less. However, when we jointly analyze equilateral and squeezed limits, the degeneracies start to break and the bispectrum marginalized errors start to be competitive against the power spectrum ones. Obviously, the full bispectrum comprises even more kinds of configurations that can further ease the remaining degeneracies. We can therefore conclude that, as we claimed in section 4, the reason why the bispectrum constraining power is not as sensitive as the power spectrum one to parameter degeneracies can be explained with the presence of many more modes that are all weighted differently when the value of the parameters changes.

## 6.2 Validation against Planck noise and foreground contamination

The Planck satellite has provided a full sky map of the Compton-$\gamma$ parameter [7]. This data was used to measure the tSZ power spectrum and the bispectrum, and the former was used to constrain cosmological parameters. It is interesting to understand what are the prospective results for a full fledged analysis of the bispectrum recovered from Planck
<table><caption>**Table 6.** Conditional and marginalized error forecasts for an experiment with Planck-like noise level and foregrounds, $\ell_{max} = 1000$, and $f_{sky} = 0.47$.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Conditional</th><th colspan="3">Marginalized</th></tr><tr><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.0044</td><td>0.018</td><td>0.0034</td><td>16</td><td>0.45</td><td>0.31</td></tr><tr><td>h</td><td>0.026</td><td>0.065</td><td>0.022</td><td>27</td><td>1.9</td><td>1.7</td></tr><tr><td>σ<sub>8</sub></td><td>0.0043</td><td>0.015</td><td>0.0034</td><td>8.1</td><td>0.25</td><td>0.25</td></tr><tr><td>n<sub>S</sub></td><td>0.1</td><td>0.12</td><td>0.09</td><td>16</td><td>1.9</td><td>1.5</td></tr><tr><td>w<sub>0</sub></td><td>0.28</td><td>0.88</td><td>0.24</td><td>210</td><td>2.8</td><td>2.1</td></tr><tr><td>b<sub>HSM</sub></td><td>0.013</td><td>0.046</td><td>0.01</td><td>730</td><td>3.7</td><td>3.3</td></tr><tr><td>P<sub>0</sub></td><td>0.21</td><td>0.72</td><td>0.17</td><td>12000</td><td>62</td><td>55</td></tr><tr><td>α<sub>P</sub></td><td>0.21</td><td>0.14</td><td>0.13</td><td>150</td><td>0.62</td><td>0.59</td></tr><tr><td>c<sub>500,c</sub></td><td>0.012</td><td>0.051</td><td>0.0095</td><td>9.2</td><td>0.38</td><td>0.38</td></tr><tr><td>α<sub>G</sub></td><td>0.0087</td><td>0.031</td><td>0.0069</td><td>140</td><td>0.34</td><td>0.33</td></tr><tr><td>β<sub>G</sub></td><td>0.048</td><td>0.13</td><td>0.035</td><td>4.4</td><td>0.19</td><td>0.19</td></tr><tr><td>γ<sub>G</sub></td><td>0.025</td><td>0.069</td><td>0.02</td><td>430</td><td>1.0</td><td>0.93</td></tr></tbody></table>

<table><caption>**Table 7.** Same as table 3 for an experiment with Planck-like noise level and foregrounds, $\ell_{max} = 1000$, and $f_{sky} = 0.47$.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Gas fixed</th><th colspan="3">Planck prior</th></tr><tr><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.89</td><td>0.43</td><td>0.27</td><td>P<sub>0</sub></td><td>7400</td><td>45</td><td>38</td></tr><tr><td>h</td><td>17</td><td>1.7</td><td>1.6</td><td>α<sub>P</sub></td><td>25</td><td>0.42</td><td>0.22</td></tr><tr><td>σ<sub>8</sub></td><td>4.6</td><td>0.25</td><td>0.24</td><td>c<sub>500,c</sub></td><td>6.0</td><td>0.38</td><td>0.37</td></tr><tr><td>n<sub>S</sub></td><td>5.5</td><td>1.5</td><td>1.2</td><td>α<sub>G</sub></td><td>53</td><td>0.33</td><td>0.31</td></tr><tr><td>w<sub>0</sub></td><td>18</td><td>2.2</td><td>1.5</td><td>β<sub>G</sub></td><td>3.9</td><td>0.19</td><td>0.19</td></tr><tr><td>b<sub>HSM</sub></td><td>19</td><td>1.9</td><td>1.5</td><td>γ<sub>G</sub></td><td>220</td><td>0.97</td><td>0.88</td></tr></tbody></table>

data. To do so, we consider the Planck Gaussian noise and foreground determined in [7], and we restrict our analysis in the multipole range [70, 1000] to remove the bulk of non-Gaussian noise and remaining spurious contributions. This was performed according to what has been done in [7, 28]. As discussed, future surveys will possibly have to re-evaluate if this suffices, or if other contributions will have to be accounted for in a different manner (see [63]). We show our results in tables 6 and 7. As it is well known, the parameters cannot be meaningfully constrained separately, since they are degenerate at power spectrum level. Instead, they can be combined in a parameter that controls the overall amplitude, such as $F = \sigma_8 \Omega_m^{0.4} (1-b)^{-0.4} h^{-0.21}$ [10]. In fact, we see that after marginalization the relative errors are $\mathcal{O}(100)$, while using the full bispectrum even Planck data could lead to $\mathcal{O}(1)$ relative errors even after marginalizing over all other parameters.

Introducing the cuts in mass and redshift have a comparatively similar impact as in the previous case, as shown in table 8.

## 6.3 Forecast for future realistic surveys

To assess what are the prospects of full bispectrum analysis in the near future and in the mid-term, we use the results on section 5 in the Fisher matrix formalism to model the
<table><caption>**Table 8.** Same as table 6 (Planck-like noise and foregrounds, $\ell_{max} = 1000$, $f_{sky} = 0.47$) but with the halo redshift and mass boundaries observationally set by Planck.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Conditional</th><th colspan="3">Marginalized</th></tr><tr><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.0041</td><td>0.016</td><td>0.0028</td><td>3.8</td><td>1.7</td><td>0.27</td></tr><tr><td>h</td><td>0.027</td><td>0.063</td><td>0.021</td><td>91</td><td>5.1</td><td>5.0</td></tr><tr><td>σ<sub>8</sub></td><td>0.0042</td><td>0.015</td><td>0.003</td><td>4.7</td><td>0.64</td><td>0.4</td></tr><tr><td>n<sub>S</sub></td><td>0.18</td><td>0.15</td><td>0.15</td><td>89</td><td>8.0</td><td>3.7</td></tr><tr><td>w<sub>0</sub></td><td>0.25</td><td>0.61</td><td>0.21</td><td>47</td><td>5.7</td><td>2.4</td></tr><tr><td>b<sub>HSM</sub></td><td>0.012</td><td>0.037</td><td>0.0085</td><td>99</td><td>7.7</td><td>6.6</td></tr><tr><td>P<sub>0</sub></td><td>0.19</td><td>0.57</td><td>0.14</td><td>3600</td><td>130</td><td>110</td></tr><tr><td>α<sub>P</sub></td><td>0.33</td><td>0.15</td><td>0.13</td><td>36</td><td>1.3</td><td>0.85</td></tr><tr><td>c<sub>500,c</sub></td><td>0.011</td><td>0.041</td><td>0.0078</td><td>5.3</td><td>0.96</td><td>0.92</td></tr><tr><td>α<sub>G</sub></td><td>0.008</td><td>0.024</td><td>0.0059</td><td>51</td><td>0.85</td><td>0.78</td></tr><tr><td>β<sub>G</sub></td><td>0.043</td><td>0.22</td><td>0.029</td><td>31</td><td>4.0</td><td>3.4</td></tr><tr><td>γ<sub>G</sub></td><td>0.024</td><td>0.063</td><td>0.018</td><td>96</td><td>1.8</td><td>1.4</td></tr></tbody></table>

impact of Gaussian noise and foreground. As already discussed, the bulk of non-Gaussian contaminations can be removed enforcing a cut on the lowest multipole, which has little to no impact on the forecasted constraints.

In particular we consider two specific surveys, that are representative of what can be achieved in the next decade and in the next 30 years: Simons Observatory and the proposed Voyage 2050 Spectro-Polarimeter, respectively.

Simons Observatory (SO) [71] is an observational facility currently being built in the Atacama desert, that will be devoted to the measurement of gravitational lensing, tSZ effect, CMB temperature and polarization on very small scales. It comprises one 6 m Large Aperture Telescope (LAT) and three 0.5 m Small Aperture Telescopes. Here we just consider the LAT instrument, whose instrumental and atmospheric noise profiles are stated in tables 1 and 2 of [71]. The observation field will amount to 40% of the sky but, since we expect that some masking will be needed in any case, we conservatively set $f_{sky} = 0.30$ in the relative forecast. The results are reported in table 9. Despite the lower sky coverage, the higher raw sensitivity and beam size will allow SO to greatly improve constraints over Planck.

The Voyage 2050 Spectro-Polarimeter (V-SP) [72] is an L-class mission concept that has been put forward for the ESA call Voyage 2050. V-SP will map the whole sky, and therefore we conservatively assume it will use a mask very similar to the one employed by Planck to generate the $y$ map. Hence, when forecasting its performance we will take $f_{sky} = 0.47$. We use the noise profile reported in table 1 of [72]. Tables 11 and 12 show the results in this scenario; the performance incrementally increases over the SO forecast.

It is interesting to notice that V-SP will already be close to saturating the cosmic variance limited constraints. Since from the technological point of view we could already have the raw sensitivity to get close to the cosmic variance, this reinforces the need for developing new better methods of foreground removal, that will allow us to make use of this potential.

For a visual comparison between a CVL survey, Planck, SO and V-SP, we refer to appendix D.
<table><caption>**Table 9.** Conditional and marginalized error forecasts for an experiment with SO-like noise level and foregrounds, $\ell_{max} = 5000$, and $f_{sky} = 0.30$.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Conditional</th><th colspan="3">Marginalized</th></tr><tr><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.00082</td><td>0.0035</td><td>0.00059</td><td>7.7</td><td>0.096</td><td>0.065</td></tr><tr><td>h</td><td>0.0075</td><td>0.016</td><td>0.0063</td><td>45</td><td>0.25</td><td>0.23</td></tr><tr><td>σ<sub>8</sub></td><td>0.00076</td><td>0.0026</td><td>0.00057</td><td>2.6</td><td>0.032</td><td>0.03</td></tr><tr><td>n<sub>S</sub></td><td>0.0097</td><td>0.051</td><td>0.0058</td><td>4.3</td><td>0.21</td><td>0.16</td></tr><tr><td>w<sub>0</sub></td><td>0.0098</td><td>0.028</td><td>0.0076</td><td>90</td><td>0.31</td><td>0.29</td></tr><tr><td>b<sub>HSM</sub></td><td>0.0024</td><td>0.0087</td><td>0.0017</td><td>700</td><td>0.80</td><td>0.57</td></tr><tr><td>P<sub>0</sub></td><td>0.034</td><td>0.12</td><td>0.026</td><td>15000</td><td>10</td><td>8.3</td></tr><tr><td>α<sub>P</sub></td><td>0.0098</td><td>0.043</td><td>0.0059</td><td>8.6</td><td>0.21</td><td>0.15</td></tr><tr><td>c<sub>500,c</sub></td><td>0.0026</td><td>0.012</td><td>0.0018</td><td>6.1</td><td>0.046</td><td>0.045</td></tr><tr><td>α<sub>G</sub></td><td>0.0015</td><td>0.0052</td><td>0.0011</td><td>51</td><td>0.042</td><td>0.04</td></tr><tr><td>β<sub>G</sub></td><td>0.012</td><td>0.018</td><td>0.0073</td><td>2.1</td><td>0.022</td><td>0.022</td></tr><tr><td>γ<sub>G</sub></td><td>0.0032</td><td>0.009</td><td>0.0025</td><td>3700</td><td>0.12</td><td>0.10</td></tr></tbody></table>

<table><caption>**Table 10.** Same as table 3 for an experiment with SO-like noise level and foregrounds, $\ell_{max} = 5000$, $f_{sky} = 0.30$.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Gas fixed</th><th rowspan="2"></th><th colspan="3">Planck prior</th></tr><tr><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.18</td><td>0.067</td><td>0.036</td><td>P<sub>0</sub></td><td>3700</td><td>6.8</td><td>4.8</td></tr><tr><td>h</td><td>2.0</td><td>0.21</td><td>0.18</td><td>α<sub>P</sub></td><td>4.0</td><td>0.13</td><td>0.022</td></tr><tr><td>σ<sub>8</sub></td><td>0.63</td><td>0.031</td><td>0.028</td><td>c<sub>500,c</sub></td><td>1.2</td><td>0.041</td><td>0.04</td></tr><tr><td>n<sub>S</sub></td><td>1.1</td><td>0.2</td><td>0.12</td><td>α<sub>G</sub></td><td>22</td><td>0.039</td><td>0.037</td></tr><tr><td>w<sub>0</sub></td><td>2.7</td><td>0.16</td><td>0.087</td><td>β<sub>G</sub></td><td>0.40</td><td>0.021</td><td>0.021</td></tr><tr><td>b<sub>HSM</sub></td><td>3.2</td><td>0.24</td><td>0.17</td><td>γ<sub>G</sub></td><td>120</td><td>0.11</td><td>0.091</td></tr></tbody></table>

# 7 Discussion and conclusions

In this work, we have carried out a joint Fisher analysis of the thermal Sunyaev-Zeldovich power spectrum and bispectrum. In our analysis, we have gone beyond similar studies by significantly increasing the level of accuracy and realism of our forecasts. This has been done in several ways. First of all, we have evaluated the full power spectrum and bispectrum (auto- and cross-) covariance in the analysis, including all NG contributions. This turns out to be important, since we have found that the covariance is dominated by the connected 6-point component. Furthermore, we have considered a multi-component energy spectrum scenario and we have modeled the effects of component separation in our forecasts, via an effective noise term, which was obtained from a simple ILC procedure. Finally, rather than focusing just on $\sigma_8$ or on a small number of parameters, we have considered an extended set of both cosmological and ICM parameters, with the aim to accurately assess correlations and degeneracies and how these are dealt with by both the power spectrum and bispectrum.

In the end, we find out that the tSZ bispectrum is a very powerful observable, able to produce even stronger constraints than the tSZ power spectrum, after marginalization in a multi-parameter analysis. This is shown in our main results, summarized in figure 8 and table 4.
<table><caption>**Table 11.** Conditional and marginalized error forecasts for an experiment with V-SP-like noise level and foregrounds, $\ell_{max} = 5000$, and $f_{sky} = 0.47$.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Conditional</th><th colspan="3">Marginalized</th></tr><tr><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.00064</td><td>0.0025</td><td>0.00045</td><td>1.3</td><td>0.037</td><td>0.026</td></tr><tr><td>h</td><td>0.0059</td><td>0.012</td><td>0.0049</td><td>8.3</td><td>0.12</td><td>0.11</td></tr><tr><td>σ<sub>8</sub></td><td>0.00059</td><td>0.0019</td><td>0.00043</td><td>0.63</td><td>0.016</td><td>0.015</td></tr><tr><td>n<sub>S</sub></td><td>0.0073</td><td>0.037</td><td>0.0043</td><td>1.1</td><td>0.12</td><td>0.093</td></tr><tr><td>w<sub>0</sub></td><td>0.0076</td><td>0.021</td><td>0.0057</td><td>16</td><td>0.15</td><td>0.14</td></tr><tr><td>b<sub>HSM</sub></td><td>0.0018</td><td>0.0064</td><td>0.0013</td><td>110</td><td>0.24</td><td>0.23</td></tr><tr><td>P<sub>0</sub></td><td>0.027</td><td>0.088</td><td>0.02</td><td>2600</td><td>4.0</td><td>3.7</td></tr><tr><td>α<sub>P</sub></td><td>0.0074</td><td>0.032</td><td>0.0044</td><td>4.1</td><td>0.059</td><td>0.054</td></tr><tr><td>c<sub>500,c</sub></td><td>0.002</td><td>0.0079</td><td>0.0014</td><td>1.3</td><td>0.023</td><td>0.022</td></tr><tr><td>α<sub>G</sub></td><td>0.0012</td><td>0.0037</td><td>0.00086</td><td>11</td><td>0.022</td><td>0.021</td></tr><tr><td>β<sub>G</sub></td><td>0.0091</td><td>0.0099</td><td>0.0053</td><td>0.46</td><td>0.011</td><td>0.011</td></tr><tr><td>γ<sub>G</sub></td><td>0.0025</td><td>0.0066</td><td>0.0019</td><td>73</td><td>0.058</td><td>0.054</td></tr></tbody></table>

<table><caption>**Table 12.** Same as table 3 for an experiment with V-SP-like noise level and foregrounds, $\ell_{max} = 5000$, $f_{sky} = 0.47$.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Gas fixed</th><th rowspan="2"></th><th colspan="3">Planck prior</th></tr><tr><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th><th>PS 1σ</th><th>BS 1σ</th><th>PS⊕BS 1σ</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.11</td><td>0.033</td><td>0.023</td><td>P<sub>0</sub></td><td>1800</td><td>3.1</td><td>2.9</td></tr><tr><td>h</td><td>1.3</td><td>0.11</td><td>0.1</td><td>α<sub>P</sub></td><td>2.3</td><td>0.042</td><td>0.016</td></tr><tr><td>σ<sub>8</sub></td><td>0.38</td><td>0.015</td><td>0.015</td><td>c<sub>500,c</sub></td><td>0.57</td><td>0.022</td><td>0.022</td></tr><tr><td>n<sub>S</sub></td><td>0.66</td><td>0.1</td><td>0.069</td><td>α<sub>G</sub></td><td>10</td><td>0.021</td><td>0.02</td></tr><tr><td>w<sub>0</sub></td><td>1.6</td><td>0.073</td><td>0.051</td><td>β<sub>G</sub></td><td>0.22</td><td>0.011</td><td>0.011</td></tr><tr><td>b<sub>HSM</sub></td><td>2</td><td>0.13</td><td>0.1</td><td>γ<sub>G</sub></td><td>60</td><td>0.057</td><td>0.052</td></tr></tbody></table>

Several reasons have been already pointed out in previous studies, which explain why the bispectrum is so useful in this type of analysis. For example it has been observed that the bispectrum is less affected by uncertainties in ICM parameters, because its contributions come from (better understood) low-redshift, high-mass clusters [14]. Another important point is that the bispectrum can break degeneracies that are present in a power spectrum-only analysis, through a different amplitude scaling with parameters [16, 18].

On top of these previously known aspects, our main finding is that the bispectrum is extremely efficient at breaking degeneracies not only via amplitude scaling, but also — and mostly — due to the fact that different triangle shapes (e.g. equilateral and squeezed triangles) are affected differently by parameter changes. This is a somewhat counter-intuitive result: tSZ statistics are dominated by the one-halo term, which leads to a weak shape dependence on cosmology. Yet, we see that this is already strong enough in the bispectrum to produce large improvements in the final results. Furthermore, the triangle shape dependence on ICM parameters is of course stronger than on cosmology, significantly lowering the impact of astrophysical uncertainties in the cosmological analysis (as well as allowing for a precise measurement of astrophysical parameters themselves, possibly complementing [24]). Given
its importance, we have studied this effect in detail. To this purpose, we have isolated two specific types of triangles, namely equilateral and squeezed. We have then verified that including only one configuration type in the bispectrum analysis leads to only modest improvements in the final results, whereas including both at the same time produces significantly better forecasts; we have checked that is precisely due to a slightly different response to changes in parameters, in the squeezed and equilateral limit. To further investigate such behaviour, we have isolated the regions in the M-z plane which mostly contribute to the bispectrum in different limits, showing that such regions do not fully overlap. In other words, when we compute equilateral and squeezed bispectrum configurations, we effectively integrate the halo mass function over slightly different intervals in mass and redshift. The effect of this on parameters is fairly small for a single triangle. However, this adds up in a significant way when we produce the final forecasts, by integrating over the very large number of available configurations.

The results obtained in our study clearly suggest that a joint power spectrum-bispectrum analysis of, e.g., Planck data, or a complete likelihood study of this kind, using mock datasets for future experiments, is clearly worth pursuing. This is the object of ongoing work. Another interesting subject for future investigation consists in accounting for spectral corrections to tSZ distortions, arising from relativistic speeds of electrons in clusters [73] (the so called relativistic SZ). Such corrections are mass and redshift dependent [74, 75] through the cluster temperature scaling [76] and therefore could provide further help to break parameter degeneracies if better sensitivity is achieved. Further synergies, motivated by the recent interest in the bispectrum of other cosmological probes [53, 77–81], could be achieved considering the cross-bispectrum of tSZ with other tracers of the low-redshift Universe matter distribution.

## Acknowledgments

We thank the referee, whose insightful comments helped improving the manuscript. We also thank Boris Bolliet, Jens Chluba and Filippo Oppizzi for useful discussion, and Aditya Rotti for extensive comments on an early version of the draft. AR was supported by the ERC Consolidator Grant CMBSPEC (No. 725456) as part of the European Union’s Horizon 2020 research and innovation program. ML was supported by the project “Combining Cosmic Microwave Background and Large Scale Structure data: an Integrated Approach for Addressing Fundamental Questions in Cosmology”, funded by the MIUR Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN) Bando 2017 — grant 2017YJYZAH. ML also acknowledges support by the University of Padova under the STARS Grants programme CoGITO, Cosmology beyond Gaussianity, Inference, Theory and Observations. FL was supported by a postdoctoral grant from Centre National d’Études Spatiales (CNES). This work has made use of the Horizon Cluster hosted by Institut d’Astrophysique de Paris. We thank Stephane Rouberol for maintenance and running of this computing cluster.

# A Projection of the 3D tSZ field

## A.1 Correlation functions for the projected Compton-y field

The total Compton-y field in a given direction $n$ of the sky is defined via the following line-of-sight integral [8]

$$y(\mathbf{n}) = \int dt \, y_{3D}(\chi(t)\mathbf{n}) = \int d\chi \, a(\chi) y_{3D}(\chi(t)\mathbf{n}). \qquad (A.1)$$
The three-dimensional Compton-y field is directly related to the electron pressure profile of a single halo $P_e$ via the following re-scaling

$$y_{3D}(\mathbf{r}) = \frac{\sigma_T}{m_e c^2} P_e(\mathbf{r}), \qquad (A.2)$$

$\mathbf{r}$ being the comoving distance from the center of the halo. Under the assumption of small enough scales, we can focus on the two-dimensional Fourier transform of the field (A.1) employing the flat-sky approximation

$$\tilde{y}_{\ell} = \int d\theta_n y(\theta_n) e^{i\ell \cdot \theta_n}, \quad y(\theta_n) = \int d\chi a(\chi) y_{3D}(\chi, \theta_n \chi). \qquad (A.3)$$

In Fourier space, we define the connected part of the n-point correlation function of the two-dimensional field $\tilde{y}_\ell$

$$\langle \tilde{y}_{\ell_1} \dots \tilde{y}_{\ell_2} \rangle_c \equiv (2\pi)^2 P_n(\ell_1, \dots, \ell_n) \delta_D(\ell_1 + \dots + \ell_n). \qquad (A.4)$$

We call the quantity $P_n(\ell_1, \dots, \ell_n)$ polyspectrum of order $n$. In eq. (A.4) we described the field on a flat sky, i.e. we approximate the full spherical harmonics decomposition of the real field with a simple two-dimensional Fourier transform. By replacing eq. (A.3) into eq. (A.4), we can derive the expression for the flat-sky polyspectra. In general, the redshift integration appearing in eq. (A.3) would naturally translate into a complex n-dimensional one. To simplify this calculation, we make use of the Limber approximation [82]: we assume that the three-dimensional matter polyspectra have a weak dependence on the momenta component corresponding to the line-of-sight direction. Consequently, the projection collapses into a simple one-dimension redshift integration and we can relate the angular multipoles $\ell$ to the three-dimensional momenta $\mathbf{k}$ via the well known Limber relation $\mathbf{k}(\ell, z) \approx \{\ell/\chi(z), 0\}$. Finally, the general n-point polyspectrum $P_n(\ell_1, \dots, \ell_n)$ relates to the same order one $P_n^{y_{3D}}(\mathbf{k}_1, \dots, \mathbf{k}_n)$ for the three-dimensional Compton-y field (A.2) via

$$\begin{aligned} P_n(\ell_1, \dots, \ell_n) &= \int_0^\infty d\chi \chi^{2-2n} a^n(\chi) P_n^{y_{3D}}(\mathbf{k}(\ell_1, \chi), \dots, \mathbf{k}(\ell_n, \chi)) \\ &= \int_0^\infty dz Q^{(n)}(z) P_n^{y_{3D}}(\mathbf{k}(\ell_1, z), \dots, \mathbf{k}(\ell_n, \chi)), \end{aligned} \qquad (A.5)$$

where we define the kernel

$$Q^{(n)}(z) \equiv \chi^{2-2n}(z) a^n(z) \frac{d\chi}{dz} = \frac{\chi^{2-2n}(z)}{H(z)} a^n(z). \qquad (A.6)$$

For sake of completeness, let us recall the definition for the quantity $P_n^{y_{3D}}(\mathbf{k}_1, \dots, \mathbf{k}_n)$

$$\langle \tilde{y}_{3D}(\mathbf{k}_1) \dots \tilde{y}_{3D}(\mathbf{k}_n) \rangle_c \equiv (2\pi)^3 P_n^{y_{3D}}(\mathbf{k}_1, \dots, \mathbf{k}_n) \delta_D(\mathbf{k}_1 + \dots + \mathbf{k}_n), \qquad (A.7)$$

$$\tilde{y}_{3D}(\mathbf{k}) = \int d^3\mathbf{x} y_{3D}(\mathbf{x}) e^{+i\mathbf{k}\cdot\mathbf{x}}. \qquad (A.8)$$

As we do rely on the halo model for the matter clustering, we can write the electron pressure profile $y_{3D}(\mathbf{x}) \propto P_e(\mathbf{x})$ (A.2) (and its Fourier transform) as the sum of contributions from different halos, the mass of the i-th halo being $m_i$, centered at position $\mathbf{x}_i$

$$\tilde{y}_{3D}(\mathbf{k}) = \sum_{i}^{n^{\circ} \text{ halos}} \int dm \, d^3\mathbf{x} \, \tilde{y}_{3D}(\mathbf{k}, m) e^{+i\mathbf{x}\cdot\mathbf{k}} \delta(m - m_i) \delta(\mathbf{x} - \mathbf{x}_i). \qquad (A.9)$$
The contribution $\tilde{y}_{3D}(\mathbf{k}, m)$ from the single halo of mass $m_{500}^c$ can be obtained as

$$ 
\begin{aligned}
\tilde{y}_{3D}(\mathbf{k}, m) &= 4\pi \frac{\sigma_{\mathrm{T}}}{m_e c^2} \int dr \, r^2 j_0(kr) P_e(m, r) \\
&= 4\pi r_{500} \frac{\sigma_{\mathrm{T}}}{m_e c^2} \int dx \, x^2 j_0(kx r_{500}) P_e(m_{500}^c, x),
\end{aligned}
\quad (A.10) $$

where we employed the mass definition $m_{500}^c$ and the variable $x$ as requested for employing the parametrisation (2.11). By replacing eq. (A.9) within eq. (A.7), it is possible to obtain any desired order of correlation for the three-dimensional Compton-$y$ parameter. Furthermore, by splitting the sum in eq. (A.9) into contributions from different multi-halo configurations, we do obtain the well known hierarchy of the different halo terms. We can finally obtain the observable of interest via eq. (A.4). For consistency with the literature, we call the 2-, the 3- and the 4-point polyspectrum power spectrum, bispectrum and trispectrum, respectively:

$$ C_{\ell} \equiv P_2(\ell), \quad (A.11) $$

$$ b_{\ell_1 \ell_2 \ell_3} \equiv P_3(\ell_1, \ell_2, \ell_3), \quad (A.12) $$

$$ T(\ell_1, \ell_2, \ell_3, \ell_4) \equiv P_4(\ell_1, \ell_2, \ell_3, \ell_4). \quad (A.13) $$

We underline that the assumption of an isotropic and homogeneous Universe (Cosmological Principle) allows us to reduce the actual dependencies of the polyspectra. The power spectrum is expressed as function of the module $\ell$ of the momentum $\ell$ (we define $\ell \equiv |\ell|$) and the bispectrum has a dependence on just 3 degrees of freedom, i.e. the edges of the associated triangular configuration [83]. In the following, we will show in details the equations required for our implementation as obtained from the formalism above.

## A.2 Data vector

As far as the data vector of our analyses is concerned, we can write the power spectrum as the sum of the one- and two-halo term

$$ C_{\ell} = C_{\ell}^{1h} + C_{\ell}^{2h}, \quad (A.14) $$

$$ C_{\ell}^{1h} = \int dz \, Q^{(2)}(z) \int d\nu \, f(\nu, z) \frac{m_{\nu}}{\rho_{\mathrm{m},0}} |\tilde{y}_{3D}(k_{\ell}^{z}, m_{\nu})|^2, \quad (A.15) $$

$$ C_{\ell}^{2h}(\ell) = \int dz \, Q^{(2)}(z) P^{\mathrm{lin.}}(k_{\ell}^{z}, z) \left[ \int d\nu \, f(\nu, z) \frac{m_{\nu}}{\rho_{\mathrm{m},0}} b^{(1)}(z, m_{\nu}) \tilde{y}_{3D}(k_{\ell}^{z}, m_{\nu}) \right]^2, \quad (A.16) $$

where we introduced the abbreviation $k_{\ell}^z \equiv |\mathbf{k}(\ell, \chi(z))|$ and the mass $m_{\nu}$ is related to the parameter $\nu$ via $\sigma(m_{\nu}) = \delta/\nu$. To simplify the expression for the general matter polyspectrum, we can introduce the following quantity

$$ I_{\mu}^{\beta}(k_1, \dots, k_{\mu}; z) = \int d\nu \, f(\nu, z) \frac{m_{\nu}}{\rho_{\mathrm{m},0}} b^{(\beta)}(z, m_{\nu}) \prod_{i=1}^{\mu} \tilde{y}_{3D}(k_i, m_{\nu}), \quad (A.17) $$

where $b^{(0)} \equiv 1$ and we omitted the redshift dependence. Then, the Compton-$y$ power spectrum (A.14) can be written in a more synthetic way as

$$ C_{\ell} = \int dz \, Q^{(2)}(z) \left[ I_{2}^{0}(k_{\ell}^{z}, k_{\ell}^{z}, z) + \left[ I_{1}^{1}(k_{\ell}^{z}, z) \right]^{2} P^{\mathrm{lin.}}(k_{\ell}^{z}, z) \right]. \quad (A.18) $$
Borrowing the above notation, we can write the bispectrum used in this work as the sum of the one- and two-halo term

$$b_{\ell_1 \ell_2 \ell_3} = b_{\ell_1 \ell_2 \ell_3}^{1h} + b_{\ell_1 \ell_2 \ell_3}^{2h}, \qquad (A.19)$$

$$b_{\ell_1 \ell_2 \ell_3}^{1h} = \int dz \, Q^{(3)}(z) I_3^0(k_{\ell_1}^z, k_{\ell_2}^z, k_{\ell_3}^z; z), \qquad (A.20)$$

$$b_{\ell_1 \ell_2 \ell_3}^{2h} = \int dz \, Q^{(3)}(z) \left( I_1^1(k_{\ell_1}^z) I_2^1(k_{\ell_2}^z, k_{\ell_3}^z) P^{\text{lin.}}(k_{\ell_1}^z, z) \right. \\ \left. + I_1^1(k_{\ell_3}^z) I_2^1(k_{\ell_1}^z, k_{\ell_2}^z) P^{\text{lin.}}(k_{\ell_3}^z, z) + I_1^1(k_{\ell_2}^z) I_2^1(k_{\ell_3}^z, k_{\ell_1}^z) P^{\text{lin.}}(k_{\ell_2}^z, z) \right). \qquad (A.21)$$

As already mentioned in the main text, we do not include the three-halo term as it would require the second order halo bias for which a fit has not yet been performed.

### A.3 Covariance matrix

Moving to the covariance matrix, we employ the full expressions (A.14) and (A.19) whenever power spectra and bispectra are required. For higher order correlation functions entering the covariance computation, we rely on their respective one-halo component, which can be written in a general fashion as

$$P_n^{1h}(\ell_1, \dots, \ell_n) = \int dz \, Q^{(n)}(z) I_n^0(k_{\ell_1}^z, \dots, k_{\ell_n}^z; z). \qquad (A.22)$$

## B Three-halo term

The bispectrum three-halo term reads [28]

$$b_{\ell_1 \ell_2 \ell_3}^{3h} = \int dz \frac{d^2 V}{dz d\Omega} D_+^4(z) \int dM_1 dM_2 dM_3 \prod_{i=1}^3 \left[ \frac{dn}{dM_i}(z, M_i) \tilde{y}_{\ell_i}(z, \tilde{M}_i) \right] B_{\text{halo}}, \qquad (B.1)$$

where $B_{\text{halo}} = B_{\text{halo}}^{2PT} + B_{\text{halo}}^{b_2}$, and

$$B_{\text{halo}}^{2PT} = b(M_1) b(M_2) b(M_3) F(k_{\ell_1}, k_{\ell_2}, \theta_{12}) P_m(k_{\ell_1}) P_m(k_{\ell_2}) + 2 \text{Perm}, \qquad (B.2)$$

$$F(k_{\ell_i}, k_{\ell_j}, \theta_{ij}) = \frac{5}{7} + \frac{1}{2} \cos(\theta_{ij}) \left( \frac{k_{\ell_i}}{k_{\ell_j}} + \frac{k_{\ell_j}}{k_{\ell_i}} \right) + \frac{2}{7} \cos^2(\theta_{ij}), \\ B_{\text{halo}}^{b_2} = b(M_1) b(M_2) b_2(M_3) P_m(k_{\ell_1}) P_m(k_{\ell_2}) + 2 \text{Perm.} \qquad (B.3)$$

where the $\theta_{ij} = \theta_{ij}(k_{\ell_i}, k_{\ell_j}, k_{\ell_l})$ is the angle between the wavevectors $k_{\ell_i}$ and $k_{\ell_j}$ and we omitted the dependence on redshift, masses and wavevectors in all the halo bispectrum terms. The linear bias is defined in eq. (2.9), while the second order bias for this HMF parametrization has been derived using the peak-background split technique in [84]$^7$

$$b_2(\nu) = \frac{2\phi(42\gamma\nu^2 + 8\delta_c - 84\eta + 42\phi - 21)}{21\delta_c^2[(\beta\nu)^{2\phi} + 1]} + \frac{\gamma\nu^2(21\gamma\nu^2 + 8\delta_c - 84\eta - 63)}{21\delta_c^2} + \frac{-16\delta_c\eta - 8\delta_c + 84\eta^2 + 42\eta}{21\delta_c^2}. \qquad (B.4)$$

<footnote>$^7$Notice that the condition $\int d\nu f(\nu) b_2(\nu) = 0$ is satisfied within numerical accuracy for the parametrization they provide.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Comparison of the one-, two- and three-halo contributions to the bispectrum. The bispectrum configurations are ordered such that $\ell_1 \le \ell_2 \le \ell_3$. In the left panel crests correspond to equilateral configurations and troughs to squeezed configurations. In the right panel we multiply the bispectrum by the perimeter $P = \ell_1 + \ell_2 + \ell_3$ to the third power to flatten it.</figcaption></figure>

The resulting one-, two- and three-halo components of the bispectrum are shown in figure 9. While the two-halo term on squeezed configuration is of the same order as the one-halo term, we verified that the three-halo term is always $\le 10\%$ of the two-halo term. Therefore we do not account for it in the rest of the analysis.

## C Noise and foreground spectral shape

Here we enumerate the energy and angular dependence of the spectral components we used in the forecasts for next generation surveys. For convenience, we express the SED in terms of thermodynamic temperature rather than antenna temperature dividing them by the black-body derivative

$$ \mathcal{G}(x) = \frac{x^2 e^x}{(e^x - 1)^2}, \qquad (C.1) $$

where $x = h\nu k_B^{-1} T_{CMB}^{-1}$.

### Instrumental noise.

$$ \Theta^{\mathrm{N}}(\nu) = \Delta T_{\nu}^{2}, \quad \xi^{\mathrm{N}} \to 0, \quad C_{\ell}^{\mathrm{N}} = 8 \ln 2 e^{\ell(\ell+1)} \theta_{\nu}^{2} \theta_{\nu}^{2}, \qquad (C.2) $$

with $\theta_{\nu} = (1/\sqrt{8\ln 2})\pi/(180 \times 60)$ FWHM$_{\nu}$, where $\Delta T_{\nu}$ and FWHM$_{\nu}$ are the instrument thermodynamic-temperature-error and Full Width at Half Maximum (FWHM) measured in $T_{CMB}$ units and arcseconds respectively, for the channel $\nu$.

### CMB.

$$ \Theta^{\mathrm{CMB}}(\nu) = 1, \quad \xi^{\mathrm{CMB}} \to \infty, \quad C_{\ell}^{\mathrm{CMB}} = C_{\ell}^{\mathrm{TT}}, \qquad (C.3) $$

where $C_{\ell}^{TT}$ is the dimensionless CMB temperature power spectrum.
**Free free.**

$$ \Theta^{\text{FF}}(\nu) \propto \left\{ 1 + \ln \left[ 1 + \left( \frac{\nu_{\text{FF}}}{\nu} \right)^{\sqrt{3}/\pi} \right] \right\} \mathcal{G}^{-1}(x), \quad \nu_{\text{FF}} = 255.33 \, \text{GHz} \left( \frac{7000 \, \text{K}}{1000 \, \text{K}} \right)^{3/2}, \qquad (C.4) \\ \xi^{\text{FF}} = 0.02, \quad C_{\ell}^{\text{FF}} = \frac{(70 \, \mu\text{K})^2}{[\Theta^{\text{FF}}(31 \, \text{GHz})]^2} \ell^{-3}. $$

**Thermal dust.**

$$ \Theta^{\text{Dust}}(\nu) \propto x_D^{1.53} \frac{x^3}{e^x - 1} \mathcal{G}^{-1}(x), \quad x_D = \frac{h\nu}{k_B 21 \, \text{K}}, \qquad (C.5) \\ \xi^{\text{Dust}} = 0.3, \quad C_{\ell}^{\text{Dust}} = \frac{(24 \, \mu\text{K})^2}{[\Theta^{\text{Dust}}(90 \, \text{GHz})]^2} \ell^{-3}. $$

**Synchrotron.**

$$ \Theta^{\text{Synch}}(\nu) \propto \left( \frac{\nu}{100 \, \text{GHz}} \right)^{-0.82} \left[ 1 + \frac{1}{2} 0.2 \ln^2 \left( \frac{\nu}{100 \, \text{GHz}} \right) \right] \mathcal{G}^{-1}(x) \qquad (C.6) \\ \xi^{\text{Synch}} = 0.15, \quad C_{\ell}^{\text{Synch}} = \frac{(101 \, \mu\text{K})^2}{[\Theta^{\text{Synch}}(19 \, \text{GHz})]^2} \ell^{-2.4}. $$

**Radio point sources.**

$$ \Theta^{\text{Radio}}(\nu) \propto \nu^{-0.5} \mathcal{G}^{-1}(x), \quad \xi^{\text{Radio}} = 0.5, \quad C_{\ell}^{\text{Radio}} = \frac{(\sqrt{3} \, \mu\text{K})^2}{[\Theta^{\text{Radio}}(31 \, \text{GHz})]^2} \frac{2\pi}{\ell(\ell+1)} \frac{\ell^2}{3000^2}. \qquad (C.7) $$

**Infrared point sources.**

$$ \Theta^{\text{IR}}(\nu) \propto x_{\text{IR}}^{0.86} \frac{x^3}{e^x - 1} \mathcal{G}^{-1}(x), \quad x_{\text{IR}} = \frac{h\nu}{k_B 18.8 \, \text{K}}, \qquad (C.8) \\ \xi^{\text{IR}} = 0.3, \quad C_{\ell}^{\text{IR}} = \frac{2\pi}{\ell(\ell+1)} \frac{1}{[\Theta^{\text{IR}}(31 \, \text{GHz})]^2} \left[ (7 \, \mu\text{K})^2 \frac{\ell^2}{3000^2} + (5.7 \, \mu\text{K})^2 \frac{\ell^{2-1.2}}{3000^{2-1.2}} \right]. $$

# D Full triangle plots

For completeness we report here the triangle plot for all the parameters. In figure 10 we compare the power spectrum and bispectrum covariance ellipses for a noiseless survey with $\ell_{max} = 5000$, $f_{sky} = 1$, and no foreground contamination, assuming the full integration domain in redshift and masses from [31]. A comparison to the case in which the integration boundaries have been set according to the Planck observation is drawn in figure 11.

Figure 12 and figure 13 show the triangle plots for actual surveys, taking into account both instrumental noise and foreground contamination; Planck in figure 12, and SO and V-SP in 13. For reference, we also show there the results for the CVL survey with $\ell_{max} = 1000$ and 5000 respectively.
<figure><img src="image_11.png" /><figcaption>Figure 10. Triangle plots for a Cosmic variance limited experiment with perfect foreground separation, $l_{max} = 5000$, and $f_{sky} = 1$. Notice that the power spectrum $1\sigma$ ellipses have been rescaled by a factor to fit in the same graph. The grey bands are the power spectrum conditional errors. Figure 7 is a sub-sample of this one.</figcaption></figure>

## E Binning

As there is not a general recipe to choose a priori the best binning scheme, we validated our choice both calculating the correlation between the binned and unbinned bispectrum, and by comparing different schemes.

In general, one can evaluate the correlation between two signals, and verify if they can be distinguished from each other using a given set of data, calculating theirs scalar product using the covariance matrix as metric [85–87]

$$ \text{Corr}(b, b') \equiv \frac{1}{N} \vec{b} \cdot \text{Cov}[b, b]\vec{b}', \qquad N = \sqrt{\vec{b} \cdot \text{Cov}[b, b]\vec{b}} \sqrt{\vec{b}' \cdot \text{Cov}[b, b]\vec{b}'} . \quad (\text{E.1}) $$
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Triangle plot for the joint power spectrum and bispectrum analysis with a CVL survey with $l_{max}$ = 5000 and $f_{sky}$ = 1. The figure proposes a comparison of the two models detailed in section 6.1. We also try to fix $n_s$ and $\alpha_P$ to show the impact of the most degenerate parameters, that could be determined using external data.</figcaption></figure>

In [85–87], this method was employed in the case of the CMB primary bispectrum. That relatively simpler case allowed them to carry an exact calculation of the correlation. In our case this is not possible, and we have to resort to two approximations. First, in principle, the covariance would not depend on the tested models. In the case of the primary bispectrum, the weakly non-Gaussian limit can be employed, so that the covariance only depends on the $C_\ell$ and instrumental noise, both of which are independent of the bispectrum. For an actual survey, the covariance could instead be constructed by the collected data and, as such, it would again be independent from the theoretical model. Instead, in our case, we assume that the covariance is the one calculated using our fiducial model, and we do not vary it while
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Triangle plot for the joint power spectrum and bispectrum analysis with a CVL survey with $\ell_{max} = 1000$ and $f_{sky} = 1$, compared with Planck.</figcaption></figure>

we compare another model to the fiducial one. Second, the theoretically sound way to test
the binning would be to compare the binned bispectrum to the full bispectrum calculated on
every multipole, just repeating the representative value of each bin for every configuration
in the bin. However, this leads to the problem that the covariance matrix quickly becomes
ill-conditioned when the binning is made finer and finer. This problem was again not present
in the case of the primary bispectrum since in the weakly non-Gaussian limit the covariance
matrix is diagonal and, as such, trivial to invert. To overcome this problem, we exploit
the fact that the bispectrum is a monotonous function of $(\ell_1, \ell_2, \ell_3)$. The values of the
bispectrum in the bin that deviate the most from the bin representative will therefore be the
two calculated at the extreme bin boundaries, e.g. (100,100,100) and (150,150,150) for the
bin [100, 150] × [100, 150] × [100, 150]. Therefore, we compare the fiducial bispectrum $b$ with
<figure><img src="image_14.png" /><figcaption>Figure 13. Triangle plot for the joint power spectrum and bispectrum analysis with a CVL survey with $l_{max}$ = 5000 and $f_{sky}$ = 1, compared with SO ($f_{sky}$ = 0.30) and V-SP ($f_{sky}$ = 0.47).</figcaption></figure>

the bispectrum calculated on the higher boundary $b^{hi}$ and the lower one $b^{low}$. We obtain

$$ \text{Corr}(b, b^{\text{hi}}) = 99\%, \quad \text{Corr}(b, b^{\text{low}}) = 97\%. \qquad (\text{E.2}) $$

We have verified that changing the maximum and minimum multipoles ([10, 5000], [70, 5000], [10, 3000]), the number of bins (27, 25, 32) has little effect ($\lesssim 10\%$) on the signal to noise ratio in the case of logarithmically spaced bins. The same applies when using (according to [28]) a combination of linearly and logarithmically spaced bins (linear with $\Delta l = 64$ in [32, 1440], 10 logarithmically spaced ones in [1440, 5000]).
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Effect of parameter variations on the bispectrum kernel for two isoperimetric configurations: (10,2995,2995) and (2000,2000,2000). In the top row we assess the change due to an increase in $\Omega_m$, and in the bottom row an increase in $\beta_G$. Left and centre) contour plot of the bispectrum kernel in the redshift-mass plane. Solid lines refer to the kernel calculated with the fiducial set of parameters, dashed lines to the kernel calculated with $\Omega_m$ or $\beta_G$ increased by 10%. The two left panels display a squeezed bispectrum configuration, while the two central panels an equilateral one. The contour lines are traced at 1/10th and half maximum. Right) finite differences, as defined in eq. (F.1), of the kernel integrated over all masses.</figcaption></figure>

## F Bispectrum derivatives and breaking the degeneracies

To partially motivate the difference among the derivatives of squeezed and equilateral bispectrum, we investigated the impact of parameter changes on the bispectrum kernel. Not to have to deal with numerical convergence of derivatives in each point we consider the bispectrum kernel evaluated with the fiducial parameters, compared with the same evaluated with one parameter increased by 10%, and finite differences defined as

$$ \Delta_{\theta} \frac{\partial b_{\ell_1 \ell_2 \ell_3}}{\partial \ln z} = \frac{\partial b_{\ell_1 \ell_2 \ell_3}}{\partial \ln z}(\theta, \tilde{\theta}) - \frac{\partial b_{\ell_1 \ell_2 \ell_3}}{\partial \ln z}(1.1\theta, \tilde{\theta}) \quad (\text{F.1}) $$

where $\theta$ is the parameter being varied and $\tilde{\theta}$ is the vector of all the others which have been kept fixed. We already argued that isoperimetric bispectrum configurations probe quite similar regions in the $z-M$ plane, but here we need to focus on the subtle differences that are still present. In the left four panels of figure 14 we therefore compare the bispectrum kernel of the isoperimetric (10, 2995, 2995) and (2000, 2000, 2000) configurations, calculated for different combination of parameters. Different parameters modify the kernel in different ways (first vs second row), but also different configurations react differently (left vs central column). To see this, we show contour lines for the kernel calculated with the parameter
fiducial value (solid) and with one of the parameters ($\Omega_m$ on top, $\beta_G$ at the bottom) increased by 10%. This can be better appreciated in the right column, where we carried out the mass integral and plotted the finite differences as defined in eq. (F.1). The bispectrum derivatives (which are in a sense related to the integral of the curves plotted there) are not only different between configurations (squeezed or equilateral) but different parameters affect the various configurations differently (the ratio of the curves in the top panel is different from the ratio of the two curves in the bottom panel).

References

[1] Y.B. Zeldovich and R.A. Sunyaev, *The interaction of matter and radiation in a hot-model universe*, Astrophys. Space Sci. **4** (1969) 301.

[2] R.A. Sunyaev and Y.B. Zeldovich, *The interaction of matter and radiation in the hot model of the universe*, II, Astrophys. Space Sci. **7** (1970) 20.

[3] R.A. Sunyaev and Y.B. Zeldovich, *The Observations of relic radiation as a test of the nature of X-Ray radiation from the clusters of galaxies*, Comments Astrophys. Space Phys. **4** (1972) 173 [INSPIRE].

[4] T. Mroczkowski et al., *Astrophysics with the Spatially and Spectrally Resolved Sunyaev-Zeldovich Effects: A Millimetre/Submillimetre Probe of the Warm and Hot Universe*, Space Sci. Rev. **215** (2019) 17 [arXiv:1811.02310] [INSPIRE].

[5] E. Komatsu and U. Seljak, *The Sunyaev-Zel'dovich angular power spectrum as a probe of cosmological parameters*, Mon. Not. Roy. Astron. Soc. **336** (2002) 1256 [astro-ph/0205468] [INSPIRE].

[6] PLANCK collaboration, *Planck 2013 results. XXI. Power spectrum and high-order statistics of the Planck all-sky Compton parameter map*, Astron. Astrophys. **571** (2014) A21 [arXiv:1303.5081] [INSPIRE].

[7] PLANCK collaboration, *Planck 2015 results. XXII. A map of the thermal Sunyaev-Zeldovich effect*, Astron. Astrophys. **594** (2016) A22 [arXiv:1502.01596] [INSPIRE].

[8] J.C. Hill and E. Pajer, *Cosmology from the thermal Sunyaev-Zel'dovich power spectrum: Primordial non-Gaussianity and massive neutrinos*, Phys. Rev. D **88** (2013) 063526 [arXiv:1303.4726] [INSPIRE].

[9] M. Roncarelli, C. Carbone and L. Moscardini, *The effect of massive neutrinos on the Sunyaev-Zel'dovich and X-ray observables of galaxy clusters*, Mon. Not. Roy. Astron. Soc. **447** (2015) 1761 [arXiv:1409.4285] [INSPIRE].

[10] B. Bolliet, B. Comis, E. Komatsu and J.F. Macías-Pérez, *Dark energy constraints from the thermal Sunyaev-Zeldovich power spectrum*, Mon. Not. Roy. Astron. Soc. **477** (2018) 4957 [arXiv:1712.00788] [INSPIRE].

[11] I.G. Mccarthy, S. Bird, J. Schaye, J. Harnois-Deraps, A.S. Font and L. Van Waerbeke, *The BAHAMAS project: the CMB-large-scale structure tension and the roles of massive neutrinos and galaxy formation*, Mon. Not. Roy. Astron. Soc. **476** (2018) 2999 [arXiv:1712.02411] [INSPIRE].

[12] B. Bolliet, T. Brinckmann, J. Chluba and J. Lesgourgues, *Including massive neutrinos in thermal Sunyaev Zeldovich power spectrum and cluster counts analyses*, Mon. Not. Roy. Astron. Soc. **497** (2020) 1332 [arXiv:1906.10359] [INSPIRE].

[13] A. Cooray, W. Hu and M. Tegmark, *Large scale Sunyaev-Zel'dovich effect: Measuring statistical properties with multifrequency maps*, Astrophys. J. **540** (2000) 1 [astro-ph/0002238] [INSPIRE].
[14] S. Bhattacharya, D. Nagai, L. Shaw, T. Crawford and G.P. Holder, *Bispectrum of the Sunyaev-Zel'dovich Effect*, Astrophys. J. 760 (2012) 5 [arXiv:1203.6368] [INSPIRE].

[15] M.J. Wilson et al., *The Atacama Cosmology Telescope: A Measurement of the Thermal Sunyaev-Zel'dovich Effect Using the Skewness of the CMB Temperature Distribution*, Phys. Rev. D 86 (2012) 122005 [arXiv:1203.6633] [INSPIRE].

[16] J.C. Hill and B.D. Sherwin, *Cosmological constraints from moments of the thermal Sunyaev-Zel'dovich effect*, Phys. Rev. D 87 (2013) 023527 [arXiv:1205.5794] [INSPIRE].

[17] T.M. Crawford et al., *A Measurement of the Secondary-CMB and Millimeter-wave-foreground Bispectrum using 800 deg² of South Pole Telescope Data*, Astrophys. J. 784 (2014) 143 [arXiv:1303.3535] [INSPIRE].

[18] G. Hurier and F. Lacasa, *Combined analysis of galaxy cluster number count, thermal Sunyaev-Zel'dovich power spectrum, and bispectrum*, Astron. Astrophys. 604 (2017) A71 [arXiv:1701.09067] [INSPIRE].

[19] I. Zubeldia and A. Challinor, *Cosmological constraints from Planck galaxy clusters with CMB lensing mass bias calibration*, Mon. Not. Roy. Astron. Soc. 489 (2019) 401 [arXiv:1904.07887] [INSPIRE].

[20] R. Makiya, S. Ando and E. Komatsu, *Joint analysis of the thermal Sunyaev-Zel'dovich effect and 2MASS galaxies: probing gas physics in the local Universe and beyond*, Mon. Not. Roy. Astron. Soc. 480 (2018) 3928 [arXiv:1804.05008] [INSPIRE].

[21] R. Makiya, C. Hikage and E. Komatsu, *New constraints on the mass bias of galaxy clusters from the power spectra of the thermal Sunyaev-Zel'dovich effect and cosmic shear*, Publ. Astron. Soc. Jap. 72 (2020) 26 [arXiv:1907.07870] [INSPIRE].

[22] L. Salvati, M. Douspis, A. Ritz, N. Aghanim and A. Babul, *Mass bias evolution in tSZ cluster cosmology*, Astron. Astrophys. 626 (2019) A27 [arXiv:1901.03096] [INSPIRE].

[23] F. Ruppin, F. Mayet, J.F. Macías-Pérez and L. Perotto, *Impact of the mean pressure profile of galaxy clusters on the cosmological constraints from the Planck tSZ power spectrum*, Mon. Not. Roy. Astron. Soc. 490 (2019) 784 [arXiv:1905.05129] [INSPIRE].

[24] S. Pandey, E.J. Baxter and J.C. Hill, *Constraining the properties of gaseous halos via cross-correlations of upcoming galaxy surveys and thermal Sunyaev-Zel'dovich maps*, Phys. Rev. D 101 (2020) 043525 [arXiv:1909.00405] [INSPIRE].

[25] ACT collaboration, *The Atacama Cosmology Telescope: A Measurement of the Thermal Sunyaev-Zel'dovich One-Point PDF*, arXiv:1411.8004 [INSPIRE].

[26] L. Thiele, J.C. Hill and K.M. Smith, *Accurate analytic model for the thermal Sunyaev-Zel'dovich one-point probability distribution function*, Phys. Rev. D 99 (2019) 103511 [arXiv:1812.05584] [INSPIRE].

[27] A. Cooray and R.K. Sheth, *Halo Models of Large Scale Structure*, Phys. Rept. 372 (2002) 1 [astro-ph/0206508] [INSPIRE].

[28] F. Lacasa, *Non-Gaussianity and extragalactic foregrounds to the Cosmic Microwave Background*, Ph.D. Thesis, Orsay, IAS (2013) [arXiv:1406.0441] [INSPIRE].

[29] L. Salvati, M. Douspis and N. Aghanim, *Constraints from thermal Sunyaev-Zel'dovich cluster counts and power spectrum combined with CMB*, Astron. Astrophys. 614 (2018) A13 [arXiv:1708.00697] [INSPIRE].

[30] N. Taburet, C. Hernandez-Monteagudo, N. Aghanim, M. Douspis and R.A. Sunyaev, *The ISW-tSZ cross correlation: ISW extraction out of pure CMB data*, Mon. Not. Roy. Astron. Soc. 418 (2011) 2207 [arXiv:1012.5036] [INSPIRE].

[31] J.L. Tinker et al., *The Large Scale Bias of Dark Matter Halos: Numerical Calibration and Model Tests*, Astrophys. J. 724 (2010) 878 [arXiv:1001.3162] [INSPIRE].
[32] N. Kaiser, *On the Spatial correlations of Abell clusters*, Astrophys. J. Lett. **284** (1984) L9
[INSPIRE].

[33] J.M. Bardeen, J.R. Bond, N. Kaiser and A.S. Szalay, *The Statistics of Peaks of Gaussian Random Fields*, Astrophys. J. **304** (1986) 15 [INSPIRE].

[34] S. Cole and N. Kaiser, *Biased clustering in the cold dark matter cosmogony*, Mon. Not. Roy. Astron. Soc. **237** (1989) 1127 [INSPIRE].

[35] H.J. Mo and S.D.M. White, *An Analytic model for the spatial clustering of dark matter halos*, Mon. Not. Roy. Astron. Soc. **282** (1996) 347 [astro-ph/9512127] [INSPIRE].

[36] PLANCK collaboration, *Planck 2015 results. XIII. Cosmological parameters*, Astron. Astrophys. **594** (2016) A13 [arXiv:1502.01589] [INSPIRE].

[37] PLANCK collaboration, *Planck 2018 results. VI. Cosmological parameters*, Astron. Astrophys. **641** (2020) A6 [arXiv:1807.06209] [INSPIRE].

[38] A.J. Mead, T. Tröster, C. Heymans, L. Van Waerbeke and I.G. McCarthy, *A hydrodynamical halo model for weak-lensing cross correlations*, Astron. Astrophys. **641** (2020) A130 [arXiv:2005.00009] [INSPIRE].

[39] N. Battaglia, J.R. Bond, C. Pfrommer and J.L. Sievers, *On the Cluster Physics of Sunyaev-Zel'dovich Surveys II: Deconstructing the Thermal SZ Power Spectrum*, Astrophys. J. **758** (2012) 75 [arXiv:1109.3711] [INSPIRE].

[40] M. Arnaud, G.W. Pratt, R. Piffaretti, H. Boehringer, J.H. Croston and E. Pointecouteau, *The universal galaxy cluster pressure profile from a representative sample of nearby systems (REXCESS) and the $Y_{SZ}-M_{500}$ relation*, Astron. Astrophys. **517** (2010) A92 [arXiv:0910.1234] [INSPIRE].

[41] A.J.S. Hamilton, C.D. Rimes and R. Scoccimarro, *On measuring the covariance matrix of the nonlinear power spectrum from simulations*, Mon. Not. Roy. Astron. Soc. **371** (2006) 1188 [astro-ph/0511416] [INSPIRE].

[42] M. Takada and S. Bridle, *Probing dark energy with cluster counts and cosmic shear power spectra: including the full covariance*, New J. Phys. **9** (2007) 446 [arXiv:0705.0163] [INSPIRE].

[43] R. de Putter, C. Wagner, O. Mena, L. Verde and W. Percival, *Thinking Outside the Box: Effects of Modes Larger than the Survey on Matter Power Spectrum Covariance*, JCAP **04** (2012) 019 [arXiv:1111.6596] [INSPIRE].

[44] B. Joachimi, P. Schneider and T. Eifler, *Analysis of two-point statistics of cosmic shear*. 3. Covariances of shear measures made easy, Astron. Astrophys. **477** (2008) 43 [arXiv:0708.0387] [INSPIRE].

[45] B. Joachimi, X. Shi and P. Schneider, *Bispectrum covariance in the flat-sky limit*, Astron. Astrophys. **508** (2009) 1193 [arXiv:0907.2901] [INSPIRE].

[46] I. Kayo, M. Takada and B. Jain, *Information content of weak lensing power spectrum and bispectrum: including the non-Gaussian error covariance matrix*, Mon. Not. Roy. Astron. Soc. **429** (2013) 344 [arXiv:1207.6322] [INSPIRE].

[47] R. Scoccimarro, S. Colombi, J.N. Fry, J.A. Frieman, E. Hivon and A. Melott, *Nonlinear evolution of the bispectrum of cosmological perturbations*, Astrophys. J. **496** (1998) 586 [astro-ph/9704075] [INSPIRE].

[48] R. Scoccimarro, E. Sefusatti and M. Zaldarriaga, *Probing primordial non-Gaussianity with large-scale structure*, Phys. Rev. D **69** (2004) 103513 [astro-ph/0312286] [INSPIRE].

[49] M. Bucher, B. Racine and B. van Tent, *The binned bispectrum estimator: template-based and non-parametric CMB non-Gaussianity searches*, JCAP **05** (2016) 055 [arXiv:1509.08107] [INSPIRE].
[50] K.C. Chan, A. Moradinezhad Dizgah and J. Noreña, *Bispectrum Supersample Covariance*, *Phys. Rev. D* **97** (2018) 043532 [arXiv:1709.02473] [INSPIRE].

[51] A. Barreira, E. Krause and F. Schmidt, *Complete super-sample lensing covariance in the response approach*, *JCAP* **06** (2018) 015 [arXiv:1711.07467] [INSPIRE].

[52] A. Barreira, E. Krause and F. Schmidt, *Accurate cosmic shear errors: do we need ensembles of simulations?*, *JCAP* **10** (2018) 053 [arXiv:1807.04266] [INSPIRE].

[53] M. Rizzato, K. Benabed, F. Bernardeau and F. Lacasa, *Tomographic weak lensing bispectrum: a thorough analysis towards the next generation of galaxy surveys*, *Mon. Not. Roy. Astron. Soc.* **490** (2019) 4688 [arXiv:1812.07437] [INSPIRE].

[54] N.S. Sugiyama, S. Saito, F. Beutler and H.-J. Seo, *Perturbation theory approach to predict the covariance matrices of the galaxy power spectrum and bispectrum in redshift space*, *Mon. Not. Roy. Astron. Soc.* **497** (2020) 1684 [arXiv:1908.06234] [INSPIRE].

[55] I. Kayo and M. Takada, *Cosmological parameters from weak lensing power spectrum and bispectrum tomography: including the non-Gaussian errors*, arXiv:1306.4684 [INSPIRE].

[56] A. Barreira, *The squeezed matter bispectrum covariance with responses*, *JCAP* **03** (2019) 008 [arXiv:1901.01243] [INSPIRE].

[57] PLANCK collaboration, *Planck 2015 results. XVII. Constraints on primordial non-Gaussianity*, *Astron. Astrophys.* **594** (2016) A17 [arXiv:1502.01592] [INSPIRE].

[58] F. Lacasa, N. Aghanim, M. Kunz and M. Frommert, *Characterisation of the non-Gaussianity of radio and IR point-sources at CMB frequencies*, *Mon. Not. Roy. Astron. Soc.* **421** (2012) 1982 [arXiv:1107.2251] [INSPIRE].

[59] J. Chluba, J.C. Hill and M.H. Abitbol, *Rethinking CMB foregrounds: systematic extension of foreground parametrizations*, *Mon. Not. Roy. Astron. Soc.* **472** (2017) 1195 [arXiv:1701.00274] [INSPIRE].

[60] M.H. Abitbol, J. Chluba, J.C. Hill and B.R. Johnson, *Prospects for Measuring Cosmic Microwave Background Spectral Distortions in the Presence of Foregrounds*, *Mon. Not. Roy. Astron. Soc.* **471** (2017) 1126 [arXiv:1705.01534] [INSPIRE].

[61] A. Rotti and J. Chluba, *Combining ILC and moment expansion techniques for extracting average-sky signals and CMB anisotropies*, *Mon. Not. Roy. Astron. Soc.* **500** (2020) 976 [arXiv:2006.02458] [INSPIRE].

[62] PLANCK collaboration, *Planck 2015 results. XXVII. The Second Planck Catalogue of Sunyaev-Zeldovich Sources*, *Astron. Astrophys.* **594** (2016) A27 [arXiv:1502.01598] [INSPIRE].

[63] W.R. Coulton et al., *Non-Gaussianity of secondary anisotropies from ACTPol and Planck*, *JCAP* **09** (2018) 022 [arXiv:1711.07879] [INSPIRE].

[64] M. Tegmark, D.J. Eisenstein, W. Hu and A. de Oliveira-Costa, *Foregrounds and forecasts for the cosmic microwave background*, *Astrophys. J.* **530** (2000) 133 [astro-ph/9905257] [INSPIRE].

[65] J. Dunkley et al., *The Atacama Cosmology Telescope: likelihood for small-scale CMB data*, *JCAP* **07** (2013) 025 [arXiv:1301.0776] [INSPIRE].

[66] A. Heavens, *Statistical techniques in cosmology*, arXiv:0906.0664 [INSPIRE].

[67] J. Carron, *On the assumption of Gaussianity for cosmological two-point statistics and parameter dependent covariance matrices*, *Astron. Astrophys.* **551** (2013) A88 [arXiv:1204.4724] [INSPIRE].

[68] E. Sellentin, M. Quartin and L. Amendola, *Breaking the spell of Gaussianity: forecasting with higher order Fisher matrices*, *Mon. Not. Roy. Astron. Soc.* **441** (2014) 1831 [arXiv:1401.6892] [INSPIRE].

[69] E. Sellentin, *A fast, always positive definite and normalizable approximation of non-Gaussian likelihoods*, *Mon. Not. Roy. Astron. Soc.* **453** (2015) 893 [arXiv:1506.04866] [INSPIRE].
[70] PLANCK collaboration, *Planck Intermediate Results. V. Pressure profiles of galaxy clusters from the Sunyaev-Zeldovich effect*, Astron. Astrophys. **550** (2013) A131 [arXiv:1207.4061] [INSPIRE].

[71] SIMONS OBSERVATORY collaboration, *The Simons Observatory: Science goals and forecasts*, JCAP **02** (2019) 056 [arXiv:1808.07445] [INSPIRE].

[72] J. Delabrouille et al., *Microwave Spectro-Polarimetry of Matter and Radiation across Space and Time*, arXiv:1909.01591 [INSPIRE].

[73] J. Chluba, D. Nagai, S. Sazonov and K. Nelson, *A fast and accurate method for computing the Sunyaev-Zeldovich signal of hot galaxy clusters*, Mon. Not. Roy. Astron. Soc. **426** (2012) 510 [arXiv:1205.5778] [INSPIRE].

[74] M. Remazeilles, B. Bolliet, A. Rotti and J. Chluba, *Can we neglect relativistic temperature corrections in the Planck thermal SZ analysis?*, Mon. Not. Roy. Astron. Soc. **483** (2019) 3459 [arXiv:1809.09666] [INSPIRE].

[75] M. Remazeilles and J. Chluba, *Mapping the relativistic electron gas temperature across the sky*, Mon. Not. Roy. Astron. Soc. **494** (2020) 5734 [arXiv:1907.00916] [INSPIRE].

[76] E. Lee, J. Chluba, S.T. Kay and D.J. Barnes, *Relativistic SZ temperature scaling relations of groups and clusters derived from the BAHAMAS and MACSIS simulations*, Mon. Not. Roy. Astron. Soc. **493** (2020) 3274 [arXiv:1912.07924] [INSPIRE].

[77] V. Yankelevich and C. Porciani, *Cosmological information in the redshift-space bispectrum*, Mon. Not. Roy. Astron. Soc. **483** (2019) 2078 [arXiv:1807.07076] [INSPIRE].

[78] C. Hahn, F. Villaescusa-Navarro, E. Castorina and R. Scoccimarro, *Constraining $M_\nu$ with the bispectrum. Part I. Breaking parameter degeneracies*, JCAP **03** (2020) 040 [arXiv:1909.11107] [INSPIRE].

[79] D. Karagiannis, A. Lazanu, M. Liguori, A. Raccanelli, N. Bartolo and L. Verde, *Constraining primordial non-Gaussianity with bispectrum and power spectrum from upcoming optical and radio surveys*, Mon. Not. Roy. Astron. Soc. **478** (2018) 1341 [arXiv:1801.09280] [INSPIRE].

[80] J.C. Hill, *Foreground Biases on Primordial Non-Gaussianity Measurements from the CMB Temperature Bispectrum: Implications for Planck and Beyond*, Phys. Rev. D **98** (2018) 083542 [arXiv:1807.07324] [INSPIRE].

[81] C. Heinrich and O. Doré, *Cross-bispectra constraints on modified gravity theories from the Nancy Grace Roman Space Telescope and the Rubin Observatory Legacy Survey of Space and Time*, Phys. Rev. D **102** (2020) 123549 [arXiv:2006.03138] [INSPIRE].

[82] M. LoVerde and N. Afshordi, *Extended Limber Approximation*, Phys. Rev. D **78** (2008) 123506 [arXiv:0809.5112] [INSPIRE].

[83] O.H. Marcori and T.S. Pereira, *Two-point Correlation Functions in Inhomogeneous and Anisotropic Cosmologies*, JCAP **02** (2017) 032 [arXiv:1612.01994] [INSPIRE].

[84] A. Lazanu, T. Giannantonio, M. Schmittfull and E.P.S. Shellard, *Matter bispectrum of large-scale structure: Three-dimensional comparison between theoretical models and numerical simulations*, Phys. Rev. D **93** (2016) 083517 [arXiv:1510.04075] [INSPIRE].

[85] K.M. Smith and M. Zaldarriaga, *Algorithms for bispectra: Forecasting, optimal analysis, and simulation*, Mon. Not. Roy. Astron. Soc. **417** (2011) 2 [astro-ph/0612571] [INSPIRE].

[86] J.R. Fergusson and E.P.S. Shellard, *The shape of primordial non-Gaussianity and the CMB bispectrum*, Phys. Rev. D **80** (2009) 043510 [arXiv:0812.3413] [INSPIRE].

[87] J.R. Fergusson, M. Liguori and E.P.S. Shellard, *General CMB and Primordial Bispectrum Estimation I: Mode Expansion, Map-Making and Measures of $f_{\text{NL}}$, Phys. Rev. D **82** (2010) 023502 [arXiv:0912.5516] [INSPIRE].