The Full-Sky Angular Bispectrum in Redshift Space

TL;DR

This paper presents a comprehensive calculation of the full-sky angular bispectrum of galaxy counts in redshift space, incorporating nonlinear bias and redshift space distortions, with implications for future large-scale structure analyses.

Contribution

It introduces a novel formalism for computing the full-sky angular bispectrum including RSD and nonlinear bias, extending previous flat-sky and simplified models.

Findings

Nonlinear RSD effects are comparable to matter density perturbations in narrow redshift bins.

The formalism accurately models the bispectrum without flat-sky approximation.

Differences between galaxy number counts and CMB bispectra are highlighted.

Abstract

We compute the redshift-dependent angular bispectrum of galaxy number counts at tree-level, including nonlinear clustering bias and estimating numerically for the first time the effect of redshift space distortions (RSD). We show that for narrow redshift bins the amplitude of nonlinear RSD is comparable with the matter density perturbations. While our numerical results only include terms relevant on sub-horizon scales, the formalism can readily be extended to the full tree-level bispectrum. Our approach does not rely on the flat-sky approximation and it can be easily generalized to different sources by including the appropriate bias expansion. We test the accuracy of Limber approximation for different z-bins. We highlight the subtle but relevant differences in the angular bispectrum of galaxy number counts with respect to CMB, due to the different scale dependence of perturbations. Our formalism can also be directly applied to the angular HI intensity mapping bispectrum.

Figures

Equations

$$\Delta \left(\mathbf{n},z\right)=\frac{N\left(\mathbf{n},z\right)-⟨N⟩\left(z\right)}{⟨N⟩\left(z\right)}$$

$$d{s}^2=a^2\left[-\left(1+2\Psi \right)d{\eta }^2+\left(1-2\Phi \right)d{\mathbf{x}}^2\right],$$

$${\Delta }^{(1)}=b_1{\delta }^{(1)}+{\mathcal{H}}^{-1}{\partial }_r^2{v}^{(1)},$$

$${\Delta }^{(2)}(\mathbf{n},z)$$

$${\delta }_g=b_1\delta +\frac{1}{2}{b}_2{\delta }^2+b_s{s}^2.$$

$$s_{ij}=\frac{2}{3{\Omega }_M{\mathcal{H}}^2}{\partial }_i{\partial }_j\Phi -\frac{1}{3}{\delta }_{ij}\delta .$$

$$b_s=-\frac{2}{7}(b_1-1).$$

$$b_1\left(z\right)$$

$$b_2\left(z\right)$$

$$B({\mathbf{n}}_1,{\mathbf{n}}_2,{\mathbf{n}}_3,z_1,z_2,z_3)$$

$$B_{\mathrm{tree}}({\mathbf{n}}_1,{\mathbf{n}}_2,{\mathbf{n}}_3,z_1,z_2,z_3)$$

$$\Delta (\mathbf{n},z)=\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }{a}_{\ell m}(z)Y_{\ell m}(\mathbf{n}),a_{\ell m}(z)=\int d{\Omega }_{\mathbf{n}}\Delta (\mathbf{n},z)Y_{\ell m}^{\ast }(\mathbf{n})$$

$$B({\mathbf{n}}_1,{\mathbf{n}}_2,{\mathbf{n}}_3,z_1,z_2,z_3)=\sum _{\begin{array}{c}{\ell }_1,{\ell }_2,{\ell }_3\\ m_1,m_2,m_3\end{array}}{B}_{{\ell }_1{\ell }_2{\ell }_3}^{m_1{m}_2{m}_3}(z_1,z_2,z_3)Y_{{\ell }_1{m}_1}({\mathbf{n}}_1)Y_{{\ell }_2{m}_2}({\mathbf{n}}_2)Y_{{\ell }_3{m}_3}({\mathbf{n}}_3),$$

$$B_{{\ell }_1{\ell }_2{\ell }_3}^{m_1{m}_2{m}_3}(z_1,z_2,z_3)$$

$$B_{{\ell }_1{\ell }_2{\ell }_3}^{m_1{m}_2{m}_3}(z_1,z_2,z_3)={\mathcal{G}}_{{\ell }_1,{\ell }_2,{\ell }_3}^{m_1,m_2,m_3}{b}_{{\ell }_1{\ell }_2{\ell }_3}(z_1,z_2,z_3),$$

$${\mathcal{G}}_{{\ell }_1,{\ell }_2,{\ell }_3}^{m_1,m_2,m_3}$$

$$b_1(z_2)b_1(z_3)⟨{\mathcal{H}}^{-1}{\partial }_r^2{v}^{(2)}\left({\mathbf{n}}_1,z_1\right){\delta }^{(1)}\left({\mathbf{n}}_2,z_2\right){\delta }^{(1)}\left({\mathbf{n}}_3,z_3\right)⟩.$$

$$v^{(2)}\left({\mathbf{k}}_1,z\right)=-\frac{\mathcal{H}(z)f(z)}{{\left(2\pi \right)}^3{k}_1^2}\int d^3{k}_2{d}^3{k}_3{\delta }_D\left({\mathbf{k}}_1-{\mathbf{k}}_2-{\mathbf{k}}_3\right)G_2\left({\mathbf{k}}_2,{\mathbf{k}}_3\right)\delta \left({\mathbf{k}}_2,z\right)\delta \left({\mathbf{k}}_3,z\right),$$

$$G_2\left({\mathbf{k}}_2,{\mathbf{k}}_3\right)=\frac{3}{7}+\frac{1}{2}\frac{{\mathbf{k}}_2\cdot {\mathbf{k}}_3}{k_2{k}_3}\left(\frac{k_2}{k_3}+\frac{k_3}{k_2}\right)+\frac{4}{7}{\left(\frac{{\mathbf{k}}_2\cdot {\mathbf{k}}_3}{k_2{k}_3}\right)}^2.$$

$$G_2\left(k_1,k_2,k_3\right)$$

$${\delta }_D\left({\mathbf{k}}_1+{\mathbf{k}}_2+{\mathbf{k}}_3\right)$$

$$b_{{\ell }_1{\ell }_2{\ell }_3}^{v^{(2{)}^{\prime }}}\left(z_1,z_2,z_3\right)$$

$$j_{{\ell }_1}^{\prime \prime }\left(k_1{r}_1\right)=\int dx{\delta }_D(x-r_1)j_{{\ell }_1}^{\prime \prime }\left(k_{1}x\right)=\int dx{k}_1^{-2}{\delta }_D^{\prime \prime }(x-r_1)j_{{\ell }_1}\left(k_{1}x\right).$$

$$\left(\frac{{\partial }^2}{\partial {\chi }^2}+\frac{2}{\chi }\frac{\partial }{\partial \chi }-\frac{\ell \left(\ell +1\right)}{{\chi }^2}\right)j_{\ell }\left(k\chi \right)=-k^2{j}_{\ell }\left(k\chi \right),$$

$$\int d{k}_1{j}_{{\ell }_1}\left(k_{1}x\right)j_{{\ell }_1}\left(k_1\chi \right)=\frac{\pi }{2\left(1+2{\ell }_1\right)}\left({\chi }^{-1-{\ell }_1}{x}^{{\ell }_1}{\Theta }_H\left(\chi -x\right)+x^{-1-{\ell }_1}{\chi }^{{\ell }_1}{\Theta }_H\left(x-\chi \right)\right),$$

$$\int d{k}_1{k}_1^2{j}_{{\ell }_1}\left(k_{1}x\right)j_{{\ell }_1}\left(k_1\chi \right)=\frac{\pi }{2x^2}{\delta }_D\left(\chi -x\right),$$

$$\int d{k}_1{k}_1^4{j}_{{\ell }_1}\left(k_{1}x\right)j_{{\ell }_1}\left(k_1\chi \right)=\frac{\pi }{2x^2}\left[-\frac{{\partial }^2}{\partial {\chi }^2}-\frac{2}{\chi }\frac{\partial }{\partial \chi }+\frac{{\ell }_1\left({\ell }_1+1\right)}{{\chi }^2}\right]{\delta }_D\left(\chi -x\right).$$

$$b_{{\ell }_1{\ell }_2{\ell }_3}^{v^{(2{)}^{\prime }}}\left(z_1,z_2,z_3\right)=$$

$$D_{\ell }=-\frac{{\partial }^2}{\partial {\chi }^2}+\frac{2}{\chi }\frac{\partial }{\partial \chi }+\frac{\ell \left(\ell +1\right)-2}{{\chi }^2}.$$

$$b_{{\ell }_1{\ell }_2{\ell }_3}^{v^{(2{)}^{\prime }}}$$

$$b_1(z_i)T_{\delta }\left(k_i,z_i\right)j_{{\ell }_i}\left(k_i{r}_i\right)\to b_1(z_i)T_{\delta }\left(k_i,z_i\right)j_{{\ell }_i}\left(k_i{r}_i\right)-f(z_i)T_{\delta }\left(k_i,z_i\right)j_{\ell }^{\prime \prime }\left(k_i{r}_i\right).$$

$$B^{{\delta }^{(2)}}({\mathbf{n}}_1,{\mathbf{n}}_2,{\mathbf{n}}_3,z_1,z_2,z_3)=$$

Full text

1 Introduction

The next generation of galaxy surveys will open a new era of precision cosmology and yield new capacities to study the observable Universe via three-dimensional matter distribution. The optical/infrared surveys that are planned for Euclid [1], LSST [2] and DESI [3], will cover huge volumes of the Universe with very high galaxy numbers, and errors on standard cosmological parameters will be dominated by systematics. The new radio surveys planned for the SKA [4] will cover even larger volumes, but new challenges from systematics in the radio will need to be overcome.

These upcoming surveys are based on tremendous advances in experimental precision. In order to fully exploit their great potential, it is necessary to develop also new theoretical tools and improve theoretical precision. In this spirit, we develop a new analysis of the angular bispectrum for galaxy number counts in terms of directly observable quantities: redshift $z$ and angular position ${\bm{n}}$. In addition to the power spectrum, the bispectrum will be increasingly important for improving constraints and breaking degeneracies (see e.g. [5]). So far, most analyses have used the Cartesian Fourier-space bispectrum (see e.g. Refs. [6, 7]) or the 3-point correlation function (for a recent treatment including redshift space distortions (RSD) see e.g. [8]), which imposes a plane-parallel approximation, and is unable to incorporate the effects of lensing magnification. For current surveys, this may be a reasonable approximation. But next-generation surveys require the inclusion of wide-angle correlations, given their large sky area, and of lensing magnification effects, given their high redshift reach. The direct way to include full-sky and lensing magnification effects is to use the angular bispectrum, as explained in detail in [9] (see also [10, 11, 12]).

In the analysis developed in Ref. [9] the nonlinear clustering bias and the full effect of RSD were not included. Here we include the nonlinear bias, and for the first time, we include the full RSD in the numerical computation of the angular bispectrum (up to terms only relevant on scales close to the Hubble horizon [13] that are nonetheless straightforward, although cumbersome, to further incorporate in our formalism). Including efficiently RSD perturbations in terms of the 3-point function or angular bispectrum is a non-trivial problem, and some different approaches have been recently proposed, see Ref. [14].

We derive the bispectrum of galaxy number counts in terms of directly observable quantities – angles and redshifts. Our approach presents two advantages: first, by adopting only observable quantities we do not need to assume any cosmological model to convert angles and redshifts into distances; second, it naturally provides a full-sky description of the bispectrum. In particular it includes the description of effects imprinted at the largest scales, e.g. primordial non-Gaussianity and relativistic lightcone effects as well as wide-angle effects, which are typically of the same magnitude as these other effects (see e.g. [15, 16, 17, 18, 19]).

We also apply our formalism to describe the angular HI intensity mapping bispectrum. Indeed while HI intensity mapping foregrounds can limit the ability to perform tomography, the need for very thin redshift bins requires one to include RSD perturbations in the angular power spectrum and bispectrum of intensity mapping. At the level of bispectrum, this has never been considered previously and our approach introduces a proper description of RSD on very narrow redshift bins.

In section 2 we introduce the dominant terms contributing to galaxy number counts up to second order in redshift space, including galaxy bias. In section 3 we compute the angular bispectrum and present a novel derivation of the RSD term, well-suited for numerical estimation. In section 4 we present numerical results and discuss the accuracy of the Limber approximation. In section 5 we explain how our formalism can be readily adapted for intensity mapping studies and we discuss its limitations. We conclude in section 6. Appendix A discusses details about the RSD bispectrum relevant for numerical computations. In appendix B we provide an alternative derivation of the density bispectrum, closer to the novel approach introduced here for RSD. In appendix C we list geometrical factors relevant for the angular bispectrum. Finally, in appendix D we compute cosmic variance of the galaxy number count bispectrum at lowest order.

2 Galaxy number counts

A spectroscopic galaxy survey provides the redshift $z$ and the angular position ${\bm{n}}$ for each source. From the number of galaxies $N\left({\bm{n}},z\right)$ we can define the galaxy number count fluctuation as

[TABLE]

where $\langle(\cdots)\rangle$ denotes the angular average at fixed observed redshift $z$. Assuming Gaussian initial conditions, non-Gaussian correlations are generated by non-linear gravitational clustering. To determine the tree-level bispectrum we need to compute the fluctuations up to second order in perturbation theory. This has been done recently in terms of galaxy number counts by [20, 21, 13] (see also [22, 23] for early works).

In our work we consider a perturbed FLRW metric described by

[TABLE]

where the purely scalar metric perturbations agree (to first order) with the gauge-invariant Bardeen potentials. At linear order, considering only the dominant terms, which scale as $\left(k/\mathcal{H}\right)^{2}\Phi^{(1)}\sim\delta^{(1)}$, i.e. the terms which dominate on sub-Hubble scales, we have222See [24, 25] for the full relativistic expression.

[TABLE]

where $\mathcal{H}=\dot{a}/a$ is the comoving Hubble parameter, $r(z)$ is the comoving distance to redshift $z$, $\delta^{(1)}$ is the first-order matter density perturbation in the comoving gauge and $v^{(1)}$ is the first-order velocity perturbation in the longitudinal gauge.

We neglect here the lensing term which is parametrically of the same order but usually, for $z<2$, with equal redshifts and thin z-bins, the lensing term is much smaller than density and redshift space distortion. Also, the novelty of the present work focuses on RSD. The lensing contributions to the tree-level bispectrum have been already computed in Ref. [9]. It is worth pointing out that lensing and RSD dominate in opposite regimes. Indeed while lensing dominates the cosmological signal for large radial separation, or at high redshifts in a signal averaged over a wide redshift binning, RSD decay quickly for sources well separated in redshift and in wide redshift bins.

At second order, we use the convention $\Delta=\Delta^{(1)}+\Delta^{(2)}$. Again including only the dominant terms, which scale as $\left(k/\mathcal{H}\right)^{4}[\Phi^{(1)}]^{2}$, we have

[TABLE]

where $s$ is related to the clustering bias tidal field (see below). The full second-order expression is much more cumbersome, covering several pages (see [20, 21, 13] for number counts, and [26] for intensity mapping). In this ‘quasi-Newtonian’ approximation, we can set $\Phi=\Psi$.

For the nonlinear clustering bias, we assume a local bias model and neglect stochastic bias terms. Following  [27], we use the convention

[TABLE]

Here $\delta=\delta^{(1)}+\delta^{(2)}$ is the matter over-density in comoving gauge and the bias coefficients are scale-independent. The tidal bias coefficient is $b_{s}$, where $s^{2}=s_{ij}s^{ij}$ and the tidal field is given by

[TABLE]

In the simplest local bias model, there is no tidal bias at the time of galaxy formation (see e.g. Ref. [28, 27]) which leads to

[TABLE]

We adopt the bias parameters presented in Ref. [7] for a Euclid-like survey, converted to our bias convention. The different bias values are shown in Fig. 1. Here we provide also a simple polynomial fit

[TABLE]

3 The angular bispectrum

In this section we compute the tree-level bispectrum induced by the 3-point function containing one second order term,333In principle we are interested in the connected part of the bispectrum, which means we have to replace $\Delta^{(2)}({\bm{n}})$ by $\Delta^{(2)}({\bm{n}})-\langle\Delta^{(2)}({\bm{n}})\rangle$. But due to statistical isotropy $\langle\Delta^{(2)}({\bm{n}})\rangle$ only contributes to the monopole and not to the bispectrum $B_{\ell_{1}\ell_{2}\ell_{3}}^{m_{1}m_{2}m_{3}}$ for $\ell_{i}\neq 0$. For this reason we ignore this subtlety here.

[TABLE]

where $\circlearrowleft$ denotes the two additional permutations where the second order term is evaluated at $\left({\bm{n}}_{2},z_{2}\right)$ and $\left({\bm{n}}_{3},z_{3}\right)$, respectively. We expand $\Delta$ in spherical harmonics,

[TABLE]

so that

[TABLE]

where

[TABLE]

Statistical isotropy demands that $B({\bm{n}}_{1},{\bm{n}}_{2},{\bm{n}}_{3},z_{1},z_{2},z_{3})$ only depends on the scalar products ${\bm{n}}_{1}\cdot{\bm{n}}_{2}\,,~{}{\bm{n}}_{1}\cdot{\bm{n}}_{3}$ and ${\bm{n}}_{1}\cdot{\bm{n}}_{2}$. This dictates the $m_{i}$ dependence of the bispectrum,

[TABLE]

where ${\cal G}^{m_{1},m_{2},m_{3}}_{\ell_{1},\ell_{2},\ell_{3}}$ is the Gaunt integral which is simply related to the Clebsch–Gordan coefficients [29]

[TABLE]

In the last equation we have expressed the Gaunt integral in terms of the Wigner $3j$ symbols, see e.g. [29]. The Gaunt integral is invariant under permutations of its columns and it is non-vanishing only if $m_{1}+m_{2}+m_{3}=0$ and the triangle inequality is satisfied, i.e., $|\ell_{2}-\ell_{3}|\leq\ell_{1}\leq\ell_{2}+\ell_{3}$. Furthermore, the sum $\ell_{1}+\ell_{2}+\ell_{3}$ has to be even. The quantity $b_{\ell_{1}\ell_{2}\ell_{3}}$ is called the reduced bispectrum. It contains all non-trivial physical information.

We assume Gaussian initial conditions so that linearly evolved perturbations are Gaussian and $\langle\Delta^{(1)}({\bm{n}}_{1},z_{1})\Delta^{(1)}({\bm{n}}_{2},z_{2})\Delta^{(1)}({\bm{n}}_{3},z_{3})\rangle=0$. The first non-vanishing contribution to the bispectrum is therefore $B_{\rm{tree}}$ defined in (3.2), and we want to compute this contribution to $b_{\ell_{1}\ell_{2}\ell_{3}}(z_{1},z_{2},z_{3})$.

Below we show the explicit form for all the different contributions to $\Delta^{(2)}({\bm{n}},z)$. In order to somewhat simplify the notation, we derive the bispectrum first by setting $\Delta^{(1)}\left({\bm{n}},z\right)\sim b_{1}\delta^{(1)}\left({\bm{n}},z\right)$. Including the linear redshift space distortion is however simple and we do it at the end, see eq. (3.25).

3.1 RSD $\mathcal{H}^{-1}\partial_{r}^{2}v^{(2)}$

We start considering the contribution of the pure second order redshift-space distortion. At this point we want to stress that the RSD contribution to the bispectrum has not been computed previously in terms of directly observable quantities like the full-sky angular bispectrum. Our approach does not rely on any (geometrical) approximation, like flat-sky or Limber. Let us point out that we will consider the Newtonian kernels [30, 31] in the rest of the paper, which is in line with neglecting terms of higher order in $\mathcal{H}/k$, see [32]. Our results can be generalized to an arbitrary (separable) non-linear velocity kernel $G_{2}\left({\bm{k}}_{1},{\bm{k}}_{2}\right)$, including the full higher order GR dynamics [33, 32, 34, 35].

We follow the approach developed in Ref. [11] in order to exchange the order of integration, providing a fully analytical expression for the bispectrum. We begin with the following 3-point function

[TABLE]

In Fourier space the second order velocity potential is given by [31]

[TABLE]

where $f=d{\rm ln}D_{+}/d{\rm ln}a$ in terms of the linear growth factor $D_{+}$ and scale factor $a$,444As discussed in section 4, for numerical computations we interface our code to the class Boltzmann solver, where the growth factor enters only implicitly in perturbation equations via the velocity divergence $\Theta^{(1)}({\bm{k}},z)=-\mathcal{H}(z)f({\bm{k}},z)\delta^{(1)}({\bm{k}},z)$ [36]. The presence of radiation and neutrinos actually leads to a scale-dependent growth function, but we neglect this small effect. The second-order solutions presented here are themselves valid for an Einstein-de Sitter Universe, and in a $\Lambda$CDM cosmology, the growth factor can be simply replaced by the corresponding one, see e.g. [31] for details. and

[TABLE]

We replace the angular dependence ${\bm{k}}_{2}\cdot{\bm{k}}_{3}$ with $k^{2}_{1}=\left({\bm{k}}_{2}+{\bm{k}}_{3}\right)^{2}=k_{2}^{2}+k_{3}^{3}+2{\bm{k}}_{2}\cdot{\bm{k}}_{3}$,

[TABLE]

In this way we have separated the velocity kernel $G_{2}$ into terms proportional to different powers of $k_{1}$ which contain the angular dependence induced by the scalar product ${\bm{k}}_{2}\cdot{\bm{k}}_{3}$.

As in Ref. [9], we now rewrite the Dirac-delta distribution in eq. (3.12) as

[TABLE]

With the same manipulations as detailed in Section 3.4.1 of Ref. [9] we then obtain for the $\mathcal{H}^{-1}\partial_{r}^{2}v^{(2)}$ contribution to the bispectrum

[TABLE]

where $T_{\delta}\left(k,z\right)$ denotes the matter density transfer function, $P_{R}\left(k\right)$ is the primordial curvature power spectrum and $r_{i}=r(z_{i})$ is the comoving distance to redshift $z_{i}$. We assume transfer functions to be normalized to some initial value of curvature perturbations, as in [37]. Given the symmetry in the second and third indices of eq. (3.11), there are only two more permutations (the even ones) of the doublets $(\ell_{i},z_{i})$ to be added to the one shown here and denoted by $\circlearrowleft$.

While in principle it is possible to reduce the integral by one dimension by integrating analytically over $\chi$ as done in Ref. [38], we show below that it is more convenient to exchange the order of integration and express the result in terms of Dirac deltas. However, as pointed out in [9, 38, 11], we stress that in general care must be taken when exchanging the order of integration to avoid UV divergences. While redshift binning (necessary in practice) alleviates UV issues, divergences can be safely avoided only if redshift bins are larger than the non-linearity scale, which does not correspond to an optimal configuration for extracting the maximal amount of cosmological information [39], and it would severely degrade the redshift accuracy of spectroscopic galaxy surveys.

The main difference between the integrand of the bispectrum induced by RSD with respect to density perturbations (see Appendix B) comes from the second derivative of one spherical Bessel function. Indeed, while the angular dependence in the kernel $G_{2}$, parametrized by the variable $k_{1}$, can be expressed by a few Legendre polynomials (as for the density perturbation in Sec. 3.2), the second order derivative of the Bessel function introduces an additional factor $k_{1}^{-2}$. Because of this factor we cannot expand the part containing $G_{2}^{(0)}$ in a finite number of Legendre polynomials555See reference [14] for an alternative expansion.. To solve this problem we write

[TABLE]

This also allows us to use the identity [11]

[TABLE]

to evaluate the integral over the second derivative of the delta-Dirac distribution. With this we can integrate analytically the products666The Heaviside function is normalized such that $\Theta_{H}\left(0\right)=1/2$.

[TABLE]

Using these identities as well as (3.17) we can perform the $k_{1}$ integrations of (3.16) and we find the following expression for the RSD contribution to the bispectrum

[TABLE]

where we have introduced the dimensionless power spectrum $\mathcal{P}_{R}(k)=\frac{k^{3}}{2\pi^{2}}P_{R}(k)$ and defined the differential operator

[TABLE]

For the second equal sign we have also performed the $x$-integration and an integration by parts on the $G^{(0)}_{2}$ term. For the second and the fourth terms we could perform both the $x$ and the $\chi$ integrations. The first term still contains a triple integral. The terms proportional to $G^{(2)}_{2}$ and $G^{(4)}_{2}$ are simply sums of products of 1-dimensional integrals.

Summing up the previous expression we obtain

[TABLE]

In order to add also the RSD in the linear term, we simply replace

[TABLE]

The important point for numerical efficiency when evaluating $b^{v^{(2)^{\prime}}}_{\ell_{1}\ell_{2}\ell_{3}}$ is that the kernels $G_{2}^{(i)}\left(k_{2},k_{3}\right)$, with $i=0,2,4$, are separable in $k_{2}$ and $k_{3}$. Because of that we can reduce the dimensionality of the integrals of eq. (3.24). In Appendix A we rewrite eq. (3.24) explicitly as a sum of products of 1-dimensional integrals, in terms of generalized spectra defined in eq. (3.31), except for the terms involving an additional integral along $\chi$. Indeed, naively exchanging the integration order for these latter terms would lead to divergences for a $\Lambda$CDM cosmology due to the powers of $k_{2}$ and $k_{3}$ in the kernel $G_{2}^{(0)}$.

3.2 Density $\delta^{(2)}$

The second term we consider is the second order density perturbation. In order to have a consistent bias expansion we consider also the local bias parameters $b_{1}$, $b_{2}$ and $b_{s}$

[TABLE]

with

[TABLE]

and

[TABLE]

Again, we follow [9], where the bispectrum is computed separately for the monopole $b^{\delta 0}_{\ell_{1}\ell_{2}\ell_{3}}$, dipole $b^{\delta 1}_{\ell_{1}\ell_{2}\ell_{3}}$ and quadrupole $b^{\delta 2}_{\ell_{1}\ell_{2}\ell_{3}}$ terms for the case $b_{1}=1$, $b_{2}=b_{s}=0$, where the multipole expansion refers to the Legendre polynomials in eq. (3.28). Indeed, it is trivial to generalize this result. We denote the term proportional to $b_{2}/2$ by $b^{\delta^{2}}_{\ell_{1}\ell_{2}\ell_{3}}$ and the tidal term proportional to $b_{s}$ by $b^{s^{2}}_{\ell_{1}\ell_{2}\ell_{3}}$. The term proportional to $b_{2}/2$ is just a monopole with pre-factor one, hence

[TABLE]

while the tidal term becomes777 In Fourier space the tidal term can be expressed as

$s^{2}\left({\bm{k}},z\right)=\frac{1}{(2\pi)^{3}}\int d^{3}k_{1}d^{3}k_{2}\delta_{D}\left({\bm{k}}-{\bm{k}}_{1}-{\bm{k}}_{2}\right)S_{2}\left({\bm{k}}_{1},{\bm{k}}_{2}\right)\delta\left({\bm{k}}_{1},z\right)\delta\left({\bm{k}}_{2},z\right)$

where

$S_{2}\left({\bm{k}}_{1},{\bm{k}}_{2}\right)=-\frac{1}{3}+\frac{\left({\bm{k}}_{1}\cdot{\bm{k}}_{2}\right)^{2}}{k_{1}^{2}k_{2}^{2}}=\frac{2}{3}P_{2}\left(\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{2}\right)$

Comparing with the quadrupole of eq. (3.28) we obtain eq. (3.30).

[TABLE]

We define generalized angular power spectra as in [9]:

[TABLE]

Here $\Delta^{A}_{\ell}\left(k,r\right)$ is the angular transfer function related to the perturbation $A$. In particular we need (using the same notation as [37, 9, 40, 38])

[TABLE]

where $T_{\delta}(k,z)$ denotes the linear transfer function for the density $\delta$ and for the velocity we have $T_{V}(k,z)=-(\mathcal{H}(z)/k)f(z)T_{\delta}(k,z)$. For $n=0$ we shall drop this pre-superscript in $\prescript{n}{}{c}^{AB}_{\ell\ \ell^{\prime}}$ and if $\ell^{\prime}=\ell$ we just indicate it by one subscript $\ell$ so that e.g. $\prescript{0}{}{c}^{AB}_{\ell\ \ell}\equiv{c}^{AB}_{\ell}$.

With this notation we can write

[TABLE]

where we introduce

• •

Monopole:

[TABLE]

• •

Dipole:

[TABLE]

The geometrical factors $g_{\ell_{1}\ell_{2}\ell_{3}}$ and $Q_{\ell\ \ell^{\prime}\ell^{\prime\prime}}^{\ell_{1}\ell_{2}\ell_{3}}$ are defined in Appendix C. The quantity $Q_{1\ \ell^{\prime}\ell^{\prime\prime}}^{\ell_{1}\ell_{2}\ell_{3}}$ is zero unless $\ell^{\prime}=\ell_{2}\pm 1$ and $\ell^{\prime\prime}=\ell_{1}\pm 1$ so that $i^{\ell^{\prime}+\ell^{\prime\prime}}(-i)^{\ell_{1}+\ell_{2}}=\pm 1$.

• •

Quadrupole:

[TABLE]

where $Q_{2\ \ell^{\prime}\ell^{\prime\prime}}^{\ell_{1}\ell_{2}\ell_{3}}$ is zero unless $\ell^{\prime}=\ell_{2}\pm 2,\ell_{2}$ and $\ell^{\prime\prime}=\ell_{1}\pm 2,\ell_{1}$. We then have again $i^{\ell^{\prime}+\ell^{\prime\prime}}(-i)^{\ell_{1}+\ell_{2}}=\pm 1$.

In Appendix B we also present an alternative derivation of the density bispectrum using the approach adopted in the previous section for RSD.

3.3 Products of linear terms

Using Wick’s theorem, one finds that all the terms in the second line of eq. (2.4) induce a bispectrum which is the sum of products of certain power spectra. See [9, 38] for an explicit computation. In detail we find

[TABLE]

The two additional permutations are the even ones of the doublets $(\ell_{i},z_{i})$.

3.4 The full Bispectrum

Adding all the terms computed in the previous sections, we obtain the leading galaxy bispectrum as

[TABLE]

4 Numerical results

4.1 The different contributions to the bispectrum

In this section we show the numerical redshift dependent angular bispectrum for different configurations. All the results are obtained with the following cosmological parameters: $h=0.67$, $\Omega_{b}=0.05$, $\Omega_{cdm}=0.27$, and consistently with the previous section vanishing curvature. The amplitude of the primordial curvature power spectrum is set to $A_{s}=2.3\times 10^{-9}$, the pivot scale is $k_{\rm pivot}=0.05{\rm Mpc}^{-1}$, the spectral index is $n_{s}=0.962$ and we assume no running. As discussed in Section 3, bispectra are computed correlating a second-order term with two first-order terms, where the latter include density and RSD. Clustering bias (figure 1) is always included.

Our numerical computations rely on byspectrum,888The code is found at https://gitlab.com/montanari/byspectrum. a Python code making use of a low-level C++ library for the computation of generalized spectra that wraps the class code [36, 37] to retrieve the transfer functions.

In figures 2-4, we plot the ratio between the bispectrum contributions normalized with respect to the square of the power spectra (which include density and RSD perturbations), i.e.

[TABLE]

We set $z_{1}=z_{2}=z_{3}=1$ and show the dependence on the three multipoles. Considering equal redshifts has the advantage of allowing a straightforward relation between multipoles and modes in Fourier space $k\approx\ell/r(z)$, given the line-of-sight comoving distance $r(z)$. The choice $\ell_{1}=200$ then corresponds to a comoving wavenumber $k_{1}\approx 0.06/$Mpc (or equivalently a 100 Mpc comoving scale), which is quite linear at $z=1$. The leading terms are typically of the same order, $\sim 0.1$ to $0.2$ for the different configurations. The density dipole term amounts to about 10 to 30% of the monopole while the quadrupole is 3 to 10%. The terms shown in the second row of fig. 4 are significantly smaller. Indeed the latter, even if computed at the same redshift, contain integrals of spherical Bessel function which are out of phase leading to a strong suppression of their amplitudes. All velocity terms have the largest power in the squeezed limit (top-left corner in the figures), whereas the density monopole somewhat prefers the equilateral shape. Note, however, that the bispectra are nearly constant with little variation over the range shown in the figures.

In figure 5 we show the different contributions to the tree-level bispectrum in the configuration $\ell_{1}=\ell_{2}=\ell_{3}/2$ at equal redshifts $z_{1}=z_{2}=z_{3}=1$ as a function of $\ell_{1}$. Clustering bias (figure 1) is very important and due to a negative $b_{2}(z)$ the term $b^{v^{\prime}\delta}_{\ell_{1}\ell_{2}\ell_{3}}$ dominates at large scales. Consistently with the previous plots, contributions involving the velocity terms $b^{\delta^{\prime}v}_{\ell_{1}\ell_{2}\ell_{3}}$ and $b^{v^{\prime\prime}v}_{\ell_{1}\ell_{2}\ell_{3}}$ are significantly suppressed with respect to the other terms. The physical oscillations in these velocity terms are roughly in phase.

In figure 6 we consider the density term as a function of redshift difference for $z_{1}=0.7$ and $z_{2}=z_{3}=z_{1}+\Delta z$ in the following configurations999We refer to a given shape considering only multipole triangles. Of course, a given shape at equal or at different redshifts corresponds to a different configuration in comoving space. In particular, care should be taken when comparing to results obtained in Fourier space (at equal redshifts). (always such that $\ell_{1}+\ell_{2}+\ell_{3}$ is even):

• •

Equilateral, $\ell_{1}=\ell_{2}=\ell_{3}$.

• •

Squeezed, $\ell_{1}\ll\ell_{2}=\ell_{3}$.

• •

Folded, $\ell_{1}=\ell_{2}=\ell_{3}/2$.

The density bispectrum is suppressed by several orders of magnitude already at relatively small redshift differences $\Delta z$. In all cases the bispectrum reaches a negative minimum and then tends to zero. In the squeezed case the trough is wider and the suppression is somewhat less severe.

We also compare the density bispectrum to its cosmic variance. We use eq. (C.2) to write the angle-averaged bispectrum as

[TABLE]

and we define the signal-to-noise ratio for a fixed multipole configuration as

[TABLE]

The variance $\sigma_{B_{\ell_{1}\ell_{2}\ell_{3}}}$ is calculated in Appendix D. Correlations at different redshifts are subdominant compared to equal redshift contributions in the computation of the cosmic variance $\sigma_{B_{\ell_{1}\ell_{2}\ell_{3}}}^{2}$ given by equation (D). Hence, we neglect terms $c_{\ell}(z_{i},z_{j})$ with $z_{i}\neq z_{j}$ in the variance.

For all shapes101010While strictly speaking the squeezed limit corresponds to $\ell_{1}\ll\ell_{2}=\ell_{3}$, here we loosen this definition. we find a $S/N\gtrsim 0.01$ at $z_{1}=z_{2}=z_{3}=0.7$, leading to a presumably large cumulative $S/N$ when summing over all triangles. The results also show that, compared to cosmic variance, the density contribution is more sensitive to the folded and squeezed cases and the equilateral case is more strongly suppressed. We checked that the folded configuration also yields better results than the case $\ell_{1}\neq\ell_{2}\neq\ell_{3}$ (not shown in the figure). In the squeezed case, we verified that increasing the long mode $\ell_{1}$ decreases the $S/N$. Note the turnover in $S/N$ in the equilateral configuration, due to the fact that the signal decreases faster than the variance at large multipoles. In all cases the $S/N$ decreases with redshift, as non-Gaussianities get weaker. However, for different redshifts with $\Delta z\gtrsim 0.1$, for high redshifts, $z>2$ and wide redshift bins, lensing may not be negligible and hence should be included both in the signal and in the variance; see Ref. [9] for numerical results including lensing.

We stress that the $S/N$ estimates presented in this section are only meant to study the effect of cosmic variance on different bispectrum configurations. A realistic estimate must include the effect of finite radial selection functions and of shot-noise, the latter being particularly relevant at small scales.111111The validity of a tree-level bispectrum should also be better assessed at the largest multipoles and smallest redshifts here considered.

4.2 Limber approximation for density and RSD

In this section we introduce a redshift binning. To reduce the numerical effort we adopt the Limber approximation [41, 42]

[TABLE]

For the density contribution we consider the monopole of the reduced bispectrum and we compare the Limber approximation with the exact solution. We also neglect RSD contributions to the linear perturbations. Introducing redshift binning with a normalized window function centered at $z$, $W(z,z^{\prime})$, we define the following z-binned bispectrum (denoted with an over-bar)

[TABLE]

For the density bispectrum we apply the Limber approximation directly on eq. (B.6)

[TABLE]

Here we have used the same binning for each variable $z_{1},z_{2},z_{3}$ and $B_{\delta\delta\delta}\left(k_{1},k_{2},k_{3}\right)$ is the density bispectrum in Fourier space, see e.g. Ref. [31]. In fig. 8 we show the accuracy of Limber approximation for different z-binning. It is evident that the Limber approximation is not accurate to describe redshift bins more narrow than about $\Delta z\sim 0.1$ for $\ell<500$ if we require an accuracy of 10% or better.

For the RSD term we can not apply Limber approximation on all the Bessel functions, due to its second derivative appearing in eq. (3.16). Therefore, we consider

[TABLE]

such that we can apply the Limber approximation on the first term directly. For the second term we generalize the Limber approximation, which is derived from the approximation

[TABLE]

yielding

[TABLE]

Using the ordinary Limber approximation (4.4) on the first term of eq. (4.7) and the generalized approximation (4.9) on the second term, we obtain

[TABLE]

Interestingly the scale dependence of the two contributions at large $\ell$ (limit of validity of Limber approximation) are given respectively by

[TABLE]

and

[TABLE]

The two leading contributions, that scale as $\ell_{1}^{-1}$, cancel exactly leaving the z-binned RSD bispectrum scaling as $\ell^{-2}$ with respect to the density perturbations. Considering that the contribution of RSD to the z-binned power spectrum scales as $\ell^{-1}$ with respect to density perturbations and that we do expect the bispectrum to scale roughly like the square of the power spectrum, we can use the latter as a proxy to investigate the importance of RSD for different z-binning. This finding, together with the corresponding results for the power spectrum shown in Fig. 9, suggests that for binned bispectra with a wide enough binning compatible with the Limber approximation, RSD is negligibly small. Only for very slim redshift bins, as are possible for spectroscopic surveys, is RSD measurable. Furthermore, for such bins the Limber approximation cannot be trusted.

5 The HI intensity mapping angular bispectrum

Intensity mapping experiments provide maps of the emitted intensity of a given molecular or atomic line. Here we consider the $21$cm line from the spin flip (hyperfine structure) of neutral hydrogen (HI). We are interested in the low redshift, i.e. the post-reonization universe, on scales where it is well described by quasi-linear physics [43, 44]. A $21$cm intensity survey does not resolve single galaxies, but it measures the superposed diffuse emission of several sources. Therefore intensity mapping experiments do typically not provide high angular resolution. However, as they directly measure frequencies, they allow for very narrow redshift bins. Moreover, $21$cm intensity surveys face the challenge to clean the cosmological signal of foregrounds which are several orders of magnitude larger. Foreground cleaning approaches [45] are based on the different frequency dependence of foregrounds which are typically smooth while the 21cm line is very sharply peaked. This allows us to recover the 21cm intensity very well on small radial scales, but recovery becomes increasingly difficult on large scales, on which also the 21 cm intensity should have smooth correlations. This may limit for practical purposes the use of HI mapping tomography since the long range radial correlations are lost by the cleaning of the much stronger foreground. We plan in the future to include the angular formalism developed in this work in a full 3-dimensional Fourier angular spectrum [46, 47, 48] to precisely quantify and possibly overcome this limitation.

On the other hand, the HI mapping angular bispectrum, being limited to very thin z-bins, requires one to include RSD developed in our formalism. Indeed, while in galaxy surveys broad z-bins can strongly reduce the impact of RSD on the observable power spectra or higher order statistics, broad z-bins in HI mapping remain foreground dominated. Previous works, for instance [49], have not included RSD in the bispectrum analysis. Moreover they based their results on the Limber approximation, which, as we have shown, fails completely for the slim z-bins required by HI mapping.

Without the need to resolve single galaxies, intensity mapping experiments can more efficiently cover a larger fraction of the sky and go deeper in redshift, providing more information on the largest scales we can observe. A $21$cm intensity mapping survey provides the brightness temperature $T_{b}\left({\bm{n}},z\right)$ of neutral hydrogen [50, 43] as a function of sky direction ${\bm{n}}$ and observed redshift $z$. We can therefore define the brightness temperature fluctuation as

[TABLE]

where $\langle(\cdots)\rangle$ denotes the angular average at fixed observed redshift $z$. The brightness temperature in 21cm intensity is proportional to the observed neutral hydrogen number density $n_{\text{HI}}$ per unit surface and per redshift bin (see [51], eq. (22))

[TABLE]

where $d_{A}$ is the angular diameter distance.

Therefore, to any order, the HI intensity mapping fluctuation differs only by a convergence term (lensing) from the galaxy number counts. Since in our work we do not consider any term induced by the lensing potential (see Ref. [52] for HI intensity mapping lensing at second and third order) our results can be directly interpreted also in terms of the HI angular bispectrum, providing we use the appropriate bias parameters [53].

6 Conclusions

In this work we have computed the tree level angular-redshift bispectrum for number counts taking into account the contribution from density and redshift space distortions. For equal redshifts, $z_{1}=z_{2}=z_{3}$ and very narrow redshift bins as are achieved for spectroscopic number count surveys and for intensity mapping, these are the only relevant contributions for $\ell_{i}\gtrsim 10$. On large scales, $\ell\sim 10$ relativistic and wide angle effects are relevant, while for unequal redshifts and for wide redshift bins, lensing terms are relevant.

We have presented results using the bias prescription expected to hold for a Euclid-like survey. We have found that the normalized bispectrum contributions are relatively constant and of the order 0.1 to 0.2. All contributions are very similar and not strongly scale or shape dependent, apart from the $\delta^{\prime}v$ and $v^{\prime\prime}v$ terms which are suppressed by several orders of magnitude due to the fact that these terms are out of phase.

The $S/N$ of the folded configuration is largest, about a factor of 3 larger than the squeezed contribution, and up to a factor of 30 larger than the equilateral contribution where the $S/N$ of a single mode remains of order 0.01 up to $\ell=1000$, while the single mode values for the folded and squeezed configurations at $\ell\sim 1000$ become of order 0.1. In this analysis only cosmic variance is included in the noise, so it certainly over-estimates the value, but it indicates that it should be straightforward to detect this bispectrum from nonlinearities by summing over several modes.

We have considered redshift binning and shown that the Limber approximation is unreliable if the redshift bins are too slim. For $\Delta z=0.01$ we find deviations from the correct density bispectrum of up to a factor 8. This is due to the relatively high UV power of density perturbations which are a factor $k^{2}/\mathcal{H}^{2}$ larger than the gravitational potential. For sufficiently wide redshift bins which average over scales smaller than $\Delta z/\mathcal{H}(z)$, this UV power is reduced which renders the Limber approximation less inaccurate. Even though we cannot compute the binned RSD without Limber approximation we are convinced that the same is true for RSD since it has the same UV behavior as the density term. Introducing significant redshift binning reduces the RSD bispectrum by about a factor $\ell^{-2}$ faster than the density bispectrum so that in the regime where the Limber approximation is reasonably reliable, $\ell\gtrsim 100$ and $\Delta z\gtrsim 0.1$, RSD can actually be neglected.

The decomposition into spherical harmonics has always been used for the statistical description and analysis of CMB physics, but the subtle UV issues we investigated here do not affect CMB observables because the latter are proportional to the metric perturbations and not to their second derivatives, i.e., density and RSD. Therefore, their transfer functions decay much faster with $k$ and one does not encounter UV divergences.

The present study of the bispectrum is still preliminary and can be improved and probably sped up in multiple ways. For example, while in this work we express the bispectrum for all the terms except RSD as products of generalized angular power spectra, in the future we plan to replace numerical line-of-sight integrals by the more efficient FFTlog approaches advocated in [11]. This has proven to speed-up computations by a factor of order $\mathcal{O}(100)$ [54]. This numerical improvement is especially important when using redshift binning which increases the number of integrals and when we want to be efficient enough to sample cosmological parameters with MCMC methods based on a bispectrum likelihood.

Acknowledgements

It is a pleasure to thank Simone Ferraro for useful discussions. ED is supported by the Swiss National Science Foundation (No. 171494). RD is supported by the Swiss National Science Foundation. FM is supported by the Research Project FPA2015-68048-C3-3-P [MINECO-FEDER], and the Centro de Excelencia Severo Ochoa Program SEV-2016-0597. RM and OU are supported by the UK STFC grant ST/N000668/1. RM is also supported by the South African SKA Project and the National Research Foundation (Grant No. 75415).

Appendix A Details for the RSD bispectrum

The aim of this section is to provide the bispectrum for RSD, described in eq. (3.24), in terms of products of 1-dimensional integrals, useful for numerical computations. The main task is to apply the derivative to the products of Bessel functions (last two lines of eq. (3.24)) and to rewrite them in terms of spherical Bessel functions of different order. We warn the reader, that due to the recursive relation for spherical Bessel functions, the final results can be rewritten in several other equivalent forms.

Without the integrated terms along $\chi$ we can rewrite eq. (3.24) in terms of the generalized spectra (3.31), as follows

[TABLE]

The full RSD bispectrum is therefore given by eq. (A.1) plus the $\chi$-integral given below for completeness,

[TABLE]

We solve eq. (A) numerically as a 3-dimensional integral using the Suave Monte Carlo algorithm available from the Cuba library [55].121212http://www.feynarts.de/cuba/. The Suave algorithm does not need to evaluate the integrand at the integration boundaries, so the semi-infinite range of integration can be simply mapped to the unitary interval through a change of coordinates. Alternatively, note that the $\chi$-integrand is strongly peaked around $\chi\sim r_{1}$, so the semi-infinite range can be safely restricted to, e.g., $\chi\in[r_{1}/10,10\,r_{1}]$.

Figure 10 shows a striking cancellation between the separable RSD reduced bispectrum term, eq. (A.1), and the non-separable one, eq. (A), up to two orders of magnitude at the smallest multipoles at $z=0.1$. This suggests that a more convenient factorization where the sum is carried out analytically is desirable. Also note that the non-separable term becomes sub-leading at small scales ($\ell_{1}\gtrsim 100$ in the plotted configurations).

Appendix B Density $\delta^{(2)}$: Alternative derivation

We give an alternative derivation of the density bispectrum to that presented in section 3.2, following the computation presented in section 3.1 for RSD. This is useful to validate the method as we compared the analytical and numerical results for the density reduced bispectrum obtained in the two approaches. Here we neglect clustering bias for brevity, but it is clear how to re-introduce it.

Given eqs. (3.2), (3.27), we obtain the 3-point function

[TABLE]

We expand the Fourier modes as

[TABLE]

and the Dirac-delta distribution as

[TABLE]

Integrating now over the angular part of the Fourier integrals we obtain

[TABLE]

with

[TABLE]

where

[TABLE]

By expanding the integral (B.6), we note that it depends on $k_{1}$ only through

[TABLE]

To solve the integrals (B.10) and (B.11) we have used the identity (3.18). Hence we obtain

[TABLE]

This reduces the integral from a 4-dimensional integral (B.6) to a product of 1-dimensional integrals, since the kernels $F^{(i)}_{2}\left(k_{2},k_{3}\right)$ with $i=0,2,4$ are separable in $k_{2}$ and $k_{3}$. This result agrees with the analytical derivation of [9]. To introduce the linear redshift space distortions we use the substitution given in eq. (3.25).

Appendix C Generalized spectra geometrical factors

We define the geometrical factors used to compute the density bispectrum dipole, eq. (3.40), and quadrupole, eq. (3.41). For details see appendix A of [9].

The factor $g_{\ell_{1}\ell_{2}\ell_{3}}$ relates the reduced bispectrum $b_{\ell_{1}\ell_{2}\ell_{3}}$ to the angle-averaged one. It is given by

[TABLE]

One can show that

[TABLE]

Note that vanishing $g_{\ell_{1}\ell_{2}\ell_{3}}$ values correspond to unphysical multipole combinations, for which the bispectrum is undefined.

The factor $Q_{\ell\ell^{\prime}\ell^{\prime\prime}}^{\ell_{1}\ell_{2}\ell_{3}}$ can be written in terms of the Wigner 6$j$ symbol as

[TABLE]

where $I_{\ell\ell^{\prime}\ell^{\prime\prime}}^{\ell_{1}\ell_{2}\ell_{3}}$ further depends on Wigner 3$j$ symbols:

[TABLE]

Typically, only a few coefficients of $Q_{\ell\ell^{\prime}\ell^{\prime\prime}}^{\ell_{1}\ell_{2}\ell_{3}}$ entering in the equations are non-vanishing. They are determined by the fact that the $6j$ symbol is non-zero only if the triangle condition is satisfied by all the triplets

[TABLE]

The $3j$ symbols further require the following sums to be even:

[TABLE]

which also implies

[TABLE]

To compute Wigner symbols numerically we use the efficient library WIGXJPF131313http://fy.chalmers.se/subatom/wigxjpf/. [56].

Appendix D Cosmic variance for the angular bispectrum

In this section we compute the cosmic variance for an arbitrary redshift-dependent angular bispectrum, generalizing results obtained in CMB studies [57]. We introduce the bispectrum estimator (in term of the spherical harmonic expansion coefficients $a_{\ell m}(z)$) as

[TABLE]

This is related to the full bispectrum through

[TABLE]

Analogously we define the power spectrum estimator

[TABLE]

such that

[TABLE]

Assuming weak non-Gaussianities, the main contribution to the bispectrum cosmic variance comes from the dominant Gaussian part of the harmonic expansion coefficients $a_{\ell m}$. In this approximation $\langle{\hat{B}}_{\ell_{1}\ell_{2}\ell_{3}}(z_{1},z_{2},z_{3})\rangle\approx 0$ and the covariance is given by

[TABLE]

Wick’s theorem gives 15 permutations. However, those permutations involving terms such as $\langle a_{\ell_{i}m_{i}}(z_{i})a_{\ell_{j}m_{j}}(z_{j})\rangle=(-1)^{m_{j}}\langle a_{\ell_{i}m_{i}}(z_{i})a^{*}_{\ell_{j}-m_{j}}(z_{j})\rangle$ (or, similarly, $\langle a^{*}_{\ell_{i}m_{i}}(z_{i})a^{*}_{\ell_{j}m_{j}}(z_{j})\rangle$) vanish. Recalling that $\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}\\ m_{1}&m_{2}&m_{3}\end{pmatrix}$ is zero unless $m_{1}+m_{2}+m_{3}=0$, it is easy to show that such terms give a contribution to the covariance proportional to

[TABLE]

which vanishes if all multipoles are larger than zero, as in our case.141414Bispectrum multipoles $\ell\leq 2$ involve non-linear terms at the observer, not treatable within cosmological perturbation theory. There are 6 remaining non-vanishing terms from Wick’s theorem (those that only involve contractions of coefficients with their complex conjugate $\langle a_{\ell_{i}m_{i}}(z_{i})a^{*}_{\ell_{j}m_{j}}(z_{j})\rangle$). Using the definitions given above and the orthogonality relation

[TABLE]

together with a more compact notation for power spectra $c_{\ell}^{ij^{\prime}}\equiv c_{\ell}(z_{i},z^{\prime}_{j})$, the bispectrum covariance reads

[TABLE]

where $\delta_{\ell_{i}\ell_{j}\ell_{k}}^{\ell_{p}\ell_{q}\ell_{r}}\equiv\delta_{\ell_{i}\ell_{p}}\delta_{\ell_{j}\ell_{q}}\delta_{\ell_{k}\ell_{r}}$. Hence, in the Gaussian approximation the covariance is diagonal in multipole space. It is useful to consider the variance for $\ell_{1}+\ell_{2}+\ell_{3}=$ even:

[TABLE]

Hence, the variance strongly depends on whether none, two or all the multipoles are equal.

At equal redshifts eq. (D) and eq. (D) reduce to the usual CMB expressions [57]

[TABLE]

and

[TABLE]

respectively.