© ESO, 2014

1. Introduction

The cosmic microwave background (CMB) provides the most direct experimental access to the statistics of the primordial curvature perturbations. In standard single-field, slow roll inflation, these perturbations are Gaussian to a very good approximation (Maldacena 2003). However, more complicated models of inflation often predict detectable amounts of non-Gaussianities of various shapes. In particular, the primordial bispectrum B(k₁,k₂,k₃) arising from the three-point correlations of the curvature field ⟨ΦΦΦ⟩ is a sensitive probe to discriminate among models of inflation (see e.g. the reviews Liguori et al. 2010; Komatsu 2010; Yadav & Wandelt 2010). High precision measurements of the CMB recently provided by the Planck experiment (Planck Collaboration XXIV 2014) have been consistent with Gaussianity and set stronger limits on primordial non-Gaussianities. However, the availability of such limits depend on the specific shape of the bispectra under consideration.

A class of bispectra that has attracted attention in recent years are oscillating shapes. Such bispectrum oscillations can arise in a variety of theoretical models. The authors of Chen et al. (2008) calculated the primordial bispectrum in the presence of features in the inflaton potential of standard single field inflation. They provided two analytical bispectrum shapes that approximate their results. The feature model oscillates linearly with the scale and is induced by sharp features in the inflaton potential. The resonance model includes oscillations with the logarithm of the scale and is induced by periodic features in the inflaton potential. More recently, the authors of Bartolo et al. (2010) used the effective field theory of inflation to examine the influence of sharp features in the inflaton potential, which also provide oscillating bispectrum solutions. Periodical modulations of the inflaton potential appear, for example, in axion monodromy inflation models (Flauger et al. 2010; Flauger & Pajer 2011). Certain bispectrum shapes motivated by Non-Bunch-Davis vacua also include oscillations (e.g. Chen et al. 2007; Meerburg et al. 2009), as do cascade inflation models (Ashoorioon & Krause 2006). A transient reduction in the speed of sound also leads to oscillations (Achucarro et al. 2014b,a). Oscillations in the bispectrum are usually accompanied by oscillations in the primordial power spectrum. Recent searches for power spectrum oscillation were presented in Pahud et al. (2009); Peiris et al. (2013); Planck Collaboration XXII (2014); Easther & Flauger (2014); Meerburg et al. (2014); Meerburg & Spergel (2014), but no statistically significant result has been found. Combining power spectrum and bispectrum measurements can improve the sensitivity.

In the present work, we focus on the feature model shape, since it is simple and approximates some more complicated oscillating shapes. In particular, the feature model has the important property of separability. The Planck paper on non-Gaussianities (Planck Collaboration XXIV 2014) already included a targeted search for the feature and resonant shapes with the modal expansion method. In this methodology, the bispectrum under consideration is expanded into a basis of separable shapes (Fergusson et al. 2010), which allows an efficient numerical estimation by the Komatsu Spergel Wandelt (KSW) estimator (Komatsu et al. 2005). However, the separable basis functions used by Planck did not allow to represent high frequency oscillations, limiting the frequency range, which could be searched for oscillations. The situation can be improved by using a set of oscillating basis functions for the modal expansion (Meerburg 2010). In this work, we follow a more direct approach by writing the feature model bispectrum in separable form, which makes it possible to search for oscillations of arbitrary frequency, which are only limited by the resolution of the maps. The gain in computational time with respect to the modal expansion is of the order of the number of modes that would be necessary to approximate the shape. We study the properties of the estimator in detail, including the precision with which phase and frequency of the oscillation can be determined. We also give an illuminating interpretation of the KSW estimator for oscillations in position space.

2. Oscillations in the primordial bispectrum

2.1. Bispectrum shape and experimental constraints

The general translation and rotation invariant primordial bispectrum of the curvature potential Φ can be written as $\mathrm{⟨}\mathrm{\Phi }\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\mathrm{\Phi }\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\mathrm{\Phi }\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\mathrm{\left(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}\mathit{\delta }\mathrm{\right(}{{k}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}\mathrm{\left)}{\mathit{f}}_{\mathrm{NL}}{\mathit{B}}_{\mathrm{\Phi }}\mathrm{\right(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{)}\mathit{,}$(1)where the bispectrum B_Φ is a function of the magnitude of the wave numbers k_i and f_NL is the amplitude of the bispectrum. The feature model bispectrum that we are primarily interested in is ${\mathit{B}}_{\mathrm{\Phi }}^{\mathrm{feat}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{6}{\mathrm{\Delta }}_{\mathrm{\Phi }}^{\mathrm{2}}{\mathit{f}}_{\mathrm{NL}}}{\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{2}}}\sin \left(\frac{\mathrm{2}\mathit{\pi K}}{\mathrm{3}{\mathit{k}}_{\mathrm{c}}}\mathrm{+}\mathit{\phi }\right)\mathit{.}$(2)It oscillates linearly with the mean, $\mathit{K}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\mathrm{(}{\mathit{k}}_{\mathrm{1}}\mathrm{+}{\mathit{k}}_{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{3}}\mathrm{)}$, of the wave numbers. Here Δ_Φ is the primordial power spectrum amplitude and the 1 /k⁶ factor compensates for the phase space factor. The bispectrum is parametrised by the amplitude f_NL by the oscillation scale k_c and by the phase φ. The oscillation scale k_c implies an efficient multipole periodicity of l_c ≃ k_c[τ − τ_rec], where τ − τ_rec is the conformal distance to recombination.

Planck has searched for the feature bispectrum shape for sample frequencies in the range 0.01 <k_c< 0.1 at four different phases φ = 0,π/ 4,π/ 2,3π/ 4Planck Collaboration XXIV (2014). Here, k_c = 0.01 corresponds to an effective multipole periodicity l_c = 140. The best fit model has k_c = 0.0185 (l_c = 260) and phase Φ = 0 with a significance of 3σ. This may correspond to weak hints for oscillation in the full bispectrum reconstruction, which were found for l< 500. However, the statistical significance becomes much lower when one takes the number of statistically uncorrelated feature models into account that were searched. We note that the range of the oscillation frequency was constrained because of the limitations of the analysis method. With the present work, we target the unexplored range k_c< 0.01 in particular.

2.2. Separability

For convenience, we rewrite the feature model (2) as a sum of sine and cosine contributions as ${\mathit{B}}_{\mathrm{\Phi }}^{\mathrm{feat}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{6}{\mathrm{\Delta }}_{\mathrm{\Phi }}^{\mathrm{2}}}{\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{2}}}\left[{\mathit{f}}_{\mathrm{1}}\sin \left(\frac{\mathrm{2}\mathit{\pi K}}{\mathrm{3}{\mathit{k}}_{\mathrm{c}}}\right)\mathrm{+}{\mathit{f}}_{\mathrm{2}}\cos \left(\frac{\mathrm{2}\mathit{\pi K}}{\mathrm{3}{\mathit{k}}_{\mathrm{c}}}\right)\right]\mathit{,}$(3)where ${\mathit{f}}_{\mathrm{NL}}\mathrm{=}\sqrt{{\mathit{f}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{f}}_{\mathrm{2}}^{\mathrm{2}}}$ and $\mathrm{\Phi }\mathrm{=}\arctan \mathrm{\left(}\frac{{\mathit{f}}_{\mathrm{2}}}{{\mathit{f}}_{\mathrm{1}}}\mathrm{\right)}$. This can be written in separated form as $\begin{array}{ccc}{\mathit{B}}_{\mathrm{\Phi }}^{\mathrm{feat}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}& \mathrm{=}& \mathrm{6}{\mathit{f}}_{\mathrm{1}}[\mathrm{-}\mathit{X}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\mathit{X}\mathrm{\left(}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}\mathit{X}\mathrm{\left(}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{+}(\mathit{X}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\mathit{Y}\mathrm{\left(}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}\mathit{Y}\mathrm{\left(}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\\ & & \mathrm{+}\mathrm{2}\mathrm{perm}\mathrm{.})]\mathrm{+}\mathrm{6}{\mathit{f}}_{\mathrm{2}}[\mathit{Y}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\mathit{Y}\mathrm{\left(}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}\mathit{Y}\mathrm{\left(}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\\ & & \mathrm{-}(\mathit{X}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\mathit{X}\mathrm{\left(}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}\mathit{Y}\mathrm{\left(}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{+}\mathrm{2}\mathrm{perm}\mathrm{.})]\mathit{,}\end{array}$(4)where we have defined $\mathit{X}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\Delta }}_{\mathrm{\Phi }}^{\mathrm{2}\mathit{/}\mathrm{3}}}{{\mathit{k}}^{\mathrm{2}}}\sin \left(\frac{\mathrm{2}\mathit{\pi k}}{\mathrm{3}{\mathit{k}}_{\mathrm{c}}}\right)\mathrm{and}\mathit{Y}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\Delta }}_{\mathrm{\Phi }}^{\mathrm{2}\mathit{/}\mathrm{3}}}{{\mathit{k}}^{\mathrm{2}}}\cos \left(\frac{\mathrm{2}\mathit{\pi k}}{\mathrm{3}{\mathit{k}}_{\mathrm{c}}}\right)\mathrm{·}$(5)This separability property allows efficient computation and estimation of the bispectrum. One may also include an exponential decay factor of the form $\exp \left(\mathrm{-}\frac{\mathit{K}}{\mathit{\mu }}\right)$, while retaining separability with identical formulas up to trivial replacements.

3. Oscillations in the CMB bispectrum

From the separable expression (4) for the primordial bispectrum, one can calculate the CMB bispectrum with the standard line-of-sight integration method Seljak & Zaldarriaga (1996). The reduced bispectrum is then $\begin{array}{ccc}{\mathit{b}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}& \mathrm{=}& \mathrm{6}{\mathit{f}}_{\mathrm{1}}\mathrm{\int }\mathrm{d}\mathit{r}{\mathit{r}}^{\mathrm{2}}[\mathrm{-}{\mathit{X}}_{{\mathit{l}}_{\mathrm{1}}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{X}}_{{\mathit{l}}_{\mathrm{2}}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{X}}_{{\mathit{l}}_{\mathrm{3}}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}\mathrm{\left(}{\mathit{X}}_{{\mathit{l}}_{\mathrm{1}}}\mathrm{\right(}\mathit{r}\mathrm{\left)}{\mathit{Y}}_{{\mathit{l}}_{\mathrm{2}}}\mathrm{\right(}\mathit{r}\mathrm{\left)}{\mathit{Y}}_{{\mathit{l}}_{\mathrm{3}}}\mathrm{\right(}\mathit{r}\mathrm{)}\\ & & \mathrm{+}\mathrm{2}\mathrm{perm}\mathrm{.}\mathrm{)}]\mathrm{+}\mathrm{6}{\mathit{f}}_{\mathrm{2}}\mathrm{\int }\mathrm{d}\mathit{r}{\mathit{r}}^{\mathrm{2}}[{\mathit{Y}}_{{\mathit{l}}_{\mathrm{1}}}\mathrm{(}\mathit{r}\mathrm{\left)}{\mathit{Y}}_{{\mathit{l}}_{\mathrm{2}}}\mathrm{\right(}\mathit{r}\mathrm{\left)}{\mathit{Y}}_{{\mathit{l}}_{\mathrm{3}}}\mathrm{\right(}\mathit{r}\mathrm{)}\\ & & \mathrm{-}\mathrm{\left(}{\mathit{X}}_{{\mathit{l}}_{\mathrm{1}}}\mathrm{\right(}\mathit{r}\mathrm{\left)}{\mathit{X}}_{{\mathit{l}}_{\mathrm{2}}}\mathrm{\right(}\mathit{r}\mathrm{\left)}{\mathit{Y}}_{{\mathit{l}}_{\mathrm{3}}}\mathrm{\right(}\mathit{r}\mathrm{)}\mathrm{+}\mathrm{2}\mathrm{perm}\mathrm{.}\mathrm{)}]\mathit{.}\end{array}$(6)where we have defined the functions $\begin{array}{ccc}{\mathit{X}}_{\mathit{l}}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}\frac{\mathrm{2}}{\mathit{\pi }}\mathrm{\int }\mathrm{d}\mathit{k}{\mathit{k}}^{\mathrm{2}}\mathit{X}\mathrm{\left(}\mathit{k}\mathrm{\right)}{\mathit{j}}_{\mathit{l}}\mathrm{\left(}\mathit{kr}\mathrm{\right)}{\mathrm{\Delta }}_{\mathit{l}}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \\ & \mathrm{=}\frac{\mathrm{2}}{\mathit{\pi }}\mathrm{\int }\mathrm{d}\mathit{k}{\mathrm{\Delta }}_{\mathrm{\Phi }}^{\mathrm{2}\mathit{/}\mathrm{3}}\sin \mathrm{(}{\frac{\mathrm{2}\mathit{\pi k}}{\mathrm{3}{\mathit{k}}_{\mathrm{c}}}}^{\mathrm{)}}{\mathit{j}}_{\mathit{l}}\mathrm{\left(}\mathit{kr}\mathrm{\right)}{\mathrm{\Delta }}_{\mathit{l}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}& \\ {\mathit{Y}}_{\mathit{l}}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}\frac{\mathrm{2}}{\mathit{\pi }}\mathrm{\int }\mathrm{d}\mathit{k}{\mathit{k}}^{\mathrm{2}}\mathit{Y}\mathrm{\left(}\mathit{k}\mathrm{\right)}{\mathit{j}}_{\mathit{l}}\mathrm{\left(}\mathit{kr}\mathrm{\right)}{\mathrm{\Delta }}_{\mathit{l}}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \\ & \mathrm{=}\frac{\mathrm{2}}{\mathit{\pi }}\mathrm{\int }\mathrm{d}\mathit{k}{\mathrm{\Delta }}_{\mathrm{\Phi }}^{\mathrm{2}\mathit{/}\mathrm{3}}\cos \mathrm{(}{\frac{\mathrm{2}\mathit{\pi k}}{\mathrm{3}{\mathit{k}}_{\mathrm{c}}}}^{\mathrm{)}}{\mathit{j}}_{\mathit{l}}\mathrm{\left(}\mathit{kr}\mathrm{\right)}{\mathrm{\Delta }}_{\mathit{l}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}& \end{array}$and where j_l are spherical Bessel functions and Δ_l are the CMB transfer functions that we evaluate numerically with CAMB Lewis et al. (2000). We note that these functions depend on k_c and have to be evaluated for each oscillation frequency of interest. Both the transfer functions Δ_l(k) (in k space) and the Bessel functions j_l are highly oscillatory integrals, which must be evaluated with sufficient sampling. Examples of the X(r) and Y(r) functions are given in Fig. 1.

Fig. 1 Function X and Y of the feature model for kc = 0.01 (top) and kc = 0.001 (bottom) as a function of l for r = τ0 − τrec. X and Y have units of Mpc-1.

To calculate the resulting CMB bispectrum, we perform the r integral in Eq. (6). We choose a quadrature of about 2000 points in r with a higher sampling in the range of recombination. The contribution of different values of r is examined in more detail in Sect. 6. To plot the bispectrum, it is convenient to normalise by the constant bispectrum with natural k^-6 scaling, as proposed in Fergusson & Shellard (2009). The constant primordial bispectrum is given by ${\mathit{B}}_{\mathrm{\Phi }}^{\mathrm{const}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}$(9)and its large angle Sach-Wolfe CMB solution is Fergusson & Shellard (2009)$\begin{array}{ccc}{\mathit{b}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}^{\mathrm{const}}& \mathrm{=}& {\left(\frac{\mathrm{1}}{\mathrm{3}}\right)}^{\mathrm{3}}\frac{\mathrm{1}}{\mathrm{(}\mathrm{2}{\mathit{l}}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{l}}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{l}}_{\mathrm{3}}\mathrm{+}\mathrm{1}\mathrm{)}}\\ & & \mathrm{\times }\left[\frac{\mathrm{1}}{{\mathit{l}}_{\mathrm{1}}\mathrm{+}{\mathit{l}}_{\mathrm{2}}\mathrm{+}{\mathit{l}}_{\mathrm{3}}\mathrm{+}\mathrm{3}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{l}}_{\mathrm{1}}\mathrm{+}{\mathit{l}}_{\mathrm{2}}\mathrm{+}{\mathit{l}}_{\mathrm{3}}}\right]\mathrm{·}\end{array}$(10)In Fig. 2, we show a simple 1-dimensional visualisation, the equal l bispectrum b_lll normalised by ${\mathit{b}}_{\mathit{lll}}^{\mathrm{const}}$ for different frequencies and phases of the feature model.

Fig. 2 Equal l reduced bispectrum normalised by the large angle solution of the constant bispectrum for feature models with different parameters.

4. The KSW estimator for the oscillating bispectrum

4.1. KSW estimator

The optimal KSW estimator for a sum of bispectra bⁱ in the presence of a mask and non-uniform noise is (Komatsu et al. 2005; Babich 2005; Komatsu 2010) ${\mathit{f}}_{\mathit{i}}\mathrm{=}\sum _{\mathit{j}}{\mathrm{\left(}{\mathit{F}}^{-1}\mathrm{\right)}}_{\mathit{ij}}{\mathit{S}}_{\mathit{j}}\mathit{,}$(11)with $\begin{array}{ccc}{\mathit{S}}_{\mathit{i}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{6}}\sum _{\mathit{lm}}{\mathit{G}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}^{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}{\mathit{b}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}^{\mathit{i}}[\mathrm{\left(}{\mathit{C}}^{-1}\mathit{a}{\mathrm{)}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}\mathrm{\right(}{\mathit{C}}^{-1}\mathit{a}{\mathrm{)}}_{{\mathit{l}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\mathrm{(}{\mathit{C}}^{-1}\mathit{a}{\mathrm{)}}_{{\mathit{l}}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}\\ & & \mathrm{-}\mathrm{3}{\mathrm{(}{\mathit{C}}^{-1}\mathrm{)}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}\mathit{,}{\mathit{l}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\mathrm{(}{\mathit{C}}^{-1}\mathit{a}{\mathrm{)}}_{{\mathit{l}}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}]\end{array}$(12)and with the Fisher matrix given by $\begin{array}{ccc}{\mathit{F}}_{\mathit{ij}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{6}}\sum _{\mathit{lm}}\sum _{{\mathit{l}}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}{\mathit{G}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}^{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}{\mathit{b}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}^{\mathit{i}}\mathrm{\left(}{\mathit{C}}^{-1}{\mathrm{)}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}\mathit{,}{\mathit{l}}_{\mathrm{1}}^{\mathrm{\prime }}{\mathit{m}}_{\mathrm{1}}^{\mathrm{\prime }}}\mathrm{\right(}{\mathit{C}}^{-1}{\mathrm{)}}_{{\mathit{l}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}\mathit{,}{\mathit{l}}_{\mathrm{2}}^{\mathrm{\prime }}{\mathit{m}}_{\mathrm{2}}^{\mathrm{\prime }}}\\ & & \mathrm{\times }\mathrm{(}{\mathit{C}}^{-1}{\mathrm{)}}_{{\mathit{l}}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}\mathit{,}{\mathit{l}}_{\mathrm{3}}^{\mathrm{\prime }}{\mathit{m}}_{\mathrm{3}}^{\mathrm{\prime }}}{\mathit{b}}_{{\mathit{l}}_{\mathrm{1}}^{\mathrm{\prime }}{\mathit{l}}_{\mathrm{2}}^{\mathrm{\prime }}{\mathit{l}}_{\mathrm{3}}^{\mathrm{\prime }}}^{\mathit{j}}{\mathit{G}}_{{\mathit{l}}_{\mathrm{1}}^{\mathrm{\prime }}{\mathit{l}}_{\mathrm{2}}^{\mathrm{\prime }}{\mathit{l}}_{\mathrm{3}}^{\mathrm{\prime }}}^{{\mathit{m}}_{\mathrm{1}}^{\mathrm{\prime }}{\mathit{m}}_{\mathrm{2}}^{\mathrm{\prime }}{\mathit{m}}_{\mathrm{3}}^{\mathrm{\prime }}}\mathit{.}\end{array}$(13)In the present case of an oscillation with a single scale k_c, the bispectrum is a sum b = f₁b¹ + f₂b² of sine and cosine components, as given by Eq. (6).

To numerically evaluate the KSW estimator terms S_i we need the weighted maps $\begin{array}{ccc}{\mathit{M}}_{\mathit{X}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{lm}}\mathrm{\left(}{\mathit{C}}^{-1}\mathit{a}{\mathrm{)}}_{\mathit{lm}}{\mathit{X}}_{\mathit{l}}\mathrm{\right(}\mathit{r}\mathrm{\left)}{\mathit{Y}}_{\mathit{lm}}\mathrm{\right(}\mathit{n̂}\mathrm{)}\mathit{,}\\ {\mathit{M}}_{\mathit{Y}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{lm}}\mathrm{\left(}{\mathit{C}}^{-1}\mathit{a}{\mathrm{)}}_{\mathit{lm}}{\mathit{Y}}_{\mathit{l}}\mathrm{\right(}\mathit{r}\mathrm{\left)}{\mathit{Y}}_{\mathit{lm}}\mathrm{\right(}\mathit{n̂}\mathrm{)}\mathit{.}\end{array}$(14)The cubic KSW estimator, which is exact in the case of a full sky observation, is given in terms of these maps by $\begin{array}{ccc}{\mathit{S}}_{\mathrm{1}}^{\mathrm{cub}}& \mathrm{=}& \mathrm{\int }{\mathit{r}}^{\mathrm{2}}\mathrm{d}\mathit{r}\mathrm{\int }\mathrm{d\Omega }\mathrm{[}\mathrm{-}{\mathit{M}}_{\mathit{X}}^{\mathrm{3}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}\mathrm{+}{\left(\mathrm{3}{\mathit{M}}_{\mathit{X}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}{\mathit{M}}_{\mathit{Y}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}{\mathit{M}}_{\mathit{Y}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}\right)}^{\mathrm{]}}\mathit{,}\\ {\mathit{S}}_{\mathrm{2}}^{\mathrm{cub}}& \mathrm{=}& \mathrm{\int }{\mathit{r}}^{\mathrm{2}}\mathrm{d}\mathit{r}\mathrm{\int }\mathrm{d\Omega }\mathrm{[}{\mathit{M}}_{\mathit{Y}}^{\mathrm{3}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}\mathrm{-}{\left(\mathrm{3}{\mathit{M}}_{\mathit{X}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}{\mathit{M}}_{\mathit{X}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}{\mathit{M}}_{\mathit{Y}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}\right)}^{\mathrm{]}}\mathit{.}\end{array}$(15)Partial sky coverage can be taken into account by incorporating the linear term of the KSW estimator, so that ${\mathit{S}}_{\mathit{i}}\mathrm{=}{\mathit{S}}_{\mathit{i}}^{\mathrm{cub}}\mathrm{+}{\mathit{S}}_{\mathit{i}}^{\mathrm{lin}}$ with $\begin{array}{ccc}{\mathit{S}}_{\mathrm{1}}^{\mathrm{lin}}& \mathrm{=}& \mathrm{-}\mathrm{3}\mathrm{\int }{\mathit{r}}^{\mathrm{2}}\mathrm{d}\mathit{r}\mathrm{\int }\mathrm{d\Omega }[\mathrm{-}{\mathit{M}}_{\mathit{X}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}\mathrm{⟨}{\mathit{M}}_{\mathit{X}}^{\mathrm{2}}\mathrm{(}\mathit{r,}\mathit{n̂}{\mathrm{)}}^{\mathrm{⟩}}\\ & & \mathrm{+}{\mathit{M}}_{\mathit{X}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}\mathrm{⟨}{\mathit{M}}_{\mathit{Y}}^{\mathrm{2}}\mathrm{(}\mathit{r,}\mathit{n̂}{\mathrm{)}}^{\mathrm{⟩}}\mathrm{+}\mathrm{2}{\mathit{M}}_{\mathit{Y}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}⟨{\mathit{M}}_{\mathit{X}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}{\mathit{M}}_{\mathit{Y}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}⟩]\mathit{,}\\ {\mathit{S}}_{\mathrm{2}}^{\mathrm{lin}}& \mathrm{=}& \mathrm{-}\mathrm{3}\mathrm{\int }{\mathit{r}}^{\mathrm{2}}\mathrm{d}\mathit{r}\mathrm{\int }\mathrm{d\Omega }[{\mathit{M}}_{\mathit{Y}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}\mathrm{⟨}{\mathit{M}}_{\mathit{Y}}^{\mathrm{2}}\mathrm{(}\mathit{r,}\mathit{n̂}{\mathrm{)}}^{\mathrm{⟩}}\\ & & \mathrm{-}{\mathit{M}}_{\mathit{Y}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}\mathrm{⟨}{\mathit{M}}_{\mathit{X}}^{\mathrm{2}}\mathrm{(}\mathit{r,}\mathit{n̂}{\mathrm{)}}^{\mathrm{⟩}}\mathrm{-}\mathrm{2}{\mathit{M}}_{\mathit{X}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}⟨{\mathit{M}}_{\mathit{X}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}{\mathit{M}}_{\mathit{Y}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}⟩]\mathit{,}\end{array}$(16)where the expectation values have to be evaluated by Monte Carlo averaging over Gaussian realisations drawn with the same beam, mask, and noise properties as expected in the data. To make the KSW estimator optimal for a non-uniform sky coverage, it is necessary to perform an inverse covariance weighting with the non-diagonal covariance matrix. This is a computationally challenging problem (Smith et al. 2009; Elsner & Wandelt 2012). It was noted in Planck Collaboration XXIV (2014) that one can also achieve excellent results by assuming a diagonal covariance matrix Ĉ_l = C_l + N_l, where N_l assumes homogeneous noise, and by using a diffusive inpainting on the masked areas. In this approximation, the Fisher matrix scales proportionally to the visible fraction of the sky f_sky.

For the remainder of this paper, we assume full sky coverage so that $\begin{array}{ccc}{\mathit{S}}_{\mathit{i}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{6}}\sum _{\mathit{lm}}\frac{{\mathit{G}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}^{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}{\mathit{b}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}^{\mathit{i}}}{{\mathit{C}}_{{\mathit{l}}_{\mathrm{1}}}{\mathit{C}}_{{\mathit{l}}_{\mathrm{2}}}{\mathit{C}}_{{\mathit{l}}_{\mathrm{3}}}}{\mathit{a}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}{\mathit{a}}_{{\mathit{l}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}{\mathit{a}}_{{\mathit{l}}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}\mathit{,}\\ {\mathit{F}}_{\mathit{ij}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{6}}\sum _{\mathit{l}}{\mathit{I}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}\frac{{\mathit{b}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}^{\mathit{i}}{\mathit{b}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}^{\mathit{j}}}{{\mathit{C}}_{{\mathit{l}}_{\mathrm{1}}}{\mathit{C}}_{{\mathit{l}}_{\mathrm{2}}}{\mathit{C}}_{{\mathit{l}}_{\mathrm{3}}}}\mathit{,}\end{array}$(17)where ${\mathit{I}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}\mathrm{=}{\sum }_{\mathrm{all}\mathit{m}}{{}^{\mathrm{\right(}}{\mathit{G}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}^{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}^{\mathrm{\left)}}}^{\mathrm{2}}$.

From the inverse Fisher matrix (i.e. the covariance), one obtains the correlation between the sine and cosine terms. For example, for l = 1000 and k_c = 0.01, the correlation matrix is $\mathrm{corr}\mathrm{\left(}{\mathit{b}}^{\mathit{i}}\mathit{,}{\mathit{b}}^{\mathit{j}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{F}}_{\mathit{ij}}^{-1}}{\sqrt{{\mathit{F}}_{\mathit{ii}}^{-1}{\mathit{F}}_{\mathit{jj}}^{-1}}}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{1}& -0.04\\ -0.04& \mathrm{1}\end{array}\right)$(18)which shows a weak correlation as expected. From f₁ and f₂ of Eq. (6), we can calculate the amplitude and phase of the oscillation. The variance of the quantities f_NL and Φ can be calculated by error propagation from f₁ and f₂. For example ${\mathit{\sigma }}_{\mathit{f}}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{f}}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{f}}^{\mathrm{2}}}{\mathit{\sigma }}_{{\mathit{f}}_{\mathrm{1}}}^{\mathrm{2}}\mathrm{+}\frac{{\mathit{f}}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{f}}^{\mathrm{2}}}{\mathit{\sigma }}_{{\mathit{f}}_{\mathrm{2}}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\frac{{\mathit{f}}_{\mathrm{1}}{\mathit{f}}_{\mathrm{2}}}{{\mathit{f}}^{\mathrm{2}}}\mathit{cov}\mathrm{\left(}{\mathit{f}}_{\mathrm{1}}\mathit{,}{\mathit{f}}_{\mathrm{2}}\mathrm{\right)}\mathit{.}$(19)In the approximation of a diagonal covariance matrix and σ_f1 = σ_f2, this gives σ_f = σ_f1. In this approximation the sensitivity on the phase depends on f as ${\mathit{\sigma }}_{\mathrm{\Phi }}\mathrm{=}\frac{{\mathit{\sigma }}_{{\mathit{f}}_{\mathrm{1}}}}{\sqrt{\mathit{f}}}$.

The Fisher matrix allows us to forecast the precision that can be obtained on the bispectrum parameters f₁,f₂. With the approximation that the f₁,f₂ covariance matrix is a multiple of the unit matrix, the Fisher forecast on f₁,f₂ equals the forecast on f_NL. The precision on f_NL is then given by ${\mathit{\sigma }}_{{\mathit{f}}_{\mathrm{NL}}}\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathit{F}}}$. For a noiseless full-sky experiment, the Fisher forecast for the feature model is shown in Fig. 3 for a number of different oscillation frequencies.

Fig. 3 Fisher forecast of σfNL for the feature model for different kc, assuming a noiseless full sky experiment.

4.2. Estimating the frequency

The estimator described above explicitly estimates the amplitude f_NL and phase φ of the oscillation for a fixed frequency k_c. To estimate the oscillation frequency k_c, it is necessary to sample the frequency space with the KSW estimator and search for peaks in the significance of the estimated amplitude. We assume the primordial bispectrum to be given by a single oscillation frequency and not a spectrum of contributions. We consider only the sine component of the bispectrum first, meaning that the phase is φ = 0 (see below for the generalisation to the phase). In this case, the estimator for a frequency k_i is $\mathit{f̂}\mathit{i}\mathrm{=}\mathrm{(}{\mathit{F}}^{-1}{\mathrm{)}}_{\mathit{ii}}\mathit{S}i$(20)with covariance (in the usual Gaussian approximation) $\begin{array}{ccc}\mathrm{cov}\mathrm{\left(}\mathit{f̂}\mathit{i}\mathit{,}\mathit{f̂}\mathit{j}\mathrm{\right)}& \mathrm{=}& ⟨\mathit{f̂}\mathit{i}\mathit{f̂}\mathit{j}⟩\mathrm{-}⟨\mathit{f̂}\mathit{i}⟩⟨\mathit{f̂}\mathit{j}⟩\\ & \mathrm{=}& \frac{{\mathit{F}}_{\mathit{ij}}}{{\mathit{F}}_{\mathit{ii}}{\mathit{F}}_{\mathit{jj}}}\mathrm{·}\end{array}$The Fisher matrix is given by Eq. (13), where the index i now goes over frequency sampling points. An example of the Fisher matrix F_ij is shown in Fig. 4 (top). The corresponding correlation matrix is $\mathrm{corr}\mathrm{\left(}\mathit{f̂}{}_{\mathit{i}}\mathit{,}\mathit{f̂}{}_{\mathit{j}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{F}}_{\mathit{ij}}}{\sqrt{{\mathit{F}}_{\mathit{ii}}{\mathit{F}}_{\mathit{jj}}}}$, shown in Fig. 4 (middle).

Fig. 4 Top: fisher matrix Fij for 100 logarithmically spaced frequencies between kc = 0.002 and kc = 0.01. Middle: corresponding correlation matrix $\mathrm{corr}\mathrm{\left(}{\mathit{f}}_{\mathit{i}}^{\mathrm{\prime }}\mathit{,}{\mathit{f}}_{\mathit{j}}^{\mathrm{\prime }}\mathrm{\right)}$. Bottom: correlation matrix one-dimensional slice corr(fi,fj) for ki = 0.005.

A one-dimensional slice of the correlation matrix is shown in Fig. 4 (bottom) for k_c = 0.005. The plot shows strong anti-peaks to both sides of the maximum and several small secondary peaks. This is not an artefact of the chosen estimator, but the physical overlap of the CMB bispectra induced by different primordial oscillation frequencies. The frequency sampling must be at least sufficient to resolve the peak structure of this plot. However, the width of the central peak does not directly limit the precision σ_k with which the primordial frequency can be determined.

To evaluate the precision with which k_c can be estimated, we note that the correlation matrix is identical to the bispectrum correlator, $\mathit{C}\mathrm{\left(}\mathit{B,}{\mathit{B}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{𝒩}\mathrm{\left(}\mathit{B}\mathrm{\right)}\mathrm{𝒩}\mathrm{\left(}{\mathit{B}}^{\mathrm{\prime }}\mathrm{\right)}}\sum _{\mathit{l}}\frac{{\mathit{B}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}^{\mathit{i}}{\mathit{B}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}^{\mathit{j}}}{{\mathit{C}}_{{\mathit{l}}_{\mathrm{1}}}{\mathit{C}}_{{\mathit{l}}_{\mathrm{2}}}{\mathit{C}}_{{\mathit{l}}_{\mathrm{3}}}}\mathit{,}$(23)which is the usual measure to discriminate bispectra (see e.g. Fergusson & Shellard 2009). It gives the estimated proportion of f that is recovered when estimating a spectrum B′ when the true underlying spectrum is B. If the underlying bispectum is f_kcB_kc, the expectation value of the estimated amplitude at a different frequency k is thus ). The variance at each data point is independent of the signal. Thus, the correlation matrix approximates the estimates in the k range where the signal dominates the variance (compare Fig. 4 with the estimation example in Fig. 6).

Knowing the means f(k) and the covariance matrix Cov(k₁,k₂) and using the property that each estimator is Gaussian distributed, we can write the continuum likelihood for the estimated amplitudes , if the true bispectrum is f_kcB_kc: $\begin{array}{ccc}-2& & \ln \mathrm{ℒ}\mathrm{\left(}\mathit{f̂}\mathrm{\right(}\mathit{k}\mathrm{\left)}\mathrm{\right|}{\mathit{f}}_{{\mathit{k}}_{\mathrm{c}}}\mathit{,}{\mathit{k}}_{\mathrm{c}}\mathrm{)}\mathrm{=}\ln \mathrm{|}\mathrm{Cov}\mathrm{|}\\ & & \mathrm{+}\left[\mathrm{\int }\mathrm{d}{\mathit{k}}_{\mathrm{1}}\mathrm{d}{\mathit{k}}_{\mathrm{2}}\mathrm{\left(}\mathit{f̂}\mathrm{\right(}{\mathit{k}}_{\mathrm{1}}\mathrm{)}\mathrm{-}\mathit{\mu }\mathrm{(}{\mathit{k}}_{\mathrm{1}}\mathrm{\left)}\mathrm{\right)}{\mathrm{Cov}}^{-1}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}\mathrm{\left(}\mathit{f̂}\mathrm{\right(}{\mathit{k}}_{\mathrm{2}}\mathrm{)}\mathrm{-}\mathit{\mu }\mathrm{(}{\mathit{k}}_{\mathrm{2}}\mathrm{\left)}\mathrm{\right)}\right]\mathit{,}\end{array}$(24)with mean $\mathit{\mu }\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\mathit{f}}_{{\mathit{k}}_{\mathrm{c}}}\mathrm{Corr}\mathrm{\left(}{\mathit{B}}_{{\mathit{k}}_{\mathrm{c}}}\mathit{,}{\mathit{B}}_{\mathit{k}}\mathrm{\right)}\mathit{.}$(25)This likelihood could be explored by Monte Carlo to find maximum likelihood estimates of f_kc and k_c if a significant peak of f_NL is found in the spectrum. The Fisher matrix to forecast optimal precision on the estimated parameters is ${\mathit{F}}_{\mathit{\theta }{\mathit{\theta }}^{\mathrm{\prime }}}\mathrm{=}⟨\frac{\mathit{\partial }\mathrm{ℒ}}{\mathit{\partial \theta \partial }{\mathit{\theta }}^{\mathrm{\prime }}}⟩\mathrm{=}\mathrm{\int }\mathrm{d}{\mathit{k}}_{\mathrm{1}}\mathrm{d}{\mathit{k}}_{\mathrm{2}}\frac{\mathit{\partial \mu }\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}}{\mathit{\partial \theta }}{\mathrm{Cov}}^{-1}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}\frac{\mathit{\partial \mu }\mathrm{\left(}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}}{\mathit{\partial }{\mathit{\theta }}^{\mathrm{\prime }}}\mathit{,}$(26)where θ ∈ { f_kc,k_c }. This can be integrated numerically for any given fiducial parameters.

In the above discussion, we assumed that the bispectrum only has a sine contribution. The generalisation of Eq. (24) to a free phase is $\begin{array}{ccc}& & \mathrm{-}\mathrm{2}\ln \mathrm{ℒ}\mathrm{\left(}\mathit{f̂}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right(}\mathit{k}\mathrm{\left)}\mathrm{\right|}{\mathit{f}}_{\mathrm{1}}\mathit{,}{\mathit{f}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{c}}\mathrm{)}\mathrm{=}\ln \mathrm{|}\mathrm{Cov}\mathrm{|}\\ & & \mathrm{+}\begin{array}{c}\sum \\ \mathit{i}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{2}\mathit{j}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{2}\end{array}\mathrm{\int }\mathrm{d}{\mathit{k}}_{\mathrm{1}}\mathrm{d}{\mathit{k}}_{\mathrm{2}}\mathrm{\left(}\mathit{f̂}\mathit{i}\mathrm{\right(}{\mathit{k}}_{\mathrm{1}}\mathrm{)}\mathrm{-}{\mathit{\mu }}^{\mathit{i}}\mathrm{(}{\mathit{k}}_{\mathrm{1}}\mathrm{\left)}\mathrm{\right)}{\mathrm{Cov}}_{\mathit{ij}}^{-1}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}\mathrm{\left(}\mathit{f̂}\mathit{j}\mathrm{\right(}{\mathit{k}}_{\mathrm{2}}\mathrm{)}\mathrm{-}{\mathit{\mu }}^{\mathit{j}}\mathrm{(}{\mathit{k}}_{\mathrm{2}}\mathrm{\left)}\mathrm{\right)}\mathit{,}\end{array}$(27)with the mean ${\mathit{\mu }}^{\mathit{i}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\mathit{f}}_{{\mathit{k}}_{\mathrm{c}}}\mathrm{Corr}\mathrm{(}{\mathit{f}}_{\mathrm{1}}{\mathit{B}}_{{\mathit{k}}_{\mathrm{c}}}^{\mathrm{1}}\mathrm{+}{\mathit{f}}_{\mathrm{2}}{\mathit{B}}_{{\mathit{k}}_{\mathrm{c}}}^{\mathrm{2}}\mathit{,}{{\mathit{B}}_{\mathit{k}}^{\mathit{i}}}^{\mathrm{)}}\mathit{.}$(28)The correlation matrix can be split into terms of f₁ and f₂ for efficient evaluation of the likelihood.

As we have seen, the estimation of the frequency requires to run a large number of estimators (around 100 to cover the frequency interval k_c = 0.001 to k_c = 0.01). For each of these estimators, it is necessary to calculate the linear term in Eq. (16) via a Monte Carlo averaging procedure over Gaussian map realisations with the same mask and noise properties as present in the experiment. Convergence of the linear term is usually achieved with 100 Monte Carlo realisations, although several hundreds can be used to improve accuracy. After the optimisation of the conformal distance integral, that is reviewed in Sect. 6, one can expect, with an angular resolution of l_max = 2000, a calculation time of about one day on a single cpu for a single frequency, including the Monte Carlo averaging. On a computation grid, one can thus easily cover the frequency range of interest.

Fig. 5 Parameter estimation histograms for 100 maps with linear model non-Gaussianity with kc = 0.01 and lmax = 1000. The maps were created with fNL = 500 and φ = 0° (first two rows) and φ = 45° (rows three and four). Top: reconstructed f1 (green) and f2 (blue) amplitudes for φ = 0°. Second row: corresponding reconstructed amplitudes fNL and phases φ. Third and fourth row: same as above but with phase φ = 45°.

5. Map making

To verify the implementation and unbiased nature of the estimator, it is useful to be able to generate maps with the bispectrum signature of interest. The authors of Smith & Zaldarriaga (2011) introduced an algorithm that generates maps for arbitrary bispectra in the weak non-Gaussian limit. A map is constructed from a linear combination of a Gaussian and a non-Gaussian contribution as ${\mathit{a}}_{\mathit{lm}}\mathrm{=}{\mathit{a}}_{\mathit{lm}}^{\mathit{L}}\mathrm{+}{\mathit{f}}_{\mathrm{NL}}{\mathit{a}}_{\mathit{lm}}^{\mathrm{NL}}$. The straightforward implementation of the algorithm gives a non-Gaussian contribution of the form $\begin{array}{ccc}{\mathit{a}}_{\mathit{lm}}^{\mathrm{NL}}& \mathrm{=}& {\mathit{f}}_{\mathrm{1}}\mathrm{\int }{\mathit{r}}^{\mathrm{2}}\mathrm{d}\mathit{r}[\mathrm{-}{\mathit{X}}_{\mathit{l}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{\int }\mathrm{d\Omega }{\mathit{Y}}_{\mathit{lm}}^{\mathrm{\ast }}{\mathit{M}}_{\mathit{X}}^{\mathrm{2}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}\\ & & \mathrm{+}\left({\mathit{X}}_{\mathit{l}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{\int }\mathrm{d\Omega }{\mathit{Y}}_{\mathit{lm}}^{\mathrm{\ast }}{\mathit{M}}_{\mathit{Y}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}{\mathit{M}}_{\mathit{Y}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}\mathrm{+}\mathrm{2}\mathrm{perm}\mathit{.}\right)]\\ & & \mathrm{+}{\mathit{f}}_{\mathrm{2}}\mathrm{\int }{\mathit{r}}^{\mathrm{2}}\mathrm{d}\mathit{r}[{\mathit{Y}}_{\mathit{l}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{\int }\mathrm{d\Omega }{\mathit{Y}}_{\mathit{lm}}^{\mathrm{\ast }}{\mathit{M}}_{\mathit{Y}}^{\mathrm{2}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}\\ & & \mathrm{-}\left({\mathit{Y}}_{\mathit{l}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{\int }\mathrm{d\Omega }{\mathit{Y}}_{\mathit{lm}}^{\mathrm{\ast }}{\mathit{M}}_{\mathit{X}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}{\mathit{M}}_{\mathit{X}}\mathrm{\left(}\mathit{r,}\mathit{n̂}\mathrm{\right)}\mathrm{+}\mathrm{2}\mathrm{perm}\mathit{.}\right)]\mathit{.}\end{array}$(29)It is thus easy to create maps of arbitrary phase from the sine and cosine terms.

An example of two sets of 100 maps, created and then estimated with the algorithms presented here, is shown in Fig. 5. The means and variances in these histograms are compatible with their Fisher forecast. We note that the distribution of f_NL is not Gaussian but follows a Rayleigh distribution, since f_NL represents the magnitude of the vector of the two directional components f₁ and f₂.

An example of a frequency sweep can be found in Fig. 6 for the phase ϕ = 0°. It shows the secondary peaks that are expected from the correlation of different frequencies.

Fig. 6 Left: frequency sweep over a map with kc = 0.005 and fNL = 2000. Right: as left but in units of σfNL.

6. Computational speed improvement by optimisation of the r-integral

Due to the large number of frequencies that have to be sampled with the corresponding estimator, the search for oscillations is computationally challenging. The time critical steps are the calculation of the Fisher matrix and the necessity of calculating a large number of r dependent KSW filtered maps. The latter problem becomes even more severe if one has to estimate many Monte Carlo generated maps for the calculation of the linear term. The situation can be improved by an analysis of the Fisher matrix, as shown in Smith & Zaldarriaga (2011).

For a separable shape, the Fisher matrix can also be expressed by a sum over contributions of different r sampling points that arise when numerically evaluating the bispectrum integral in Eq. (6). The total Fisher matrix is then given as a sum $\mathit{F}\mathrm{=}{\sum }_{\mathit{i,j}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{fact}}}{\mathit{F}}_{\mathit{ij}}$, where F_ij is the Fisher matrix element between the sampling points i and j and N_fact is the number of sampling points. The Fisher matrix elements are then given by $\begin{array}{ccc}{\mathit{F}}_{\mathit{ij}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{6}}\sum _{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}\frac{\mathrm{(}\mathrm{2}{\mathit{l}}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{l}}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{l}}_{\mathrm{3}}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{4}\mathit{\pi }}\\ & & \mathrm{\times }{\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)}^{\mathrm{2}}\frac{{\mathit{b}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}^{\mathit{i}}{\mathit{b}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}^{\mathit{j}}}{{\mathit{C}}_{{\mathit{l}}_{\mathrm{1}}}{\mathit{C}}_{{\mathit{l}}_{\mathrm{2}}}{\mathit{C}}_{{\mathit{l}}_{\mathrm{3}}}}\mathrm{·}\end{array}$(30)We now explicitly consider the sine term of the linear model (the cosine term is analogous), where the contribution of a distance r_i is given by $\begin{array}{ccc}{\mathit{b}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}^{\mathit{i}}& \mathrm{=}& \mathrm{\left(}\mathrm{\Delta }{\mathit{r}}_{\mathit{i}}\mathrm{\right)}{\mathit{r}}_{\mathit{i}}^{\mathrm{2}}[\mathrm{-}{\mathit{X}}_{{\mathit{l}}_{\mathrm{1}}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathrm{\right)}{\mathit{X}}_{{\mathit{l}}_{\mathrm{2}}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathrm{\right)}{\mathit{X}}_{{\mathit{l}}_{\mathrm{3}}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathrm{\right)}\\ & & \mathrm{+}\left[{\mathit{X}}_{{\mathit{l}}_{\mathrm{1}}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathrm{\right)}{\mathit{Y}}_{{\mathit{l}}_{\mathrm{2}}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathrm{\right)}{\mathit{Y}}_{{\mathit{l}}_{\mathrm{3}}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathrm{\right)}\mathrm{+}\mathrm{2}\mathrm{perm}\mathrm{.}\right]]\mathit{.}\end{array}$(31)To give an impression of the bispectrum contribution of different distances r_i, we plot the diagonal elements F_ii in Fig. 7. The plots show the contribution of recombination (r = 14 000 Mpc), reionisation (r ≃ 10 500 Mpc), and ISW (r> 5000 Mpc). As expected, the dominant contribution comes from the time around recombination, and one can get a good approximation to the integral by sampling only a window around recombination. This is particularly useful to quickly scan a parameter space, in the present case the oscillation frequency, k_c.

In Smith & Zaldarriaga (2011), it was shown that one can go further and optimise the r sampling points to find a quadrature with surprisingly few sampling points that give an almost identical estimator. Their algorithm constructs a new bispectrum B′ from the original bispectrum B by choosing a subsample of points and weighting them so that the Fisher distance between the two is minimised: $\mathit{F}\mathrm{\left(}\mathit{B,}{\mathit{B}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{6}}\sum _{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}\frac{\mathrm{(}{\mathit{B}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}\mathrm{-}{\mathit{B}}_{{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{3}}}^{\mathrm{\prime }}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{C}}_{{\mathit{l}}_{\mathrm{1}}}{\mathit{C}}_{{\mathit{l}}_{\mathrm{2}}}{\mathit{C}}_{{\mathit{l}}_{\mathrm{3}}}}\mathrm{·}$(32)This means that bispectrum values with a small signal-to-noise are allowed to be very different. Using this algorithm, as an example, we obtain an approximate bispectrum B′ consisting of 30 sample points leading to a separability between B and B′ of 0.1σ assuming k_c = 0.01 and f_NL = 1000.

The most computationally demanding task in the estimation pipeline remains the calculation of the Fisher matrix in Eq. (30). It was also shown in Smith & Zaldarriaga (2011) that one can factorise this equation by inserting the integral representation of the Wigner symbol ${\left(\begin{array}{c}\\ {\mathit{l}}_{\mathrm{1}}& {\mathit{l}}_{\mathrm{2}}& {\mathit{l}}_{\mathrm{3}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{1}}\mathrm{d}\mathit{z}{\mathit{P}}_{{\mathit{l}}_{\mathrm{1}}}\mathrm{\left(}\mathit{z}\mathrm{\right)}{\mathit{P}}_{{\mathit{l}}_{\mathrm{2}}}\mathrm{\left(}\mathit{z}\mathrm{\right)}{\mathit{P}}_{{\mathit{l}}_{\mathrm{3}}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{.}$(33)The integral over z can be computed efficiently by Gauss Legendre integration. However, even with this expression, the calculation of F_ij needs many CPU hours depending on the chosen initial quadrature point number. If one does not want to calculate an optimised quadrature at every frequency points, but only wants to know the normalisation F(k_c) of the estimators, one can calculate this normalisation on a much wider frequency spacing and interpolate in between. This can be seen from the frequency dependent Fisher matrix in Fig. 4, where the diagonal elements vary slowly.

Fig. 7 Fii as a function of conformal distance ri. Top: kc = 0.01, bottom: kc = 0.005.

7. A position space interpretation of the KSW estimator for oscillations

In Appendix A, we show that the three-point function of the feature model in position space peaks for configurations that lie on a circumcircle of radius sin(θ) = 1 / (2k₀η), where $\mathrm{2}{\mathit{k}}_{\mathrm{0}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}\mathit{\pi }}{\mathit{k}}_{\mathrm{c}}$ in the convention of Eq. (2). This suggests searching for bispectrum oscillations in real space by convoluting the CMB map with a ring kernel of varying radius.

For a radially symmetric kernel, the convolution can be done efficiently in harmonic space as s_lm = K_lr_lm, where the kernel is given by the Legendre transformation, ${\mathit{K}}_{\mathit{l}}\mathrm{=}\mathrm{2}\mathit{\pi }{\mathrm{\int }}_{-1}^{\mathrm{1}}\mathit{K}\mathrm{\left(}\mathit{z}\mathrm{\right)}{\mathit{P}}_{\mathit{l}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{d}\mathit{z,}$(34)and K(z) is a narrow window function in z = cos(θ). The estimate ℰ^ring(z) is given by the sum over the pixels of the cube of the convoluted map, ${\mathrm{ℰ}}^{\mathrm{ring}}\mathrm{\left(}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{\int }\mathrm{d\Omega }{\left[\sum _{\mathit{lm}}{\mathit{K}}_{\mathit{l}}{\mathit{a}}_{\mathit{lm}}{\mathit{Y}}_{\mathit{lm}}\right]}^{\mathrm{3}}\mathit{.}$(35)Figure 8 shows an example of this estimator for a map that was simulated with the algorithm of the preceding section. The plots show ${\mathrm{ℰ}}_{\mathrm{NG}}^{\mathrm{ring}}\mathrm{-}{\mathrm{ℰ}}_{\mathrm{G}}^{\mathrm{ring}}$, which means that the estimate from the Gaussian map was subtracted. The red line shows the predicted position of the maximum. We note that this maximum would be harder to locate in real data where the Gaussian contribution cannot simply be subtracted.

Fig. 8 Ring estimate ${\mathrm{ℰ}}_{\mathrm{NG}}^{\mathrm{ring}}\mathrm{-}{\mathrm{ℰ}}_{\mathrm{G}}^{\mathrm{ring}}$. Top left: kc = 0.01. Top right: kc = 0.005. Bottom: kc = 0.001. The kernel width in all plots was Δθ = 0.005, fNL was 10−20 times the optimal Fisher error. The red line shows the expected position of the maximum.

It is interesting to compare the “intuitive” ring kernel K(θ) with the KSW kernel that is known to give optimal results. The KSW estimator is of form, ${\mathrm{ℰ}}^{\mathrm{KSW}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{\propto }\mathrm{\int }{\mathit{r}}^{\mathrm{2}}\mathrm{d}\mathit{r}\mathrm{\int }\mathrm{d\Omega }{\left[\sum _{\mathit{lm}}\frac{{\mathit{X}}_{\mathit{l}}{\mathit{a}}_{\mathit{lm}}}{{\mathit{C}}_{\mathit{l}}}{\mathit{Y}}_{\mathit{lm}}\right]}^{\mathrm{3}}\mathrm{+}\mathit{....}$(36)The largest contribution to this integral comes from decoupling at r_rec. We plot the KSW kernel function at decoupling X_rec(θ) in Fig. 9. It shows the maximum at the expected position and has roughly the expected window shape.

Fig. 9 KSW kernel X(θ) at decoupling radius rrec. Top left: kc = 0.01. Top right: kc = 0.005. Bottom: kc = 0.001. The red line is the predicted maximum of the kernel.

8. Conclusion

In this paper, we have presented and extensively studied the KSW estimator for linear bispectrum oscillations. The main motivation for this approach is that the oscillating bispectrum shapes are difficult to represent with a modal expansion and, thus, have not yet been constrained at high oscillation frequency. We have provided the equations for estimation and map making and validated them with Monte Carlo simulations. Unlike many of the well know bispectum shapes, oscillations have two free parameters in addition to the common amplitude parameter f_NL. We have developed the methodology to estimate and constrain the oscillation phase φ and the frequency k_c. Our work will therefore allow one to explore a parameter space that was not previously accessible for a theoretically well-motivated bispectrum shape. Finally, we have found an interesting position space interpretation of the KSW estimator for oscillations, based on an approximate evaluation of the corresponding three-point function.

Acknowledgments

The authors acknowledge support from NSF Grant NSF AST 09-08693 ARRA. B.D.W. acknowledges funding from an ANR Chaire d’Excellence (ANR-10-CEXC-004-01), the UPMC Chaire Internationale in Theoretical Cosmology, and NSF grants AST-0908 902 and AST-0708849. This work made in the ILP LABEX (under reference ANR-10-LABX-63) was supported by French state funds managed by the ANR within the Investissements d’Avenir programme under reference ANR-11-IDEX-0004-02. M.M. acknowledges funding by Centre National d’Etudes Spatiales (CNES). The authors would like to thank Daan Meerburg for useful discussions and comments.

References

Appendix A: Angular correlation function of the feature model

In this appendix, we show that the linear feature model bispectrum peaks in real space for a special class of three-point function configurations. A similar analysis was presented in Adshead et al. (2012). The corresponding calculation for logarithmic oscillations can be found in Jackson et al. (2014).

In the approximation of instantaneous CMB decoupling at time η, the CMB temperature perturbation is given in terms of the potential φ as Bashinsky & Bertschinger (2001), $\frac{\mathrm{\Delta }\mathit{T}\mathrm{\left(}\hat{{n}}\mathrm{\right)}}{\mathit{T}}\mathrm{=}\mathrm{\int }\frac{{\mathrm{d}}^{\mathrm{3}}{k}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}{\mathit{\varphi }}_{{k}}\mathrm{\left(}{\mathrm{0}}^{\mathrm{-}}\mathrm{\right)}{\mathit{T}}_{\mathrm{rad}}\mathrm{\left(}\mathit{k}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\eta }{k}\mathrm{·}\hat{{n}}}\mathit{.}$(A.1)The transfer function T_rad(k) is generally a complicated function of scale and cosmological parameters. For simplicity of an analytic answer, which accounts for finite resolution, we take it to be ${\mathit{T}}_{\mathrm{rad}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{\approx }{\mathrm{e}}^{\mathrm{-}{\mathit{k}}^{\mathrm{2}}\mathit{/}{\mathit{k}}_{\mathrm{D}}^{\mathrm{2}}}$.

The position-space primordial bispectrum is given by $\begin{array}{ccc}⟨\frac{\mathrm{\Delta }\mathit{T}\mathrm{\left(}{\hat{{n}}}_{\mathrm{1}}\mathrm{\right)}}{\mathit{T}}& & \frac{\mathrm{\Delta }\mathit{T}\mathrm{\left(}{\hat{{n}}}_{\mathrm{2}}\mathrm{\right)}}{\mathit{T}}\frac{\mathrm{\Delta }\mathit{T}\mathrm{\left(}{\hat{{n}}}_{\mathrm{3}}\mathrm{\right)}}{\mathit{T}}⟩\mathrm{=}\\ & & \mathrm{\int }\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}\frac{{\mathrm{d}}^{\mathrm{3}}{{k}}_{\mathit{i}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}{\mathit{T}}_{\mathrm{rad}}\mathrm{\left(}{\mathit{k}}_{\mathit{i}}\mathrm{\right)}\exp \left(\mathrm{-}\mathrm{i}\mathit{\eta }{{k}}_{\mathit{i}}\mathrm{·}{\hat{{n}}}_{\mathit{i}}\right)\mathrm{⟨}{\mathit{\varphi }}_{{{k}}_{\mathrm{1}}}{\mathit{\varphi }}_{{{k}}_{\mathrm{2}}}{\mathit{\varphi }}_{{{k}}_{\mathrm{3}}}\mathrm{⟩}\mathit{,}\end{array}$(A.2)where the correlation is of the form,

$\mathrm{⟨}{\mathit{\varphi }}_{{{k}}_{\mathrm{1}}}{\mathit{\varphi }}_{{{k}}_{\mathrm{2}}}{\mathit{\varphi }}_{{{k}}_{\mathrm{3}}}\mathrm{⟩}\mathrm{\equiv }{\mathit{B}}_{\mathit{\varphi }}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathit{,}{{k}}_{\mathrm{2}}\mathit{,}{{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}{\mathit{\delta }}^{\mathrm{3}}\left(\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{3}}{{k}}_{\mathit{i}}\right)\mathit{.}$Since k₁ + k₂ + k₃ = 0, the k-space correlations are categorised by the triangle formed by the k_i.

We now examine this integral for the linear model,

${\mathit{B}}_{\mathit{\varphi }}\mathrm{=}\frac{{\mathit{B}}_{\mathrm{0}}}{\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{2}}}\sin \left(\frac{{\mathit{k}}_{\mathrm{1}}\mathrm{+}{\mathit{k}}_{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{3}}}{\mathrm{2}{\mathit{k}}_{\mathrm{0}}}\right)\mathrm{·}$The delta function that couples the momenta can be written as

$\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}{\mathit{\delta }}^{\mathrm{3}}\left({{k}}_{\mathrm{1}}\mathrm{+}{{k}}_{\mathrm{2}}\mathrm{+}{{k}}_{\mathrm{3}}\right)\mathrm{=}{\mathit{\eta }}^{\mathrm{3}}\mathrm{\int }{\mathrm{d}}^{\mathrm{3}}{w}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\eta }{w}\mathrm{·}\mathrm{(}{{k}}_{\mathrm{1}}\mathrm{+}{{k}}_{\mathrm{2}}\mathrm{+}{{k}}_{\mathrm{3}}\mathrm{)}}\mathit{,}$where the factor of η has been included for future convenience. Writing the sine as exponentials and performing the integral over angles, we obtain $\begin{array}{ccc}\left(\frac{\mathrm{\Delta }\mathit{T}}{\mathit{T}}\right)\begin{array}{c}\mathrm{3}\\ \mathrm{osc}\end{array}& \mathrm{=}& \frac{{\mathit{B}}_{\mathrm{0}}{\mathit{\eta }}^{\mathrm{3}}}{\mathrm{2}\mathrm{i}}\mathrm{\int }{\mathrm{d}}^{\mathrm{3}}{w}\sum _{\mathrm{\pm }}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}\mathrm{\int }\frac{{\mathrm{d}}^{\mathrm{3}}{{k}}_{\mathit{i}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}{\mathit{k}}_{\mathit{i}}^{\mathrm{2}}}\\ & & \mathrm{\times }\mathrm{\pm }{\mathrm{e}}^{\mathrm{\pm }\mathrm{i}{\mathit{k}}_{\mathit{i}}\mathit{/}\mathrm{2}{\mathit{k}}_{\mathrm{0}}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\eta }{{k}}_{\mathit{i}}\mathrm{·}\mathrm{(}{w}\mathrm{+}\hat{{n}}\mathit{i}\mathrm{)}}{\mathrm{e}}^{\mathrm{-}{\mathit{k}}_{\mathit{i}}^{\mathrm{2}}\mathit{/}{\mathit{k}}_{\mathrm{D}}^{\mathrm{2}}}\\ & & \mathrm{=}\frac{{\mathit{B}}_{\mathrm{0}}{\mathit{\eta }}^{\mathrm{3}}}{\mathrm{2}\mathrm{i}}\mathrm{\int }{\mathrm{d}}^{\mathrm{3}}{w}\sum _{\mathrm{\pm }}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}\frac{\mathrm{\pm }\mathrm{1}}{\mathrm{\left(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}\mathrm{i}\mathit{\eta }\mathrm{\right|}{w}\mathrm{+}\hat{{n}}\mathit{i}\mathrm{|}}\\ & & \mathrm{\times }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{d}{\mathit{k}}_{\mathit{i}}}{{\mathit{k}}_{\mathit{i}}}{\mathrm{e}}^{\mathrm{\pm }\mathrm{i}{\mathit{k}}_{\mathit{i}}\mathit{/}\mathrm{2}{\mathit{k}}_{\mathrm{0}}}\mathrm{(}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{\mathit{k}}_{\mathrm{i}}\mathit{\eta }\mathrm{|}{w}\mathrm{+}\hat{{n}}\mathit{i}\mathrm{|}}\mathrm{-}{\mathrm{e}}^{\mathrm{i}{\mathit{k}}_{\mathrm{i}}\mathit{\eta }\mathrm{|}{w}\mathrm{+}\hat{{n}}\mathit{i}\mathrm{|}}\mathrm{)}{\mathrm{e}}^{\mathrm{-}{\mathit{k}}_{\mathit{i}}^{\mathrm{2}}\mathit{/}{\mathit{k}}_{\mathrm{D}}^{\mathrm{2}}}\mathit{.}\end{array}$Defining the dimensionless parameter x_i ≡ kη, the momentum integral is $\begin{array}{ccc}\mathrm{2}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{d}{\mathit{x}}_{\mathit{i}}}{{\mathit{x}}_{\mathit{i}}}& & (\cos [\mathrm{(}\mathit{\eta }{\mathit{k}}_{\mathrm{0}}{\mathrm{)}}^{-1}\mathrm{+}\mathrm{|}{w}\mathrm{+}\hat{{n}}\mathit{i}\mathrm{|}]{\mathit{x}}_{\mathit{i}}\\ & & \mathrm{-}\cos \mathrm{[}\mathrm{(}\mathit{\eta }{\mathit{k}}_{\mathrm{0}}{\mathrm{)}}^{-1}\mathrm{-}\mathrm{|}{w}\mathrm{+}\hat{{n}}\mathit{i}{\mathrm{|}}^{\mathrm{]}}{\mathit{x}}_{\mathit{i}}){\mathrm{e}}^{\mathrm{-}\mathrm{2}{\mathit{x}}_{\mathit{i}}^{\mathrm{2}}\mathit{/}{\mathit{\eta }}^{\mathrm{2}}{\mathit{k}}_{\mathrm{D}}^{\mathrm{2}}}\mathit{.}\end{array}$This integral can be evaluated exactly in terms of a hypergeometric function, but there is a simplifying limit we can take. The low-momentum, long-distance approximation allows

$\cos \mathrm{\left[}\mathrm{\right(}{\mathit{k}}_{\mathrm{0}}\mathit{\eta }{\mathrm{)}}^{-1}\mathrm{\pm }\mathrm{|}{w}\mathrm{+}\hat{{n}}\mathit{i}{\mathrm{\left|}}^{\mathrm{\right]}}\mathit{x}\mathrm{\approx }{\mathrm{e}}^{\mathrm{-}{\left[\mathrm{(}{\mathit{k}}_{\mathrm{0}}\mathit{\eta }{\mathrm{)}}^{-1}\mathrm{\pm }\mathrm{|}{w}\mathrm{+}\hat{{n}}\mathit{i}\mathrm{|}\right]}^{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}\mathit{/}\mathrm{2}}\mathit{.}$The integral then has the simple analytic solution $\begin{array}{ccc}\left(\frac{\mathrm{\Delta }\mathit{T}}{\mathit{T}}\right)\begin{array}{c}\mathrm{3}\\ \mathrm{osc}\end{array}\mathrm{\approx }& & \frac{{\mathit{B}}_{\mathrm{0}}{\mathit{\eta }}^{\mathrm{3}}}{\mathrm{2}\mathrm{i}}\mathrm{\int }{\mathrm{d}}^{\mathrm{3}}{w}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}\frac{\mathrm{1}}{\mathrm{\left(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}\mathit{i\eta }\mathrm{\right|}{w}\mathrm{+}\hat{{n}}\mathit{i}\mathrm{|}}\\ & & \mathrm{\times }\ln \left(\frac{{\mathrm{[}\mathrm{(}\mathrm{2}{\mathit{k}}_{\mathrm{0}}\mathit{\eta }{\mathrm{)}}^{-1}\mathrm{+}\mathrm{|}{w}\mathrm{+}\hat{{n}}\mathit{i}{\mathrm{|}}^{\mathrm{]}}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathrm{(}\mathit{\eta }{\mathit{k}}_{\mathrm{D}}{\mathrm{)}}^{-2}}{{\left[\mathrm{(}\mathrm{2}{\mathit{k}}_{\mathrm{0}}\mathit{\eta }{\mathrm{)}}^{-1}\mathrm{-}\mathrm{|}{w}\mathrm{+}\hat{{n}}\mathit{i}\mathrm{|}\right]}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathrm{(}\mathit{\eta }{\mathit{k}}_{\mathrm{D}}{\mathrm{)}}^{-2}}\right)\mathrm{·}\end{array}$In the ηk_D ≫ 1 limit, this maximally peaks when all three products peak near

$\mathrm{(}\mathrm{2}\mathit{\eta }{\mathit{k}}_{\mathrm{0}}{\mathrm{)}}^{-1}\mathrm{-}\mathrm{|}{{w}}_{\mathrm{0}}\mathrm{+}\hat{{n}}\mathit{i}\mathrm{|}\mathrm{=}\mathrm{0}\mathit{,i}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3.}$Squaring then subtracting, we obtain

${{w}}_{\mathrm{0}}\mathrm{·}{{N}}_{\mathit{ij}}\mathrm{=}\mathrm{0}\mathit{,}{{N}}_{\mathit{ij}}\mathrm{\equiv }\hat{{n}}\mathit{i}\mathrm{-}\hat{{n}}\mathit{j}\mathit{,i}\ne \mathit{j}\mathit{.}$We can take the three vectors N₁₂, N₂₃, and N₃₁ and arrange them in the $\hat{{x}}\mathrm{-}\hat{{y}}$ plane, so that $\begin{array}{ccc}\hat{{n}}\mathit{i}& \mathrm{=}& \sin \mathit{\theta }\cos {\mathrm{\Theta }}_{\mathit{i}}\hat{{x}}\mathrm{+}\sin \mathit{\theta }\sin {\mathrm{\Theta }}_{\mathit{i}}\hat{{y}}\mathrm{+}\cos \mathit{\theta }\hat{{z}}\mathit{,}\\ {{w}}_{\mathrm{0}}& \mathrm{=}& \mathrm{(}\mathrm{\pm }\mathit{\rho }\mathrm{-}\cos \mathit{\theta }\mathrm{)}\hat{{z}}\mathit{,}\end{array}$where

$\mathit{\rho }\mathrm{=}\sqrt{\mathrm{(}\mathrm{2}\mathit{\eta }{\mathit{k}}_{\mathrm{0}}{\mathrm{)}}^{-2}\mathrm{-}{\sin }^{\mathrm{2}}\mathit{\theta }}\mathit{.}$We now have the value of w for which the three factors are in resonance for a given configuration (n₁,n₂,n₃). The value of the factor at the resonance point is largest when $\mathrm{|}{w}\mathrm{+}\hat{{n}}{}_{\mathit{i}}\mathrm{|}$ is at its minimum, which is the case for ρ = 0. We conclude that the three-point function peaks for triangle configurations (n₁,n₂,n₃) that lie on a circle with radius given by $\sin \mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}{\mathit{k}}_{\mathrm{0}}\mathit{\eta }}$.

All Figures

	Fig. 1 Function X and Y of the feature model for kc = 0.01 (top) and kc = 0.001 (bottom) as a function of l for r = τ0 − τrec. X and Y have units of Mpc-1.
In the text

	Fig. 2 Equal l reduced bispectrum normalised by the large angle solution of the constant bispectrum for feature models with different parameters.
In the text

	Fig. 3 Fisher forecast of σfNL for the feature model for different kc, assuming a noiseless full sky experiment.
In the text

	Fig. 4 Top: fisher matrix Fij for 100 logarithmically spaced frequencies between kc = 0.002 and kc = 0.01. Middle: corresponding correlation matrix $\mathrm{corr}\mathrm{\left(}{\mathit{f}}_{\mathit{i}}^{\mathrm{\prime }}\mathit{,}{\mathit{f}}_{\mathit{j}}^{\mathrm{\prime }}\mathrm{\right)}$. Bottom: correlation matrix one-dimensional slice corr(fi,fj) for ki = 0.005.
In the text

Fig. 5 Parameter estimation histograms for 100 maps with linear model non-Gaussianity with kc = 0.01 and lmax = 1000. The maps were created with fNL = 500 and φ = 0° (first two rows) and φ = 45° (rows three and four). Top: reconstructed f1 (green) and f2 (blue) amplitudes for φ = 0°. Second row: corresponding reconstructed amplitudes fNL and phases φ. Third and fourth row: same as above but with phase φ = 45°.

In the text

	Fig. 6 Left: frequency sweep over a map with kc = 0.005 and fNL = 2000. Right: as left but in units of σfNL.
In the text

	Fig. 7 Fii as a function of conformal distance ri. Top: kc = 0.01, bottom: kc = 0.005.
In the text

	Fig. 8 Ring estimate ${\mathrm{ℰ}}_{\mathrm{NG}}^{\mathrm{ring}}\mathrm{-}{\mathrm{ℰ}}_{\mathrm{G}}^{\mathrm{ring}}$. Top left: kc = 0.01. Top right: kc = 0.005. Bottom: kc = 0.001. The kernel width in all plots was Δθ = 0.005, fNL was 10−20 times the optimal Fisher error. The red line shows the expected position of the maximum.
In the text

	Fig. 9 KSW kernel X(θ) at decoupling radius rrec. Top left: kc = 0.01. Top right: kc = 0.005. Bottom: kc = 0.001. The red line is the predicted maximum of the kernel.
In the text