An Accurate Analytic Model for the Thermal Sunyaev-Zel'dovich One-Point PDF

TL;DR

This paper introduces a refined analytic halo model for the thermal Sunyaev-Zel'dovich one-point PDF, incorporating halo overlaps and clustering effects, validated against simulations, to enhance cosmological analyses with upcoming CMB data.

Contribution

It presents a new, more accurate analytic model for the tSZ PDF that accounts for halo overlaps and clustering, improving upon previous simplified models.

Findings

Model matches numerical simulations accurately.

Necessary for analysis in upcoming CMB experiments.

Potential applicability to other cosmological observables.

Abstract

Non-Gaussian statistics of late-time cosmological fields contain information beyond that captured in the power spectrum. Here we focus on one such example: the one-point probability distribution function (PDF) of the thermal Sunyaev-Zel'dovich (tSZ) signal in maps of the cosmic microwave background (CMB). It has been argued that the one-point PDF is a near-optimal statistic for cosmological constraints from the tSZ signal, as most of the constraining power in tSZ $N$-point functions is contained in their amplitudes (rather than their shapes), which probe differently-weighted integrals over the halo mass function. In this paper, we develop a new analytic halo model for the tSZ PDF, discarding simplifying assumptions made in earlier versions of this approach. In particular, we account for effects due to overlaps of the tSZ profiles of different halos, as well as effects due to the clustering of halos. We verify the accuracy of our analytic model via comparison to numerical simulations. We demonstrate that this more accurate model is necessary for the analysis of the tSZ PDF in upcoming CMB experiments. The novel formalism developed here may be useful in modeling the one-point PDF of other cosmological observables, such as the weak lensing convergence field.

Figures

Equations

$$y(\mathbf{n})=\frac{{\sigma }_{\mathrm{T}}}{m_e{c}^2}{\int }_{\mathrm{LOS}}dr{n}_e(r,\mathbf{n})k_{\mathrm{B}}{T}_e(r,\mathbf{n}),$$

$$g(\nu )=x\coth (x/2)-4.$$

$$\frac{\Delta T_{\mathrm{tSZ}}(\mathbf{n},\nu )}{T_{\mathrm{CMB}}}=g(\nu )y(\mathbf{n}).$$

$$p_i={\int }_{y_i}^{y_{i+1}}dyP(y).$$

$$p_i=\int dzdM\frac{{\chi }^2}{H}\frac{dn}{dM}\pi \left({\theta }^2(y_i)-{\theta }^2(y_{i+1})\right)+{\delta }_i(1-F_{\mathrm{clust}}),$$

$$P_{200\mathrm{c}}=\frac{200G{M}_{200\mathrm{c}}{\rho }_{\mathrm{c}}(z){\Omega }_{\mathrm{B}}}{2{\Omega }_{\mathrm{M}}{r}_{200\mathrm{c}}},$$

$$\frac{P_{\mathrm{th}}(x)}{P_{200\mathrm{c}}}={\Pi }_0\frac{(x/x_{\mathrm{c}}{)}^{\gamma }}{(1+(x/x_{\mathrm{c}}{)}^{\alpha }{)}^{\beta }}.$$

$${\Pi }_0(\tilde{M},z)$$

$$x_{\mathrm{c}}(\tilde{M},z)$$

$$\beta (\tilde{M},z)$$

$$\tilde{P}(\lambda )\equiv \int dyP(y)e^{i\lambda y}.$$

$$\frac{dn}{d\Omega }=\frac{{\chi }^2(z)}{H(z)}\frac{dn(M,z)}{dM}dMdz$$

$$\tilde{P}(\lambda)=\big{\langle}\exp(i\lambda y(\mathbf{n}))\big{\rangle},$$

$$\tilde{P}(\lambda)=\Big{(}1-\frac{dn}{d\Omega}\pi\theta_{\rm max}^{2}\Big{)}+\frac{dn}{d\Omega}\int_{0}^{\theta_{\rm max}}d\theta\,2\pi\theta\,e^{i\lambda y(\theta)},$$

$$\tilde{Y}(M,z,\lambda )\equiv \int d\theta 2\pi \theta \left(e^{i\lambda y_0(M,z,\theta )}-1\right),$$

$$\tilde{P}(\lambda )=1+\frac{{\chi }^2(z)}{H(z)}\frac{dn(M,z)}{dM}\tilde{Y}(M,z,\lambda )dMdz.$$

$$\tilde{P}(\lambda )=\exp \left(\frac{{\chi }^2(z)}{H(z)}\frac{dn(M,z)}{dM}\tilde{Y}(M,z,\lambda )dMdz\right)$$

$$P(y)=\underset{N\to \infty }{\lim }{P}_{M_1,z_1}(y)\otimes \cdots \otimes P_{M_N,z_N}(y).$$

$$\tilde{P}(\lambda )$$

$$=\prod _{i=1}^{\infty }\exp \left(\frac{{\chi }^2(z_i)}{H(z_i)}\frac{dn(M_i,z_i)}{d{M}_i}\tilde{Y}(M_i,z_i,\lambda )d{M}_{i}d{z}_i\right)$$

$$=\exp \left(\sum _{i=1}^{\infty }\frac{{\chi }^2(z_i)}{H(z_i)}\frac{dn(M_i,z_i)}{d{M}_i}\tilde{Y}(M_i,z_i,\lambda )d{M}_{i}d{z}_i\right)$$

$$=\exp \left(\int \frac{{\chi }^2(z)}{H(z)}\frac{dn(M,z)}{dM}\tilde{Y}(M,z,\lambda )dMdz\right)$$

$${\tilde{P}}_{\delta }(\lambda ,\mathbf{n})$$

$$\displaystyle\quad\times\frac{dn(M,z)}{dM}[1+b(M,z)\delta_{\rm lin}(\mathbf{n},z)]\Big{\}}.$$

$${\tilde{P}}_{\mathrm{cl}}(\lambda )=⟨{\tilde{P}}_{\delta }(\lambda ,\mathbf{n}){⟩}_{\delta (\mathbf{n},z)}.$$

$$A(\lambda ,\mathbf{n})$$

$$\alpha (z,\lambda )$$

$${\tilde{P}}_{\delta }(\lambda ,\mathbf{n})={\tilde{P}}_u(\lambda )e^{A(\lambda ,\mathbf{n})}.$$

$${\tilde{P}}_{\mathrm{cl}}(\lambda )={\tilde{P}}_{\mathrm{u}}(\lambda )\exp \frac{1}{2}⟨A^2(\lambda ,\mathbf{n}){⟩}_{\delta }.$$

$$⟨A^2(\lambda ,\mathbf{n}){⟩}_{\delta }=\int dzH(z)D^2(z){\alpha }^2(z,\lambda )\times \int \frac{kdk}{2\pi }{P}_{\mathrm{lin}}(k).$$

$${\tilde{P}}_{\mathrm{cl}}(\lambda )$$

$$\displaystyle\hskip 71.13188pt\times\int\frac{k\,dk}{2\pi}P_{\rm lin}(k)\bigg{)}$$

$$F_{\ell }=C_{\ell }^{yy}/(C_{\ell }^{yy}+N_{\ell }^{yy}),$$

$${\sigma }_{yy}^2=\sum _{\ell }\frac{2\ell +1}{4\pi }{N}_{\ell }^{yy}{F}_{\ell }^2{p}_{\ell }^2,$$

Full text

I Introduction

Cosmological inference has traditionally focused on measurements of the power spectrum (or its real-space analogue, the two-point correlation function). For a Gaussian random field, this approach is sensible, as the power spectrum contains all statistical information in the data. The primary cosmic microwave background (CMB) temperature and polarization anisotropies are canonical examples of such Gaussian fields Planck2018params ; Planck2015NG ; WMAP9params . However, although the initial conditions for cosmic structure formation captured in the CMB are consistent with Gaussianity, non-Gaussian features inevitably develop in the late-time universe, due to non-linear gravitational evolution and complex baryonic physics on small scales. Thus, the information content of late-time cosmological data sets, e.g., weak gravitational lensing maps or maps of the galaxy distribution, is not completely captured by the power spectrum. For highly non-Gaussian fields, the amount of additional information in higher-order statistics can be significant.

In Ref. Hill2014ACTPDF (hereafter H14), it was argued that the one-point probability distribution function (PDF) is a near-optimal non-Gaussian statistic for cosmological inference from maps of the thermal Sunyaev-Zel’dovich (tSZ) effect. The tSZ effect is the up-scattering of CMB photons to higher energies due to Thomson scattering off hot, free electrons, which produces a unique distortion in the energy spectrum of the CMB ZS1969 ; SZ1970 . The tSZ effect probes the integrated pressure of free electrons along the line-of-sight (LOS); thus it is a biased tracer of free electrons, due to its dependence on the product of the electron number density and temperature. In particular, because of this temperature dependence, the tSZ signal is predominantly sourced by electrons in galaxy groups and clusters, where electrons are virialized to high temperatures. As the distribution of such objects in our Hubble volume is nearly Poissonian, maps of the tSZ effect are extremely non-Gaussian, dominated by individual rare, bright sources. H14 proposed the tSZ one-point PDF as an efficient statistic with which to extract the information in this non-Gaussianity.

Non-Gaussian properties of the tSZ signal have long been used for cosmological constraints, in the more familiar guise of inferring parameters via measurements of the halo mass function, in which candidate clusters are identified via the tSZ effect, their existence is confirmed and redshifts are estimated via multi-wavelength follow-up observations, and their masses are estimated via follow-up observations and/or scaling relations (e.g., Hasselfield2013 ; Planck2015clusters ; Bocquet2018 ). However, this approach utilizes only the brightest tSZ objects in the sky, e.g., those with signal-to-noise ($S/N$) greater than some threshold in the map. An alternative strategy has thus been developed in the past two decades, in which “indirect” statistics of the tSZ signal are directly utilized as cosmological probes, such as the power spectrum (e.g., Komatsu-Seljak2002 ; Hill-Pajer2013 ; McCarthy2014 ; Horowitz-Seljak2017 ; Sievers2013 ; George2015 ; Planck2013ymap ; Planck2015ymap ), bispectrum/skewness Wilson2012 ; Bhattacharya2013 ; Crawford2014 ), and cross-power spectra with gravitational lensing maps (e.g., HS2014 ; vanWaerbeke2014 ; Hojjatti2017 ). In such applications, no individual galaxy clusters are identified; instead, these statistics utilize information in objects below the $S/N$ threshold for individual detection, modeling their properties at the ensemble level. In particular, it has been shown that tSZ statistics beyond the power spectrum contain a significant amount of cosmological information beyond that contained in cluster counts (for current $S/N$ thresholds) or the power spectrum Wilson2012 ; Bhattacharya2013 ; Hill-Sherwin2013 ; Hill2014ACTPDF .

A unifying feature of these statistical tSZ analyses is that their cosmological constraining power arises almost entirely from the one-halo term. In other words, these statistics are just indirect methods of counting halos, weighted in different ways.222Thermal SZ – gravitational lensing cross-correlations are an exception to this statement, but even these statistics have only moderate sensitivity to the two-halo term in current data HS2014 ; Hojjatti2017 . For example, the tSZ power spectrum is dominated by the one-halo term at all $\ell\gtrsim 50$ Komatsu-Kitayama1999 ; Komatsu-Seljak2002 ; Hill-Pajer2013 . Thus, rather than deriving cosmological constraints from spatial clustering information, these statistics do so via the halo mass function. In particular, due to the tSZ signal’s bias toward electrons in high-mass (i.e., high-temperature) halos, these statistics probe the exponential tail of the mass function, which makes them very sensitive probes of the amplitude of fluctuations, $\sigma_{8}$ (e.g., Komatsu-Seljak2002 ; Wilson2012 ; Hill-Sherwin2013 ; Bhattacharya2013 ). Importantly, this sensitivity to $\sigma_{8}$ is almost entirely encoded in the amplitude of these statistics, rather than their shape; due to the one-halo term’s dominance (as mentioned above), the shape encodes information about intracluster medium (ICM) physics (and weak dependence on non-$\sigma_{8}$ cosmological parameters), while the amplitude is directly connected to integrals over the halo mass function, and hence $\sigma_{8}$. This underlies the argument presented in H14 that the one-point PDF is an optimal statistic for cosmological inference from the tSZ signal. The one-point PDF effectively captures the information in the amplitude of all $N$-point functions (or zero-lag moments), at the expense of information contained in the shape of the $N$-point functions (we anticipate that some shape information could be restored by considering the PDF on multiple smoothing scales, but such issues are not the focus of this paper). While not optimal for constraining ICM parameters, the tSZ PDF does allow for the breaking of degeneracies between these parameters and $\sigma_{8}$, due to the different dependence of each moment on these parameters and $\sigma_{8}$ (this is a generalization of the argument for parameter degeneracy-breaking using the tSZ two- and three-point functions presented in Ref. Hill-Sherwin2013 ).

The improved cosmological constraining power of the tSZ PDF over other tSZ statistics was demonstrated in practice in H14, which analyzed data from the Atacama Cosmology Telescope (ACT) at 148 GHz (using some 218 GHz data for foreground control as well). Compared to an analysis of the same data using the tSZ skewness Wilson2012 , the error bar on $\sigma_{8}$ was decreased by a factor of two in H14 (yielding $\sigma_{8}=0.793\pm 0.018$), simply due to the improved cosmological sensitivity of the tSZ PDF over the skewness alone. However, the $S/N$ of the measurement was not high enough to simultaneously constrain cosmological and ICM parameters (though the latter were marginalized over in obtaining the final $\sigma_{8}$ constraint). The Planck Collaboration subsequently applied the PDF statistic to an analysis of their component-separated tSZ map, obtaining $\sigma_{8}=0.77\pm 0.02$ Planck2015ymap .

However, the theoretical modeling that was developed for the tSZ PDF in H14 made several simplifying assumptions, limiting its utility in upcoming measurements (its sufficiency for the analysis of the ACT data in H14 was verified explicitly using end-to-end simulations). The analytic halo model of H14 assumed that clusters were Poisson-distributed on the sky and did not overlap, allowing the tSZ PDF to be computed via a simple integral over the mass function, given a model for the tSZ profile of each halo. The primary goal of this paper is to remove these simplifying assumptions and generalize the analytic model from H14, thereby allowing its use in analyses of the tSZ PDF from ongoing and upcoming CMB experiments (e.g., Advanced ACT Henderson2016 , SPT-3G Benson2014 , Simons Observatory SO2018 , and CMB-S4 CMBS4ScienceBook ). The assumptions are related to the distribution of halos sourcing the tSZ signal. In H14, it was assumed that these halos were sufficiently rare that they never overlapped on the sky. For massive clusters, this assumption is valid, but as the tSZ signal of progressively lower mass halos is included in the PDF, this assumption breaks down. For an experiment with relatively high noise levels (e.g., $\gtrsim 20\,\mu{\rm K}$-arcmin with $\sim$arcmin-scale beams), the assumption is valid, since low-mass clusters are subsumed into the noise. However, current and upcoming CMB experiments have noise levels well below this threshold, necessitating an improved model.

The other assumption from H14 that we will discard in this analysis is the neglect of halo clustering effects. These effects are relevant due to the LOS projection inherent in the tSZ signal. Note that the one-point PDF of 3D cosmic fields (i.e., the “voxel” PDF) does not receive any clustering contributions: as long as halo exclusion is enforced, only the one-halo term is necessary to compute the 3D PDF in the formalism used in this paper. For projected 2D fields, however, this is not true. Due to halo clustering, there is an excess probability for two halos to overlap along the LOS, compared to the Poisson expectation. For the tSZ field, we will find (Fig. 1 below) that the clustering effect is relatively weak, but for future extensions of this formalism to the weak lensing convergence field, we expect that it will be significant.333See, e.g., Refs. Liu2016 ; Patton2017 ; Liu-Madhavacheril2018 for simulation-based analyses of the weak lensing one-point PDF. Note that another important difference between the tSZ and weak lensing cases, which will require further theoretical work, is the existence of negative-signal regions in the latter (cosmic voids), whereas the tSZ effect is strictly positive.

The remainder of this paper is organized as follows. In §II, we review the halo model formalism of H14 and the associated assumptions, before proceeding to generalize the model in §III and compare the results to those obtained in the simpler approach. In §IV, we compare results from our analytic halo model to those obtained from numerical simulations, demonstrating its validity and accuracy. In §V, we use the analytic model to investigate the physical origin of the tSZ PDF signal. §VI presents the cosmological and ICM parameter dependence of the tSZ PDF. We then include noise and non-tSZ foregrounds in §VII, and demonstrate the sufficiency of our new model for the analysis of upcoming, high-precision CMB data sets. We conclude in §VIII.

Our fiducial cosmology is flat $\Lambda$CDM with dimensionless Hubble constant $h=0.7$, matter density $\Omega_{\rm M}=0.25$, spectral index $n_{s}=0.96$, $\sigma_{8}=0.8$, baryon density $\Omega_{\rm B}=0.043$, sum of the neutrino masses $\Sigma m_{\nu}=0$ eV, and CMB temperature $T_{\rm CMB}=2.726\,{\rm K}$.

II Background and Previous Model

H14 introduced a novel, simple analytic model for the tSZ one-point PDF. Here we review this model and its assumptions, laying the groundwork for the more accurate model derived in the following section.

The tSZ signal is quantified by the Compton-$y$ parameter, which measures the integrated pressure of free electrons along the LOS:

[TABLE]

where $\sigma_{\rm T}$ is the Thomson cross-section, $m_{e}c^{2}$ is the electron rest-mass energy, $r$ is physical distance along the LOS, $n_{e}$ is the electron number density, $k_{\rm B}$ is Boltzmann’s constant, and $T_{e}$ is the electron temperature.

The tSZ effect produces a non-blackbody distortion in the energy spectrum of the CMB, which is negative (positive) at frequencies below (above) $\approx 218$ GHz. Defining the dimensionless frequency $x\equiv h_{\rm Pl}\nu/(k_{\rm B}T_{\rm CMB})$, where $h_{\rm Pl}$ is Planck’s constant and $\nu$ is the photon frequency, the tSZ spectral function is given by:

[TABLE]

The tSZ-induced fluctuations in the CMB temperature field are then given by:

[TABLE]

In this work, we neglect relativistic corrections to the tSZ effect (e.g., Nozawa2006 ), which are important for massive, hot clusters; these effects must be included in an actual data analysis (as they were in H14). Throughout the rest of the paper, our results will generally be given either in terms of Compton-$y$ or in terms of the tSZ field at a reference frequency of 148 GHz, where $g(148\,{\rm GHz})=-0.97881$. We will also typically denote the tSZ-induced temperature fluctuation as simply $T\equiv\Delta T_{\rm tSZ}$.

We denote the (differential) tSZ one-point PDF as $P(y)$. The binned version of the PDF used in much of this work is given by

[TABLE]

The concept underpinning the H14 model is to note that $p_{i}$ quantifies the fraction of sky subtended by Compton-$y$ values in the range $[y_{i},y_{i+1}]$. Thus, for a single spherically symmetric halo with an azimuthally symmetric projected $y$-profile $y(\theta)$, this would correspond to the area in the annulus between $\theta(y_{i})$ and $\theta(y_{i+1})$, where $\theta(y_{i})$ is the angular distance from the center of the halo to the radius where $y(\theta)=y_{i}$. If one then makes the approximation that halos sourcing the tSZ signal are sufficiently rare that they never overlap on the sky, the final result for the full tSZ PDF is simply given by adding up such annular area contributions from all halos:

[TABLE]

where $\chi(z)$ is the comoving distance to redshift $z$, $H(z)$ is the Hubble parameter, $dn(M,z)/dM$ is the halo mass function (i.e., the number of halos of mass $M$ at redshift $z$ per unit mass and comoving volume), and $\theta(y,M,z)$ is the inverse function of $y(\theta,M,z)$ (i.e., the Compton-$y$ profile of a halo of mass $M$ at redshift $z$), $F_{\rm clust}$ is the total sky area subtended by halos (assuming some radial cutoff), $\delta_{i}$ is unity if $y=0$ lies in the bin and zero otherwise, and redshift and mass dependences have been suppressed in the equation for compactness.

Eq. (5) makes two strong assumptions: (1) halos sourcing the tSZ signal are rare enough that their projected $y$-profiles never overlap; (2) the clustering of these halos (which would make overlaps more likely) can also be neglected. These assumptions were valid for the analysis of ACT data in Hill2014ACTPDF , where the noise level was sufficiently high that the PDF could be modeled considering only halos for which these approximations are true (due to these halos being massive and hence rare). The remainder of this paper is focused on discarding these assumptions, thus yielding a more accurate and general model for analysis of the tSZ PDF in ongoing and upcoming CMB experiments. We will not focus on the modeling of non-tSZ foregrounds and noise (which were considered in detail in H14, and must be in any future analysis as well), but rather only on the modeling of the physical tSZ PDF signal.

In this work, we adopt the same models for the physical quantities underlying the tSZ PDF as used in H14. We compute the halo mass function $dn/dM$ and halo bias $b(M,z)$ using the fitting functions of Ref. Tinker2010 . We compute electron thermal pressure profiles using the fitting function given by Ref. Battaglia2012 (hereafter B12), in order to facilitate direct comparison with their hydrodynamical simulations (§IV.2). For convenience, we give the relevant formulae here. The thermal gas pressure at $r=r_{200\rm c}$ is given by

[TABLE]

with $\rho_{\rm c}(z)$ the critical density. Defining $x\equiv r/r_{200\rm c}$ and the core scale length $x_{\rm c}$, the pressure profile is parameterized as

[TABLE]

While $\alpha$ and $\gamma$ are held fixed, the remaining parameters are taken to have the following mass and redshift dependence:

[TABLE]

where $\tilde{M}\equiv M_{200\rm c}/10^{14}M_{\odot}$. We take the fiducial values $\alpha^{\rm fid}=1$, $\gamma^{\rm fid}=-0.3$, $P_{0}^{\rm fid}=18.1$, $x_{\rm c,0}^{\rm fid}=0.497$, and $\beta_{0}^{\rm fid}=4.35$. In §VI we consider the effect of changes in $\alpha$, $P_{0}$, and $\beta_{0}$ on the one-point PDF. All masses in the remainder of this work are given in terms of $M\equiv M_{200\rm m}$ unless otherwise stated. When needed, we convert to $M_{200\rm c}$ using the NFW profile NFW1997 and the concentration-mass relation of Ref. Duffy2008 , making use of the Colossus package Colossus . We note that this concentration model is calibrated in a mass range narrower than our integration boundaries, but since the contributions from the low- and high-mass ends are generally small as shown in §V, we do not expect this to be a significant source of error. In our fiducial result we apply a radial cutoff to the $y$-profiles at $r_{\rm out}=2\,r_{\rm vir}$, where we define the virial radius with a redshift dependent overdensity approximated according to Ref. BryanNorman1998 . We choose integration boundaries such that $10^{11}\leq M\,[h^{-1}M_{\odot}]\leq 10^{16}$, and $0.005\leq z\leq 6$. With these boundaries all integrals are converged.

III New Model

In this section, we derive the main result of this paper: an analytic model for the one-point PDF $P(y)$, making no simplifying assumptions about non-overlapping or non-clustering properties of halos. We will do this in three steps. First, we calculate the PDF for a narrow bin in mass $M$ and redshift $z$, such that overlaps can be neglected (§III.1). Second, we combine the contributions from all $(M,z)$-values accounting for overlaps, but with large-scale clustering of halos neglected (§III.2). Finally, we show how to include halo clustering (§III.3).

It is convenient to work in conjugate space by introducing the Fourier transform (FT) of the one-point PDF:

[TABLE]

Note that a calculation of $\tilde{P}(\lambda)$ is equivalent to a calculation of $P(y)$.444This approach bears similarities to the traditional $P(D)$ analysis used in radio point source studies Scheuer1957 ; Condon1974 ; however, our method accounts for tSZ sources’ non-trivial profiles, and the overlaps (and clustering) associated with these. We also note related work focused on the PDF of the tSZ power spectrum bandpowers Zhang-Sheth2007 .

III.1 PDF for a narrow mass-redshift bin

Let us consider now a narrow bin in halo mass and redshift of width $dM\,dz$ centered around mass $M$ and redshift $z$, in which the number of halos is sufficiently small that halo overlaps can be neglected (by “overlaps”, we mean overlaps of the projected $y$-profiles of these objects on the sky). We will calculate $\tilde{P}(\lambda)$ considering only halos in this bin.

We define the angular halo density in the narrow bin:

[TABLE]

and let $y_{0}(M,z,\theta)$ denote the $y$-profile of a halo with mass $M$ and redshift $z$, where $\theta$ is the angular distance from the halo center. We assume that the profile has a finite radius $\theta_{\rm max}$, i.e., $y_{0}(M,z,\theta)=0$ for $\theta>\theta_{\rm max}$. This is simply the projection of the radial cutoff defined in §II.

We write $\tilde{P}(\lambda)$ as an expectation value:

[TABLE]

where the sky location $\mathbf{n}$ is fixed, and the expectation value runs over random halo placements. Since we are neglecting overlaps, we can trade this expectation value with an integral over the angular profile:

[TABLE]

where the first term is the probability that no halo overlaps sky location ${\mathbf{n}}$, and the second term integrates over overlap locations. Introducing the auxiliary quantity

[TABLE]

we rewrite Eq. (14) in the form:

[TABLE]

For reasons that will be apparent in the next section, it will be convenient to write this as:

[TABLE]

which is equivalent for a differential bin $(dM\,dz)$.

In Appendix A, we show how to obtain the formalism used in H14 from Eq. (16), under the assumption that halo overlaps can be neglected.

III.2 PDF neglecting halo clustering

Now we calculate the one-point PDF from all masses and redshifts, accounting for halo overlaps, but neglecting large-scale halo clustering.

Suppose we define a large number of narrow mass-redshift bins $(M_{1},z_{1}),\cdots,(M_{N},z_{N})$. If halo clustering is neglected, then the total $y$-signal is the sum of independent contributions from each bin. Therefore, the complete PDF from all $N$ mass-redshift bins is then obtained by convolution:

[TABLE]

Taking Fourier transforms, the convolution becomes multiplication, and simplifies as follows:

[TABLE]

where we have used Eq. (17) in the second line.

III.3 Including Halo Clustering

Finally, we include the effect of halo clustering. In this subsection, we will denote the “unclustered” PDF found in Eq. (19) by $\tilde{P}_{\rm u}(\lambda)$, and denote the “clustered” PDF by $\tilde{P}_{\rm cl}(\lambda)$.

The derivation proceeds in two steps. First, we compute the one-point PDF $\tilde{P}_{\delta}(\lambda,\mathbf{n})$ in a fixed realization of the linear density field $\delta_{\rm lin}(\mathbf{n},z)$. Second, we average over realizations of the Gaussian field $\delta_{\rm lin}$ to obtain $\tilde{P}_{\rm cl}(\lambda)$. (Note that $\tilde{P}_{\delta}$ depends on $\mathbf{n}$, since translation invariance is broken by a particular realization of $\delta_{\rm lin}$, but $\tilde{P}_{\rm cl}$ is a translation-invariant PDF as usual.)

The quantity $\tilde{P}_{\delta}(\lambda,\mathbf{n})$ can be obtained from Eq. (19) for the unclustered PDF, by biasing the halo mass function with the halo bias $b(M,z)$, i.e.,

[TABLE]

We note that this expression is only meaningful as long as we can define a sufficiently large environment around $\mathbf{n}$ in which $\delta_{\rm lin}(\mathbf{n}+\mathbf{n}^{\prime},z)\simeq\delta_{\rm lin}(\mathbf{n},z)$. This assumption is justified because $\delta_{\rm lin}(\mathbf{n},z)$ varies slowly in comparison to the typical cluster radius.

The PDF including clustering is obtained by averaging over realizations of the linear density field $\delta_{\rm lin}(\mathbf{n},z)$ in Eq. (III.3):

[TABLE]

Introducing

[TABLE]

we write Eq. (III.3) as:

[TABLE]

Using the identity $\langle e^{x}\rangle=e^{\langle x^{2}\rangle/2}$ for a Gaussian random variable $x$, we obtain:

[TABLE]

The expectation value $\langle A^{2}\rangle$ can be evaluated using the Limber approximation and the LOS integral representation in Eq. (22). The result is:

[TABLE]

where $P_{\rm lin}(k)$ is the linear matter power spectrum at $z=0$ and $D(z)$ is the growth factor with $D(0)=1$. This gives our final expression for the one-point PDF:

[TABLE]

where $\tilde{P}_{\rm u}(\lambda)$ is the “unclustered” PDF in Eq. (19).

III.4 Quantifying the Effects of Overlaps and Clustering

In the following, we denote the PDF integrated over a given temperature bin by $p_{T}$ and evaluate the temperature decrement corresponding to a specific $y$-signal at $\nu=148\,{\rm GHz}$. Unless otherwise stated, we bin the PDF into temperature bins of width $1\,\mu{\rm K}$. The linear power spectrum $P_{\rm lin}(k)$ and growth factor $D(z)$ are computed using CAMB CAMB . We briefly mention two numerical properties of our analytic method. If the integration boundaries (in mass and cluster radius) are chosen too narrow, the Fourier transform $\tilde{P}(\lambda)$ does not vanish for $\lambda\rightarrow\infty$. This gives rise to ringing in the PDF. Furthermore, the PDF at high $|T|$ is only properly converged if the $y$-profiles are evaluated on a very fine angular grid, which is due to their rapid variation for angles close to the clusters’ center.

We now turn to concrete results from our analytic model. In Fig. 1 we show the effects of clustering and overlaps. The fiducial result is obtained using Eq. (27), the result neglecting halo clustering is calculated from Eq. (19), and the result neglecting both overlaps and clustering follows from the formalism developed in H14 (Eq. 5). Overlaps have a much larger impact on the PDF, in particular for low $|T|$ values. These temperature bins are dominated by numerous low-mass halos, which have a larger probability to overlap. Clustering increases the PDF for almost all $|T|$ values, but it is less important than overlaps. This is explained by the fact that it only plays a role if the gravitational interaction aligns two nearby clusters along the line of sight, which has subdominant probability in comparison to random alignments on arbitrarily large scales. Clustering slightly decreases the PDF for $T>-4\,\mu\rm K$. These bins are dominated by very numerous low mass halos, which have a relatively high probability to gravitationally interact and produce an alignment, which would push their contribution to higher $|T|$ values. The lowest $|T|$ bin, on the other hand, sees an increase in the PDF due to clustering, which is clear because clustering increases the clear sky fraction $(1-F_{\rm clust})$. Note that the unphysical divergence of the H14 model in the lowest $|T|$ bin (arising from the neglect of overlaps) is removed by the new Fourier-based approach developed here.

IV Comparison to Simulations

We check the validity of our analytic approach by comparing to the results of two different simulation methods. First, we produce “simplified”random maps, in which unclustered halos are randomly distributed on a simulated flat-sky map and assigned Compton-$y$ profiles computed with the B12 pressure profile. We then measure the average tSZ PDF from these maps. By construction, the one-point PDF of our simplified simulations should agree perfectly with our “unclustered” analytic result in Eq. (19), but verifying the agreement is a strong check on details of the implementation, which are nontrivial (see §III.4).

Second, we measure the tSZ PDF from Compton-$y$ maps constructed directly from the cosmological hydrodynamics simulations of B12.555We thank N. Battaglia for sharing these maps. The simulated maps in this case only include tSZ signal from halos at $z<1$, and thus in this section we set the upper redshift integration boundary to $z=1$ in our analytic calculations, in order to enable direct comparison with the simulated maps. Note that our fiducial cosmology is identical to that used in B12.

IV.1 Simplified Simulations

The simulations described in this section are produced as follows. We construct individual maps of area $\Omega=10\times 10\,{\rm deg}^{2}$, with square pixels of side-length $3\,{\rm arcsec}$. We then consider discrete, narrow bins in mass and redshift of size $dM\,dz$ centered around $M,z$, for which we compute the average number of such halos in the map via $n(M,z)=\Omega\,dn(M,z)/d\Omega$. For each such mass-redshift bin, we populate the map with $y$-profiles (using the B12 pressure profile), whose number is given by the probability distribution $w(\lceil\!n\!\rceil)=n-\lfloor\!n\!\rfloor$, and $w(\lfloor\!n\!\rfloor)=\lceil\!n\!\rceil-n$. Since this distribution reproduces the correct mean, it is valid to use it to find the average tSZ PDF computed from many maps. We find that sampling the number of halos from the physically more realistic Poisson distribution leads to relatively slow convergence of the average PDF, but we have confirmed that it yields a consistent result with the more rapid approach. Note that we do not include halo clustering in these maps.

Fig. 2 shows the average tSZ PDF computed from 507 simplified simulations (dashed red curve). As the maps do not include clustering of halos, we compare the simulation-derived PDF to the analytic result from §III.2, in which clustering effects are not included (but halo overlaps are). The discrepancy with our analytic result is on average $\approx 0.2\%$ and decreases as more maps are added to the average. This confirms the validity of our analytic formalism, in the limit where halo clustering effects can be neglected.

IV.2 Cosmological Hydrodynamics Simulations

We now compare the results of the full analytic calculation presented in §III.3 (including halo clustering) to measurements of the tSZ PDF from cosmological hydrodynamics simulations. We use Compton-$y$ maps constructed by direct LOS integration (to $z=1$) of randomly rotated and translated simulated volumes from Ref. Battaglia2010 . These are the same simulations from which the B12 pressure profile fitting function was extracted; thus, the comparison here is a direct test of our analytic formalism for the tSZ PDF, with no additional tuning of ICM parameters required. The simulations were performed using the smoothed particle hydrodynamics code GADGET-2, with a custom implementation of a sub-grid prescription for feedback from active galactic nuclei. The pressure profile model extracted from the simulations has subsequently been found to agree with a wide range of tSZ and X-ray measurements (e.g., Sun2011 ; Planck2013stack ; Hajian2013 ; HS2014 ; Greco2015 ; Hilletal2018 ). A full description of the simulations can be found in Ref. Battaglia2010 .

We consider 390 Compton-$y$ maps extracted from the simulations, each of area $\Omega=4.09\times 4.09\,{\rm deg}^{2}$, with square pixels of side-length $6\,{\rm arcsec}$. The average tSZ PDF measured from this suite of $y$-maps is shown in Fig. 2 (blue curve with error bars). As mentioned above, an upper redshift cut at $z=1$ is applied in the construction of the maps, so as to minimize correlations arising from high-redshift objects common to multiple maps. We apply this redshift cut in the analytic calculation in this subsection (and only this subsection) for consistency with the simulated maps.

The agreement in Fig. 2 between the tSZ PDF measured from the hydrodynamical simulations and our analytic result is much worse than the agreement for the simplified simulations seen in the previous subsection. However, this difference can be traced back to the discrepancy between the halo mass function found in B12 and the mass function Tinker2010 used in this work. We explicitly confirm that this halo mass function difference can indeed give rise to discrepancies in the one-point PDF as observed in Fig. 2.

In Fig. 3 we show how different interpolations of the halo mass function (HMF) given in B12 (see their Fig. 11) affect the PDF. This is compared to the discrepancy found between our fiducial result and the PDF measured from the B12 simulations. Since only information for redshift $z=0$ is available in Fig. 11 of B12, we take the ratio between the interpolated halo mass function and the fitting function Tinker2008 as constant across redshift. Although this is likely a gross oversimplification, it can be seen from the figure that reasonable interpolations can already reproduce the observed discrepancy very well. We thus conclude that our analytic formalism passes this check, although future comparisons with additional hydrodynamical simulations will also be useful.

V Origin of the Signal

We now turn to the different contributions to the PDF. First, we consider cluster mass and redshift. In Fig. 4 we plot the absolute fractional deviation of the PDF as a function of the maximum mass and redshift included in the calculation. These results were obtained using the simplified simulations generated as described in §IV. It should be noted that care must be taken in interpreting these plots, since overlaps can shift part of the contribution from a certain $(M,z)$-bin to higher $|T|$ values. Furthermore, we note that direct comparison with the analogous results given in H14 is not possible, because their results included instrumental noise, non-tSZ foregrounds, and beam convolution.

Regarding the mass contributions, we note two general trends in Fig. 4. First, as $|T|$ increases, the transition region, i.e., the range of relevant masses for the specific temperature bin, gets smaller (ignoring the temperature bin that contains the clear sky fraction for now). At high $|T|$, the PDF in a given temperature bin is dominated by clusters in a narrow mass range. On the other hand, for low $|T|$ a variety of sources contributes. Second, the relevant masses are higher for high $|T|$, which is expected.

Regarding the redshift contributions, the most relevant interval broadens and shifts to larger redshift as $|T|$ decreases. For the temperature bin containing the clear sky fraction the behavior is drastically different, with a very small interval in mass and redshift being the dominant contribution.

Now we consider the effect of the radial cutoff. In our fiducial computation we considered the $y$-profile up to $r_{\rm out}=2\,r_{\rm vir}$. In Fig. 5 we show how the choice of this outer radius affects the PDF. The effect is largest for the low-$|T|$ regime. These temperature bins receive significant contributions from the outskirts of clusters. Furthermore, as discussed above, overlaps have their largest impact on these bins, and reducing the outer radius decreases the amount of overlaps. The clear sky fraction is increased as $r_{\rm out}$ is decreased, which is intuitively clear. We note that the choice of $r_{\rm out}=2\,r_{\rm vir}$ does not correspond to convergence, since there is still a significant discrepancy to the result obtained with $r_{\rm out}=3\,r_{\rm vir}$. However, it is physically not justified to suppose the validity of the pressure profile fitting function up to infinite radius. Physically, the virial shock will lead to a sharp decline in the pressure profile at $r\approx 2-2.5r_{200c}$. Upon convolution with instrumental noise (as described in §VII) the precise choice of radial cutoff becomes irrelevant to the PDF prediction, as the extremely small $y$ values in the outskirts are subsumed into the noise.

VI Parameter Dependence

We now turn to the dependence of the PDF on cosmological and pressure profile parameters. In Fig. 6 we plot the effect of changing the cosmological parameters $\sigma_{8}$, $\Omega_{\rm M}$, and $\Sigma m_{\nu}$, as well as the pressure profile parameters $P_{0}$, $\beta_{0}$, and $\alpha$ as defined in §II. When considering non-zero neutrino masses we take $N_{\rm eff}=3.046$.

The impact of changing the cosmology is as follows. The neutrino mass sum has its largest effect on high $|T|$ bins, consistent with the fact that it changes the matter power spectrum most on small spatial scales, which here leads to a suppression in the number of clusters (however, to hold $\sigma_{8}$ fixed, the initial scalar amplitude $A_{s}$ must be increased, thus leading to an overall increase in the number of clusters; this is simply an artifact of this choice of cosmological parameters). The impact of changing $\sigma_{8}$ is similar, yielding a mild degeneracy with $\Sigma m_{\nu}$. On the other hand, changing $\Omega_{\rm M}$ has largest impact on low $|T|$ bins. Thus, degeneracy between $\sigma_{8}$ and $\Omega_{\rm M}$ could be broken by a measurement of the tSZ PDF over a sufficiently large $|T|$ range.

The pressure profile parameters affect the PDF as follows. $P_{0}$ changes the PDF relatively constantly across temperature, which is explained by the fact that it is an overall normalization of the pressure profile. Changing the logarithmic slope at large radii, $\beta_{0}$, produces a similar effect, although with the opposite sign. Increasing the logarithmic slope at intermediate radii, $\alpha$, decreases the PDF in low $|T|$ bins and increases it in large $|T|$ bins.

VII Including noise and foregrounds

Here we consider a representative upcoming CMB experiment and demonstrate the necessity of our improved analytic approach for modeling the tSZ PDF sufficiently accurately to perform unbiased inference from this observable. In lieu of full parameter forecasts, which require treatment of the covariance matrix and likelihood function, we simply compute the predictions from our full analytic model (Eq. 27) and from the H14 model (Eq. 5), and compare these to predicted error bars on a measurement of the tSZ PDF.

Specifically, we consider a measurement with the Simons Observatory (SO), which will cover $f_{\rm sky}=0.4$ with six frequency channels, reaching a depth and resolution at 145 GHz of $6\,\mu$K-arcmin and FWHM = 1.4 arcmin, respectively SO2018 .666The sensitivity values considered here are for the “goal” configuration presented in Ref. SO2018 . We include the effects of non-tSZ foregrounds via the post-component-separation Compton-$y$ noise power spectra that are publicly available via the SO Collaboration.777https://simonsobservatory.org/assets/supplements/20180822_SO_Noise_Public.tgz In particular, we use the $y$ noise power spectrum, $N_{\ell}^{yy}$, derived via the “standard internal linear combination”, without additional deprojection constraints (see SO2018 for further details about the foreground modeling and component separation). In practice, deprojection options will have to be explored in the foreground cleaning as well, but this is beyond the scope of our work here.

We construct a Wiener filter to optimally weight the harmonic-space tSZ signal via

[TABLE]

where $C_{\ell}^{yy}$ is the tSZ power spectrum, which we compute using the fiducial model of Hill-Pajer2013 . The filter is smoothly tapered to zero at the boundaries of the multipole range provided in the $N_{\ell}^{yy}$ data file ($\ell_{\rm min}\approx 40$ and $\ell_{\rm max}\approx 8000$). We apply this filter to the $y$-profiles of all halos in our analytic calculation, which captures the suppression of modes lost due to foregrounds and noise. Note that the multifrequency information of SO (and Planck, which is also used) allows large-scale tSZ modes to be included (because the CMB can be removed using spectral information), which were lost due to CMB “noise” in the single-frequency ACT analysis of H14. Thus the filter extends to lower multipoles than in H14 (which used a filter originally constructed in Wilson2012 ).

After calculating the analytic prediction for the filtered tSZ PDF, we convolve the result with a Gaussian noise (+residual foreground) PDF, whose variance is computed via

[TABLE]

where $p_{\ell}$ is the pixel window function (here assumed to be 0.5 arcmin circular pixels, although this has negligible effect). We then rescale the results from Compton-$y$ to 148 GHz temperature to match those shown elsewhere in the paper (although the application of the filter and noise convolution mean the temperature values are not comparable to those in earlier plots). We bin the results into bins of width $5\,\mu$K.

The results are shown in Fig. 7. The solid blue curve shows the prediction from the FT-based formalism developed in this paper (Eq. 27). The dashed green curve shows the prediction from the H14 model (Eq. 5). The error bars shown on the blue curve are computed using the diagonal elements of the covariance matrix $C_{ij}$ estimated from the simplified simulations described in §IV.1, but with Poisson distributed cluster numbers, rather than the modified distribution described earlier (which achieved faster convergence at the expense of only capturing the mean correctly), so that the variance is correctly captured. Specifically, we compute the error on the value of the PDF in the $i$th bin as

[TABLE]

where $f_{\rm sky}$ are the sky fractions and $d_{\rm pixel}$ the pixel side lengths. We take $f_{\rm sky}^{\rm SO}=0.4$ and $d_{\rm pixel}^{\rm SO}=0.5\,{\rm arcmin}$; the simulation parameters are $f_{\rm sky}^{\rm maps}=6.21\times 10^{-3}$ and $d_{\rm pixel}^{\rm maps}=0.1\,{\rm arcmin}$. For a survey with the properties described above, the difference between the no-overlaps case and our fiducial result is considerably larger than the errors for essentially all bins plotted. We note that the error bars themselves should not be taken at face value, because significant bin-to-bin correlations exist as discussed in H14. Nevertheless, we conclude that if the earlier model of H14 were used in an analysis of the tSZ PDF from SO, cosmological and ICM parameter inference would clearly be biased. With our accurate model in hand, we plan to pursue full parameter forecasts for ongoing and upcoming CMB experiments in future work.

VIII Outlook

In this paper we have presented a new analytic model for the tSZ one-point PDF, building upon and substantially improving the model first developed in H14. In particular, by working in Fourier conjugate space, we have shown how to account for effects due to overlaps in the tSZ profiles of halos on the sky, as well as contributions due to the clustering of halos (which arise because of the LOS projection). For the tSZ PDF, the effects due to overlaps are non-negligible, but the clustering effects are rather small. We have verified the accuracy of the model via comparison to numerical simulations, both simplified simulations containing randomly distributed clusters and full-scale cosmological hydrodynamics simulations. However, issues related to the halo mass function in the latter simulations rendered a precise test of the clustering effects challenging; future simulation comparisons will thus be useful. Finally, we have demonstrated that the use of this more accurate analytic model will be necessary in analyses of the tSZ PDF in upcoming, high-sensitivity CMB data sets.

We anticipate a number of interesting next steps in this line of research. An obvious first step is to compute the covariance matrix in this formalism and the likelihood function associated with the PDF observable. Given the challenges observed in this context in H14, it may be useful to pursue novel approaches such as likelihood-free inference (although this could render the analytic model redundant) Alsing2018 . We expect that the forecast cosmological constraints using the tSZ PDF will significantly improve upon those for the tSZ power spectrum alone (e.g., as presented in Ref. SO2018 ). An optimal combination with constraints from individually detected clusters is clearly also a pressing issue, and will lead to further improvements.

Beyond the tSZ signal, the formalism developed here likely has applications to other cosmological fields. An obvious candidate is the one-point PDF of the weak lensing convergence field, which has already been investigated in simulations Liu2016 ; Patton2017 ; Liu-Madhavacheril2018 . We expect that the clustering effects computed in this paper will be more important for this application than for the tSZ field. In addition, further development to treat negative-convergence regions (voids) will be necessary. Nevertheless, a full, non-perturbative model for the one-point PDF of the projected density field is clearly a goal worth pursuing.

Acknowledgements.

We thank Nick Battaglia and David Spergel for helpful conversations. JCH is supported by the Friends of the Institute for Advanced Study. LT acknowledges support by the Studienstiftung des deutschen Volkes. KMS was supported by an NSERC Discovery Grant and an Ontario Early Researcher Award. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Research & Innovation. This is not an official SO Collaboration paper.

Appendix A Equivalence to Formalism of H14

In this appendix we show that our analytic model is equivalent to the formalism used in H14 under the assumption that no overlaps occur. Denote $\int dM\int dz(\chi^{2}/H)(dn/dM)$ by $\int_{M,z}$ for brevity. The arguments $M$ and $z$ are understood for $y_{0}(M,z,\theta)$, etc. Integrating Eq. (16) (which is equivalent to a first order expansion of Eq. 19) over mass and redshift, we obtain

[TABLE]

The one-point PDF in $y$-space is given by:

[TABLE]

where we denote the inverse function to $y_{0}(\theta)$ by $\theta_{0}(y)$. The PDF binned into $y_{i}\leq y\leq y_{i+1}$ is then found as

[TABLE]

where $\delta_{i}$ equals one if $y=0$ is contained in the integration interval and zero otherwise. $\theta_{\rm max}=\theta_{\rm max}(M,z)$ corresponds to the radial cutoff, so that $1-F_{\rm clust}$ is the clear-sky fraction. This is the expression used in H14, i.e., Eq. (5) presented earlier.