{'fraction_non_alphanumeric': 0.03781560006317642, 'fraction_numerical': 0.01538797090876272, 'mean_word_length': 4.417070686494195, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 26, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cosmology from the Thermal Sunyaev-Zel'dovich Power Spectrum: Primordial non-Gaussianity and Massive Neutrinos

We carry out a comprehensive analysis of the possible constraints on cosmological and astrophysical parameters achievable with measurements of the thermal Sunyaev-Zel'dovich (tSZ) power spectrum from upcoming full-sky CMB observations, with a particular focus on one-parameter extensions to the LCDM standard model involving local primordial non-Gaussianity (described by fNL) and massive neutrinos (described by Mnu). We include all of the relevant physical effects due to these additional parameters, including the change to the halo mass function and the scale-dependent halo bias induced by local primordial non-Gaussianity. [...] We compute forecasts for Planck, PIXIE, and a cosmic variance (CV)-limited experiment, using multifrequency subtraction to remove foregrounds and implementing two masking criteria based on the ROSAT and eROSITA cluster catalogs to reduce the significant CV errors at low multipoles. We find that Planck can detect the tSZ power spectrum with>30-sigma significance, regardless of the masking scenario. However, neither Planck or PIXIE is likely to provide competitive constraints on fNL from the tSZ power spectrum due to CV noise at low-ell overwhelming the unique signature of the scale-dependent bias. A future CV-limited experiment could provide a 3-sigma detection of fNL~37, which is the WMAP9 maximum-likelihood result. The outlook for neutrino masses is more optimistic: Planck can reach levels comparable to the current upper bounds ~<0.3 eV with conservative assumptions about the intracluster medium (ICM); stronger ICM priors could allow Planck to provide 1-2-sigma evidence for massive neutrinos from the tSZ power spectrum, depending on the true value of the sum of the neutrino masses. We also forecast a ~<10% constraint on the outer slope of the ICM pressure profile using the unmasked Planck tSZ power spectrum.


I. INTRODUCTION
The thermal Sunyaev-Zel'dovich (tSZ) effect is a spectral distortion of the cosmic microwave background (CMB) that arises due to the inverse Compton scattering of CMB photons off hot electrons that lie between our vantage point and the surface of last scattering [1]. The vast majority of these hot electrons are located in the intracluster medium (ICM) of galaxy clusters, and thus the tSZ signal is dominated by contributions from these massive objects. The tSZ effect has been used for many years to study individual clusters in pointed observations (e.g., [2][3][4][5]) and in recent years has been used as a method with which to find and characterize massive clusters in blind millimeter-wave surveys [6][7][8][9][10]. Moreover, recent years have brought the first detections of the angular power spectrum of the tSZ effect through its contribution to the power spectrum in arcminute-resolution maps of the microwave sky made by the Atacama Cosmology Telescope (ACT) 1 and the South Pole Telescope (SPT) 2 [11][12][13][14]. In addition, three-point statistics of the tSZ signal have been detected within the past year: first, the real-space skewness was detected in ACT data [15] using methods first anticipated by [16], and second, the Fourier-space bispectrum was very recently detected in SPT data [17]. The amplitudes of these measurements were shown to be consistent in the SPT analysis, despite observing different regions of sky and using different analysis methods. Note that the tSZ signal is highly non-Gaussian since it is dominated by contributions from massive collapsed objects in the late-time density field; thus, higher-order tSZ statistics contain significant information beyond that found in the power spectrum. Furthermore, the combination of multiple different N -point tSZ statistics provides an avenue to extract tighter constraints on cosmological parameters and the astrophysics of the ICM than the use of the power spectrum alone, through the breaking of degeneracies between ICM and cosmological parameters [21,22]. The recent SPT bispectrum detection used such methods in order to reduce the error bar on the tSZ power spectrum amplitude by a factor of two [17].
Thus far, tSZ power spectrum detections have been limited to measurements or constraints on the power at a single multipole (typically = 3000) because ACT and SPT do not have sufficient frequency coverage to fully separate the tSZ signal from other components in the microwave sky using its unique spectral signature. However, this situation will shortly change with the imminent release of full-sky temperature maps from the Planck satellite 3 . Planck has nine frequency channels that span the spectral region around the tSZ null frequency ≈ 218 GHz. Thus, it should be possible to separate the tSZ signal from other components in the sky maps to high accuracy, allowing for a measurement of the tSZ power spectrum over a wide range of multipoles, possibly ∼ 100 < ∼ < ∼ 1500, as we demonstrate in this paper. The proposed Primordial Inflation Explorer (PIXIE) experiment [24] will also be able to detect the tSZ power spectrum at high significance, as its wide spectral coverage and high spectral resolution will allow for very accurate extraction of the tSZ signal. However, its angular resolution is much lower than Planck's, and thus the tSZ power spectrum will be measured over a much smaller range of multipoles ( < ∼ 200). However, the PIXIE data (after masking using X-ray cluster catalogs -see below) will permit tSZ measurements on large angular scales ( < ∼ 100) that are essentially inaccessible to Planck due to its noise levels; these multipoles are precisely where one would expect the signature of scale-dependent bias induced by primordial non-Gaussianity to arise in the tSZ power spectrum. Assessing the amplitude and detectability of this signature is a primary motivation for this paper.
The tSZ power spectrum has been suggested as a potential cosmological probe by a number of authors over the past few decades (e.g., [18][19][20]). Nearly all studies in the last decade have focused on the small-scale tSZ power spectrum ( > ∼ 1000) due to its role as a foreground in high angular resolution CMB measurements, and likely because without multi-frequency information, the tSZ signal has only been able to be isolated by looking for its effects on small scales (e.g., using an -space filter to upweight tSZ-dominated small angular scales in CMB maps). Much of this work was driven by the realization that the tSZ power spectrum is a very sensitive probe of the amplitude of matter density fluctuations, σ 8 [19]. The advent of multi-frequency data promises measurements of the large-scale tSZ power spectrum very shortly, and thus we believe it is timely to reassess its value as a cosmological probe, including parameters beyond σ 8 and including a realistic treatment of the uncertainties due to modeling of the ICM. We build on the work of [20] to compute the full angular power spectrum of the tSZ effect, including both the one-and two-halo terms, and moving beyond the Limber/flat-sky approximations where necessary.
Our primary interest is in assessing constraints from the tSZ power spectrum on currently unknown parameters beyond the ΛCDM standard model: the amplitude of local primordial non-Gaussianity, f NL , and the sum of the neutrino masses, M ν ≡ m ν . The values of these parameters are currently unknown, and determining their values is a key goal of modern cosmology.
Primordial non-Gaussianity is one of the few known probes of the physics of inflation. Models of single-field, minimally-coupled slow-roll inflation predict negligibly small deviations from Gaussianity in the initial curvature perturbations [25,26]. In particular, a detection of a non-zero bispectrum amplitude in the so-called "squeezed" limit (k 1 k 2 , k 3 ) would falsify essentially all single-field models of inflation [26,27]. This type of non-Gaussianity can be parametrized using the "local" model, in which f NL describes the lowest-order deviation from Gaussianity [28][29][30]: where Φ is the primordial potential and Φ G is a Gaussian field. Note that Φ = 3 5 ζ, where ζ is the initial adiabatic curvature perturbation. Local-type non-Gaussianity can be generated in multi-field inflationary scenarios, such as the curvaton model [31][32][33], or by non-inflationary models for the generation of perturbations, such as the new ekpyrotic/cyclic scenario [34][35][36]. It is perhaps most interesting when viewed as a method with which to rule out single-field inflation, however. Current constraints on f NL are consistent with zero [37,38], but the errors will shrink significantly very soon with the imminent CMB results from Planck. We review the effects of f NL = 0 on the large-scale structure of the universe in Section II.
In contrast to primordial non-Gaussianity, massive neutrinos are certain to exist at a level that will be detectable within the next decade or so; neutrino oscillation experiments have precisely measured the differences between the squared masses of the three known species, leading to a lower bound of ≈ 0.05 eV on the total summed mass [39]. The remaining questions surround their detailed properties, especially their absolute mass scale. The contribution of massive neutrinos to the total energy density of the universe today can be expressed as where M ν is the sum of the masses of the three known neutrino species. Although this contribution appears to be small, massive neutrinos can have a significant influence on the small-scale matter power spectrum. Due to their large thermal velocities, neutrinos free-stream out of gravitational potential wells on scales below their free-streaming scale [40,41]. This suppresses power on scales below the free-streaming scale. Current upper bounds from various cosmological probes assuming a flat ΛCDM+M ν model are in the range M ν < ∼ 0.3 eV [12,42,43,101], although a 3σ detection near this mass scale was recently claimed in [44]. Should the true total mass turn out to be near 0.3 eV, its effect on the tSZ power spectrum may be marginally detectable using the Planck data even for fairly conservative assumptions about the ICM physics, as we show in this paper. With stronger ICM priors, Planck could achieve a ∼ 2 − 3σ detection for masses at this scale, using only the primordial CMB temperature power spectrum and the tSZ power spectrum. We review the effects of massive neutrinos on the large-scale structure of the universe in Section II.
In addition to the effects of both known and currently unknown cosmological parameters, we also model the effects of the physics of the ICM on the tSZ power spectrum. This subject has attracted intense scrutiny in recent years after the early measurements of tSZ power from ACT [11] and SPT [45] were significantly lower than the values predicted from existing ICM pressure profile models (e.g., [122]) in combination with WMAP5 cosmological parameters. Subsequent ICM modeling efforts have ranged from fully analytic approaches (e.g., [123]) to cosmological hydrodynamics simulations (e.g., [116,117]), with other authors adopting semi-analytic approaches between these extremes (e.g., [46,118]), in which dark matter-only N -body simulations are post-processed to include baryonic physics according to various prescriptions. In addition, recent SZ and X-ray observations have continued to further constrain the ICM pressure profile from data, although these results are generally limited to fairly massive, nearby systems (e.g., [119,120]). We choose to adopt a parametrized form of the ICM pressure profile known as the Generalized NFW (GNFW) profile, with our fiducial parameter values chosen to match the constrained pressure profile fit from hydrodynamical simulations in [116]. In order to account for uncertainty in the ICM physics, we free two of the parameters in the pressure profile (the overall normalization and the outer logarithmic slope) and treat them as additional parameters in our model. This approach is discussed in detail in Section III B. We use this profile to compute the tSZ power spectrum following the halo model approach, for which we provide a complete derivation in Appendix A.
In addition to a model for the tSZ signal, we must compute the tSZ power spectrum covariance matrix in order to forecast parameter constraints. There are two important issues that must be considered in computing the expected errors or signal-to-noise ratio (SNR) for a measurement of the tSZ power spectrum. First, we must assess how well the tSZ signal can actually be separated from the other components in maps of the microwave sky, including the primordial CMB, thermal dust, point sources, and so on. Following [47] and [48], we choose to implement a multi-frequency subtraction technique that takes advantage of the unique spectral signature of the tSZ effect, and also takes advantage of the current state of knowledge about the frequency-and multipole-dependence of the foregrounds. Although there are other approaches to this problem, such as using internal linear combination techniques to construct a Compton-y map from the individual frequency maps in a given experiment (e.g., [50][51][52]), we find this method to be fairly simple and robust. We describe these calculations in detail in Section IV.
Second, we must account for the extreme cosmic variance induced in the large-angle tSZ power spectrum by massive clusters at low redshifts. The one-halo term from these objects dominates the angular trispectrum of the tSZ signal, even down to very low multipoles [59]. The trispectrum represents a large non-Gaussian contribution to the covariance matrix of the tSZ power spectrum [19], which is especially problematic at low multipoles. However, the trispectrum can be greatly suppressed by masking massive, low-redshift clusters using existing X-ray, optical, or SZ catalogs [20]. This procedure can greatly increase the SNR for the tSZ power spectrum at low multipoles. For constraints on f NL it also has the advantage of enhancing the relative importance of the two-halo term compared to the one-halo, thus showing greater sensitivity to the scale-dependent bias at low-. Moreover, even in a Gaussian cosmology, the inclusion of the two-halo term slightly changes the shape of the tSZ power spectrum, which likely helps break degeneracies amongst the several parameters which effectively only change the overall amplitude of the one-halo term; the relative enhancement of the two-halo term due to masking should help further in this regard. We consider two masking scenarios motivated by the flux limits of the cluster catalogs from all-sky surveys performed with the ROSAT 4 X-ray telescope and the upcoming eROSITA 5 X-ray telescope. These scenarios are detailed in Section V 2; by default all calculations and figures are computed for the unmasked scenario unless they are labelled otherwise.
Earlier studies have investigated the consequences of primordial non-Gaussianity for the tSZ power spectrum [53,54], though we are not aware of any calculations including the two-halo term (and hence the scale-dependent bias) or detailed parameter constraint forecasts. We are also not aware of any previous work investigating constraints on massive neutrinos from the tSZ power spectrum, although previous authors have computed their signature [55].
Other studies have investigated detailed constraints on the primary ΛCDM parameters from the combination of CMB and tSZ power spectrum measurements [57]. Many authors have investigated constraints on f NL and M ν from cluster counts, though the results depend somewhat on the cluster selection technique and mass estimation method. Considering SZ cluster count studies only, [55] and [56] investigated constraints on M ν from a Planck-derived catalog of SZ clusters (in combination with CMB temperature power spectrum data). The earlier paper found a 1σ uncertainty of ∆M ν ≈ 0.28 eV while the later paper found ∆M ν ≈ 0.06 − 0.12 eV; the authors state that the use of highly degenerate nuisance parameters degraded the results in the former study. In either case, the result is highly sensitive to uncertainties in the halo mass function, as the clusters included are deep in the exponential tail of the mass function. We expect that our results using the tSZ power spectrum should be less sensitive to uncertainties in the tail of the mass function, as the power spectrum is dominated at most angular scales by somewhat less massive objects (10 13 − 10 14 M /h) [19]. Finally, a very recent independent study [58] found ∆M ν ≈ 0.3 − 0.4 eV for Planck SZ cluster counts (with CMB temperature power spectrum information added), although they estimated that this bound could be improved to ∆M ν ≈ 0.08 eV with the inclusion of stronger priors on the ICM physics.
Our primary findings are as follows: • The tSZ power spectrum can be detected with a total SNR > 30 using the imminent Planck data up to = 3000, regardless of masking; • The tSZ power spectrum can be detected with a total SNR between ≈ 6 and 22 using the future PIXIE data up to = 300, with the result being sensitive to the level of masking applied to remove massive, nearby clusters; • Adding the tSZ power spectrum information to the forecasted constraints from the Planck CMB temperature power spectrum and existing H 0 data is unlikely to significantly improve constraints on the primary cosmological parameters, but may give interesting constraints on the extensions we consider: -If the true value of f NL is near the WMAP9 ML value of ≈ 37, a future CV-limited experiment combined with eROSITA-masking could provide a 3σ detection, completely independent of the primordial CMB temperature bispectrum; alternatively, PIXIE could give 1 − 2σ evidence for such a value of f NL with this level of masking; -If the true value of M ν is near 0.1 eV, the Planck tSZ power spectrum with eROSITA masking can provide upper limits competitive with the current upper bounds on M ν ; with stronger external constraints on the ICM physics, Planck with eROSITA masking could provide 1 − 2σ evidence for massive neutrinos from the tSZ power spectrum, depending on the true neutrino mass; • Regardless of the cosmological constraints, Planck will allow for a very tight constraint on the logarithmic slope of the ICM pressure profile in the outskirts of galaxy clusters, and may also provide some information on the overall normalization of the pressure profile (which sets the zero point of the Y − M relation).
The remainder of this paper is organized as follows. In Section II, we describe our models for the halo mass function and halo bias, as well as the effects of primordial non-Gaussianity and massive neutrinos on large-scale structure. In Section III, we describe our halo model-based calculation of the tSZ power spectrum, including the relevant ICM physics. We also demonstrate the different effects of each parameter in our model on the tSZ power spectrum. In Section IV, we consider the extraction of the tSZ power spectrum from the other components in microwave sky maps via multifrequency subtraction techniques. Having determined the experimental noise levels, in Section V we detail our calculation of the covariance matrix of the tSZ power spectrum, and discuss the role of masking massive nearby clusters in reducing the low-cosmic variance. In Section VI, we use our tSZ results to forecast constraints on cosmological and astrophysical parameters from a variety of experimental set-ups and masking choices. We also compute the expected SNR of the tSZ power spectrum detection for each possible scenario. We discuss our results and conclude in Section VII. Finally, in Appendix B, we provide a brief comparison between our forecasts and the Planck tSZ power spectrum results that were publicly released while this manuscript was under review [137].
The WMAP9+eCMB+BAO+H 0 maximum-likelihood parameters [69] define our fiducial model (see Section III C for details). All masses are quoted in units of M /h, where h ≡ H 0 /(100 km s −1 Mpc −1 ) and H 0 is the Hubble parameter today. All distances and wavenumbers are in comoving units of Mpc/h. All tSZ observables are computed at ν = 150 GHz, since ACT and SPT have observed the tSZ signal at (or very near) this frequency, where the tSZ effect leads to a temperature decrement in the CMB along the line-of-sight (LOS) to a galaxy cluster.

II. MODELING LARGE-SCALE STRUCTURE
In order to compute statistics of the tSZ signal, we need to model the comoving number density of halos as a function of mass and redshift (the halo mass function) and the bias of halos with respect to the underlying matter density field as a function of mass and redshift. Moreover, in order to extract constraints on f NL and M ν from the tSZ power spectrum, we must include the effects of these parameters on large-scale structure. We describe our approach to these computations in the following.

A. Halo Mass Function
We define the mass of a dark matter halo by the spherical overdensity (SO) criterion: M δ,c (M δ,d ) is the mass enclosed within a sphere of radius r δ,c (r δ,d ) such that the enclosed density is δ times the critical (mean matter) density at redshift z. To be clear, c subscripts refer to masses referenced to the critical density at redshift z, ρ cr (z) = 3H 2 (z)/8πG with H(z) the Hubble parameter at redshift z, whereas d subscripts refer to masses referenced to the mean matter density at redshift z,ρ m (z) ≡ρ m (this quantity is constant in comoving units).
We will generally work in terms of a particular SO mass, the virial mass, which we denote as M . The virial mass is the mass enclosed within a radius r vir [70]: where ∆ cr (z) = 18π 2 +82(Ω(z)−1)−39(Ω(z)−1) 2 and Ω(z) = Ω m (1+z) 3 /(Ω m (1+z) 3 +Ω Λ ). For many calculations, we need to convert between M and various other SO masses (e.g., M 200c or M 200d ). We use the NFW density profile [66] and the concentration-mass relation from [67] in order to do these conversions, which require solving the following non-linear equation for r δ,c (or r δ,d ): where c vir ≡ r vir /r N F W is the concentration parameter (r N F W is the NFW scale radius) and we replace the critical density ρ cr (z) with the mean matter densityρ m in this equation in order to obtain r δ,d instead of r δ,c . After solving Eq. (4) to find r δ,c , we calculate M δ,c via M δ,c = 4 3 πr 3 δ,c ρ cr (z)δ. The halo mass function, dn(M, z)/dM describes the comoving number density of halos per unit mass as a function of redshift. We employ the approach developed from early work by Press and Schechter [71] and subsequently refined by many other authors (e.g., [72][73][74][75]): where σ 2 (M, z) is the variance of the linear matter density field smoothed with a (real space) top-hat filter on a scale at redshift z: where P lin (k, z) is the linear theory matter power spectrum at wavenumber k and redshift z. Note that the window function W (k, R) is a top-hat filter in real space, which in Fourier space is given by where x ≡ kR. In Eq. (5), the function f (σ(M, z)) is known as the halo multiplicity function. It has been measured to increasingly high precision from large N -body simulations over the past decade [74][75][76][77]. However, many of these calibrated mass functions are specified in terms of the friends-of-friends (FOF) mass rather than the SO mass, hindering their use in analytic calculations such as ours. For this reason, we use the parametrization and calibration from [75], where computations are performed in terms of the SO mass with respect to the mean matter density, M δ,d , for a variety of overdensities. The halo multiplicity function in this model is parametrized by where {A, a, b, c} are (redshift-and overdensity-dependent) parameters fit from simulations. We use the values of these parameters appropriate for the M 200,d halo mass function from [75] with the redshift-dependent parameters given in their Eqs. (5)-(8); we will hereafter refer to this as the Tinker mass function. Note that the authors of that study caution against extrapolating their parameters beyond the highest redshift measured in their simulations (z = 2.5) and recommend setting the parameters equal to their z = 2.5 values at higher redshifts; we adopt this recommendation in our calculations. Also, note that our tSZ power spectrum calculations in Section III are phrased in terms of the virial mass M , and thus we compute the Jacobian dM 200,d /dM using the procedure described in Eq. (4)  dM . We compute the smoothed matter density field in Eq. (6) by first obtaining the linear theory matter power spectrum from CAMB 6 at z in = 30 and subsequently rescaling σ 2 (M, z) by D 2 (z), where D(z) is the linear growth factor. We normalize D(z) by requiring that D(z) → 1/(1 + z) deep in the matter-dominated era (e.g., at z in ). The resulting σ 2 (M, z) is then used to compute the mass function in Eq. (5).
Note that we assume the mass function to be known to high enough precision that the parameters describing it can be fixed; in other words, we do not consider {A, a, b, c} to be free parameters in our model. These parameters are certainly better constrained at present than those describing the ICM pressure profile (see Section III B), and thus this assumption seems reasonable for now. However, precision cosmological constraints based on the mass function should in principle consider variations in the mass function parameters in order to obtain robust results, as has been done in some recent X-ray cluster cosmology analyses [78]. However, we leave the implications of these uncertainties for tSZ statistics as a topic for future work.

Effect of Primordial non-Gaussianity
The influence of primordial non-Gaussianity on the halo mass function has been studied by many authors over the past two decades using a variety of approaches (e.g., [79][80][81][82][83][84][85][86][87]). The physical consequences of the model specified in Eq. (1) are fairly simple to understand for the halo mass function, especially in the exponential tail of the mass function where massive clusters are found. Intuitively, the number of clusters provides information about the tail of the probability distribution function of the primordial density field, since these are the rarest objects in the universe, which have only collapsed recently. For positive skewness in the primordial density field (f NL > 0), one obtains an increased number of massive clusters at late times relative to the f NL = 0 case, because more regions of the smoothed density field have δ > δ c , the collapse threshold (≈ 1.686 in the spherical collapse model). Conversely, for negative skewness in the primordial density field (f NL < 0), one obtains fewer massive clusters at late times relative to the f NL = 0 case, because fewer regions of the smoothed density field are above the collapse threshold. As illustrated in recent analytic calculations and simulation measurements [88][89][90][91]111], these changes can be quite significant for the number of extremely massive halos (∼ 10 15 M /h) in the late-time universe; for example, the z = 0 abundance of such halos for f NL ≈ 250 can be ≈ 1.5 − 2 times larger than in a Gaussian cosmology. These results have been used as a basis for recent studies constraining f NL by looking for extremely massive outliers in the cluster distribution (e.g., [92][93][94][95][96][97]).
We model the effect of f NL on the halo mass function by multiplying the Tinker mass function by a non-Gaussian correction factor: where dn/dM is given by Eq. (5). We use the model for R N G (M, z, f NL ) given by Eq. (35) in [88] (the "log-Edgeworth" mass function). In this approach, the density field is approximated via an Edgeworth expansion, which captures small deviations from Gaussianity. The Press-Schechter approach is then applied to the Edgeworth-expanded density field to obtain an expression for the halo mass function in terms of cumulants of the non-Gaussian density field. The results of [88] include numerical fitting functions for these cumulants obtained from N -body simulations. We use both the expression for R N G (M, z, f NL ) and the cumulant fitting functions from [88] to compute the non-Gaussian correction to the mass function. This prescription was shown to accurately reproduce the non-Gaussian halo mass function correction factor measured directly from N -body simulations in [88], and in particular improves upon the similar prescription derived in [86] (the "Edgeworth" mass function).
Note that we apply the non-Gaussian correction factor R N G (M, z, f NL ) to the Tinker mass function in Eq. (9), which is an SO mass function, as mentioned above. The prescription for computing R N G (M, z, f NL ) makes no assumption about whether M is an FOF or SO mass, so there is no logical flaw in this procedure. However, the comparisons to N -body results in [88] were performed using FOF halos. Thus, without having tested the results of Eq. (9) on SO mass functions from simulations, our calculation assumes that the change in the mass function due to non-Gaussianity is quasi-universal, even if the underlying Gaussian mass function itself is not. This assumption was tested in [91] for the non-Gaussian correction factor from [86] (see Fig. 9 in [91]) and found to be valid; thus, we choose to adopt it here. We will refer to the non-Gaussian mass function computed via Eq. (9) using the prescription from [88] as the LVS mass function.

Effect of Massive Neutrinos
It has long been known that massive neutrinos suppress the amplitude of the matter power spectrum on scales below their free-streaming scale, k f s [41]: Neutrinos do not cluster on scales much smaller than this scale (i.e., k > k f s ), as they are able to free-stream out of small-scale gravitational potential wells. This effect leads to a characteristic decrease in the small-scale matter power spectrum of order ∆P/P ≈ −8Ω ν /Ω m in linear perturbation theory [40,41]. Nonlinear corrections increase this suppression to ∆P/P ≈ −10Ω ν /Ω m for modes with wavenumbers k ∼ 0.5 − 1 Mpc/h [40]. The neutrino-induced suppression of the small-scale matter power spectrum leads one to expect that the number of massive halos in the low-redshift universe should also be decreased. Several papers in recent years have attempted to precisely model this change in the halo mass function using both N -body simulations and analytic theory [98][99][100]. In [98], N -body simulations are used to show that massive neutrinos do indeed suppress the halo mass function, especially for the largest, latest-forming halos (i.e., galaxy clusters). Moreover, the suppression is found to arise primarily from the suppression of the initial transfer function in the linear regime, and not due to neutrino clustering effects in the N -body simulations. This finding suggests that an analytic approach similar to the Press-Schecter theory should work for massive neutrino cosmologies as well, and the authors subsequently show that a modified Sheth-Tormen formalism [72] gives a good fit to their simulation results. Similar N -body simulations are examined in [99], who find generally similar results to those in [98], but also point out that the effect of M ν > 0 on the mass function cannot be adequately represented by simply rescaling σ 8 to a lower value in an analytic calculation without massive neutrinos. Finally, [100] study the effect of massive neutrinos on the mass function using analytic calculations with the spherical collapse model. Their results suggest that an accurate approximation is to simply input the M νsuppressed linear theory (cold+baryonic-only) matter power spectrum computed at z in to a ΛCDM-calibrated mass function fit (note that a similar procedure was used in some recent X-ray cluster-based constraints on M ν [101]). The net result of this suppression can be quite significant at the high-mass end of the mass function; for example, M ν = 0.1 eV leads to a factor of ∼ 2 decrease in the abundance of 10 15 M /h halos at z = 1 as compared to a massless-neutrino cosmology [100]. We follow the procedure used in [100] in our work, although we input the suppressed linear theory matter power spectrum to the Tinker mass function rather than that of [77], as was done in [100]. We will refer to the M ν -suppressed mass function computed with this prescription as the IT mass function.

B. Halo Bias
Dark matter halos are known to cluster more strongly than the underlying matter density field; they are thus biased tracers. This bias can depend on scale, mass, and redshift (e.g., [102][103][104]). We define the halo bias b(k, M, z) by where P hh (k, M, z) is the power spectrum of the halo density field and P (k, z) is the power spectrum of the matter density field. Knowledge of the halo bias is necessary to model and extract cosmological information from the clustering of galaxies and galaxy clusters. For our purposes, it will be needed to compute the two-halo term in the tSZ power spectrum, which requires knowledge of P hh (k, M, z).
In a Gaussian cosmology, the halo bias depends on mass and redshift but is independent of scale for k < ∼ 0.05 Mpc/h, i.e. on large scales (e.g., [105]). We compute this linear Gaussian bias, b G (M, z), using the fitting function in Eq. (6) of [105] with the parameters appropriate for M 200,d SO masses (see Table 2 in [105]). This fit was determined from the results of many large-volume N -body simulations with a variety of cosmological parameters and found to be quite accurate. We will refer to this prescription as the Tinker bias model.
Although the bias becomes scale-dependent on small scales even in a Gaussian cosmology, it becomes scaledependent on large scales in the presence of local primordial non-Gaussianity, as first shown in [106]. The scaledependence arises due to the coupling of long-and short-wavelength density fluctuations induced by local f NL = 0. We model this effect as a correction to the Gaussian bias described in the preceding paragraph: where the non-Gaussian correction is given by [106] ∆b Here, δ c = 1.686 (the spherical collapse threshold) and relates the linear density field to the primordial potential via δ(k, z) = α(k, z)Φ(k). Note that T (k) is the linear matter transfer function, which we compute using CAMB. Since the original derivation in [106], the results in Eqs. (13) and (14) have subsequently been confirmed by other authors [107][108][109] and tested extensively on N -body simulations (e.g., [90,106,110,112]). The overall effect is a steep increase in the large-scale bias of massive halos, which is even larger for highly biased tracers like galaxy clusters. We will refer to this effect simply as the scale-dependent halo bias. The influence of massive neutrinos on the halo bias has been studied far less thoroughly than that of primordial non-Gaussianity. Recent N -body simulations analyzed in [99] indicate that massive neutrinos lead to a nearly scaleindependent increase in the large-scale halo bias. This effect arises because of the mass function suppression discussed in Section II A 2: halos of a given mass are rarer in an M ν > 0 cosmology than in a massless neutrino cosmology (for fixed A s ), and thus they are more highly biased relative to the matter density field. However, the amplitude of this change is far smaller than that induced by f NL = 0, especially on very large scales. For example, the results of [99] indicate an overall increase of ∼ 10% in the mean bias of massive halos at z = 1 for M ν = 0.3 eV as compared to M ν = 0. Our implementation of the scale-dependent bias due to local f NL yields a factor of ∼ 100 − 1000 increase in the large-scale (k ∼ 10 −4 h/Mpc) bias of objects in the same mass range at z = 1 for f NL = 50. Clearly, the effect of f NL is much larger than that of massive neutrinos, simply because it is so strongly scale-dependent, while M ν only leads to a small scale-independent change (at least on large scales; the small-scale behavior may be more complicated). Moreover, the change in bias due to M ν is larger at higher redshifts (z > ∼ 1), whereas most of the tSZ signal originates at lower redshifts. Lastly, due to the smallness of the two-halo term in the tSZ power spectrum compared to the one-halo term (see Section III), small variations in the Gaussian bias cause essentially no change in the total signal. For all of these reasons, we choose to neglect the effect of massive neutrinos on the halo bias in our calculations.

III. THERMAL SZ POWER SPECTRUM
The tSZ effect results in a frequency-dependent shift in the CMB temperature observed in the direction of a galaxy group or cluster. The temperature shift ∆T at angular position θ with respect to the center of a cluster of mass M at redshift z is given by [1] ∆T ( θ, M, z) T CMB = g ν y( θ, M, z) where g ν = x coth(x/2)−4 is the tSZ spectral function with x ≡ hν/k B T CMB , y is the Compton-y parameter, σ T is the Thomson scattering cross-section, m e is the electron mass, and P e ( r) is the ICM electron pressure at location r with respect to the cluster center. We have neglected relativistic corrections in Eq. (15) (e.g., [113]), as these effects are relevant only for the most massive clusters in the universe ( > ∼ 10 15 M /h). Such clusters contribute non-negligibly to the tSZ power spectrum at low-, and thus our results in unmasked calculations may be slightly inaccurate; however, the optimal forecasts for cosmological constraints arise from calculations in which such nearby, massive clusters are masked (see Section VI), and thus these corrections will not be relevant. Therefore, we do not include them in our calculations.
Note that we only consider spherically symmetric pressure profiles in this work, i.e. P e ( r) = P e (r) in Eq. (15). The integral in Eq. (15) is computed along the LOS such that r 2 = l 2 + d A (z) 2 θ 2 , where d A (z) is the angular diameter distance to redshift z and θ ≡ | θ| is the angular distance between θ and the cluster center in the plane of the sky (note that this formalism assumes the flat-sky approximation is valid; we provide exact full-sky results for the tSZ power spectrum in Appendix A). In the flat-sky limit, a spherically symmetric pressure profile implies that the temperature decrement (or Compton-y) profile is azimuthally symmetric in the plane of the sky, i.e., ∆T ( θ, M, z) = ∆T (θ, M, z). Finally, note that the electron pressure P e ( r) is related to the thermal gas pressure via P th = P e (5X H +3)/2(X H +1) = 1.932P e , where X H = 0.76 is the primordial hydrogen mass fraction. We calculate all tSZ power spectra in this paper at ν = 150 GHz, where the tSZ effect is observed as a decrement in the CMB temperature (g 150 GHz = −0.9537). We make this choice simply because recent tSZ measurements have been performed at this frequency using ACT and SPT (e.g., [12,14,15,17]), and thus the temperature values in this regime are perhaps more familiar and intuitive. All of our calculations can be phrased in a frequency-independent manner in terms of the Compton-y parameter, and we will often use "y" as a label for tSZ quantities, although they are calculated numerically at ν = 150 GHz.
In the remainder of this section, we outline the halo model-based calculations used to compute the tSZ power spectrum, discuss our model for the gas physics of the ICM, and explain the physical effects of each cosmological and astrophysical parameter on the tSZ power spectrum.

A. Halo Model Formalism
We compute the tSZ power spectrum using the halo model approach (see [114] for a review). We provide complete derivations of all the relevant expressions in Appendix A, first obtaining completely general full-sky results and then specializing to the flat-sky/Limber-approximated case. Here, we simply quote the necessary results and refer the interested reader to Appendix A for the derivations. Note that we will work in terms of the Compton-y parameter; the results can easily be multiplied by the necessary g ν factors to obtain results at any frequency.
The tSZ power spectrum, C y , is given by the sum of the one-halo and two-halo terms: The exact expression for the one-halo term is given by Eq. (A17): where χ(z) is the comoving distance to redshift z, d 2 V /dzdΩ is the comoving volume element per steradian, dn/dM is the halo mass function discussed in Section II A,ỹ 3D (k; M, z) is given in Eq. (A9), and J +1/2 (x) is a Bessel function of the first kind. In the flat-sky limit, the one-halo term simplifies to the following widely-used expression (given in e.g. Eq. (1) of [19]), which we derive in Eq. (A22): whereỹ Here, r s is a characteristic scale radius (not the NFW scale radius) of the y 3D profile given by y 3D ( r) = σ T mec 2 P e ( r) and s = a(z)χ(z)/r s = d A (z)/r s is the multipole moment associated with the scale radius. For the pressure profile from [116] used in our calculations, the natural scale radius is r 200,c . In our calculations, we choose to implement the flat-sky result for the one-halo term at all -see Appendix A for a justification of this decision and an assessment of the associated error at low-(the only regime where this correction would be relevant).
The exact expression for the two-halo term is given by Eq. (A23): where P lin (k, z in ) is the linear theory matter power spectrum at z in (which we choose to set equal to 30), b(k, M, z) is the halo bias discussed in Section II B, andỹ kχ(z) (M, z) refers to the expression forỹ (M, z) given in Eq. (19) evaluated with + 1/2 = kχ. This notation is simply a mathematical convenience; no flat-sky or Limber approximation was used in deriving Eq. (A23), and no appears inỹ kχ (M, z). In the Limber approximation [115], the two-halo term simplifies to the result given in [20], which we derive in Eq. (A25): We investigate the validity of the Limber approximation in detail in Appendix A. We find that it is necessary to compute the exact expression in Eq. (20) in order to obtain sufficiently accurate results at low-, where the signature of the scale-dependent bias induced by f NL is present (looking for this signature is our primary motivation for computing the two-halo term to begin with). In particular, we compute the exact expression in Eq. (20) for < 50, while we use the Limber-approximated result in Eq. (21) at higher multipoles. The fiducial integration limits in our calculations are 0.005 < z < 4 for all redshift integrals, 5 × 10 11 M /h < M < 5 × 10 15 M /h for all mass integrals, and 10 −4 h/Mpc < k < 3 h/Mpc for all wavenumber integrals. We check that extending the wavenumber upper limit further into the nonlinear regime does not affect our results. Note that the upper limit in the mass integral becomes redshift-dependent in the masked calculations that we discuss below, in which the most massive clusters at low redshifts are removed from the computation.
We use the halo mass functions discussed in Section II A (Tinker, LVS, and IT) and the bias models discussed in Section II B (Tinker and scale-dependent bias) in Eqs. (18), (20), and (21). The only remaining ingredient needed to complete the tSZ power spectrum calculation is a prescription for the ICM electron pressure profile as a function of mass and redshift. Note that this approach to the tSZ power spectrum calculation separates the cosmology-dependent component (the mass function and bias) from the ICM-dependent component (the pressure profile). This separation arises from the fact that the small-scale baryonic physics that determines the structure of the ICM pressure profile effectively decouples from the large-scale physics described by the background cosmology and linear perturbation theory. Thus, it is a standard procedure to constrain the ICM pressure profile from cosmological hydrodynamics simulations (e.g., [116,118]) or actual observations of galaxy clusters (e.g., [119,120], which are obtained for a fixed cosmology in either case (at present, it is prohibitively computationally expensive to run many large hydrodynamical simulations with varying cosmological parameters). Of course, it is also possible to model the ICM analytically and obtain a pressure profile (e.g., [122,123]. Regardless of its origin (observations/simulations/theory), the derived ICM pressure profile can then be applied to different background cosmologies by using the halo mass function and bias model appropriate for that cosmology in the tSZ power spectrum calculations. We follow this approach.
Note that because the tSZ signal is heavily dominated by contributions from collapsed objects, the halo model approximation gives very accurate results when compared to direct LOS integrations of numerical simulation boxes (see Figs. 7 and 8 in [116] for direct comparisons). In particular, the halo model agrees very well with the simulation results for < ∼ 1000, which is predominantly the regime we are interested in for this paper (on smaller angular scales effects due to asphericity and substructure become important, which are not captured in the halo model approach). These results imply that contributions from the intergalactic medium, filaments, and other diffuse structures are unlikely to be large enough to significantly impact the calculations and forecasts in the remainder of the paper. Contamination from the Galaxy is a separate issue, which we assume can be minimized to a sufficient level through sky cuts and foreground subtraction (see Section IV).

B. Modeling the ICM
We adopt the parametrized ICM pressure profile fit from [116] as our fiducial model. This profile is derived from cosmological hydrodynamics simulations described in [117]. These simulations include (sub-grid) prescriptions for radiative cooling, star formation, supernova feedback, and feedback from active galactic nuclei (AGN). Taken together, these feedback processes typically decrease the gas fraction in low-mass groups and clusters, as the injection of energy into the ICM blows gas out of the cluster potential. In addition, the smoothed particle hydrodynamics used in these simulations naturally captures the effects of non-thermal pressure support due to bulk motions and turbulence, which must be modeled in order to accurately characterize the cluster pressure profile in the outskirts. . fNL values of this magnitude are highly disfavored by current constraints, but we plot them here to clarify the influence of primordial non-Gaussianity on the tSZ power spectrum. The effect of fNL on the one-halo term is simply an overall amplitude shift due to the corresponding increase or decrease in the number of massive clusters in the universe, as described in Section II A 1. The effect of fNL on the two-halo term includes not only an amplitude shift due to the change in the mass function, but also a steep upturn at low-due to the influence of the scale-dependent halo bias, as described in Section II B. Note that for fNL < 0 the two effects cancel for ≈ 4−5. The relative smallness of the two-halo term (compared to the one-halo term) makes the scale-dependent bias signature subdominant for all values except < ∼ 7 − 8. However, masking of nearby massive clusters suppresses the low-one-halo term in the power spectrum (in addition to decreasing the cosmic variance, as discussed in Section V), which increases the relative importance of the two-halo term and thus the dependence of the total signal on fNL at low-.
The ICM thermal pressure profile in this model is parametrized by a dimensionless GNFW form, which has been found to be a useful parametrization by many observational and numerical studies (e.g., [3,119,120,125]): where P th (x) = 1.932P e (x) is the thermal pressure profile, x is the dimensionless distance from the cluster center, x c is a core scale length, P 0 is a dimensionless amplitude, α, β, and γ describe the logarithmic slope of the profile at intermediate (x ∼ x c ), large (x x c ), and small (x x c ) radii, respectively, and P 200,c is the self-similar amplitude for pressure at r 200,c given by [126,127]: In [116] this parametrization is fit to the stacked pressure profiles of clusters extracted from the simulations described above. Note that due to degeneracies the parameters α and γ are not varied in the fit; they are fixed to α = 1.0 and . Mν values of this magnitude are at the lower bound allowed by neutrino oscillation measurements, and thus an effect of at least this magnitude is expected in our universe. The effect of Mν on both the one-and two-halo terms is effectively an overall amplitude shift, although the effect tapers off very slightly at high-. It may be puzzling at first to see an increase in the tSZ power when Mν > 0, but the key fact is that we are holding σ8 constant when varying Mν (indeed, we hold all of the other parameters constant). In order to keep σ8 fixed despite the late-time suppression of structure growth due to Mν > 0, we must increase the primordial amplitude of scalar perturbations, As. Thus, the net effect is an increase in the tSZ power spectrum amplitude, counterintuitive though it may be. If our ΛCDM parameter set included As rather than σ8, and we held As constant while Mν > 0, we would indeed find a corresponding decrease in the tSZ power (and in σ8, of course). γ = −0.3, which agree with many other studies (e.g., [3,119,120,125]. In addition to constraining the amplitude of the remaining parameters, [116] also fit power-law mass and redshift dependences, with the following results: Note that the denominator of the mass-dependent factor has units of M rather than M /h as used elsewhere in this paper. The mass and redshift dependence of these parameters captures deviations from simple self-similar cluster pressure profiles. These deviations arise from non-gravitational energy injections due to AGN and supernova feedback, star formation in the ICM, and non-thermal processes such as turbulence and bulk motions [116,121]. Eqs. (22)- (26) completely specify the ICM electron pressure profile as a function of mass and redshift, and provide the remaining ingredient needed for the halo model calculations of the tSZ power spectrum described in Section III A, in addition to the halo mass function and halo bias. We will refer to this model of the ICM pressure profile as the Battaglia model. Although it is derived solely from numerical simulations, we note that the Battaglia pressure profile is in good agreement with a number of observations of cluster pressure profiles, including those based on the REXCESS X-ray sample of massive, z < 0.3 clusters [119], independent studies of low-mass groups at z < 0.12 with Chandra [128], and early Planck measurements of the stacked pressure profile of z < 0.5 clusters [120].
We allow for a realistic degree of uncertainty in the ICM pressure profile by freeing the amplitude of the parameters that describe the overall normalization (P 0 ) and the outer logarithmic slope (β). To be clear, we do not free the mass and redshift dependences for these parameters given in Eqs. (24) and (26), only the overall amplitudes in those expressions. The outer slope β is known to be highly degenerate with the scale radius x c (e.g., [3,116]), and thus it is only feasible to free one of these parameters. The other slope parameters in Eq. (22) are fixed to their Battaglia values, which match the standard values in the literature. We parametrize the freedom in P 0 and β by introducing new parameters C P0 and C β defined by: These parameters thus describe multiplicative overall changes to the amplitudes of the P 0 and β parameters. The fiducial Battaglia profile corresponds to {C P0 , C β } = {1, 1}. We discuss our priors for these parameters in Section VI.

C. Parameter Dependences
Including both cosmological and astrophysical parameters, our model is specified by the following quantities: As a reminder, the ΛCDM parameters are (in order of their appearance in Eq. (29)) the physical baryon density, the physical cold dark matter density, the vacuum energy density, the rms matter density fluctuation on comoving scales of 8 Mpc/h at z = 0, and the scalar spectral index. The ICM physics parameters C P0 and C β are defined in Eqs. (27) and (28), respectively, f NL is defined by Eq. (1), and M ν is the sum of the neutrino masses in units of eV.
We have placed f NL and M ν in parentheses in Eq. (29) in order to make it clear that we only consider scenarios in which these parameters are varied separately: for all cosmologies that we consider with f NL = 0, we set M ν = 0, and for all cosmologies that we consider with M ν > 0, we set f NL = 0. In other words, we only investigate one-parameter extensions of the ΛCDM concordance model. For the primary ΛCDM cosmological parameters, we use the parametrization adopted by the WMAP team (e.g., [69]), as the primordial CMB data best constrain this set. The only exception to this convention is our use of σ 8 , which stands in place of the primordial amplitude of scalar perturbations, A s . We use σ 8 both because it is conventional in the tSZ power spectrum literature and because it is a direct measure of the low-redshift amplitude of matter density perturbations, which is physically related more closely to the tSZ signal than A s . However, this choice leads to slightly counterintuitive results when considering cosmologies with M ν > 0, because in order to keep σ 8 fixed for such scenarios we must increase A s (to compensate for the suppression induced by M ν in the matter power spectrum).
For the fiducial model specified by the values in Eq. (30), we find that the tSZ power spectrum amplitude at = 3000 is ( + 1)C y /2π = 7.21 µK 2 at ν = 150 GHz. This corresponds to 7.59 µK 2 at ν = 148 GHz (the relevant ACT frequency) and 6.66 µK 2 at ν = 152.9 GHz (the relevant SPT frequency). The most recent measurements from ACT and SPT find corresponding constraints at these frequencies of 3.4 ± 1.4 µK 2 [12] and 3.09 ± 0.60 µK 2 [17] (using their more conservative error estimate). Note that the SPT constraint includes information from the tSZ bispectrum, which reduces the error by a factor of ∼2. Although it appears that our fiducial model predicts a level of tSZ power too high to be consistent with these observations, the results are highly dependent on the true value of σ 8 , due to the steep dependence of the tSZ power spectrum on this parameter. For example, recomputing our model predictions for σ 8 = 0.79 gives 5.52 µK 2 at ν = 148 GHz and 4.84 µK 2 at ν = 152.9 GHz, which are consistent at 3σ with the corresponding ACT and SPT constraints. Given that σ 8 = 0.79 is within the 2σ error bar for WMAP9 [69], it is difficult to assess the extent to which our fiducial model may be discrepant with the ACT and SPT results. The difference can easily be explained by small changes in σ 8 and is also sensitive to variations in the ICM physics, which we have kept fixed in these calculations. We conclude that our model is not in significant tension with current tSZ measurements (or other cosmological parameter constraints), and is thus a reasonable fiducial case around which to consider variations. Figs. 1 and 2 show the tSZ power spectra for our fiducial model and several variations around it, including the individual contributions of the one-and two-halo terms. In the fiducial case, the two-halo term is essentially negligible for > ∼ 300, as found by earlier studies [20], and it only overtakes the one-halo term at very low-( < ∼ 4). However, for f NL = 0, the influence of the two-halo term is greatly enhanced due to the scale-dependent bias described in Section II B, which leads to a characteristic upturn in the tSZ power spectrum at low-. In addition, f NL induces an overall amplitude change in both the one-and two-halo terms due to its effect on the halo mass function described in Section II A 1. While this amplitude change is degenerate with the effects of other parameters on the tSZ power spectrum (e.g., σ 8 ), the low-upturn caused by the scale-dependent bias is a unique signature of primordial non-Gaussianity, which motivates our assessment of forecasts on f NL using this observable later in the paper. Fig. 2 shows the results of similar calculations for M ν > 0. In this case, the effect is simply an overall amplitude shift in the one-and two-halo terms, and hence the total tSZ power spectrum. The amplitude shift is caused by the change in the halo mass function described in Section II A 2. Note that the sign of the amplitude change is somewhat counterintuitive, but arises due to our choice of σ 8 as a fundamental parameter instead of A s , as mentioned above. In order to keep σ 8 fixed while increasing M ν , we must increase A s , which leads to an increase in the tSZ power spectrum amplitude. Although this effect is degenerate with that of σ 8 and other parameters, the change in the tSZ power spectrum amplitude is rather large even for small neutrino masses (≈ 12% for M ν = 0.1 eV, which is larger than the amplitude change caused by f NL = 100). This sensitivity suggests that the tSZ power spectrum may be a useful observable for constraints on the neutrino mass sum.
We demonstrate the physical effects of each parameter in our model on the tSZ power spectrum in Figs. 3-11, including the effects on both the one-and two-halo terms individually. Note that the limits on the vertical axis in each plot differ, so care must be taken in assessing the amplitude of the change caused by each parameter. Except for f NL and M ν , the figures show ±1% variations in each of the parameters, which facilitates easier comparisons between their relative influences on the tSZ power spectrum. On large angular scales ( < ∼ 100), the most important parameters (neglecting f NL and M ν ) are σ 8 , Ω Λ , and C β . On very large angular scales ( < 10), the effect of f NL is highly significant, but its relative importance is difficult to assess, since the true value of f NL may be unmeasurably small. Note, however, that M ν is important over the entire range we consider, even if its true value is as small as 0.1 eV. Comparison of Figs. 4 and 8 indicates that the amplitude change induced by M ν = 0.1 eV (for fixed σ 8 ) is actually slightly larger than that caused by a 1% change in σ 8 around its fiducial value.
We now provide physical interpretations of the effects shown in Figs. 3-11: • f NL (3): The change to the halo mass function discussed in Section II A 1 leads to an overall increase (decrease) in the amplitude of the one-halo term for f NL > 0 (< 0). This increase or decrease is essentially -independent, is also seen at > 100 in the two-halo term, and is ±10% for f NL = ±100. More significantly, the influence of the scale-dependent halo bias induced by f NL = 0 is clearly seen in the dramatic increase of the two-halo term at low-. This increase is significant enough to be seen in the total power spectrum despite the typical smallness of the two-halo term relative to the one-halo term for a Gaussian cosmology.
• M ν (4): The presence of massive neutrinos leads to a decrease in the number of galaxy clusters at late times, as discussed in Section II A 2. This decrease would lead one to expect a corresponding decrease in the amplitude of the tSZ signal, but Fig. 4 shows an increase. This increase is a result of our choice of parameters -we hold σ 8 constant while increasing M ν , which means that we must simultaneously increase A s , the initial amplitude of scalar fluctuations. This increase in A s (for fixed σ 8 ) leads to the increase in the tSZ power spectrum amplitude seen in Fig. 4. The effect appears to be essentially -independent, although it tapers off slightly at very high-.
• Ω b h 2 (5): Increasing (decreasing) the amount of baryons in the universe leads to a corresponding increase (decrease) in the amount of gas in galaxy clusters, and thus a straightforward overall amplitude shift in the tSZ power spectrum (which goes like f 2 gas ).
• Ω c h 2 (6): In principle, one would expect that changing Ω c h 2 should change the tSZ power spectrum, but it turns out to have very little effect, as pointed out in [19], who argue that the effect of increasing (decreasing) Ω c h 2 on the halo mass function is cancelled in the tSZ power spectrum by the associated decrease (increase) in the comoving volume to a given redshift. We suspect that the small increase (decrease) seen in Fig. (6) when decreasing (increasing) Ω c h 2 is due to the fact that we hold Ω Λ constant when varying Ω c h 2 . Thus, Ω m ≡ 1 − Ω Λ is also held constant, and thus Ω b = Ω m − Ω c is decreased (increased) when Ω c is increased (decreased). This decrease (increase) in the baryon fraction leads to a corresponding decrease (increase) in the tSZ power spectrum amplitude, as discussed in the previous item. The slight -dependence of the Ω c h 2 variations may be due to the associated change in h required to keep Ω m constant, which leads to a change in the angular diameter distance to each cluster, and hence a change in the angular scale associated with a given physical scale. Increasing (decreasing) Ω c h 2 requires increasing (decreasing) h in order to leave Ω m unchanged, which decreases (increases) the distance to each cluster, shifting a given physical scale in the spectrum to lower (higher) multipoles. However, it is hard to completely disentangle all of the effects described here, and in any case the overall influence of Ω c h 2 is quite small.
• Ω Λ (7): An increase (decrease) in Ω Λ has several effects which all tend to decrease (increase) the amplitude of the tSZ power spectrum. First, Ω m is decreased (increased), which leads to fewer (more) halos, although this effect is compensated by the change in the comoving volume as described above. Second, for fixed Ω b /Ω c , this decrease (increase) in Ω m leads to fewer (more) baryons in clusters, and thus less (more) tSZ power. Third, more (less) vacuum energy leads to more (less) suppression of late-time structure formation due to the decaying of gravitational potentials, and thus less (more) tSZ power. All of these effects combine coherently to produce the fairly large changes caused by Ω Λ seen in Fig. 7. The slight -dependence may be due to the associated change in h required to keep Ω c h 2 and Ω b h 2 constant, similar (though in the opposite direction) to that discussed in the Ω c h 2 case above. Regardless, this effect is clearly subdominant to the amplitude shift caused by Ω Λ , which is only slightly smaller on large angular scales than that caused by σ 8 (for a 1% change in either parameter).
• σ 8 (8): Increasing (decreasing) σ 8 leads to a significant overall increase (decrease) in the amplitude of the tSZ power spectrum, as has been known for many years (e.g. [19]). The effect is essentially -independent and appears in both the one-and two-halo terms.
• n s (9): An increase (decrease) in n s leads to more (less) power in the primordial spectrum at wavenumbers above (below) the pivot, which we set at the WMAP value k piv = 0.002 Mpc −1 (no h). Since the halo mass function on cluster scales probes much smaller scales than the pivot (i.e., much higher wavenumbers k ∼ 0.1 − 1 h/Mpc), an increase (decrease) in n s should lead to more (fewer) halos at late times. However, since we require σ 8 to remain constant while increasing (decreasing) n s , we must decrease (increase) A s in order to compensate for the change in power on small scales. This is similar to the situation for M ν described above. Thus, an increase (decrease) in n s actually leads to a small decrease (increase) in the tSZ power spectrum on most scales, at least for the one-halo term. The cross-over in the two-halo term is likely related to the pivot scale after it is weighted by the kernel in Eq. (21), but this is somewhat non-trivial to estimate. Regardless, the overall effect of n s on the tSZ power spectrum is quite small.
• C P0 (10): Since C P0 sets the overall normalization of the ICM pressure profile (or, equivalently, the zero-point of the Y − M relation), the tSZ power spectrum simply goes like C 2 P0 . • C β (11): Since C β sets the logarithmic slope of the ICM pressure profile at large radii (see Eq. (22)), it significantly influences the total integrated thermal energy of each cluster, and thus the large-angular-scale behavior of the tSZ power spectrum. An increase (decrease) in C β leads to a decrease (increase) in the pressure profile at large radii, and therefore a corresponding decrease (increase) in the tSZ power spectrum on angular scales corresponding to the cluster outskirts and beyond. On smaller angular scales, the effect should eventually vanish, since the pressure profile on small scales is determined by the other slope parameters in the pressure profile. This trend is indeed seen at high-in Fig. 11. Note that a 1% change in C β leads to a much larger change in the tSZ power spectrum at nearly all angular scales than a 1% change in C P0 , suggesting that simply determining the zero-point of the Y − M relation may not provide sufficient knowledge of the ICM physics to break the long-standing ICM-cosmology degeneracy in tSZ power spectrum measurements. It appears that constraints on the shape of the pressure profile itself will be necessary.

IV. EXPERIMENTAL CONSIDERATIONS
In this section we estimate the noise in the measurement of the tSZ power spectrum. The first ingredient is instrumental noise. We describe it for the Planck experiment and for an experiment with the same specifications as the proposed PIXIE satellite [24]. The second ingredient is foregrounds 7 . We try to give a rather complete account of   all these signals and study how they can be handled using multifrequency subtraction. Our final results are in Fig. 12. Because of the several frequency channels, Planck and to a much larger extent PIXIE can remove all foregrounds and have a sensitivity to the tSZ power spectrum mostly determined by instrumental noise.

A. Multifrequency Subtraction
We discuss and implement multifrequency subtraction 8 along the lines of [47,48]. The main idea is to find a particular combination of frequency channels that minimize the variance of some desired signal, in our case the tSZ power spectrum. We hence start fromâ where a SZ refers to our estimator for the tSZ signal at 150 GHz (the conversion to a different frequency is straightforward), ν i are the different frequency channels relevant for a given experiment, w i are the weights for each channel, a m (ν i ) are spherical harmonic coefficients of the total measured temperature anisotropies at each frequency and finally g νi is the tSZ spectral function defined in Section III, allowing us to convert from Compton-y to ∆T . We can decompose the total signal according to a m = a SZ m + f a f m with f enumerating all other contributions. We will assume that a f m a f m ∝ δ f f , i.e. different contributions are uncorrelated with each other. Dropping for the moment the and m indices, the variance ofâ SZ is then found to be where C SZ = C SZ is the tSZ power spectrum at 150 GHz as given in Eq. (16) and C f (ν i , ν j ) (again the index is implicit) is the cross-correlation at different frequencies of the a f m of each foreground component (we will enumerate and describe these contributions shortly). To simplify the notation, in the following we will use We now want to minimize (â SZ ) 2 with the constraint that the weights describe a unit response to a tSZ signal, i.e., i w i = 1. This can be done using a Langrange multiplier λ and solving the system Because of the constraint i w i = 1, the C SZ term in (â SZ ) 2 is independent of w i (alternatively one can keep this term and see that it drops out at the end of the computation). Then the solution of the first equation can be written as where e j = 1 is just a vector with all ones and (C −1 ) ij is the inverse of C ij in Eq. (33). This solution can then be plugged back into the constraint i w i = 1 to give which is our final solution for the minimum-variance weights. From Eq. (32) we see that the total noise in eachâ SZ m after multifrequency subtraction is and the partial contributions to N from each foreground can be obtained by substituting C with C f (recall that there is an implicit index on C ij ). Notice that averaging over all m's for each and assuming that a given experiment covers only a fraction f sky of the sky, the final noise in each is N /(f sky (2 + 1)).

B. Foregrounds
We will consider the following sources of noise: instrumental noise (N ), CMB (CM B), synchrotron (Synch), freefree (f f ), radio and IR point sources (Radio and IR) and thermal dust (Dust). We now discuss each of them in turn.
We assume that the noise is Gaussian with a covariance matrix diagonal in l-space and uncorrelated between different frequencies. Then [129] where the beam size in radians at each frequency is θ(ν) = FWHM(ν)(8 ln 2) −1/2 × π/(180 × 60). The frequency channels ν, the FWHM(ν) (Full Width at Half Maximum) and ∆T (ν) depend on the experiment. In the following we consider the Planck satellite with specifications given in Table I and the proposed PIXIE satellite [24]. The latter is a fourth generation CMB satellite targeting primordial tensor modes through the polarization of the CMB. PIXIE will cover frequencies between 30 GHz and 6 THz with an angular resolution of 1 • .6 Gaussian FWHM corresponding to max ≡ θ −1 84. The frequency coverage will be divided into 400 frequency channels each with a typical sensitivity of ∆I = 4 × 10 −24 W m −2 sr −1 Hz −1 in each of 49152 sky pixels. In order to get ∆T we can use Planck's law with respect to CMB temperature where we used the numerical value of fundamental constants and T CM B = 2.725 K to write hν/(k B T CM B ) = ν/(56.8GHz). For example one finds ∆T (150GHz) 1.00 µK. For all foregrounds except point sources we use the models and parameters discussed in [47]. We assume that different components are uncorrelated and for each component f we define where Θ f (ν) encodes the frequency dependence, C f provides the -dependence and normalization at some fiducial frequency ν 0 (which will be different for different components) and finally R(ν i , ν j ) accounts for the frequency coherence. The latter ingredient was used in [47,48] and discussed in [49]. The general picture is that the auto-correlation of some contribution f at two different frequencies might not be perfect. Instrumental noise is an extreme case of this in which two different frequency channels have completely uncorrelated noise, i.e. R(ν i , ν j ) = δ ij . The CMB sits at the opposite extreme in that it follows a blackbody spectrum to very high accuracy, hence being perfectly coherent between any two frequencies: R(ν i , ν j ) = 1 for any i, j. All other foregrounds lie in between these two extrema, having an R(ν i , ν j ) that starts at unity for i = j and goes to zero as the frequencies are taken apart from each other. To model this Tegmark [49] proposed using where ξ f depends on the foreground and can be estimated as ξ f ∼ √ 2∆α −1 with ∆α being the variance across the sky of the spectral index of that particular component f . In the following we will write the frequency covariance as R(ν i , ν j , ∆α), e.g. for the we CMB we will have R(ν i , ν j , 0) while for instrumental noise R(ν i , ν j , ∞).
We will parameterize the frequency dependences of the various components as where we used Eq. (39) and Adding the information about the angular scale dependence we get For the CMB, C CM B is obtained using CAMB with the parameters of our fiducial cosmology. The parameters in the free-free, synchrotron and thermal dust components have been taken from the Middle Of the Road values in [47] (their Table 2 and text). The parameterization of the IR and Radio point sources follows [130].

C. Noise After Multifrequency Subtraction
Using the formulae in the last two sections we can estimate what the total variance inâ SZ will be after multifrequency subtraction. We denote the final result by N for the total noise and by N f for each foreground component (see around Eq. (37)). Then we plot [ ( + 1)N /2π] 1/2 and ( + 1)N f /2π in units of µK for Planck and PIXIE in Fig. 12.
With this choice we can compare directly with the results of [48] and see that they agree for Planck once one accounts for the fact that we are constraining tSZ at 150 GHz while there the tSZ in the Rayleigh-Jeans tail is considered, which brings a factor of about 4 difference in N .
The results for PIXIE are new. The proposed PIXIE design features 400 logarithmically-spaced frequency channels. This would require one to work with a very large multifrequency matrix which quickly becomes computationally expensive. Also, since C is very close to a singular matrix, the numerical inversion introduces some unavoidable error that becomes too large for matrices larger than about 35 × 35. For this reason, we decide to perform the computation binning the initial 400 channel into a smaller more manageable number. As we decrease the number of bins (i.e. bin more and more channels together) we expect two main effects to influence the final result. First, when the number of channels become comparable with the number of foregrounds that we want to subtract, the multifrequency subtraction will become very inefficient. Since we stay well away from this limit of very heavy binning, this is not an issue for us. Second, as we decrease the number of bins, the separation in frequency between adjacent bins grows larger. Because of the frequency decoherence (see discussion around Eq. (41)), when the bins are very far apart, they are contaminated by uncorrelated foregrounds and again the subtraction becomes inefficient. For a rough estimate of when this happens we take b 400 log 6000GHz 30GHz where b is the number of channels that we put in a bin, (6000/30) 1/400 is the logarithmic spacing and for ∆α we take the largest one appearing in the foregrounds, i.e. ∆α ∼ .5 for radio sources (dust and IR point sources have a comparable value). Then one finds that Eq. (55) starts being violated around b ∼ 8, which is what we take in our analysis. The last point is that if we want to cover the frequency most relevant for tSZ with 35 bins each containing 8 channels, we cannot start from the lowest frequency covered by PIXIE, namely 30 GHz. We decide instead to start at 45 GHz, since the signal at lower frequencies is swamped anyhow by synchrotron and free-free radiation. Summarizing, we take 35 logarithmically spaced frequency channels between 45 GHZ and 1836 GHz and use them for the multifrequency subtraction. Given the arguments above we do not expect that using more channels will improve the final noise appreciably. The final noise for PIXIE after multifrequency subtraction is shown in the bottom panel of Fig. 12. We show the weights w i at = 30 (as an example) for Planck and PIXIE in Fig. 13. As expected in both cases the weights are close to zero at 217 GHz, which is the null of the tSZ signal. Also in both cases, very low and very high frequencies have very small weights. Finally in Fig. 14 we compare the total noises for Planck and PIXIE after summing over m's and for a partial sky coverage f sky = 0.7, i.e. N [f sky (2 + 1)] −1/2 , with the expected tSZ signal. PIXIE leads to an improvement of more than two orders of magnitude at low 's, which is where the effect of the scale-dependent bias due to primordial non-Gaussianity arises. Since the PIXIE beam corresponds to max ∼ 84, the PIXIE noise becomes very large beyond ∼ few hundred where Planck is still expected to have signal-to-noise greater than one.

V. COVARIANCE MATRIX OF THE TSZ POWER SPECTRUM
In order to forecast parameter constraints and the detection SNR of the tSZ power spectrum, we must compute its covariance matrix. The covariance matrix contains a Gaussian contribution from the total (signal+noise) tSZ power spectrum observed in a given experiment, as well as a non-Gaussian cosmic variance contribution from the tSZ angular trispectrum. We compute the covariance matrix for three different experiments: Planck, PIXIE, and a future cosmic variance (CV)-limited experiment. The experimental noise after foreground subtraction is computed for Planck and PIXIE using the methods described in Section IV. For PIXIE, we assume a maximum multipole max = 300, while for Planck and the CV-limited experiment we assume a maximum multipole max = 3000. In the PIXIE and Planck cases, these values are well into the noise-dominated regime, so there is no reason to go to higher multipoles. For the CV-limited experiment, one can clearly compute up to as high a multipole as desired; however, it is unrealistic to imagine a satellite experiment being launched in the foreseeable future with noise levels better than PIXIE and angular resolution better than Planck, so we choose to adopt the semi-realistic value of max = 3000 for the CV-limited experiment. In all cases, we assume that the total available sky fraction used in the analysis is f sky = 0.7, i.e., 30% of the sky is masked due to unavoidable contamination from foregrounds in our Galaxy.
In the remainder of this section, we outline the halo model-based calculations used to compute the tSZ power spectrum covariance matrix and then discuss in detail the different masking scenarios that we consider to reduce the level of cosmic variance error in the results.

Halo Model Formalism
We compute the tSZ power spectrum covariance matrix using the halo model approach, as was used for the power spectrum itself in Section III A. We provide additional background on these calculations in Appendix A. The total tSZ power spectrum covariance matrix, M y , is given by Eq. (A28): where C y is the tSZ power spectrum given by Eq. (16), N is the noise power spectrum after foreground removal given by Eq. (37), and T y is the tSZ angular trispectrum. Note that we have neglected an additional term in the The amplitude of the one-halo term is suppressed relative to the two-halo term, leading to an increase in the multipole where their contributions are equal. In addition there is an overall decrease in the total signal at < ∼ 200, as one would expect after masking the massive, nearby clusters that dominate the one-halo term in this regime.
covariance matrix that arises from the so-called "halo sample variance" (HSV) effect (e.g., see Eq. (18) in [138] although that result is for the weak lensing power spectrum, the tSZ result is directly analogous). The HSV term becomes negligible in the limit of a full-sky survey, which is all we consider in this paper; thus, we do not expect this approximation to affect our results. Furthermore, we approximate the trispectrum contribution in Eq. (56) by the one-halo term only, which has been shown to dominate the trispectrum on nearly all angular scales [59]. We also restrict ourselves to the flat-sky limit, although the exact result is given in Appendix A. The tSZ trispectrum is thus given by Eq. (A27): Justifications for our approximations are given in Appendix A. We compute Eq. (56) for our fiducial WMAP9 cosmology using each of the masking scenarios discussed in the following section. These results are then combined with the parameter variations discussed in Section III C in order to compute Fisher matrix forecasts in Section VI.

Masking
For all of the experiments that we consider, the tSZ power spectrum covariance matrix in Eq. (56) is dominated by the trispectrum contribution over at least part of the multipole range relevant to that experiment. This issue is worse for PIXIE than for Planck, due to its lower noise levels, but even the Planck tSZ covariance is dominated by the trispectrum at some multipoles (10 < ∼ < ∼ 100). Fortunately, the trispectrum contribution can be significantly reduced by simply masking the massive, nearby galaxy clusters that dominate the signal, especially at low- [20].
However, there are some complications in this procedure. It is not obvious a priori that one should mask as many clusters as possible, since at some point the power spectrum signal itself will begin to decrease enough that the overall SNR decreases, despite the decrease in the noise. In addition, one must take care to use a very pure and complete cluster sample with a well-known selection function to do the masking; otherwise, it will be extremely difficult to properly account for the masking in the corresponding theoretical computations of the tSZ power spectrum. Finally, there is the unavoidable problem of scatter in the cluster mass-observable relation (e.g., L X -M or Y -M ), which will introduce additional uncertainty in the theoretical calculation of the masked tSZ power spectrum. Moreover, this uncertainty may be hard to precisely quantify.
In this paper, we consider a set of masking scenarios that approximately correspond to catalogs from existing and future all-sky X-ray surveys. At present, X-ray cluster surveys likely possess the most well-understood selection functions, as compared to those derived from optical, SZ, or weak lensing data (e.g., [60,64]). Furthermore, the scatter in the mass-observable relation for the X-ray quantity Y X (a measure of the integrated gas pressure) is believed to be quite small, possibly < 10% [101,124]. Although existing X-ray catalogs are only complete for fairly high masses and low redshifts (due to the steep decrease in X-ray surface brightness with redshift), these clusters are exactly the ones that need to be masked to suppress the tSZ trispectrum. Furthermore, the future catalogs from eROSITA, an upcoming X-ray satellite designed for an all-sky survey, will be complete to much lower masses and very high redshifts. Overall, it seems that X-ray-based masking is the most robust option at present.
In order to simplify our calculations and avoid the need to specify the details of any individual X-ray survey, we compute the effects of masking by removing all clusters in the tSZ calculations that lie above a mass threshold M mask and below a redshift cutoff z mask . This method also circumvents the issue of modeling scatter in the mass-observable relation. Clearly a more sophisticated approach will be needed for the analysis of real data, but these choices allow us to explore several possibilities fairly quickly.
First, we consider a masking criterion based on the results of the ROSAT All-Sky Survey (RASS) [61]. The catalogs derived from RASS (e.g., BCS/eBCS [62,63], REFLEX [64]) are nominally complete to a flux limit f X > 3 × 10 −12 erg/cm 2 /s (0.5-2.0 keV band). At z = 0.05, this flux corresponds to a luminosity L X ≈ 9 × 10 42 h 2 erg/s. Using the scaling relations from [65] (for either optical or X-ray masses) to convert from L X to mass, this luminosity corresponds to M ≈ 10 14 M /h, where M is the virial mass defined in Section II A (we have also used the NFW profile [66] and the concentration-mass relation from [67] in this conversion). Thus, the catalogs derived from RASS should be ∼ 100% complete for clusters with M > ∼ few × 10 14 M /h at z < 0.05. To be conservative, we set M mask = 5 × 10 14 M /h, which corresponds to a luminosity and flux of L X ≈ 7 × 10 44 h 2 erg/s and f X ≈ 2 × 10 −11 erg/cm 2 /s at z = 0.05, well above the flux limit given above. We will refer to this masking choice (M mask = 5 × 10 14 M /h, z mask = 0.05) as "ROSAT-masked" in the remainder of the paper.
Second, we consider a masking criterion based on the upcoming all-sky survey conducted by eROSITA. The eROSITA cluster catalogs are expected to be nominally complete to a flux limit f X > 4 × 10 −14 erg/cm 2 /s (0.5-2.0 keV band) [68]. This limit is low enough that it will be possible to essentially mask clusters arbitrarily at low redshifts. At z = 0.05, this flux corresponds to L X ≈ 10 41 h 2 erg/s, far below the emission from even a ∼ 10 13 M /h group, for which L X ∼ 10 42 h 2 erg/s (extrapolating the scaling relations from [65]). Thus, there is a wide range of possible masking choices based on the eROSITA catalogs. In principle, it would be best to fully explore the (M mask , z mask ) parameter space and find the values that optimize the SNR for the tSZ power spectrum, or perhaps optimize the constraints on some particular parameter, such as f NL . This optimization is beyond the scope of this paper. We anticipate that masking heavily above z ∼ 0.1-0.2 will eventually begin to decrease the tSZ power spectrum signal too significantly, so we make a reasonably conservative cut at z mask = 0.15. As mentioned, at these redshifts eROSITA will detect essentially all clusters (at z = 0.15, the flux limit corresponds to L X ≈ 1.2 × 10 42 h 2 erg/s, which roughly scales to M ≈ 2 × 10 13 M /h using the same scaling relations as above). Thus, we can choose the mass threshold at essentially any value. Our final values are (M mask = 2 × 10 14 M /h, z mask = 0.15); we will refer to this masking choice as "eROSITA-masked" in the remainder of the paper.
Finally, we also consider a completely unmasked calculation, both due to its theoretical simplicity and as a means to assess how much the choice of masking can improve the SNR on the tSZ power spectrum, as well as how the forecasted parameter constraints change. Throughout the paper, unless figures or tables are labelled with a masking choice, they have been calculated in the unmasked case.
In Fig. 15 we demonstrate the effects of masking on the tSZ power spectrum computed using our fiducial parameters. The main effects are a decrease in the amplitude of the one-halo term relative to the two-halo term at low-and an overall decrease in the amplitude of the power spectrum at low-. In the unmasked case, the one-halo term and two-halo term are roughly equal at 4, while this cross-over multipole increases to 14 for the ROSAT-masked case and 56 for the eROSITA-masked case. Note that the enhancement of the two-halo term relative to the one-halo term increases the sensitivity of the power spectrum to f NL through the effect of the scale-dependent bias shown in Fig. 1. Although the amplitude of the total signal is decreased by masking, the suppression of the cosmic variance error due to the trispectrum is much larger, as seen in Figs. 16-21. In addition, we show the 1σ fractional errors on the tSZ power spectrum for each of the three experiments we consider, computed from the square root of the diagonal of the covariance matrix. In this unmasked calculation, the overwhelming influence of cosmic variance due to the tSZ trispectrum is clearly seen at low-.
Figs. 16-18 demonstrate the significant reduction in error due to masking. These figures show the fractional difference with respect to our fiducial WMAP9 results for the tSZ power spectrum computed in three different cosmologies (f NL = 100 (or 50), f NL = −100 (or −50), and M ν = 0.10 eV). Note that these power spectra are computed using the same masking criterion as used in the covariance matrix. The figures also show the 1σ fractional errors on the tSZ power spectrum for each of the three experiments we consider. These fractional errors have been computed from the diagonal of the covariance matrix, although we emphasize that the entire covariance matrix is used in all calculations presented in the paper. The overall implication of Figs. 16-18 is that masking is in general highly beneficial for tSZ power spectrum measurements. Although the improvements for Planck are marginal due to the fact that it is instrumental noise-dominated (see the subsequent figures), the improvements for PIXIE and the CV-limited case are dramatic. We note, however, that the reduction in the signal for the eROSITA-masked case is large enough at low-that the fractional errors for PIXIE are actually larger than in the ROSAT-masked case, because the signal at low-is becoming smaller than the PIXIE instrumental noise. Nonetheless, at higher multipoles the PIXIE fractional errors are in fact slightly smaller for the eROSITA-masked case than ROSAT-masked. This trend suggests that the optimal choice of masking for a given experiment depends on the noise levels of the experiment and the multipole range that one would like to measure with the highest SNR. For the CV-limited case, it is perhaps unsurprising that the fractional errors continue to decrease with heavier masking, as the trispectrum term is further and further suppressed relative to the Gaussian term in Eq. (56). It is unclear how far one can push this masking before the SNR starts to decrease due to the reduction in the signal, but Fig. 21 (discussed in the next paragraph) suggests that our eROSITA masking scenario is not too far from this limit.
Figs. [19][20][21] show the contributions to (and total of) the diagonal of the covariance matrix from the Gaussian and non-Gaussian terms in Eq. (56) for each of the experimental and masking scenarios. From Fig. 19, one can see that Planck is dominated by instrumental noise over most of its multipole range, the exception being a window  Fig. 16 except that it has been computed for the ROSAT-masked scenario (see Section V 2). Compared to Fig. 16, it is clear that the errors at low-for PIXIE and the CV-limited case have been dramatically reduced by the masking procedure. The errors from Planck are less affected because it is dominated by instrumental noise over most of its multipole range. from 10 < ∼ < ∼ 100 where the trispectrum contribution dominates. Fig. 19 also demonstrates that PIXIE is near the CV-limited case over essentially its entire multipole range. This explains the large decrease in the fractional errors for PIXIE after applying the ROSAT-masking in Fig. 17. Fig. 20 shows the covariance matrix contributions in the ROSAT-masked case, where it is clear that Planck is now completely dominated by instrumental noise at all multipoles. Thus, further masking is not beneficial for Planck. PIXIE is still near the CV-limited case for the ROSAT-masked scenario, except at low-where its instrumental noise starts to become significant. This suggests that further masking will decrease the SNR at low multipoles for PIXIE, as indeed can be seen in the eROSITA-masked scenario of Fig. 18. Finally, Fig. 21 shows the covariance matrix contributions in the eROSITA-masked case, where one can see that the trispectrum term has now been suppressed below the Gaussian term at all multipoles, even for the CV-limited experiment. This trend likely suggests that further masking will begin to lead to a reduction in the SNR even for the CV-limited case, although we have not computed this precisely, as we have no masking cases beyond eROSITA. Both PIXIE and Planck are clearly dominated by instrumental noise for the eROSITA-masked case, suggesting that further masking is unlikely to lead to significant improvements in SNR for these experiments.

VI. FORECASTED CONSTRAINTS
Having described our model for the signal and its covariance matrix, we now use the Fisher matrix formalism [129,[131][132][133] to forecast the expected constraints on the parameters listed in Eq. (29) for a total of nine different experimental specifications (Planck, PIXIE and CV-limited) and masking choices (unmasked, ROSAT-masked, eROSITA-masked). We compute derivatives of the tSZ power spectrum signal with respect to each parameter around  Fig. 16 except that it has been computed for the eROSITA-masked scenario (see Section V 2). Compared to Fig. 16, it is clear that the errors have been significantly reduced at low-, though again the effect is minimal for Planck. The comparison to the ROSAT-masked case in Fig. 17 is more complicated: it is clear that the errors are further reduced for the CV-limited case, but the PIXIE fractional errors actually increase at low-as a result of the reduction in the amplitude of the signal there due to the masking. However, the PIXIE fractional errors at higher multipoles are in fact slightly reduced compared to the ROSAT-masked case, although it is hard to see by eye in the plot.
the fiducial values given in Eq. (30). The Fisher matrix is given by where p i is the i th parameter in Eq. (29), C y is the tSZ power spectrum given by Eq. (16), and M −1 is the inverse of the covariance matrix given by Eq. (56). Note that we only consider f NL = 0 and M ν > 0 cosmologies separately; in the interest of simplicity, we only seek to constrain minimal one-parameter extensions of the ΛCDM standard model. Thus, our Fisher matrices are eight-by-eight, containing the five relevant ΛCDM parameters, the two ICM physics parameters (C P0 and C β ) and either f NL or M ν . The unmarginalized 1σ error on parameter p i is given by 1/ √ F ii , and it describes the best possible error one can obtain when all other parameters are known exactly. The marginalized 1σ error is given by (F −1 ) ii , and it describes the case in which all other parameters are constrained from the same set of data. In the following sections, we quote the unmarginalized 1σ errors on each of the parameters from the tSZ power spectrum alone, as well as marginalized 1σ errors computed by the following methods. The results are summarized in Fig. 22. In addition, we provide the complete Fisher matrices for these calculations online 9 .
In order to compute marginalized constraints, we find that it is necessary to include external data (in addition to basic priors, which are described below), because the degeneracies between the various parameters are too strong for the tSZ signal alone to provide meaningful constraints. For this purpose, we include a Fisher matrix computed for the imminent results from the Planck measurements of the primordial CMB temperature fluctuation power spectrum 10 . This Fisher matrix includes our five primary cosmological parameters Ω b h 2 , Ω c h 2 , Ω Λ , σ 8 , n s , as well as the optical depth to reionization, τ (note that it does not include any non-ΛCDM parameters, i.e., f NL and M ν are not included). We marginalize over τ , since the tSZ signal is insensitive to this parameter. In addition, we include a prior on H 0 from the results of [134], who find a ≈ 3.3% constraint on this parameter from cosmological distance ladder measurements. Although their central value is slightly discrepant from our fiducial value, we simply include the statistical power of this measurement as a representative measure of current constraints on h. To include this prior, we transform to a parameter set in which h lies on the diagonal, add a Gaussian prior of width ∆h = 0.033 × 0.697 = 0.023, where h = 0.697 is the value for our fiducial cosmology, and then transform back to our original parameters. We add the resulting "Planck CMB + H 0 " Fisher matrix to the appropriate five-dimensional sub-matrix of our tSZ Fisher matrices in order to break degeneracies between the various parameters, and investigate to what extent the tSZ signal can improve on the upcoming Planck CMB measurements of the primary cosmological parameters. In addition to the Planck CMB prior matrix for the primary cosmological parameters, we also place simple priors on the ICM physics parameters C P0 and C β . We adopt 1σ Gaussian priors of ∆C P0 = 0.2 = ∆C β (recall that these parameters are normalized versions of the GNFW pressure profile amplitude and outer logarithmic slope). These values correspond roughly to the variances determined for these parameters due to the scatter between cluster  Fig. 19 except that it has been computed for the ROSAT-masked scenario (see Section V 2). Results are shown for the total in each experiment, as well as the different Gaussian terms in each experiment and the trispectrum contribution that is identical for all three. See the text for discussion. pressure profiles in the simulations from which the fiducial model was obtained 11 . They also encompass values obtained from various X-ray and SZ observations [3,119,120], although we note that our use of P 200,c in Eq. (22) makes direct comparisons of the P 0 parameter somewhat nontrivial between our model and others in the literature. Finally, these priors guarantee that the unphysical values C P0 = 0 and C β = 0 are highly disfavored (> 5σ). We find that the Fisher forecasts computed below are often strong enough to provide constraints on the ICM parameters well below the width of these priors (at least for C β , which is currently the less well-constrained of the two), and thus we conclude that the resulting errors are robust.
We place no priors on f NL or M ν as our goal is to assess the detectability of these parameters using information in the tSZ power spectrum, with no external constraints (apart from those necessary to break degeneracies amongst the primary cosmological parameters, as implemented in our Planck CMB+H 0 prior matrix).
In the following subsections we describe our main parameter forecast results, summarized in Fig. 22, for the two separate cases in which the standard cosmological and astrophysical parameters (see Eq. (29)) have been supplemented by either f NL or M ν . Fig. 19 except that it has been computed for the eROSITA-masked scenario (see Section V 2). Results are shown for the total in each experiment, as well as the different Gaussian terms in each experiment and the trispectrum contribution that is identical for all three. See the text for discussion.

A. fNL
We discuss local primordial non-Gaussianity and its effects on the tSZ power spectrum in Sections II and III C. The current strongest bounds come from WMAP9 and correspond to −3 < f NL < 77 at 95% CL [37], but Planck will likely improve on these limits by a factor of ≈ 3. From the top table in Fig. 22, one can see that at the unmarginalized level the tSZ power spectrum is quite sensitive to f NL , with an unmarginalized CV limit of ∆f NL ≈ 7. On the other hand, f NL is very degenerate with all the other parameters and the marginalized bounds are far weaker than the current bounds (see the central right table in Fig. 22). Adding the tSZ power spectrum constraint to the current bounds in quadrature improves the overall constraint beyond the current bounds by at most a few percent, even in the most optimistic scenario (CV-limited experiment with eROSITA masking). We have checked that no single parameter is driving the marginalized error on f NL away from the unmarginalized value, which would be comparable with constraints from the Planck primordial temperature bispectrum measurements. On the contrary, in order to get closer to the unmarginalized bound, all other parameters, both cosmological and astrophysical, would need to be constrained much better. We point out that if one ever hopes to constrain f NL using the tSZ power spectrum, it is crucial to obtain high SNR on the lowest possible multipoles, where the signature of the scale-dependent bias breaks the degeneracy between f NL and other parameters that influence the overall amplitude of the power spectrum. This fact explains why PIXIE would achieve much tighter marginalized constraints on f NL using the tSZ power spectrum than Planck would (as seen in the central right table in Fig. 22), despite Planck's larger multipole range. Nevertheless, even the PIXIE constraints are unlikely to be competitive with those from other probes of primordial non-Gaussianity.
Thus far we have only considered constraints around our fiducial cosmology with f NL = 0. If instead we assume a fiducial f NL = 37, corresponding to the central value of WMAP9 [37], we find as before that Planck and PIXIE would have less than a 2σ detection. However, for a CV-limited experiment using eROSITA masking, we find a marginalized  ) are derived using only the tSZ power spectrum. The marginalized errors (two central tables) are derived by adding as an external prior the forecasted Planck constraints from the CMB T T power spectrum plus a prior on H0 from [134]. Since adding the tSZ power spectrum leads to effectively no improvement in the errors on the ΛCDM parameters with respect to the Planck T T priors, we do not show those numbers in the two central tables (marginalized tSZ plus Planck CMB+H0 priors). They are equal to the numbers in the bottom table, which show the marginalized Planck CMB+H0 priors, to ≈ 1% precision. Note that all constraints on Mν are in units of eV, while the other parameters are dimensionless.
error ∆f NL 13. Thus, if f NL turns out to be of this magnitude, one could obtain a 3σ detection using methods completely independent of the primordial CMB temperature bispectrum. The full results of these calculations are given in Fig. 23.

B. Mν
We discuss massive neutrinos and their effects on the tSZ power spectrum in Sections II and III C. The current strongest bounds are in the range M ν < ∼ 0.3 eV [12,42,43,101]. However, a 3σ detection near this mass scale was recently claimed in [44]. From the top table in Fig. 22, it is clear that the unmarginalized constraints on M ν from the tSZ power spectrum alone are quite strong. The Planck tSZ power spectrum unmarginalized error ∆M ν 0.03 eV is slightly smaller than the lower bounds from neutrino oscillations, and the unmarginalized error for a CV-limited experiment would lead to a very robust detection of M ν . The bounds from PIXIE are much weaker due to its lower angular resolution; the change induced by M ν in the tSZ power spectrum is effectively an overall amplitude shift,  Fig. 22, but for a fiducial Mν = 0.1 eV, similar to the minimum allowed value in the inverted neutrino hierarchy [39]. Note that all constraints on Mν are in units of eV.
as seen in Fig. 4, and thus one can gain much more leverage on this parameter by going to higher multipoles. For PIXIE, one must mask heavily in order to reduce the CV errors to a sufficient level at the low multipoles where it observes in order to measure the effect of small neutrino masses. Masking makes much less of a difference in the Planck and CV-limited constraints on M ν because most of their constraining power comes from higher multipoles where the masking procedure does not significantly reduce the error. After adding the external Planck CMB+H 0 priors on the ΛCDM parameters (which are summarized in the bottom table of Fig. 22), the marginalized error for Planck with eROSITA masking is ∆M ν 0.37 eV, comparable with current bounds. This could be useful to strengthen current bounds and confirm or reject a detection at this level. For a CV-limited experiment, the error is ∆M ν 0.20 eV, a factor of two better than Planck. Note that we have derived these bounds after fully marginalizing over the primary ΛCDM parameters and the ICM gas physics parameters, so they are unlikely to be overly optimistic. In addition, we have only used the information on M ν contained the tSZ power spectrum; we have not used any information from the Planck CMB temperature power spectrum regarding M ν . Including primordial CMB constraints could significantly tighten the forecasted Planck errors to a level well below the current constraints. Given the imminent release of the Planck sky maps, we leave this as an avenue to be pursued with data.
Thus far we have only considered constraints around our fiducial cosmology with M ν = 0. If instead we assume a fiducial M ν = 0.1 eV, roughly corresponding to the minimum value in the inverted neutrino hierarchy [39], we find that the marginalized Planck bound after eROSITA masking becomes slightly tighter than the current upper limits in the literature 0.3 eV. Note that the forecasted errors depend on the fiducial neutrino mass assumed; if the actual neutrino mass is near the current upper limits, the eROSITA-masked Planck 1σ error would lie below the actual mass, suggesting a possible marginal detection. However, we have focused on forecasts for M ν near the minimum allowed values in order to be conservative.
More interestingly, we find that the primary degeneracy of M ν is with the ICM physics parameters C P0 and C β . These degeneracies are comparable to but much stronger than that of M ν with any of the ΛCDM parameters (after Gaussian prior ∆CP 0 = ∆C β ≡ ∆gas placed on the gas physics parameters (our fiducial value is 0.2). Note that we vary the priors simultaneously, keeping them fixed to the same value. The plots demonstrate that even modest improvements in the external priors on CP 0 and C β could lead to significant decreases in the expected error on Mν from tSZ power spectrum measurements.
imposing the Planck CMB+H 0 prior), which suggests that strong external constraints on the ICM pressure profile would be of great use in constraining M ν using tSZ power spectrum measurements. We investigate this issue in detail in Fig. 25, which shows the forecasted 1σ error on M ν as a function of the width of the Gaussian prior placed on C P0 and C β , which is 0.2 in our standard case. We provide plots for both the fiducial M ν = 0 case and the case where the fiducial M ν = 0.1 eV. The latter plot indicates that strengthening the prior on C P0 and C β by a factor of two would lead to a 40% decrease in the forecasted error on M ν for the eROSITA-masked Planck experiment. Although the improvement eventually saturates as the degeneracies with other parameters become important for very tight priors on C P0 and C β , these results demonstrate the importance of strong external constraints on the ICM physics in tightening the bounds on cosmological parameters from tSZ measurements. Fortunately, it will be possible to derive such external constraints from detailed studies of the ICM pressure profile with Planck (e.g., [120]) and eROSITA themselves, which portends a bright future for tSZ-based cosmological constraints. We conclude that the tSZ power spectrum is a promising probe for the neutrino masses and could lead to interesting results already using Planck data. We find ∆M ν 0.3 eV around a fiducial M ν = 0 using the tSZ power spectrum alone to probe M ν (including other data in the M ν constraint would tighten this bound). Considering larger values of the fiducial M ν and/or imposing stronger priors on the ICM gas physics parameters C P0 and C β can shrink this bound considerably, perhaps to the level of a 2σ measurement using Planck data alone.

C. Other Parameters
Although our focus in this paper has been on constraining currently unknown extensions to the ΛCDM standard model, we also compute the forecasted constraints on all of the standard parameters in our model. Unfortunately, we find that adding the tSZ power spectrum data to the priors from the Planck CMB+H 0 Fisher matrix yields essentially no improvement in the forecasted errors, which are already very small using Planck+H 0 alone. Even for a CV-limited experiment with eROSITA masking, the marginalized constraints on the five ΛCDM parameters in our model only improve by ≈ 1 − 2% beyond the marginalized priors from Planck+H 0 (which are given in the bottom table in Fig. 22). The primary cosmological utility of tSZ measurements is in providing a low-redshift probe of the amplitude of fluctuations, which permits a constraint on M ν when combined with a high-redshift probe of this amplitude from the CMB.
In contrast to the results for the cosmological parameters, the forecasted constraints on the ICM physics parameters C P0 and C β are more encouraging, as seen in the central tables in Fig. 22, as well as Figs. 23 and 24. This is especially true of the outer slope parametrized by C β , for which we forecast a ≈ 6−8% constraint using the unmasked Planck data, which decreases to nearly the percent level after masking with eROSITA. Most current observational measurements of this parameter do not have reported error bars, but the analysis of SPT stacked SZ profiles in [3] found a ∼ 40% uncertainty in this parameter after marginalization, which is probably a representative value. We thus expect this error to decrease dramatically very shortly. However, it will be difficult for Planck or PIXIE to constrain the overall normalization of the ICM pressure profile using the tSZ power spectrum alone due to its strong degeneracy with the cosmological parameters (e.g., σ 8 ). This is reflected in the fact that the marginalized constraints on C P0 in Figs. [22][23][24] are often near the bound imposed by our prior, the exceptions being either of the masked Planck measurements or any of the CV-limited measurements. The best approach to obtaining constraints on C P0 is likely using cross-correlations between the tSZ signal and lensing maps of the dark matter distribution. However, it is also worth noting that we have not considered the minimal ΛCDM case in these calculations, in which both f NL and M ν are fixed to zero. In such a scenario, the bounds on the ICM physics parameters would be stronger than those quoted here. However, given that we know M ν > 0 in our universe, it is perhaps most reasonable to look at the marginalized constraints in that case as an example of future constraints on C P0 and C β .

D. Forecasted SNR
The cumulative SNR on the tSZ power spectrum using multipoles < max is given by where (M y ) −1 refers to the inverse of the max -bymax submatrix of the full covariance matrix. The cumulative SNR provides a simple way to assess the constraining power of a given experimental and masking choice on the tSZ power spectrum, without regard to constraints on particular parameters. We compute Eq. (59) for our fiducial cosmology using each of the experiment (Planck, PIXIE, CV-limited) and masking options (unmasked, ROSATmasked, eROSITA-masked) considered in the previous sections. Note that the covariance matrix in Eq. (59) includes all contributions from cosmic variance (Gaussian and non-Gaussian) and experimental noise after foreground removal, as discussed in Section V. Note that our approach includes the trispectrum (or sample variance) contribution to the covariance matrix in calculating the SNR, which is perhaps more conservative than an approach in which only the Gaussian errors are considered in assessing the SNR. The results of these calculations are shown in Fig. 26. We find that Planck can detect the tSZ power spectrum with a cumulative SNR ≈ 35 using multipoles < 3000 (see Appendix B for a comparison with the recently-released initial tSZ results from the Planck collaboration [137]). This result is essentially independent of the masking scenario, although masking more heavily can lead to greater cumulative SNR using lower values of max , as compared to the unmasked case. However, masking leads to significant improvement in the PIXIE results: the unmasked PIXIE cumulative SNR using < 300 is ≈ 5.8, while the ROSAT-and eROSITA-masked results are ≈ 8.9 and 22, respectively. These results follow from the fact that PIXIE is nearly CV-limited for < 100, as seen in the unmasked curves shown in Fig. 26 (and discussed earlier in Section V). For the masked cases, the CV errors are reduced sufficiently that the PIXIE noise starts to become important at < ∼ 10. There is one subtlety of the masking procedure that can be understood by considering Fig. 26 in combination with Figs. 19-21. Looking at the ROSAT-masked case for Planck in Fig. 20, it appears that the Planck errors are dominated by the Gaussian instrumental noise term at all , and hence that masking further for Planck should be harmful rather than helpful; this appears to be confirmed by the fact that the Planck fractional errors in the eROSITA-masked case in Fig. 18 are indeed larger than in the ROSAT-masked case in Fig. 17. However, it is clear in Fig. 26 that the eROSITA-masked Planck case has a larger total SNR than the ROSAT-masked case. The resolution of this apparent discrepancy lies in the fact that the masking continues to suppress the off-diagonal terms in the covariance matrix, which arise solely from the trispectrum term in Eq. (56), even as the on-diagonal fractional errors begin to increase. The plots in Figs. 16-21 only show the diagonal entries in the covariance matrix, and thus one may not realize the effect of the masking on the off-diagonal terms in the covariance matrix without examining the cumulative SNR (this result is also implied by the improved parameter constraints for the eROSITA-masked Planck case given in the previous sections). These considerations imply that there is likely an optimal masking choice for a given experimental noise level and survey specifications, but obtaining the precise answer to this question lies beyond the scope of this paper.
Overall, Fig. 26 indicates that near-term data promises highly significant detections of the tSZ power spectrum on larger angular scales than have been probed thus far by ACT and SPT. Note that PIXIE is close to the CV limit over its signal-dominated multipole range. The total SNR using the imminent Planck data is ≈ 35, essentially independent of the masking option used (note that masking can increase the cumulative SNR up to lower multipoles, however, as compared to the unmasked case).

VII. DISCUSSION AND OUTLOOK
In this paper we have performed a comprehensive analysis of the possible constraints on cosmological and astrophysical parameters achievable with measurements of the tSZ power spectrum from upcoming full-sky CMB observations, with a particular focus on extensions to the ΛCDM standard model parametrized by f NL and M ν . We have included all of the important physical effects due to these additional parameters, including the change to the halo mass function and the scale-dependent halo bias induced by primordial non-Gaussianity. Our halo model calculations of the tSZ power spectrum include both the one-and two-halo terms, and we use the exact expressions where necessary to obtain accurate results on large angular scales. We model the ICM pressure profile using parameters that have been found to agree well with existing constraints, and furthermore we model the uncertainty in the ICM physics by freeing two of these parameters. We also include a realistic treatment of the instrumental noise for the Planck and PIXIE experiments, accounting for the effects of foregrounds by using a multifrequency subtraction technique. Our calculations of the covariance matrix of the tSZ power spectrum include both the Gaussian noise terms and the non-Gaussian cosmic variance term due to the tSZ trispectrum. We investigate two masking scenarios motivated by the ROSAT and eROSITA all-sky surveys, which significantly reduce the large errors that would otherwise be induced by the trispectrum term, especially at low-. Finally, we use these calculations to forecast constraints on f NL , M ν , the primary ΛCDM parameters, and two parameters describing the ICM pressure profile.
Our primary findings are as follows: • The tSZ power spectrum can be detected with a total SNR > 30 using the imminent Planck data up to = 3000, regardless of masking (see Appendix B for a comparison with the initial tSZ power spectrum results released by the Planck collaboration while this manuscript was under review [137]) ; • The tSZ power spectrum can be detected with a total SNR between ≈ 6 and 22 using the future PIXIE data up to = 300, with the result being sensitive to the level of masking applied to remove massive, nearby clusters; • Adding the tSZ power spectrum information to the forecasted constraints from the Planck CMB temperature power spectrum and existing H 0 data is unlikely to significantly improve constraints on the primary cosmological parameters, but may give interesting constraints on the extensions we consider: -If the true value of f NL is near the WMAP9 ML value of ≈ 37, a future CV-limited experiment combined with eROSITA-masking could provide a 3σ detection, completely independent of the primordial CMB temperature bispectrum; alternatively, PIXIE could give 1 − 2σ evidence for such a value of f NL with this level of masking; -If the true value of M ν is near 0.1 eV, the Planck tSZ power spectrum with eROSITA masking can provide upper limits competitive with the current upper bounds on M ν ; with stronger external constraints on the ICM physics, Planck with eROSITA masking could provide 1 − 2σ evidence for massive neutrinos from the tSZ power spectrum, depending on the true neutrino mass; • Regardless of the cosmological constraints, Planck will allow for a very tight constraint on the logarithmic slope of the ICM pressure profile in the outskirts of galaxy clusters, and may also provide some information on the overall normalization of the pressure profile (which sets the zero point of the Y − M relation).
Our results are subject to a few caveats. We have made the usual Fisher matrix approximation that the likelihood function is nearly Gaussian around our fiducial parameter values, but this should be safe for small variations, which are all that we consider (in particular, σ 8 is tightly constrained by the external Planck CMB prior, and it would be most likely to have a non-Gaussian likelihood). We have also neglected any tSZ signal from the intergalactic medium, filaments, or other diffuse structures, but the comparison between simulations and halo model calculations in [116] indicates that this approximation should be quite good. We have also assumed that the mass function parameters are perfectly well known, while in reality some uncertainties remain, especially in the exponential tail. Given that our most optimistic results involve masking nearly all of the clusters that live in the exponential tail, we believe that our forecasts should be fairly robust to the mass function uncertainties, in contrast to cluster count calculations which are highly sensitive to small changes in the tail of the mass function. Finally, we have only included the flat-sky version of the one-halo term in our computations of the tSZ power spectrum covariance matrix. For the masked calculations, the flat-sky result should suffice, since massive, nearby clusters are removed; however, it is possible that the two-, three-, or four-halo terms could eventually become relevant in the masked calculations. These would be largest at low-, however, where our primary interest is in constraining f NL , which does not appear very optimistic in any case. Thus, we neglect these corrections for our purposes.
There are many future extensions of this work involving higher-order tSZ statistics and cross-correlations with other tracers of large-scale structure. Recent work on the tSZ bispectrum and skewness [15,17,21,22] indicates that significantly stronger constraints on both cosmology and the ICM physics can be obtained by using higher-order statistics. These may also be a better place to look for f NL constraints, as the additional powers of the halo bias could lead to a larger signal at low-than in the power spectrum (controlling systematics will be of paramount importance, as will masking to reduce the very large cosmic variance due to the tSZ six-point function). Determining an optimal strategy to extract the neutrino mass through combinations of tSZ statistics and cross-spectra with other tracers is also work in progress. The key factor remains breaking the degeneracy with the ICM physics, or, more optimistically, simultaneously constraining both the ICM and cosmological parameters using tSZ measurements.
In the halo model, it is assumed that all matter in the universe is bound in halos. Each halo of virial mass M is assumed to have a density profile, ρ( x; M ), and (for our purposes) an electron pressure profile, P e ( x; M ). The mass density field at position x is then given by the sum of the contributions from all halos in the universe: Similarly, the electron pressure field at position x is given by: For calculations involving the tSZ effect, it is convenient to define a "3D Compton-y" field that is simply a re-scaling of the electron pressure field: Note that the 3D Compton-y field has dimensions of inverse length (it is thus important to be careful about comoving versus physical units -in our calculations using the Battaglia pressure profile, the pressure is given in physical units, and thus so is y 3D ). The usual (2D) Compton-y field is then given by the LOS projection of y 3D ( x): where t is the age of the universe at a given epoch, χ(t) is the comoving distance to that epoch, a(χ) is the scale factor at that epoch, andn is a unit vector on the sky. We have used dt/da = 1/(aH) and dχ/da = −c/(Ha 2 ) in going from the first line to the second line, where H(a) is the Hubble parameter. Defining the projection kernel W y (χ) for the Compton-y field via y( n) = dχW y (χ)y 3D ( x(χ,n)), we thus have: Starting from the 3D Compton-y field defined in Eq. (A3), we derive the angular power spectrum of the 2D Compton-y field. First, we calculate the relevant 3D power spectrum by means of the halo model (N.B. in this expression and many others in the following, the redshift dependence will be suppressed for notational simplicity): where the one-halo term is and the two-halo term is In these equations,ỹ 3D ( k; M ) is the Fourier transform of the 3D Compton-y profile for a halo of virial mass M : where r = | r|, k = | k|, and we have assumed that y 3D ( r; M ) is spherically symmetric to obtain the second expression. Also, in Eqs. (A7) and (A8), dn(M, z)/dM is the comoving number density of halos of mass M at redshift z, b(M, z) is the bias of halos of mass M at redshift z (which we will later consider to be scale-dependent) and P lin ( k) is the linear matter power spectrum, as defined in Section II A. At this point, it is worth emphasizing that expressions analogous to Eqs. (A6)-(A8) can be written for any field defined at all points in the universe after its profile for each halo of mass M is specified. These expressions are generic consequences of the halo model. The primary assumption made is that the halo-halo power spectrum for halos of mass M 1 and M 2 is given by the linear matter power spectrum multiplied by the relevant bias parameters: (A10)

The One-Halo Term
We now compute the contribution of the one-halo term to the tSZ power spectrum. We do the exact calculation first, and then consider the flat-sky limit. Note that for the one-halo term, the notion of the "Limber approximation" is not well-defined -there are no LOS cancellations to consider, since one halo is by definition fixed at a single redshift. Thus, in going from the exact calculation to the small-angle (high-) limit, we only need to consider the projection of the electron pressure profile from 3D to 2D. Note that this generalization will only affect very massive, low-redshift clusters, which subtend a significant solid angle on the sky.
Consider a cluster of virial mass M at comoving separation χ with respect to our location. The 3D Compton-y field due to this cluster at comoving separation r with respect to the cluster center is given by: This expression is simply the inverse transform of Eq. (A9). Projecting along the LOS as in Eq. (A4) and using the Rayleigh plane wave expansion, we obtain: wheren is a unit vector on the sky andk is the direction of k. Defining the expansion coefficients y m (M ) via we can read them off from Eq. (A12): where we have again used the Rayleigh expansion and have used the orthonormality of the spherical harmonics to do the integral overk in going from the first line to the second. From this expression, we can read off the exact result for the 2D Fourier transform of the projected y-profile due to a cluster of mass M at redshift z: where we have rewritten the spherical Bessel functions in terms of Bessel functions of the first kind using j ν (x) = π 2x J ν+1/2 (x) and we have explicitly included a possible dependence of the y 3D profile on redshift (in addition to mass).
The total one-halo term in the tSZ power spectrum is then given by the sum of the individual contributions from every cluster in the universe: where d 2 V /dzdΩ = cχ 2 (z)/H(z) is the comoving volume element per steradian. Substituting Eq. (A15) into this expression, converting the χ integral to a redshift integral, and rearranging the order of the integrals then yields the final result for the exact one-halo term: (A17) Note that this expression is exact: no flat-sky approximation (or any other) has been used in deriving Eq. (A17).
In order to recover the flat-sky (i.e., small-angle) limit of Eq. (A17), we use the following → ∞ limit for the spherical Bessel functions: Applying this limit to the first line of Eq. (A15) yields Using Eq. (A9), we can simplify this expression into a familiar form: where we have performed the following change of variables in the integral over y 3D : x ≡ a(z)r/r s , where r s is a characteristic scale radius of the y 3D profile. Finally, s = a(z)χ(z)/r s = d A (z)/r s is the characteristic multipole moment associated with the scale radius, with d A (z) the angular diameter distance. Note that the change of variables involved the scale factor because we transformed from comoving coordinates to physical coordinates. Eq. (A20) is identical to the quantityỹ (M, z) defined in (for example) Eq. (2) of [19], although we have explicitly used + 1/2 rather than . This is both technically correct and reduces the error in the approximation from O( −1 ) to O( −2 ) [135]. Eq. (A20) is simply the flat-sky limit of Eq. (A15). Following the long-standing convention, we will use the same definition as that established in [19]: The flat-sky limit of the one-halo term given in Eq. (A17) is thus given by as written down in (for example) Eq. (1) of [19]. Evaluating Eq. (A17) numerically is somewhat computationally expensive, as it contains five nested integrals (including the Fourier transform to obtainỹ 3D ), two of which involve highly oscillatory Bessel functions. However, the flat-sky limit in Eq. (A22) contains only three nested integrals, and involves no oscillatory functions. Furthermore, for our fiducial cosmology, we find that the flat-sky result in Eq. (A22) only overestimates the exact result in Eq. (A17) by ≈ 13%, 5%, and 3% at = 2, 10, and 20, respectively. By = 60, the two results are identical within our numerical precision. In addition, at = 2 where the correction is largest, the one-halo term is only ≈ 67% as large as the twohalo term, and thus the total tSZ power spectrum is only overestimated by ≈ 5%. Note that for non-Gaussian cosmologies this overestimate is far smaller, because the two-halo term dominates by a much larger amount at lowthan in a Gaussian cosmology (for example, the two-halo term at = 2 is 2.2 times as large as the one-halo term for f NL = 50). Given the small size of this correction and the significant computational expense required to evaluate the exact expression, we thus use the flat-sky result in Eqs. (A21) and (A22) to compute the one-halo contribution to the tSZ power spectrum in this work.
case, and thus we neglect this small error in our results (the constraints on all other parameters are insensitive to moderate changes in the errors at the lowest few values). Moreover, in the masked cases (which present the greatest promise for cosmological constraints), the trispectrum contribution is heavily suppressed at low-(see Figs. 19,20,and 21), and the total errors are dominated by the Gaussian term. Thus, for the masked cases, the exact vs. flat-sky correction should be vanishingly small. Note that we neglect the two-halo, three-halo and four-halo contributions to the trispectrum, as it is dominated even more heavily by the Poisson term than the power spectrum is [59]. The two-halo term will contribute to some extent at low-, but it is unlikely that higher-order terms will be significant even in this regime. For the masked calculations, the two-halo term may be somewhat important, though likely only at very low-, where, as we have argued above, it seems we do not need percent-level accuracy on the errors (since the forecasts for f NL are not particularly promising, and it is the only parameter very sensitive to this region of the power spectrum).
The full covariance matrix of the tSZ power spectrum, M y is then given by [19]: M y ≡ (C y,obs − C y )(C y,obs − C y ) where the angular brackets denote an ensemble average, f sky is the sky fraction covered by a given experiment (we assume f sky = 0.7 throughout this paper), N is the power spectrum due to instrumental noise after multifrequency subtraction (computed for Planck and PIXIE in Section IV), and we approximate T y ≈ T y,1h 1 . Note that Eq. (A28) does not include the so-called "halo sample variance" term, as discussed in Section V 1, as this term is negligible for a (nearly) full-sky survey.
We can then compute the covariance matrix Cov(p i , p j ) for the cosmological and astrophysical parameters of interest p i = Ω b h 2 , Ω c h 2 , Ω Λ , σ 8 , n s , f NL , M ν , C P0 , C β : where summation over the repeated indices is implied. The Fisher matrix for these parameters is then simply given by the inverse of their covariance matrix: The Fisher matrix encodes the constraining power of the tSZ power spectrum on the cosmological and astrophysical parameters.

Appendix B: Comparison with Planck Results
While this manuscript was under review, the Planck team released its initial set of cosmological results, including the construction of a Compton-y map and estimation of the tSZ power spectrum from this map [137]. In this Appendix, we provide a brief comparison of the publicly released Planck results with the forecasts in our work. Based on the analysis presented in Fig. 16, the Planck tSZ power spectrum should be signal-dominated over roughly the multipole range 100 < ∼ < ∼ 1500. Comparing with Fig. 15 in [137], this prediction is in very good agreement. The bandpowers and associated error bars presented in Table 3 of [137] imply a detection of the tSZ power spectrum with SNR ≈ 12.3, assuming a diagonal covariance matrix (no off-diagonal terms are presented in the Planck results). Our basic Planck forecast predicts a SNR ≈ 35. There are several reasons behind the difference in SNR between our forecast and the initial Planck result: • The usable fraction of sky in the Planck analysis (f sky ≈ 0.5) is found to be somewhat lower than that used in our analysis (f sky = 0.7) -this is primarily due to heavier masking of Galactic dust contamination than we anticipated; • The number of frequency channels used in the Planck analysis (six HFI channels only) is lower than that used in our analysis (all nine of the HFI and LFI channels); • The Planck analysis explicitly accounts for uncertainties in the contributions of the relevant foreground components (clustered CIB, IR point sources, and radio point sources) to the derived tSZ power spectrum, and finds that these uncertainties dominate the overall errors on the power spectrum; in our analysis, we have used reasonable models for the foregrounds to compute their contributions to the tSZ power spectrum, but have not explicitly propagated through uncertainties in these models to the final error bars. Our choice on this issue is partly driven by the fact that it is hard to quantify these uncertainties -in the Planck analysis, simulations are used to provide an estimate of the amplitude of each residual spectrum in the derived Compton-y power spectrum, but a 50% uncertainty remains. This uncertainty dominates the derived errors on the tSZ bandpowers, which is likely the main reason the Planck analysis SNR is significantly lower than our forecast.
In addition to these differences, we also note that the Planck analysis does not consider the possibility of masking nearby, massive clusters to reduce the sample variance in the tSZ power spectrum -however, it appears that the angular trispectrum contribution to the covariance matrix may not have been included in the Planck analysis at all, in which case masking would not be relevant. Regardless, this is another difference between our forecasts and the Planck results. Ultimately, the Planck analysis constrains σ 8 (Ω m /0.28) 3.2/8.1 = 0.784 ± 0.016. These constraints are obtained in a ΛCDM framework with all other cosmological parameters fixed, using the pressure profile of [119] with a hydrostatic mass bias of 20%, and without including the Planck constraints from the primordial CMB temperature power spectrum. In addition, the amplitudes of the foreground contributions to the tSZ power spectrum are allowed to vary, and are included as nuisance parameters. Given that this framework is rather different from ours, it is difficult to compare directly our forecasted parameter constraints with those obtained in the Planck analysis. Regardless, it is clear that the tSZ power spectrum is a useful cosmological probe, especially of the low-redshift amplitude of fluctuations, provided that uncertainties related to cluster gas physics and foreground contamination are treated with care.