CMB-S4 forecast on the primordial non-Gaussianity parameter of feature models

TL;DR

This paper forecasts the constraints on primordial non-Gaussianity parameter $f_ ext{NL}$ for feature models using future CMB-S4 experiments, introducing a novel covariance orthonormalisation method to improve computation.

Contribution

It presents new forecast methods for $f_ ext{NL}$ constraints in CMB-S4, including a novel covariance matrix orthonormalisation technique, and compares sensitivities to Planck and Simons Observatory.

Findings

CMB-S4 could improve $f_ ext{NL}$ constraints by up to 2.2 times over Planck.

Baseline CMB-S4 conditions yield 1.3-1.6 times better constraints than Planck.

Temperature and E-mode polarization data are crucial for detecting oscillatory features.

Abstract

We present forecasts on the primordial non-Gaussianity parameter $f_\text{NL}$ of feature models for the future Cosmic Microwave Background Stage-4 (CMB-S4) experiments. The Fisher matrix of the bispectrum estimator was computed using noise covariances expected for preliminary CMB-S4 specifications including ones for the Simons Observatory. We introduce a novel method that improves the computation by orthonormalising the covariance matrix. The most sensitive CMB-S4 experiment with 1' beam and 1$\mu K$-arcmin noise would yield a factor of 1.7-2.2 times more stringent constraints compared to Planck. Under the Simons Observatory baseline conditions the improvement would be about 1.3-1.6 times to Planck. We also thoroughly studied the effects of various model and experimental parameters on the forecast. Detailed analysis on the constraints coming from temperature and E-mode polarisation, in particular, provided some insight into detecting oscillatory features in the CMB bispectrum.

Tables

Table 1. Table 1: Forecasts on the estimation errors of f N ​ L subscript 𝑓 𝑁 𝐿 f_{NL} for the constant model

	Planck	SO baseline + Planck	SO goal + Planck	1’ beam, 1$\mu K$ noise + Planck
$\sigma \left(f_{\text{NL}}^{\text{const}}\right)$	23.4	14.9	14.0	10.4

Table 2. Table 2: Expected improvements on estimation errors of f NL subscript 𝑓 NL f_{\text{NL}} for the CMB-S4 1’beam, 1 μ K ⋅ \mu K\cdot arcmin setup, for various bispectrum templates. The local, equilateral and orthogonal results are quoted from Abazajian ( 2016 ) .

	Local	Equilateral	Orthogonal	Constant	Feature ($\omega =200$)
${\sigma }^{\text{Planck}}/{\sigma }^{\text{CMB-S4}}$	2.5	2.1	2.4	2.3	2.0

Table 3. Table 3: Expected improvements on the estimation errors of f NL subscript 𝑓 NL f_{\text{NL}} for each combination of Planck / CMB-S4 temperature (T) and polarisation (E) data. Here the CMB-S4 assumes 1’ beam and 1 μ ​ K 𝜇 𝐾 \mu K arcmin noise. For feature model the oscillation frequency ω = 200 𝜔 200 \omega=200 and phase ϕ = 0 italic-ϕ 0 \phi=0 . The sky fraction f s ​ k ​ y = 0.4 subscript 𝑓 𝑠 𝑘 𝑦 0.4 f_{sky}=0.4 for all cases except for Planck T + Planck E.

$\sigma \left(f_{\text{NL}}^{\text{const}}\right)$ improvement		E
		Planck	CMB-S4
T	Planck	1.0	1.6
	CMB-S4	1.1	2.2

Equations

$$⟨\Phi ({\mathbf{k}}_1)\Phi ({\mathbf{k}}_2)\Phi ({\mathbf{k}}_3)⟩=(2\pi {)}^3{\delta }^{(3)}({\mathbf{k}}_1+{\mathbf{k}}_2+{\mathbf{k}}_3)B_{\Phi }(k_1,k_2,k_3),$$

$$B_{\Phi }(k_1,k_2,k_3)=f_{NL}{B}_{\Phi }^{(f_{NL}=1)}(k_1,k_2,k_3).$$

$$a_{lm}^X=4\pi (-i{)}^l\int \frac{d^3\mathbf{k}}{(2\pi {)}^3}\Phi (\mathbf{k}){\Delta }_l^X(k)Y_{lm}(\hat{\mathbf{k}}).$$

$$b_{l_1{l}_2{l}_3}^{X_1{X}_2{X}_3}={\left(\frac{2}{\pi }\right)}^3{\int }_0^{\infty }{r}^{2}dr{\int }_{{\mathcal{V}}_k}{d}^3\text{k}(k_1{k}_2{k}_3{)}^2{B}_{\Phi }(k_1,k_2,k_3)\prod _{i=1}^3\left[j_{l_i}(k_{i}r){\Delta }_{l_i}^{X_i}(k_i)\right],$$

$$B_{\Phi }^{\text{feat}}(k_1,k_2,k_3)=\frac{6A^2}{(k_1{k}_2{k}_3{)}^2}\sin (\omega K+\varphi ),$$

$$B_{\Phi }^{\text{feat}}(k_1,k_2,k_3)=\cos \varphi B_{\Phi }^{\text{sin}}(k_1,k_2,k_3)+\sin \varphi B_{\Phi }^{\text{cos}}(k_1,k_2,k_3).$$

$$b_{l_1{l}_2{l}_3}^{X_1{X}_2{X}_3,\text{feat}}$$

$$s_l^X(r)$$

$$c_l^X(r)$$

$$b_{l_1{l}_2{l}_3}^{X_1{X}_2{X}_3,\text{feat}}$$

$$S_i=\frac{1}{6}\sum _{l_j,m_j}\sum _{X_j}{\mathcal{G}}_{m_1{m}_2{m}_3}^{l_1{l}_2{l}_3}{b}_{l_1{l}_2{l}_3}^{X_1{X}_2{X}_3,(i)}(C_{l_1{m}_1,l_4{m}_4}^{-1}{)}^{X_1{X}_4}(C_{l_2{m}_2,l_5{m}_5}^{-1}{)}^{X_2{X}_5}(C_{l_3{m}_3,l_6{m}_6}^{-1}{)}^{X_3{X}_6}$$

$$\left[a_{l_4{m}_4}^{X_4}{a}_{l_5{m}_5}^{X_5}{a}_{l_6{m}_6}^{X_6}-\left(C_{l_4{m}_4,l_5{m}_5}{a}_{l_6{m}_6}^{X_6}+\text{2 cyclic}\right)\right].$$

$$\widehat{f_i}=\sum _j(F^{-1}{)}_{ij}{S}_j,$$

$$S_i=\frac{1}{6}\sum _{l_j,m_j}\sum _{X_j,X_j^{\prime }}{\mathcal{G}}_{m_1{m}_2{m}_3}^{l_1{l}_2{l}_3}{b}_{l_1{l}_2{l}_3}^{X_1{X}_2{X}_3,(i)}(C_{l_1}^{-1}{)}^{X_1{X}_1^{\prime }}(C_{l_2}^{-1}{)}^{X_2{X}_2^{\prime }}(C_{l_3}^{-1}{)}^{X_3{X}_3^{\prime }}$$

$$\left[a_{l_1{m}_1}^{X_1^{\prime }}{a}_{l_2{m}_2}^{X_2^{\prime }}{a}_{l_3{m}_3}^{X_3^{\prime }}-\left(⟨a_{l_1{m}_1}^{X_1^{\prime }}{a}_{l_2{m}_2}^{X_2^{\prime }}⟩a_{l_3{m}_3}^{X_3^{\prime }}+\text{2 cyclic}\right)\right],$$

$$F_{ij}=\frac{f_{\text{sky}}}{6}\sum _{\text{all }X,X^{\prime }}\sum _{\text{all }l}{h}_{l_1{l}_2{l}_3}^2{b}_{l_1{l}_2{l}_3}^{X_1{X}_2{X}_3,(i)}(C_{l_1}^{-1}{)}^{X_1{X}_1^{\prime }}(C_{l_2}^{-1}{)}^{X_2{X}_2^{\prime }}(C_{l_3}^{-1}{)}^{X_3{X}_3^{\prime }}{b}_{l_1{l}_2{l}_3}^{X_1^{\prime }{X}_2^{\prime }{X}_3^{\prime },(j)}.$$

$$L_l=\left(\begin{array}{cc}\frac{1}{\sqrt{C_l^{TT}}}& 0\\ \frac{-C_l^{TE}}{\sqrt{C_l^{TT}}\sqrt{C_l^{TT}{C}_l^{EE}-{\left(C_l^{TE}\right)}^2}}& \frac{C_l^{TT}}{\sqrt{C_l^{TT}}\sqrt{C_l^{TT}{C}_l^{EE}-{\left(C_l^{TE}\right)}^2}}\end{array}\right)$$

$${\tilde{\Delta }}_l^X(k)=\sum _{X^{\prime }}{L}_l^{X{X}^{\prime }}{\Delta }_l^{X^{\prime }}(k),\text{and}{\tilde{a}}_{lm}^X=\sum _{X^{\prime }}{L}_l^{X{X}^{\prime }}{a}_{lm}^{X^{\prime }}.$$

$$S_i=\frac{1}{6}\sum _{l_j,m_j}\sum _{X_j}{\mathcal{G}}_{m_1{m}_2{m}_3}^{l_1{l}_2{l}_3}{\tilde{b}}_{l_1{l}_2{l}_3}^{X_1{X}_2{X}_3,(i)}\left[{\tilde{a}}_{l_1{m}_1}^{X_1}{\tilde{a}}_{l_2{m}_2}^{X_2}{\tilde{a}}_{l_3{m}_3}^{X_3}-\left(⟨{\tilde{a}}_{l_1{m}_1}^{X_1}{\tilde{a}}_{l_2{m}_2}^{X_2}⟩{\tilde{a}}_{l_3{m}_3}^{X_3}+\text{2 cyclic}\right)\right],$$

$$F_{ij}=\frac{f_{\text{sky}}}{6}\sum _{\text{all }X}\sum _{\text{all }l}{h}_{l_1{l}_2{l}_3}^2{\tilde{b}}_{l_1{l}_2{l}_3}^{X_1{X}_2{X}_3,(i)}{\tilde{b}}_{l_1{l}_2{l}_3}^{X_1{X}_2{X}_3,(j)}.$$

$$B_{\Phi }(k_1,k_2,k_3)=f_{NL}^{\sin }{B}^{\sin }(k_1,k_2,k_3)+f_{NL}^{\cos }{B}^{\cos }(k_1,k_2,k_3),$$

$$S^{\text{cub}}={\int }_0^{\infty }{r}^{2}dr\int d^2\hat{\mathbf{n}}\left[-{M_s}^3+3M_s{M_c}^2\right]\text{and}$$

$$S^{\text{lin}}=-3{\int }_0^{\infty }{r}^{2}dr\int d^2\hat{\mathbf{n}}\left[-M_s⟨{M_s}^2⟩+M_s⟨{M_c}^2⟩+2M_c⟨M_s{M}_c⟩\right],$$

$$M_s(r,\hat{\mathbf{n}})$$

$$M_c(r,\hat{\mathbf{n}})$$

$$h_{l_1{l}_2{l}_3}^2=\frac{(2l_1+1)(2l_2+1)(2l_3+1)}{8\pi }{\int }_{-1}^{1}d\mu P_{l_1}(\mu )P_{l_2}(\mu )P_{l_3}(\mu ),$$

$$F=\frac{3}{4\pi }\int r^{2}dr\int r^{\prime 2}d{r}^{\prime }\int d\mu \left[P_{ss}^3+3P_{ss}{P}_{cc}^2-3P_{cs}^2{P}_{ss}-3P_{sc}^2{P}_{ss}+6P_{cs}{P}_{sc}{P}_{cc}\right].$$

$$P_{ss}(r,r^{\prime },\mu )$$

$$P_{sc}(r,r^{\prime },\mu )$$

$$C_l^{X_1{X}_2}\to W_l^{X_1}{W}_l^{X_2}{C}_l^{X_1{X}_2}+N_l^{X_1{X}_2},b_{l_1{l}_2{l}_3}^{X_1{X}_2{X}_3}\to W_{l_1}^{X_1}{W}_{l_2}^{X_2}{W}_{l_3}^{X_3}{b}_{l_1{l}_2{l}_3}^{X_1{X}_2{X}_3},$$

$$C_l^{X_1{X}_2}$$

$$C_l^{\text{noise},TT}=\exp \left(l(l+1){\sigma }_{\text{beam}}^2\right)N_{\text{white}},C_l^{\text{noise},EE}=2C_l^{\text{noise},TT},C_l^{\text{noise},TE}=0.$$

$$F(\omega ,\varphi )={\cos }^2\varphi F_{ss}(\omega )+{\sin }^2\varphi F_{cc}(\omega )+2\cos \varphi \sin \varphi F_{sc}(\omega ),$$

Full text

I Introduction

The Cosmic Microwave Background (CMB) radiation is one of our most valuable probes of the primordial universe. The temperature and polarisation of this ancient light contains rich statistical information both about the primordial perturbations created during inflation and also their subsequent evolution until now. This allows us to test our inflationary theories and also the history of our universe. The recent Planck CMB experiments have provided stringent tests on various models of inflation through the estimation of cosmological parameters and via primordial non-Gaussianity Planck Collaboration (2015, 2018).

The simplest model of inflation involves a single scalar field slowly rolling down a smooth potential. In this case the CMB temperature fluctuations are expected to be Gaussian distributed with only tiny deviations (e.g. Maldacena (2002)). However, many other physically well-motivated models generate larger non-Gaussian signatures at the end of inflation (see reviews of Chen (2010)). Such primordial non-Gaussianities are well constrained by three-point correlation functions of the CMB anisotropies or their Fourier transform, the CMB bispectrum. Different inflationary models predict bispectra with different momentum dependence, or ‘shapes’. We constrain these models by using an optimal estimator for their amplitude parameter, $f_{\text{NL}}$, for each specific bispectrum shape (see, e.g., Komatsu (2010); Liguori et al. (2010) for reviews).

Although all observations to date are consistent with vanishing non-Gaussianity, the models most favoured by the 2015 Planck CMB analysis were the ones with oscillations in the primordial power spectrum Planck Collaboration (2015). Among them are feature models, where the oscillations are caused by a sharp feature in either the inflationary potential Starobinsky (1992); Adams et al. (2001); Chen et al. (2007); Adshead et al. (2012); Hazra et al. (2014); Dvorkin and Hu (2010), sound speed Miranda et al. (2012); Bartolo et al. (2013), or multi-field potentials Achúcarro et al. (2011) (see Chen (2010); Chluba et al. (2015) for reviews). The primordial power spectrum then becomes scale dependent, displaying sinusoidal oscillations that are linearly spaced in momentum space. The resulting bispectrum also oscillates and is highly uncorrelated with other popular bispectrum templates Meerburg et al. (2009), therefore allowing us to constrain them independently.

Planck constrained $f_{\text{NL}}$ for feature models from CMB bispectra, but no signal above $3\sigma$ significance were found after accounting for the ‘look elsewhere effect’ as introduced in Fergusson et al. (2015a). The multi-peak statistic analysis, however, revealed some non-standard signals up to $4\sigma$ level that deserves attention Planck Collaboration (2015). There have been many other searches on signatures of oscillations. Constraints also come from the CMB power spectrum Martin and Ringeval (2004); Benetti et al. (2011); Meerburg et al. (2012, 2014a, 2014b); Fergusson et al. (2015b), the large scale structure Chantavat et al. (2011); Ballardini et al. (2016), and a combination of the two Hu and Torrado (2015); Benetti and Alcaniz (2016). We expect stronger constraints on feature models from future LSS experiments Chen et al. (2016). This paper covers the prospects of upcoming CMB experiments in constraining $f_{\text{NL}}$ for feature models.

Currently there are two implementations of the optimal estimator for constraining $f_{\text{NL}}$ for feature type models. The Planck analysis adopted the Modal estimator for which the given bispectrum is expanded using a separable basis Fergusson et al. (2012); Fergusson (2014). This method is efficient, can flexibly account for various oscillatory shapes, and is able to easily constrain all frequencies simultaneously. However, when the oscillation frequency is large the modal basis fails to converge within reasonable number of basis elements, making the method impractical. The other approach using the Komatsu-Spergel-Wandelt (KSW) estimator is viable for various shapes including the feature model Komatsu et al. (2005); Münchmeyer et al. (2014). Although this method only applies to models with separable bispectra, even highly oscillatory templates can be computed reliably. This method is however more computationally expensive as each frequency must be dealt with separately. We present further optimisations to the fast KSW estimator introduced in Yadav et al. (2007) and apply it on feature models for forecasts in this paper.

The next generation of CMB experiments, CMB Stage-4, consists of many exciting proposed experiments located at the South Pole, the Atacama Desert in Chile, and perhaps space Abazajian (2016); The Simons Observatory Collaboration (2018); The COrE collaboration (2015). One of the main goals of these experiments is to measure the polarisation signal in the CMB to the cosmic variance limit. Preliminary specifications have been released for these experiments Abazajian (2016); The Simons Observatory Collaboration (2018) and these have been used to produce some forecasts for the standard $f_{\text{NL}}$ templates but not yet for feature type models. In this paper we address this by presenting the Fisher forecasts on $f_{NL}$ for feature models based on these specifications and observe that feature type models receive larger improvements from the extra polarisation information than the standard templates, justifying this analysis.

The paper is organised as follows. First we briefly review the theory of CMB bispectrum in Section II. Bispectrum template for the feature model is defined and computed here. In Section III we formulate the bispectrum estimator and introduce a new method to further optimise its computation. The technique is applied to the case of feature model to yield equations for the Fisher forecast of $f_{NL}$. We also briefly discuss implementation details. In Section IV we present our forecast results and their dependence on model and experimental parameters. In particular, forecasts for the Simons observatory are compared with the Planck results. The results are summarised in Section V.

II Feature model bispectrum

CMB bispectrum

One of the main subjects of primordial non-Gaussianity studies is the 3-point correlation function of the primordial perturbations which is defined by;

[TABLE]

where we have assumed statistical homogeneity and isotropy. The primordial bispectrum $B_{\Phi}$ vanishes for Gaussian perturbations, but more general inflation models predict non-zero bispectra with various shapes. In order to constrain these models we re-parameterise the bispectrum into a amplitude parameter and a normalised shape part;

[TABLE]

Constraining $f_{NL}$ from the CMB measurements allows to determine how well the particular shape under consideration aligns with the data, which we can then translate into constraints on the model itself.

In order to compare the theory with measurements we first need to relate the primordial perturbations to spherical multipole modes of the late-time CMB anisotropies.

[TABLE]

Here the index $X$ is either $T$ or $E$, representing CMB temperature and E-mode polarisation, respectively. The linear CMB radiation transfer function $\Delta_{l}^{X}(k)$ can be computed from the Boltzmann solvers like CAMB Lewis et al. (2000).

Three point correlation function of $a_{lm}^{X}$’s yield the reduced bispectrum $b_{l_{1}l_{2}l_{3}}$ times a geometrical factor $\mathcal{G}^{l_{1}l_{2}l_{3}}_{m_{1}m_{2}m_{3}}$ named the Gaunt integral. After some algebraic manipulations we obtain the following useful formula for the reduced bispectrum;

[TABLE]

where $j_{l}$ is the spherical Bessel function arising from the Rayleigh expansion formula. Using this equation, we can compute the projected bispectrum from any given primordial bispectrum. Direct computation of this four-dimensional integral for every $l$ combination, however, is practically impossible. Not only is the integral in 4D but also the oscillatory integrand requires a large number of sample points in each of $k_{i}$, making the full calculation for every $l_{i}$ triple prohibitively expensive. All bispectrum estimators get around this problem by expanding $B_{\Phi}$ as a sum of separable terms. This will be explained in more detail later using the feature model template as an example.

Feature model

We follow the works of Münchmeyer et al. (2014); Fergusson et al. (2015a, b); Planck Collaboration (2015) and assume the following template for the bispectrum of feature models;

[TABLE]

where $K=k_{1}+k_{2}+k_{3}$, $A$ represents the primordial power spectrum amplitude, and $\phi$ is a phase. The oscillation ‘frequency’ $\omega$ is associated with the location and scale of feature in the inflationary potential. It is often written in terms of the oscillation scale $k_{c}$ as $\omega=2\pi/3k_{c}$. $\omega$ is measured in Mpc but we omit the unit for notational conveniences.

The feature model template has two free parameters that need to be fixed before we can constrain the model: $\omega$ and $\phi$. The phase $\phi$ can be easily dealt with by observing that

[TABLE]

Here $B_{\Phi}^{\text{sin}}$ and $B_{\Phi}^{\text{cos}}$ correspond to feature models with $\phi=0$ and $\pi/2$, respectively. Non-zero phase simply corresponds to a linear combination of the sine and cosine templates. As we will see later these two shapes are in fact highly uncorrelated. Therefore, they can be constrained independently from each other.

On the other hand, one still has a complete freedom of choice on the oscillation frequency $\omega$. Such freedom dramatically expands size of the parameter space. In practice we constrain $f_{\text{NL}}$ for each fixed value of oscillation frequency, which yields hundreds of estimates. Since there are so many estimates we are looking at, there is a good chance that we find notable signals by sheer luck. Accounting for this ‘look elsewhere effect’ has been resolved using methods in Fergusson et al. (2015a) and subsequently applied to the Planck analysis Planck Collaboration (2015); Fergusson et al. (2015b). The look-elsewhere-adjusted statistics used in the literature can be employed for the future CMB-S4 data analysis. This work, however, focusses on forecasting the ‘raw’ estimates and comparing them with those of Planck.

Separability

The bispectrum template of feature models (5) is an example of separable shape. It can be expressed as a sum of terms in the form $f(k_{1})g(k_{2})h(k_{3})$ for some functions $f$, $g$ and $h$, which dramatically simplifies the computation of reduced bispectrum $b_{l_{1}l_{2}l_{3}}$. The three-dimensional integral over the $k$ space in (4) splits into three individual one-dimensional integrals for separable shapes. Feature models for example has

[TABLE]

Here the real and imaginary parts of $b^{\text{feat}}$ correspond to the bispectra of cosine and sine feature models, respectively. Now define

[TABLE]

These are analogous to $\alpha_{l}^{X}(r)$ and $\beta_{l}^{X}(r)$ in the usual KSW estimator for local non-Gaussianity. Then (7) reduces to

[TABLE]

III Efficient computation of the estimator with polarisation

Estimator

The optimal estimator for a given bispectrum in the weak non-Gaussian limit involves computing Komatsu et al. (2005); Komatsu (2010)

[TABLE]

Here summations are over $l_{j}$, $m_{j}$, $X_{j}$ and $X^{\prime}_{j}$ for each $j=1,2,3$. The spherical multipole moments $a_{lm}^{X}$’s are computed from observations, and $b^{(i)}$ denotes the $i$th theoretical bispectrum template under consideration.

Computing this form, however, requires an inversion of the full covariance matrix $C_{lm,l^{\prime}m^{\prime}}$, which is computationally expensive. As a result we will follow the diagonal covariance approximation in Yadav et al. (2007) for the inverse covariance; $C_{l_{1}l_{4}m_{1}m_{4}}^{-1}\approx(1/C_{l_{1}})~{}\delta^{D}_{l_{1}l_{4}}\delta^{D}_{m_{1}-m_{4}}$. We also approximate the covariance in the linear term by an ensemble average over Monte Carlo simulations of Gaussian realisations; $C_{l_{4}l_{5}m_{4}m_{5}}^{X_{1}X_{2}}\approx\left<a_{l_{1}m_{1}}^{X_{1}}a_{l_{2}m_{2}}^{X_{2}}\right>$. With these simplifications the estimator takes the form

[TABLE]

where

[TABLE]

and

[TABLE]

The covariance matrix $C_{l}$ is now a $2\times 2$ matrix consisting of values $C_{l}^{TT}$, $C_{l}^{TE}$, $C_{l}^{ET}$ and $C_{l}^{EE}$. 111Note that this is equivalent to having a $2l\times 2l$ matrix with diagonal $l\times l$ block matrices $C^{TT}$, $C^{TE}$, $C^{ET}$ and $C^{EE}$ as in other literatures including Fergusson (2014). The linear terms (the second in square brackets) are required to account for anisotropies induced by masking and anisotropic noise.

$F_{ij}$ is the Fisher information matrix of the estimator. $f_{\text{sky}}$ in (14) denotes the fraction of the sky covered by the experiment, and $h_{l_{1}l_{2}l_{3}}^{2}:=\sum_{m_{j}}\left(\mathcal{G}_{m_{1}m_{2}m_{3}}^{l_{1}l_{2}l_{3}}\right){}^{2}$ is a geometric factor. Since the estimator $\hat{f}_{i}$ in (12) is nearly optimal, its 68% confidence (1$\sigma$) interval can be computed from the Fisher matrix as $\sigma_{i}:=\Delta f_{NL}^{(i)}=(F^{-1})_{ii}$.

Note that most CMB-S4 experiments are ground-based, so they can probe smaller fraction of the sky compared to Planck. Having a smaller fraction of the sky leads to increased uncertainties for the estimator. Current estimate is that the new experiments will cover 40% of the sky, significantly less than the 74% of Planck. The error bars will thus increase by a factor of 1.38 from the decrease in $f_{\text{sky}}$ alone. This may be reduced by combining Planck data for unobserved pixels in these experiments

Orthonormalising the covariance matrix

In Fergusson (2014) it was noted that orthogonalising the multipoles of temperature and polarisation maps dramatically reduces the number of terms in computation of the Modal estimators. This technique can also be applied to KSW estimators, or indeed any optimal bispectrum estimator, which is yet to be done to the authors’ knowledge.

In both (13) and (14) there are summations over indices $X$ and $X^{\prime}$ to account for correlations between the CMB temperature and E-mode polarisation. This can be simplified by essentially making a change of basis in $X$ space for each $l$ so that every $C_{l}$ becomes orthonormal. We perform a Cholesky decomposition on $C_{l}$ and invert the matrix. Then $C_{l}^{-1}=L_{l}^{T}L_{l}$, where $L_{l}$ is a lower triangular matrix given by

[TABLE]

Now let

[TABLE]

Defining $\tilde{b}_{l_{1}l_{2}l_{3}}$ to be the corresponding reduced bispectrum, (13) and (14) simplify to

[TABLE]

[TABLE]

Using this method not only makes it more mathematically concise, but also halves the number of terms involved in the summation. Linear transformations (16) only need to be done once in the beginning of the program and cost little compared to the main computation. We also found it easier to optimise the code using instruction level vectorisations after this simplification.

The only downside of this method is that we no longer can get breakdowns of signal from each of $TTT$, $TTE$, $TEE$ and $EEE$ bispectrum since our new modes are linear combinations of $T$ and $E$ modes. However, in most cases we are interested in either $T$-only or $T+E$ results, and this method works perfectly well in these cases.

Estimator for feature models

We compute the general estimator (17) and (18) for feature models. The method is similar to the one seen in Münchmeyer et al. (2014) except that now the polarisation is included and the covariance matrices are trivial thanks to the orthonormalisation process outlined above.

Consider the bispectrum shape of

[TABLE]

for a fixed value of oscillation frequency $\omega$. Here $B^{\sin}$ and $B^{\cos}$ correspond to reduced bispectra $b^{\sin}$ and $b^{\cos}$ defined in (10). The Fisher matrix $F$ is $2\times 2$ but its off-diagonal entries are 2-3 orders of magnitude smaller than diagonal ones in most cases as will be presented in the next section. Thus, the two shapes are assumed to be uncorrelated and constrained individually. Here we present detailed computations for $f_{NL}^{\sin}$ only but the cosine one can be computed similarly.

From (10) and the definition of Gaunt integral $\mathcal{G}^{l_{1}l_{2}l_{3}}_{m_{1}m_{2}m_{3}}=\int d\hat{\mathbf{n}}\;Y_{l_{1}m_{1}}(\hat{\mathbf{n}})Y_{l_{2}m_{2}}(\hat{\mathbf{n}})Y_{l_{3}m_{3}}(\hat{\mathbf{n}})$ it follows that

[TABLE]

[TABLE]

where

[TABLE]

Again, the bracket $\left<\cdot\right>$ denotes averaging over Gaussian simulations. The sum of $S^{\text{cub}}$ and $S^{\text{lin}}$ gives the final value of $S$ for sine feature model.

For efficient Fisher matrix calculation we follow Smith and Zaldarriaga (2011) and deploy the identity

[TABLE]

where $P_{l}(\mu)$ represents the Legendre polynomial. Then,

[TABLE]

where we have defined

[TABLE]

and similarly $P_{cs}$ and $P_{cc}$.

Calculations of (22) and (25) are two of the most computationally expensive steps. If we have not orthonormalised the covariance matrix, then there would be an extra summation over $X^{\prime}$ and some $2\times 2$ matrix algebra involving $(C_{l}^{-1})^{XX^{\prime}}$.

Probing beam and instrumental noise

In an ideal experiment where measurements are made on each point of the sky perfectly, the covariance matrice $C_{l}^{XX^{\prime}}$ in (13) and (14) consists purely of the signal. In reality, however, the probing beam has finite width and the sensors are noisy. These effects can be incorporated by modifying the covariance matrices and bispectra as follows.

[TABLE]

where $W_{l}^{X}$ and $N_{l}^{X_{1}X_{2}}$ represent the beam window function and the noise covariance matrix, respectively. When substituted into the KSW estimator, these changes are equivalent to modifying

[TABLE]

while keeping the bispectra same. Here we have defined the effective (beam-corrected) noise covariance matrix $C_{l}^{\text{noise}}$. Modes for which $C_{l}^{\text{noise}}$ is much larger than $C_{l}^{\text{sig}}$ contribute little to the $f_{NL}$ estimator.

For forecasting purposes we assume Gaussian beam and white uncorrelated noise until more detailed experiment specifications become available. Under these assumptions, the effective noise covariances reduce to Ng and Liu (1999)

[TABLE]

The factor of two for $EE$ mode is comes from measuring two Stokes parameters Q and U. The Gaussian beam profile is usually specified by its FWHM (full width at half maximum) in $arcmin$, which is then converted to standard deviations in radians for $\sigma_{\text{beam}}$. The noise level often comes in the units of $\mu K\cdot arcmin$. This is then divided by $T_{\text{CMB}}=2.725K$, converted to radians and squared to get $N_{\text{white}}$.

For the Planck experiment, using 5 arcmin FWHM beam and the  47 $\mu K\cdot$arcmin noise level gives good approximations to the post-component-separation noise covariances. For CMB-S4 experiments the details are not confirmed, but the beam FWHM is expected to lie between 1-5 arcmin, while the noise level will range from 1 to 9 $\mu K\cdot$arcmin. Abazajian (2016)

In real measurements there exist extra contaminations in large angular scales due to 1/f noises and the component separation process. Though most of our analysis assumes simpler form of noise covariances elaborated above, for the Simons Observatory forecasts we follow The Simons Observatory Collaboration (2018) and model 1/f noise as $N_{l}=N_{red}(l/l_{knee})^{\alpha_{knee}}+N_{white}$. The noise curves from each channel were then put together using the inverse variance method. This is a good approximation for the E mode polarisation but not for temperature, since extra degradations occur during the component separation process. Still, because dominant contributions to the feature model signal comes from polarisation data, this would be a reasonable approximation for our forecast. For Planck the full post-component-separation noise curves are available and hence used for computations.

Implementation and validation

We implemented the pipeline outlined above using the C programming language and parallelised using hybrid MPI + openMP. The code was then run in the COSMOS supercomputing system.

The transfer functions are generated from the CAMB code Lewis et al. (2000). Bessel function values were pre-computed using recursion relations and stored in a file, while the Legendre function values were computed on the fly using the GNU scientific library. The angular power spectrum data was generated from $\Lambda$CDM parameters estimated in the Planck 2015 results.

Numerical integration for variables $k$, $r$ and $r^{\prime}$ were done using simple trapezoidal methods, as they can be easily vectorised for optimisation. On the other hand, integration of $\mu$ required more care because the Legendre polynomials are highly oscillatory. We adopted the Gauss-Legendre quadrature rule with $1.5\,l_{\text{max}}+1$ points which can integrate polynomials up to order $3\,l_{\text{max}}$ exactly. The weights and nodes were computed in the beginning using the QUADPTS code Hale and Townsend (2013).

Various checks have been done to ensure that the code runs correctly. First we used the code to reproduce the Planck results, which agreed within 3% error. The code was then used to compute bispectrum for the constant model, corresponding to the case where $\omega=\phi=0$. There exists an approximate analytic form in this case Fergusson et al. (2012) which we were able to reproduce accurately. We also performed convergence tests on $r$ and $r^{\prime}$ integration by doubling the number of points for each of them. The grid was chosen to be very dense around recombination and quite dense near reionisation. We confirmed that changes in the integral are less than 0.5% for each value of $\omega$.

IV CMB-S4 forecast results

Phase dependence

We now present the CMB-S4 forecast on the error bars of primordial non-Gaussianity parameter for feature models. For notational convenience we denote the error bars for sine and cosine feature models by $\sigma_{\sin}$ and $\sigma_{\cos}$. Superscripts $T$ and $T+E$ are also put to distinguish temperature-only analysis from the full analysis including polarisation.

First of all, we check that the sine and cosine bispectrum templates defined in (5) are indeed uncorrelated and can be constrained separately. In order to do this, we see if the Fisher matrix of feature models is robust to changes in the phase for different $\omega$ values of interest. Feature model bispectra with a specific phase $\phi$ can be represented as a sum of sine and cosine ones as in (6). Hence, its Fisher matrix is given by

[TABLE]

where $F_{ss}$ is the element $F_{ij}$ of the Fisher matrix in (14) with reduced bispectra $b^{(i)}=b^{(j)}=b^{\sin}$, and so on. Correlation between sine and cosine templates can be expressed as $F_{sc}/(F_{ss}F_{cc})^{1/2}$, and this value can be learned from analysing the $\phi$ dependence of $F(\omega,\phi)$.

Figure 1 shows forecast error bars for the full phase range $[0,\pi]$ in the most sensitive experiment specification of 1’ beam and 1$\mu K\cdot arcmin$ noise. The forecast $\sigma$ varies within 1% level for every $\omega\geq 20$. In terms of the Fisher matrix, the cross term $F_{sc}$ was 2-3 orders of magnitude smaller than $F_{ss}$ and $F_{cc}$ for all cases. In other words, correlation between the sine and cosine templates was smaller than 1%. This justifies our previous choice of constraining $f_{NL}^{\sin}$ and $f_{NL}^{\cos}$ separately. We now focus our attention to $\sigma_{\sin}$ in future discussions.

For smaller values of $\omega$, the phase affects the error bar primarily through modulating the amplitude of the acoustic oscillations in the CMB itself. The radiation transfer functions are non-zero for $k$ values in $0-0.8~{}\text{Mpc}^{-1}$. The argument $\omega k$ covers less than two full periods in this $k$ range if $\omega\leq 10\;\text{Mpc}$, and phase has a direct influence on the amplitude of the acoustic peaks. In the extreme case of $\omega=0$, bispectrum vanishes completely for the $\sin$ feature model. Variations in the overall bispectrum amplitude therefore result in varying Fisher information for low frequencies.

$l_{max}$ dependence

Figure 2 shows the graph of forecast error bar $\sigma_{\sin}^{T+E}$ as we increase $l_{max}$. The forecasts were done within angular scale range $2\leq l\leq l_{max}$, the oscillation frequency $\omega$ set to 100, and assuming 1’ beam and 1$\mu K\cdot$arcmin noise. The Planck noise curves were approximated by ones for 5’ beam and 47 $\mu K\cdot$arcmin noise for this plot only, since we extend $l_{max}$ to 4000 here.

The Planck error bar essentially stalls out when $l_{max}$ reaches 2000. The forecast error bar, on the other hand, keeps decreasing until $l_{max}=4000$ thanks to the improved sensitivity in measuring small scale, or large $l$’s. Despite the information loss due to smaller sky coverage $f_{sky}$, the forecast error bar reduces to about 42% of Planck by $l_{max}=4000$. This corresponds to a factor of 2.4 times improvement to measurement precision on $f_{\text{NL}}$.

Beam and noise dependence

We explore how different beam widths and noise levels affect the forecast error bars in this section. Figure 3 shows forecast $\sigma_{\sin}^{T+E}$ for ranges of beam and noise levels. Their oscillation frequencies are also varied, but only two representatives $\omega=20$ and $2000$ are chosen here. Forecasts for the other values of $\omega$ also show similar dependences on beam width and noise level.

First of all, note that all estimated error bars in the plot are smaller than Planck, for which $\sigma_{\sin}^{T+E}=34$ when $\omega=20$ and $\sigma_{\sin}^{T+E}=610$ when $\omega=2000$. In fact even the least sensitive CMB-S4 specification of 5’ beam and 9$\mu K\cdot$arcmin noise is expected to put better bounds on feature models.

Wider beams and noisier detectors provide less signal and thus larger error bars, as expected. In this range of beam width and noise levels, noise has a bigger effect on the forecast; experiments with 1’ beam and 5$\mu K\cdot$arcmin noise yields larger error bars than the ones with 5’ beam with 1$\mu K\cdot$arcmin noise. Between the most sensitive specification of 1’ beam and 1$\mu K\cdot$arcmin and the least sensitive one with 5’ beam and 9$\mu K\cdot$arcmin, $\sigma_{\sin}$ differs by a factor of 1.6.

Oscillation frequency dependence

We now present main results of the forecast. Figure 4 summarises the $\sigma_{\sin}$ forecasts for several different CMB-S4 preliminary specifications, including the Simons Observatory (SO) baseline and goal. Note that the 1/f noise effects are incorporated in SO forecasts but not in other ones. We also provide 1$\sigma$ errors for joint estimators, for which Planck signals from the fraction of the sky not covered by CMB-S4 are combined via $\sigma_{\text{joint}}^{-2}=\sigma_{\text{CMB-S4}}^{-2}+\sigma_{\text{Planck}}^{-2}$. This method is not statistically optimal but sufficient to give an idea of the joint estimation power.

The most sensitive setup with 1’ beam and 1$\mu K\cdot$arcmin noise would yield error bars that are 47-62% of Planck, depending on the oscillation frequency $\omega$. These correspond to a factor of 1.6-2.1 improvement. Relatively smaller improvements are made for high oscillation frequencies. They correspond to smaller momentum scales $k_{\ast}=2\pi/3\omega$, or larger angular scales, which benefit less from the increased sensitivity of CMB-S4 experiments. When the results are combined with Planck the error bar further reduces to 45-57% of Planck, which is a factor of 1.7-2.2 improvement.

Forecast error bars from the SO baseline specification and the more ambitious one do not differ very much. Quoting in terms of the baseline values, $\sigma_{\sin}$ lies about 68-86% of that of Planck or equivalently, 1.2-1.5 times smaller than Planck. Numbers change to 62-74% when combined with Planck, so that the overall improvement ratio is about 1.3-1.6.

Figure 5 shows the results when only the CMB temperature data are used in the forecast. CMB-S4 would in fact be worse than Planck in terms of constraining $f_{NL}^{\text{feat}}$ for this case. The loss in information due to less sky coverage overwhelms the increased sensitivity. We see again that the real strength of CMB-S4 experiments lies in measuring the CMB polarisation.

Then how much information do we actually gain from adding E-mode polarisation? Figure 6 shows the ratio of $\sigma_{\sin}$’s between the temperature-only (T) and polarisation-included (T+E) analyses. The forecast error bars reduces up to 4.6 times smaller when the polarisation information is added, which is much larger than the corresponding Planck value of 2.2. The ratio decreases overall when the joint statistics with Planck are considered. An intriguing feature of this plot is that the ratio is maximised around $\omega=200$ before it starts dropping again.

In order to gain insight on this behaviour, we performed some simplified computations using the power spectrum. We imposed oscillations on the primordial power spectrum as $P^{\prime}(k):=P(k)(1+\sin(2\omega k+\phi))$, which is just like our feature model bispectrum template but with $\omega(k_{1}+k_{2}+k_{3})$ replaced by $\omega(k+k)$. $P^{\prime}(k)$ is then projected to the late-time harmonic space using the transfer functions;

[TABLE]

We observed that the fractional variation $(C_{l}^{\prime}-C_{l})/C_{l}$ displays some oscillations in $l$, and the largest contribution comes from a term $\propto\sin(2\omega l/\Delta\tau)$ where $\Delta\tau$ represents the conformal distance to last scattering surface. This fact can be explained by approximating the transfer function as $\Delta_{l}(k)\approx(1/3)j_{l}(k\Delta\tau)$ and noting that the spherical Bessel function has a sharp peak at $l$ for large $l$’s. The integral in (30) therefore picks up a term proportional to $\sin(2\omega l/\Delta\tau)$.

The amplitude of these ‘maximal’ oscillations in $(C_{l}^{\prime}-C_{l})/C_{l}$ were then computed using discrete Fourier transform for different values of oscillation scale $\omega$ and two different phases $\phi=0,\pi/2$ (i.e. sine and cosine). The results are shown in Figure 7. Some extra wiggles to the graph come from the phase of oscillations imposed; we indeed see that graphs of sine and cosine oscillate between each other. Some peak features near $\omega\approx 70$ and 140 arise from resonances with Baryonic Acoustic Oscillations.

We can think of the computed amplitude as a measure of information $C_{l}$’s contain about primordial oscillations. First of all, note that the amplitude in all four plots generally decreases as $\omega$ grows. Previously in Figure 4 we saw that the amount of information obtained from the CMB is smaller for larger $\omega$’s, consistent with what can be said from the amplitude analysis. Moreover, the amplitudes for the EE mode are generally larger than the TT mode ones, and their difference is the largest in the $\omega$ range of 70 to 300. This could serve as a heuristic explanation for the improvement in forecast error bars from including polarisation data being maximised around $\omega=200$, as depicted in Figure 6.

Comparison to scale invariant models

Our pipeline for forecasting $f_{NL}^{\text{feat}}$ also yields forecasts for $f_{NL}$ of the constant model. Constant models are scale invariant and have a trivial shape, so that $B(k_{1},k_{2},k_{3})\propto(k_{1}k_{2}k_{3})^{-2}$. Forecasts on $f_{NL}^{\text{const}}$ follow from out pipeline by simply setting oscillation frequency$\omega=0$ and phase $\phi=\pi/2$. Table 1 summarises the forecast results for several different CMB-S4 specifications mentioned before, using both T and E data and in combination with Planck data from the regions of the sky not covered by CMB-S4. For the 1’ beam and 1$\mu K\cdot$arcmin noise setup, the error bar is expected to be reduced by a factor of 2.3 compared to Planck.

The latest Planck constraints on $f_{\text{NL}}$ of some popular bispectrum templates are given by $f_{\text{NL}}^{\text{local}}=2.5\pm 5.7$, $f_{\text{NL}}^{\text{equil}}=-16\pm 70$, and $f_{\text{NL}}^{\text{ortho}}=-34\pm 33$ Planck Collaboration (2015). CMB-S4 experiments are expected to yield better estimates on these as well. Table 2 summarises the forecast improvement ratio given in Abazajian (2016) together with the constant and feature model ratios computed in this work.

To the authors’ surprise, the estimation error for feature models does not improve as much as other templates. Feature models benefit much more from polarisation data than other scale independent shapes; for example, $\sigma^{T}/\sigma^{T+E}=4.6$ for the feature model with $\omega=200$ in CMB-S4, while the value equals 2.8 for the constant model. Because CMB-S4 would have significantly enhanced polarisation measurement sensitivity, we originally expected the feature models to be constrained significantly better than Planck.

In order to investigate this lack of improvement, we performed a breakdown analysis on the improvements gained from CMB-S4 temperature and polarisation; we computed $\sigma(f_{\text{NL}})$ for the constant and feature models using each of the four combinations of Planck / CMB-S4 noise curves for temperature / polarisation (e.g. Planck T + CMB-S4 E). The results are summarised in Table 3.

We see that the constraints on feature models improve by a factor of 1.7 when swapping Planck polarisation noises with the CMB-S4 ones. This factor is indeed larger than that of the constant model, which equals 1.6. The difference is however not significant. It seems that the amount of feature signals in polarisation data left unexplored by Planck is not tremendously large compared to the constant model. The feature model improves less than the constant model when the temperature measurements are enhanced. In fact, for feature models the signal loss from smaller sky fraction $f_{sky}$ eclipses the signal gain from more sensitive temperature measurements. This lack of improvements from temperature causes the full CMB-S4 constraints on the feature model not to improve as much as the constant model overall.

V Conclusion

Upcoming CMB Stage-4 experiments will provide an opportunity to measure CMB temperature and polarisation with greater precision. The estimation of primordial non-Gaussianity parameters would greatly benefit from the improvement in measurement sensitivity. In this research we made forecasts on $f_{\text{NL}}$ for the feature models, which have not been done so far despite the growing interests on inflation models with primordial oscillations. For efficient forecasts we simplified the bispectrum estimator for $f_{\text{NL}}$ by orthonormalising the covariance matrix, further optimising the computation. When the most sensitive CMB Stage-4 experiment specification of 1’ beam and 1$\mu K$arcmin noise is concerned, we expect a factor of 1.7-2.2 times more stringent constraints compared to Planck. Under realistic Simons Observatory conditions the improvement would be about 1.3-1.6 to Planck.

Although this is not a massive boost in the estimation power, we can hope to verify current 4$\sigma$-level signals found in the 2015 Planck analysis. It is also worth noting that the CMB-S4 experiments would allow us to explore higher $l$ modes, especially since localised oscillations in this range are currently unconstrained. Moreover, though we have only considered linearly spaced oscillations in this work, we expect even better improvements on the models inducing log spaced oscillations. Higher $l$ modes would promote the constraining power as the oscillation slows down in small scales for this type of models. Lastly, cross-validation using these new statistically independent modes would be useful.

We also extensively studied how the forecasts depend on various parameters. Frequency dependences of the ratio between T and T+E forecasts were particularly illuminating - the improvement from adding polarisation information is maximised around $\omega=200$. Some simplified calculations were presented to heuristically address this fact. Even though the estimation power on feature models massively benefit from the polarisation data, overall expected improvements compared to Planck are quite underwhelming. Breakdown analysis on temperature and polarisation contribution revealed that the feature models would indeed improve more than other scale-independent models if only the polarisation measurement sensitivity is enhanced to the CMB-S4 standards. However, boosts in the temperature measurements affect scale-independent models more so that they gain more information overall.