CMB-S4: Forecasting Constraints on $f_{\mathrm{NL}}$ Through $\mu$-distortion Anisotropy

In the limit of purely Gaussian-distributed primordial curvature perturbations $\zeta({\vec{k}})$, the auto-correlation of these perturbations is given by

The power spectrum $P_{\zeta}(k)$ of these perturbations from single-field inflation models is given by

Here we use best-fit numbers from the Planck 2018 data release [8]: $A_{\mathrm{s}}=2.1\times 10^{-9}$; and $n_{\mathrm{s}}=0.965$, reported for a pivot scale $k_{\mathrm{p}}=0.05\ \textrm{Mpc}^{-1}$.

In the squeezed limit, where the global power spectrum of two modes on short scales $k_{\mathrm{S}}$ are modulated by a long-wavelength mode $k_{\mathrm{L}}$, the long-wavelength modes generate a local non-Gaussianity and the power spectrum becomes position-dependent

One common measure of the level of non-Gaussianity in the curvature distribution is the Fourier-domain statistic known as the curvature bispectrum $B_{\zeta}(k_{1},k_{2},k_{3})$. Analogous to the power spectrum and the auto-correlation, the bispectrum is defined through its relationship to the three-point correlation function:

In the squeezed limit, ${k_{1}/k_{3}\rightarrow 0}$, the bispectrum can be expressed as in Eq. (1), with $f_{\mathrm{NL}}$ parameterizing the amount of local non-Gaussianity coupling the power spectrum at long and short wavelengths. In the squeezed limit, $f_{\mathrm{NL}}$ obeys a consistency relationship with the primordial power spectrum [56]:

For slow-roll single-field inflation, the consistency relationship becomes

For the specific observable $\mu\times T$, however, the type of PNG produced in single-field inflation results in a vanishingly small signal, far below what is predicted by the consistency relation [61, 21]. This means any measurable signal of $\mu\times T$ would rule out single-field inflation models. We review the reasoning behind this result in Section II.3.1.

In the early Universe, energy injected into the plasma will efficiently thermalize through double Compton scattering and bremsstrahlung, producing a blackbody distribution with a new temperature. These two processes are highly efficient until a redshift $z_{i}\simeq 2\times 10^{6}$, at which point the Universe has expanded enough that the density necessary for double Compton and bremsstrahlung to frequently occur becomes too low [46]. This leaves elastic Compton scattering, which conserves photon number, as the primary method for the photon-baryon bath to reach thermal equilibrium. The inefficient thermalization introduces a chemical potential, where the mixing of different blackbody spectra at different temperatures produces a Bose-Einstein rather than Planckian distribution [80]. At small chemical potentials, the Bose-Einstein distribution can be approximated as a distorted blackbody spectrum in which

where $h$ is Planck’s constant and $k_{\mathrm{B}}$ is Boltzmann’s constant. This distortion of the CMB spectrum from a typical blackbody is known as a $\mu$-type distortion. It occurs until a redshift $z_{\mathrm{f}}\simeq 5\times 10^{4}$, where even thermalization through single Compton becomes inefficient and distortions of the $y$-type start being produced [28].

$\mu$ distortions can be generated from a variety of mechanisms. The primary contribution we will consider in this work is from diffusion damping of small-scale power. The $\mu$ distortions produced by diffusion damping can be related to the primordial curvature perturbations via

where $W(\vec{k}_{1},\vec{k}_{2})$, the window function, captures the weighted amount of dissipated modes with wavenumber $k$ [29, 21, 79]. To relate the amount of $\mu$ distortions to the primordial power spectrum, we take the ensemble average of $\mu(x)$ and ${k}_{1}={k}_{2}$

Assuming that $\mu$ distortions are generated from sub-horizon modes at the time of dissipation, the window function is

where $k_{\mathrm{D}}$ is the damping scale for energy injection from diffusion damping,

and $\mathcal{J}_{\mu}$ is the time window function for $\mu$ distortions [29, 21, 79]. The time window function is well approximated by an analytical Green’s function given in Ref. [27]

where $z_{i}\simeq 2\times 10^{6}$ is defined above as the beginning of the $\mu$-distortion era. For single-field, slow-roll inflation, the average amount of $\mu$ distortions is roughly $\langle\mu(x)\rangle\simeq 2\times 10^{-8}$ [31, 20, 28]. We find that $\mu$ distortions are tracers of primordial perturbations in the range $50\textrm{ Mpc}^{-1}\lesssim k_{\mu}\lesssim 1\times 10^{4}\textrm{ Mpc}^{-1}$ [31, 20, 28].

Anisotropies of the CMB probe the primordial curvature perturbations of inflation. CMB anisotropies are typically decomposed into spherical harmonics. For example, the real-space temperature anisotropy field $\Theta(\hat{n})=\delta T(\hat{n})/T$, can be decomposed into spherical harmonics $\Theta(\hat{n})=\sum_{\ell m}a_{\ell m}^{T}Y_{\ell m}(\hat{n})$. The spherical harmonic coefficients of the various CMB anisotropy fields are related to (Fourier-space) primordial curvature perturbations through

where $X$ denotes the type of CMB anisotropy field (we are considering only $T$ and $E$ here), and $\Delta_{\ell}^{X}(k)$ is the transfer function connecting primordial perturbations to CMB anisotropies. The coefficients can then be correlated with each other:

where $C_{\ell}^{XY}$ is the angular cross-power spectrum of CMB anisotropy fields $X$ and $Y$.

Observations of $\mu$-distortion anistropies can also be decomposed into spherical harmonics:

where $\Delta^{\mu}_{\ell}(k)$ is the transfer function of anisotropic $\mu$ distortions [21]:

and $\eta_{0}$ and $\eta_{\mathrm{rec}}$ are the conformal times corresponding to $z=0$ and the redshift of recombination $z=z_{\mathrm{rec}}$, respectively. The angular correlation of $\mu$ distortions with CMB anisotropies is then given by

where we have explicitly indicated here that, because of the very different transfer functions, this correlation probes the connection of the curvature power on very small scales through $\mu(\vec{k}_{\mathrm{S}}$) with the large-scale curvature $\zeta(\vec{k}_{\mathrm{L}})$. In other words, this measurement is sensitive to the correlation of a large-scale mode with two extremely small-scale modes, i.e., the bispectrum in the ultra-squeezed limit. The non-Gaussianities being probed are at scales of $50\textrm{ Mpc}^{-1}\lesssim k_{\mu}\lesssim 1\times 10^{4}\textrm{ Mpc}^{-1}$, much smaller than scales probed by Planck measurements of the primary CMB [9].

The ensemble average of $\left\langle\mu\left(\vec{k}_{\mathrm{S}}\right)\zeta(-\vec{k}_{\mathrm{L}})\right\rangle$ is

We use this to re-express $C_{\ell}^{\mu X}$ as

The angular power spectra $C_{\ell}^{\mu T}$ and $C_{\ell}^{\mu E}$ for $\langle\mu\rangle=2\times 10^{-8}$ and $f_{\mathrm{NL}}=1$ are plotted in Figure 1.

If we would like to know how accurately we can measure a given parameter $p_{i}$ in a data set, we can assume the likelihood $\mathcal{L}$ of measuring the parameter follows a Gaussian distribution:

where $p_{i}$ is a fiducial value of the parameter and $\hat{p}_{i}$ is the measured value. $F_{ij}$ is the Fisher matrix, which captures the covariance of measured parameters. We discuss the validity of the assumption of Gaussian likelihood in Section III.2.

We make the approximation that the amplitude of the $\mu\times T$ and $\mu\times E$ spectra are controlled by a single free parameter $f_{\mathrm{NL}}$ and write

The Fisher “matrix” in this case is a scalar:

and the 1$\sigma$ uncertainty on $f_{\mathrm{NL}}$ is

The expected noise at a given multipole $\ell$ is given by

where $f_{\mathrm{sky}}$ is the fraction of the full sky observed in the survey. Our fiducial model for forecasting is $C_{\ell}^{\mu X}=0$, and we neglect the contribution of the noise part of $C_{\ell}^{\mu X}$ to the variance, as it will always be much smaller than the product of $C_{\ell}^{\mu\mu}$ and $C_{\ell}^{XX}.$ The likelihood is thus

For the ${C_{\ell}^{\mu X}}$ cross-spectrum (the numerator or signal part of the Fisher matrix calculation), we use the formulation in Eq. (21). For our fiducial forecasts, we adopt $\langle\mu\rangle=2\times 10^{-8}$. For the denominator or noise part of the Fisher calculation, we note that the auto-power-spectra ${C_{\ell}^{XX}}$ and ${C_{\ell}^{\mu\mu}}$ can be separated into signal and noise terms

where $S$ and $N$ denote signal and noise, respectively. At CMB-S4 noise levels, measurements of both ${C_{\ell}^{TT}}$ and ${C_{\ell}^{EE}}$ are signal-dominated for $\ell\leq 2000$ (at which point our constraints on $f_{\mathrm{NL}}$ are well saturated, see Figure 6). Therefore, we neglect the effects of noise and foregrounds on our temperature and $E$-mode anisotropy maps, and can make the approximation ${C_{\ell}^{XX}}\approx{C_{\ell}^{XX,S}}$. On the other hand, measurements of ${C_{\ell}^{\mu\mu}}$ will be noise-dominated for the foreseeable future, such that ${C_{\ell}^{\mu\mu}}\approx{C_{\ell}^{\mu\mu,N}}$. Calculations for ${C_{\ell}^{XX,S}}$ are taken from CAMB¹¹1http://camb.info [55], while ${C_{\ell}^{\mu\mu,N}}$ is dependent on instrument and observation parameters, including instrumental and atmospheric noise levels. We describe how we obtain ${C_{\ell}^{\mu\mu,N}}$ from noise and foreground models in the following sections.

The constraining power of $C^{\mu T}_{\ell}$ and $C^{\mu E}_{\ell}$ for measuring $f_{\mathrm{NL}}$ are comparable to each other. Rather than having independent constraints on $f_{\mathrm{NL}}$ from $C^{\mu T}_{\ell}$ or $C^{\mu E}_{\ell}$, we can use the fact that they both probe the same underlying correlator $\langle\mu\zeta\rangle$ and perform a joint analysis. The authors of Ref. [62] and Ref. [65] demonstrated that the differing behavior of $C^{\mu T}_{\ell}$ and $C^{\mu E}_{\ell}$ with $\ell$ provides a better constraint on $f_{\mathrm{NL}}$ than an independent analysis. We follow their method and modify the likelihood function to include the correlations between $T$ and $E$ in both the signal and the covariance. The final likelihood is then

For the purposes of our forecasting, observed maps of the total intensity of the sky in direction $\hat{n}$ and frequency $\nu_{i}$, $I_{i}(\hat{n})$, can be described as a linear combination of the CMB $\mu$ and $T$ signals and a noise contribution:

where $a_{\mu,i}$ and $a_{\mathrm{T},i}$ are the $\mu$ and temperature spectral energy distributions (SEDs) at different frequency bands $i$, and $s_{\mu}(\hat{n})$ and $s_{\mathrm{CMB}}(\hat{n})$ are the true underlying $\mu$ and temperature anisotropy maps. We treat all astrophysical signals that are not CMB $\mu$ or $T$ as noise and include them in $n$. We can choose to work in spherical harmonic space, defining $I_{\ell m}$ such that $x(\hat{n})=\sum_{\ell}\sum_{m}I_{\ell m}Y_{\ell m}(\hat{n})$. We can then rewrite Eq. (33) as

Traditionally, data from CMB experiments are calibrated such that maps in all frequency bands have the same response to primary CMB temperature anisotropy—i.e., the maps are in units of CMB fluctuation temperature $\Delta T$ or fractional CMB fluctuation $\Delta T/T$. In the latter case, the CMB temperature SED $a_{\mathrm{T},i}$ is given by

for all bands. We follow Ref. [63] and approximate the $\mu$-distortion SED at frequencies of 20 GHz and above as

where

We can obtain a temperature-free $\mu$ map or its spherical harmonic transform through component separation using a constrained internal linear combination method (CILC) [64]. This method takes advantage of the known SEDs of temperature and $\mu$-distortion anisotropies and calculates weights $\bm{w}$ that, when applied to observed frequency maps, result in a $T$-free $\mu$ map and a $\mu$-free $T$ map. For example, if we assign the $\mu$ weights to the $i=0$ component of $\bm{w}_{ij}$, then the $T$-free $\mu$ map is given by

It is shown in Ref. [63] that the weights that enforce unit response to $\mu$ distortions and zero response to temperature anisotropy, and minimize total variance, are given by

where

$\mathbf{C}=\mathbf{C}^{ij}_{\ell}$ is the frequency-frequency $\ell$-space covariance matrix, and $T$ denotes transpose. The above weights can be generalized to null additional components $\bm{b}_{m}$ by generalizing $\mathrm{A}\rightarrow\left(\begin{array}[]{lllll}\bm{a}_{\mu}&\bm{a}_{\mathrm{CMB}}&\bm{b}_{1}&\ldots&\bm{b}_{m}\end{array}\right)$ and $\bm{e}^{T}\rightarrow\left(\begin{array}[]{lllll}1&0&0&\ldots&0\end{array}\right)$. We note that the weights depend on multipole number $\ell$, i.e., $\bm{w}\rightarrow\bm{w}_{\ell}$. For simplicity of notation, we leave the $\ell$ dependence implicit and continue to use $\bm{w}$.

The reduction of the full $\mathbf{C}^{ij}_{\ell m\ell^{\prime}m^{\prime}}=\left\langle I^{i}_{\ell m}I^{j}_{\ell^{\prime}m^{\prime}}\right\rangle$ to $\mathbf{C}^{ij}_{\ell}$ rests on the assumption that all sources of variance in the maps are isotropic, stationary, and Gaussian. With detector noise only and assuming no correlations between detector noise at different frequencies, the covariance matrix is diagonal ($\mathbf{C}^{\mathrm{ij}}_{\ell}=\mathbf{C}^{\mathrm{ii}}_{\ell}\delta_{ij}$). Including foregrounds and atmosphere introduces correlated fluctuations between frequency bands and requires the full frequency-frequency matrix.

The assumption of Gaussianity and statistical isotropy is reasonably satisfied by detector noise, atmospheric emission, and extragalactic foregrounds, but not especially well by Galactic foregrounds. The CMB-S4 “ultra-deep” survey, which is the main survey we present forecasts for in this work (see Section IV) is located in an area of the sky with very low Galactic foreground emission. Furthermore, non-Gaussianity in the true foreground emission not accounted for in the covariance will not bias the component separation, it will only make it slightly suboptimal. We note that the Gaussianity of the likelihood (Eq. III.1) depends on the power spectrum of the noise sources being Gaussian, not the pixel values or modes of the noise sources themselves, which is a lighter burden owing to the Central Limit Theorem.

We can then obtain ${C_{\ell}^{\mu\mu}}$ by applying the $\mu$-distortion weights to the frequency-frequency covariance matrix:

We discuss the various contributions to the frequency-frequency covariance matrix $\mathbf{C}^{\mathrm{ij}}_{\ell}$ in the following sections.

Assuming detector noise that is white (uncorrelated between time samples), uniform sky coverage, and a Gaussian instrument beam or point-spread function, the statistics of the map noise for frequency band $i$ in spherical harmonic or $\ell$ space are given by

where $N_{i}$ and $\theta_{i}$ is are the white noise level in the map and the angular resolution (beam FWHM) for frequency band $i$. The contribution to the frequency-frequency covariance matrix from this source will be

and, for detector noise that is uncorrelated between bands,

where $\delta_{ij}$ is the Kronecker delta. For white detector noise only, the $\mu$-distortion power spectrum reduces to

A major source of noise that must be considered in ground-based observations of the CMB is the emission from blobs of poorly mixed water vapor in the Earth’s atmosphere (i.e., clouds). For detailed discussions of this effect and measurements of the impact at various sites, see, e.g., Refs. [54], [19], and [58]. The spectrum of cloud sizes is such that the noise power from this source is much larger at large angular scales, and it is often modeled as a power law in $\ell$ (e.g., Ref. [2]). The total detector + atmosphere noise power in frequency band $i$ can then be parameterized with three numbers, namely the white noise level $N_{i}$, the multipole value at which the detector and atmosphere noise levels are equal $\ell_{\mathrm{knee},i}$, and the power-law index of the atmosphere noise $\alpha^{\mathrm{atmo},i}$.

Implicit in the formulation of Eq. (49) is the assumption that the atmospheric noise is uncorrelated between bands. In fact, nearly the opposite is the case, at least for detectors for which the beam patterns mostly overlap at the height of the atmospheric emission. For example, internal CMB-S4 analysis of data from the SPT-3G receiver on the South Pole Telescope (SPT, [22, 71]) found that for detectors in different frequency bands but co-located in a focal-plane pixel, the long-timescale fluctuations in the time-ordered data were over 99% correlated. To model this, we can introduce an atmospheric correlation parameter $\eta$ and rewrite the noise + atmosphere contribution to the frequency-frequency covariance matrix as

where $\pi_{ij}=\eta(1-\delta_{ij})$. In Sec. V, we will show forecasts using values of using values of $\eta$ ranging from 0 to 1. As the correlation $\eta$ increases, the component separation algorithm defined in Sec. III.2 is more effective in reducing the atmospheric contribution to the final ${C_{\ell}^{\mu\mu}}$ covariance.

In our forecasting pipeline, we also consider the effects of foregrounds. Foregrounds contribute to the frequency-frequency covariance matrix as linear additions (because they are not correlated with the other sources of covariance) such that

For this work, we approximate all foregrounds as 100% correlated across frequency bands (though we approximate each foreground source as uncorrelated with the others), such that for a given foreground type (call it “type X”)

Foreground sources can be separated into Galactic and extragalactic sources, and we treat each of these in turn below. Often in the literature, foreground amplitudes and $\ell$-space behavior are quoted in $D_{\ell}=\frac{\ell(\ell+1)}{2\pi}C_{\ell}$. When we adopt such parameterizations, we keep the description in $D_{\ell}$ but convert to $C_{\ell}$ when actually implementing the model. We show the SEDs of the primary foreground sources, along with the $T$ and $\mu$ SEDs, in Figure 2.