Cosmology with the Galaxy Bispectrum Multipoles:
Optimal Estimation and Application to BOSS Data

The galaxy bispectrum is defined as the three-point expectation of the overdensity, $\delta_{g}$:

$\displaystyle(2\pi)^{3}\delta_{\mathrm{D}}\left({\bm{k}_{123}}\right)B_{\rm ggg}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})\equiv\left\langle{\delta_{g}(\bm{k}_{1})\delta_{g}(\bm{k}_{2})\delta_{g}(\bm{k}_{3})}\right\rangle,$

(3.1)

[e.g., 95], writing $\bm{k}_{123}\equiv\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3}$ for Dirac delta $\delta_{\rm D}$. In real-space, symmetry under translations and rotations forces the bispectrum to be a function only of three variables (usually chosen to be the side lengths $k_{i}\equiv|\bm{k}_{i}|$); this implies $B_{\rm ggg}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})\to B_{\rm ggg}(k_{1},k_{2},k_{3})$. Redshift-space distortions break symmetry with respect to the line-of-sight $\hat{\bm{n}}$ (hereafter LoS), affording an additional two degrees of freedom to the bispectrum. Whilst this can be parametrize in a number of ways, a particularly well-motivated choice of variables are the angle of the triangle plane to the LoS, and the orientation of the triangle within the plane [e.g., 96, 17, 97].

In this approach, one can expand the bispectrum as a spherical harmonic series:

$\displaystyle B_{\rm ggg}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell}B_{\ell m}(k_{1},k_{2},k_{3})Y_{\ell m}(\theta_{\bm{k}},\phi_{\bm{k}}),$

(3.2)

where $\theta_{\bm{k}}$ and $\phi_{\bm{k}}$ specify the aforementioned orientation. Though this basis is complete, measuring $B_{\ell m}$ is difficult, since the spherical harmonic cannot be separably decomposed into $\bm{k}_{1}$, $\bm{k}_{2}$, and $\bm{k}_{3}$ pieces, yielding a non-factorizable estimator. This is not a problem for theoretical forecasts [e.g., 24, 22], but severely limits application to observational data. Consequently, several works [e.g., 17, 25] have considered only the $m=0$ moment (independent of $\phi$), and set $\cos\theta\equiv\hat{\bm{k}}_{3}\cdot\hat{\bm{n}}$, additionally fixing $k_{1}\leq k_{2}\leq k_{3}$. This corresponds to representing the bispectrum as a Legendre series in $\theta$:

$\displaystyle B_{\rm ggg}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})\approx\sum_{\ell=0}^{\infty}B_{\ell}(k_{1},k_{2},k_{3})\mathcal{L}_{\ell}(\hat{\bm{k}}_{3}\cdot\hat{\bm{n}}),\qquad(k_{1}\leq k_{2}\leq k_{3})$

(3.3)

where $\mathcal{L}_{\ell}$ is a Legendre polynomial and $B_{\ell}$ the corresponding coefficient.⁴⁴4Some works [e.g., 22] have instead expanded the bispectrum as a double Legendre series in the two angles. A separable choice would be to expand in, say, $\mathcal{L}_{\ell}(\hat{\bm{k}}_{2}\cdot\hat{\bm{n}})$ and $\mathcal{L}_{\ell^{\prime}}(\hat{\bm{k}}_{3}\cdot\hat{\bm{n}})$; however, the corresponding coefficients are generally difficult to estimate robustly, since the two angles are not independent once the side-lengths are specified. We note that (3.3) is not a strict equality, since the bispectrum contains higher-order moments (with $m\neq 0$) not captured within its formalism; in the below, we will instead define the multipoles directly as integrals over $\theta,\phi$.

The decomposition of (3.3) can be used to construct estimators for the bispectrum multipoles, $B_{\ell}$. For an idealized periodic-box geometry (such as an $N$-body simulation), the conventional estimator for the bispectrum multipoles is given by

	$\displaystyle\left.\widehat{B}^{abc}_{\ell}\right|_{\rm periodic}$	$\displaystyle\equiv$	$\displaystyle\frac{2\ell+1}{N_{T}^{abc}}\int_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}(2\pi)^{3}\delta_{\mathrm{D}}\left({\bm{k}_{123}}\right)\Theta^{a}(k_{1})\Theta^{b}(k_{2})\Theta^{c}(k_{3})$
			$\displaystyle\qquad\qquad\qquad\,\times\,\delta_{g}(\bm{k}_{1})\delta_{g}(\bm{k}_{2})\delta_{g}(\bm{k}_{3})\mathcal{L}_{\ell}(\hat{\bm{k}}_{3}\cdot\hat{\bm{n}}),$

where $\int_{\bm{k}}\equiv(2\pi)^{-3}\int_{\bm{k}}$ [17]. Here, $a\leq b\leq c$ specify a triplet of $k$-bins of finite radius, defined by $\Theta^{i}(k)$, which is unity if $k$ is in bin $i$, and zero else. (3.2) is simply an integral over three copies of the density field weighted by the Legendre polynomial in the longest side $\mathcal{L}_{\ell}(\hat{\bm{k}}_{3}\cdot\hat{\bm{n}})$, with translation invariance enforced by the Dirac delta. This is normalized by the isotropic bin volume, defined by

$\displaystyle N_{T}^{abc}=\int_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}(2\pi)^{3}\delta_{\mathrm{D}}\left({\bm{k}_{123}}\right)\Theta^{a}(k_{1})\Theta^{b}(k_{2})\Theta^{c}(k_{3}).$

(3.5)

In this work, we regard (3.2) as the definition of the binned bispectrum multipoles (rather than the approximate relation of 3.3).

Theoretical predictions for the bispectrum multipoles can be similarly computed from the expectation of (3.2):

	$\displaystyle\left.B^{abc}_{\ell}\right|_{\rm theory}$	$\displaystyle\equiv$	$\displaystyle\frac{2\ell+1}{N_{T}^{abc}}\int_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}(2\pi)^{3}\delta_{\mathrm{D}}\left({\bm{k}_{123}}\right)\Theta^{a}(k_{1})\Theta^{b}(k_{2})\Theta^{c}(k_{3})$
			$\displaystyle\qquad\qquad\,\times\,B_{ggg}^{\rm theory}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})\mathcal{L}_{\ell}(\hat{\bm{k}}_{3}\cdot\hat{\bm{n}}),$

for some theory model $B_{\rm theory}$ which is not yet averaged over angles. This will be discussed in §5.

In practice, we implement (3.2) by factorizing in $\bm{k}_{i}$, following Ref. [17]. This is realized by rewriting the Dirac function as an exponential, yielding the asymmetric expression

$\displaystyle\left.\widehat{B}^{abc}_{\ell}\right|_{\rm periodic}=\frac{2\ell+1}{N_{T}^{abc}}\int d\bm{x}\,F_{0}^{a}(\bm{x})F_{0}^{b}(\bm{x})F_{\ell}^{c}(\bm{x}),\quad N_{T}^{abc}=\int d\bm{x}\,D^{a}(\bm{x})D^{b}(\bm{x})D^{c}(\bm{x}),$

(3.7)

using the definitions

$\displaystyle F_{\ell}^{i}(\bm{x})\equiv\int_{\bm{k}}e^{-i\bm{k}\cdot\bm{x}}\Theta^{i}(k)\delta(\bm{k})\mathcal{L}_{\ell}(\hat{\bm{k}}\cdot\hat{\bm{n}}),\qquad D^{i}(\bm{x})\equiv\int_{\bm{k}}e^{-i\bm{k}\cdot\bm{x}}\Theta^{i}(k).$

(3.8)

Each piece can be straightforwardly evaluated using fast Fourier transforms (FFTs) with $N_{g}\log N_{g}$ complexity for $N_{g}$ grid points. If we had defined the redshift-space components using $Y_{\ell m}(\theta_{\bm{k}},\phi_{\bm{k}})$ rather than $\mathcal{L}_{\ell}(\hat{\bm{k}}_{3}\cdot\hat{\bm{n}})$, (or some other choice) the expression would not factorize in the above manner, and computation would scale as $\mathcal{O}(N_{g}^{3})$.

In realistic surveys, the LoS is not fixed, but varies depending on which galaxies are being considered.⁵⁵5Strictly, a separate line-of-sight is required for each galaxy. The effects of assuming a single line-of-sight are small for typical survey sizes however [cf. 98, 99]. In this case, we can adopt the ‘Yamamoto’ prescription [100, 17], fixing $\hat{\bm{n}}$ to the direction vector of the galaxy associated to $\bm{k}_{3}$. This corresponds to the replacement

	$\displaystyle F_{\ell}^{i}(\bm{x})$	$\displaystyle\to$	$\displaystyle\int_{\bm{k}}e^{-i\bm{k}\cdot\bm{x}}\Theta^{i}(k)\int d\bm{r}\,e^{i\bm{k}\cdot\bm{r}}\delta(\bm{r})\mathcal{L}_{\ell}(\hat{\bm{k}}\cdot\hat{\bm{r}})$
		$\displaystyle\equiv$	$\displaystyle\frac{4\pi}{2\ell+1}\sum_{m=-\ell}^{\ell}\int_{\bm{k}}e^{-i\bm{k}\cdot\bm{x}}\Theta^{i}(k)Y_{\ell m}(\hat{\bm{k}})\int d\bm{r}\,e^{i\bm{k}\cdot\bm{r}}\delta(\bm{r})Y_{\ell m}^{*}(\hat{\bm{r}}),$

with the latter equality allowing for fast estimation using the spherical harmonic addition theorem.

When applying the estimators described in §3 to observational data, we must specify the density field $\delta_{g}$. Usually, this is modelled by the pixelized field of “data-minus-randoms”; $\delta_{g}(\bm{r})\propto n_{g}(\bm{r})-\alpha\,n_{r}(\bm{r})$, where $n_{g}$ is the observed galaxy density field and $n_{r}(\bm{r})$ is the random catalog (containing $1/\alpha$ times more particles than the galaxy catalog). Since both data and randoms are multiplied by the survey mask, conventional estimators will measure only the windowed bispectrum, $B^{\rm win}_{ggg}$, rather than the true underlying statistic, $B_{ggg}$. Before bin integration, the two are related by the following convolution integral:

	$\displaystyle B_{\rm ggg}^{\rm win}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$	$\displaystyle=$	$\displaystyle\int_{\bm{p}_{1}\bm{p}_{2}\bm{p}_{3}}(2\pi)^{3}\delta_{\mathrm{D}}\left({\bm{p}_{123}}\right)$
			$\displaystyle\qquad\qquad\,\times\,W(\bm{k}_{1}-\bm{p}_{1})W(\bm{k}_{2}-\bm{p}_{2})W(\bm{k}_{3}-\bm{p}_{3})B_{\rm ggg}(\bm{p}_{1},\bm{p}_{2},\bm{p}_{3}).$

To compare theory and data, we should similarly convolve the theory model. Due to its oscillatory nature, this is a difficult and time-consuming numerical operation (though see Ref. [52] for a possible $\ell=0$ approach), thus the effect is often ignored or heavily simplified [e.g., 26, 101, 10, 102, 11, 73, 75]. This may lead to biases in data-analysis when large-scale modes (relevant to primordial non-Gaussianity studies) are included.

A major goal of this work is the estimation of unwindowed bispectrum multipoles. These are unbiased by the window function and can be robustly compared to theory models without the need to window-convolve the latter (via 4.1). Our approach follows Refs. [49, 50] for the power spectrum and $\ell=0$ monopole (as well as Ref. [103] for the higher-point CMB correlators), themselves inspired by early work on the subject in [104, 105, 106].

To define unwindowed estimators, we must first express the true bispectrum $B_{ggg}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$ in terms of the quantity of interest: the set of bispectrum coefficients $b_{\alpha}\equiv B^{abc}_{\ell}$ (using $\alpha$ to denote the radial bin indices and multipole). This relation will then be used to form an estimator for $b_{\alpha}$ via maximum-likelihood methods. As an ansatz, we will assume

$\displaystyle B_{\rm ggg}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=\sum_{\alpha}\frac{b_{\alpha}}{\Delta_{\alpha}}\left[\Theta^{a}(k_{1})\Theta^{b}(k_{2})\Theta^{c}(k_{3})\mathcal{L}_{\ell}(\hat{\bm{k}}_{3}\cdot\hat{\bm{n}})+\text{5 perms.}\right].$

(4.2)

This is similar in form to the Legendre decomposition of (3.3), but is defined for all arbitrary ordering of $\{\bm{k}_{1},\bm{k}_{2},\bm{k}_{3}\}$, with the binning functions picking out the relevant permutation, such that we can represent the full bispectrum in terms of its binned components $b_{\alpha}$ with $a\leq b\leq c$. (4.2) includes a bin-specific normalization factor $\Delta_{\alpha}$; this takes a simple form for $\ell=0$ as in Ref. [50] but is more complex in general, as we show below, due to the omitted $\phi$ integrals and exchange symmetry.

Inserting (4.2) into the expectation of our idealized estimator (3.2) gives

	$\displaystyle\left\langle{\widehat{B}^{abc}_{\ell}}\right\rangle$	$\displaystyle=$	$\displaystyle\frac{2\ell+1}{N_{T}^{abc}}\int_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}(2\pi)^{3}\delta_{\mathrm{D}}\left({\bm{k}_{123}}\right)\Theta^{a}(k_{1})\Theta^{b}(k_{2})\Theta^{c}(k_{3})\mathcal{L}_{\ell}(\hat{\bm{k}}_{3}\cdot\hat{\bm{n}})$
			$\displaystyle\,\times\,\sum_{\beta}\frac{b_{\beta}}{\Delta_{\beta}}\left[\Theta^{a^{\prime}}(k_{1})\Theta^{b^{\prime}}(k_{2})\Theta^{c^{\prime}}(k_{3})L_{\ell^{\prime}}(\hat{\bm{k}}_{3}\cdot\hat{\bm{n}})+\text{5 perms.}\right],$

where $\beta\equiv\{a^{\prime},b^{\prime},c^{\prime},\ell^{\prime}\}$. Assuming non-overlapping bins, the integral will be non-zero only when $\{a^{\prime},b^{\prime},c^{\prime}\}$ is some permutation of $\{a,b,c\}$ (again restricting to $a^{\prime}\leq b^{\prime}\leq c^{\prime}$). Invoking global rotational invariance, we can average over the LoS, making use of the relation:

$\displaystyle\int\frac{d\hat{\bm{n}}}{4\pi}\mathcal{L}_{\ell}(\hat{\bm{k}}_{i}\cdot\hat{\bm{n}})\mathcal{L}_{\ell^{\prime}}(\hat{\bm{k}}_{j}\cdot\hat{\bm{n}})=\frac{\delta_{\rm K}^{\ell\ell^{\prime}}}{2\ell+1}\mathcal{L}_{\ell}(\hat{\bm{k}}_{i}\cdot\hat{\bm{k}}_{j}).$

(4.4)

Writing out the permutations explicitly, this gives

	$\displaystyle\left\langle{\hat{B}^{abc}_{\ell}}\right\rangle$	$\displaystyle=$	$\displaystyle\frac{1}{N_{T}^{abc}}\frac{b_{\alpha}}{\Delta_{\alpha}}\int_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}(2\pi)^{3}\delta_{\mathrm{D}}\left({\bm{k}_{123}}\right)\Theta^{a}(k_{1})\Theta^{b}(k_{2})\Theta^{c}(k_{3})$
			$\displaystyle\,\times\,\left\{\left[\mathcal{L}_{\ell}(\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{3})\left[\delta_{\rm K}^{bb^{\prime}}\delta_{\rm K}^{ca^{\prime}}+\delta_{\rm K}^{ba^{\prime}}\delta_{\rm K}^{cb^{\prime}}\right]\delta_{\rm K}^{ac^{\prime}}+\mathcal{L}_{\ell}(\hat{\bm{k}}_{2}\cdot\hat{\bm{k}}_{3})\left[\delta_{\rm K}^{aa^{\prime}}\delta_{\rm K}^{cb^{\prime}}\delta_{\rm K}^{ab^{\prime}}\delta_{\rm K}^{ca^{\prime}}\right]\delta_{\rm K}^{bc^{\prime}}\right.\right.$
			$\displaystyle\qquad\qquad\left.\left.+\,\delta_{\rm K}^{aa^{\prime}}\delta_{\rm K}^{bb^{\prime}}+\delta_{\rm K}^{ab^{\prime}}\delta_{\rm K}^{ba^{\prime}}\right]\delta_{\rm K}^{cc^{\prime}}\right\}.$

The Kronecker deltas demarcate four scenarios: (1) $a\neq b\neq c$, (2) $a=b\neq c$, (3) $a\neq b=c$, (4) $a=b=c$. The latter two are more complex since they involve additional Legendre polynomials of two different $\bm{k}$ vectors. To simplify these, we define the term:

	$\displaystyle N_{\ell}^{abc}$	$\displaystyle\equiv$	$\displaystyle\int_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}(2\pi)^{3}\delta_{\mathrm{D}}\left({\bm{k}_{123}}\right)\Theta^{a}(k_{1})\Theta^{b}(k_{2})\Theta^{c}(k_{3})\mathcal{L}_{\ell}(\hat{\bm{k}}_{2}\cdot\hat{\bm{k}}_{3})$
		$\displaystyle=$	$\displaystyle\frac{4\pi}{2\ell+1}\sum_{m=-\ell}^{\ell}\int d\bm{x}\,\left[\int_{\bm{k}_{1}}e^{-i\bm{k}_{1}\cdot\bm{x}}\Theta^{a}(k_{1})\right]\left[\int_{\bm{k}_{2}}e^{-i\bm{k}_{2}\cdot\bm{x}}\Theta^{c}(k_{2})Y_{\ell m}(\hat{\bm{k}}_{2})\right]$
			$\displaystyle\qquad\,\times\,\left[\int_{\bm{k}_{3}}e^{-i\bm{k}_{3}\cdot\bm{x}}\Theta^{c}(k_{3})Y_{\ell m}^{*}(\hat{\bm{k}}_{3})\right],$

rewriting the Dirac function as an exponential in the second line, allowing expression in terms of Fourier transforms. We note that $N_{0}^{abc}$ is just the isotropic bin volume $N_{T}^{abc}$. With the above definitions, we obtain the desired result $\left\langle{\hat{B}^{abc}_{\ell}}\right\rangle=b_{\alpha}$ (i.e. an unbiased estimator) subject to the following definition:

$\displaystyle\Delta_{\alpha}$

$\displaystyle\equiv$

$\displaystyle\begin{cases}1&a\neq b\neq c\\ 2&a=b\neq c\\ \left(1+N_{\ell}^{abc}/N_{T}^{abc}\right)&a\neq b=c\\ 2\left(1+2N_{\ell}^{abc}/N_{T}^{abc}\right)&a=b=c.\end{cases}$

(4.7)

For $\ell=0$, this reduces to the symmetry factors used in [50] (1 for scalene, 2 for isosceles, 6 for equilateral). This calculation generalizes the standard bispectrum definition (3.3) to the binned bispectrum beyond the narrow bin limit (whence $a\neq b\neq c$ is guaranteed).

We now consider the estimation of bispectrum coefficients $b_{\alpha}$, given their relation to the full bispectrum $B_{\rm ggg}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$. Following Refs. [49, 50, 103], our pathway to this will be:

1. 1.

   Write down the likelihood for the observed pixellized data-minus-randoms field $\bm{d}$ in terms of the pixel correlators $\mathsf{C}_{ij}\equiv\left\langle{d_{i}d_{j}}\right\rangle$, $\mathsf{B}_{ijk}\equiv\left\langle{d_{i}d_{j}d_{k}}\right\rangle$, et cetera, where $i,j,\cdots\in[1,N_{\rm pix}]$ are pixel indices.

2. 2.

   Express the relevant correlator (here $\mathsf{B}_{ijk}$) in terms of the coefficients of interest, i.e. the binned bispectrum multipoles $b_{\alpha}$.

3. 3.

   Maximize the log-likelihood with respect to $b_{\alpha}$ forming a quasi-optimal estimator.

4. 4.

   Simplify the resulting form such that it can be efficiently implemented on data using FFTs.

In the weakly non-Gaussian regime, the likelihood of the data is given by the Edgeworth expansion [e.g., 107]

$\displaystyle-\log L[\bm{d}]=-\log L_{G}[\bm{d}]-\frac{1}{3!}\mathsf{B}^{ijk}\left(h_{i}h_{j}h_{k}-h_{i}\mathsf{C}^{-1}_{jk}-h_{j}\mathsf{C}^{-1}_{ik}-h_{k}\mathsf{C}^{-1}_{ij}\right)+\cdots$

(4.8)

where $L_{G}$ is the Gaussian piece (which we do not need here), and $h_{i}\equiv\mathsf{C}^{-1}_{ij}d^{j}$ is the Wiener-filtered data. In this formalism, the optimal estimator for $b_{\alpha}$ (which enters linearly in $\mathsf{B}^{ijk}$) is given by

$\displaystyle\widehat{b}_{\alpha}=\sum_{\beta}\left(F^{-1}\right)_{\alpha\beta}\widehat{b}_{\beta}^{\rm num},$

(4.9)

defining the numerator and normalization:

	$\displaystyle\widehat{b}_{\alpha}^{\rm num}$	$\displaystyle=$	$\displaystyle\frac{1}{6}\frac{\partial\mathsf{B}^{ijk}}{\partial b_{\alpha}}\left[h_{i}h_{j}h_{k}-\left(h_{i}\mathsf{C}^{-1}_{jk}+\text{2 perms.}\right)\right]$		(4.10)
	$\displaystyle F_{\alpha\beta}$	$\displaystyle=$	$\displaystyle\frac{1}{6}\frac{\partial\mathsf{B}^{ijk}}{\partial b_{\alpha}}\mathsf{C}^{-1}_{il}\mathsf{C}^{-1}_{jm}\mathsf{C}^{-1}_{kn}\frac{\partial\mathsf{B}^{lmn}}{\partial b_{\beta}}.$

This is just the maximum likelihood solution of (4.8).

In our case, the three-point function can be written as a Fourier-transform of the full redshift-space bispectrum $B_{ggg}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$, noting that $d_{i}\equiv n(\bm{r}_{i})\delta_{g}(\bm{r}_{i})$ for background density $n(\bm{r})$:

$\displaystyle\mathsf{B}^{ijk}=n(\bm{r}_{i})n(\bm{r}_{j})n(\bm{r}_{k})\int_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}e^{i\bm{k}_{1}\cdot\bm{r}_{i}+i\bm{k}_{2}\cdot\bm{r}_{j}+i\bm{k}_{3}\cdot\bm{r}_{k}}(2\pi)^{3}\delta_{\mathrm{D}}\left({\bm{k}_{123}}\right)B_{ggg}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3}),$

(4.11)

Inserting (4.2), we can write the cumulant derivative as

	$\displaystyle\frac{\partial\mathsf{B}^{ijk}}{\partial b_{\alpha}}$	$\displaystyle=$	$\displaystyle\frac{n(\bm{r}_{i})n(\bm{r}_{j})n(\bm{r}_{k})}{\Delta_{\alpha}}\int_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}\left[\Theta^{a}(k_{1})\Theta^{b}(k_{2})\Theta^{c}(k_{3})\mathcal{L}_{\ell}(\hat{\bm{k}}_{3}\cdot\hat{\bm{n}})+\text{5 perms.}\right]$
			$\displaystyle\qquad\qquad\qquad\qquad\,\times\,e^{i\bm{k}_{1}\cdot\bm{r}_{i}+i\bm{k}_{2}\cdot\bm{r}_{j}+i\bm{k}_{3}\cdot\bm{r}_{k}}(2\pi)^{3}\delta_{\mathrm{D}}\left({\bm{k}_{123}}\right).$

Under the Yamamoto approximation, we fix the LoS to be $\hat{\bm{n}}=\hat{\bm{r}}_{3}$, as above.

Inserting the above results into (4.10), the numerator of the bispectrum estimator is found to be:

$\displaystyle\widehat{b}^{\rm num}_{\alpha}$

$\displaystyle=$

$\displaystyle\frac{1}{\Delta_{\alpha}}\int d\bm{r}\,\bigg{[}g^{a}_{0}[\bm{d}](\bm{r})g^{b}_{0}[\bm{d}](\bm{r})g^{c}_{\ell}[\bm{d}](\bm{r})-\left(g^{a}_{0}[\bm{d}](\bm{r})\left\langle{g^{b}_{0}[\bm{a}](\bm{r})\tilde{g}^{c}_{\ell}[\bm{a}](\bm{r})}\right\rangle+\text{2 perms.}\right)\bigg{]},$

(4.13)

subject to the definitions

	$\displaystyle g^{a}_{\ell}[\bm{y}](\bm{r})=\int_{\bm{k}}e^{-i\bm{k}\cdot\bm{r}}\Theta^{a}(k)\int d\bm{r}^{\prime}e^{i\bm{k}\cdot\bm{r}^{\prime}}n(\bm{r}^{\prime})[\mathsf{H}^{-1}\bm{y}](\bm{r}^{\prime})\mathcal{L}_{\ell}(\hat{\bm{k}}\cdot\hat{\bm{r}}^{\prime})$		(4.14)
	$\displaystyle\tilde{g}^{a}_{\ell}[\bm{y}](\bm{r})=\int_{\bm{k}}e^{-i\bm{k}\cdot\bm{r}}\Theta^{a}(k)\int d\bm{r}^{\prime}e^{i\bm{k}\cdot\bm{r}^{\prime}}n(\bm{r}^{\prime})[\mathsf{A}^{-1}\bm{y}](\bm{r}^{\prime})\mathcal{L}_{\ell}(\hat{\bm{k}}\cdot\hat{\bm{r}}^{\prime}).$

$g_{0}^{a}$ is equal to the $g^{a}$ function of Ref. [50]. This is closely linked to the $F_{\ell}$ functions found in the ideal estimator (3.8), but now includes the survey mask and custom weighting functions. Two points are of note: (a) we replace the $\mathsf{C}^{-1}$ Wiener filtering by a more general weighting $\mathsf{H}^{-1}$; (b) we introduce a set of random maps $\bm{a}$ with known covariance $\mathsf{A}$ following Ref. [108]. The former allows for a simple-to-implement estimator (since the full pixel covariance is difficult to compute and harder still to invert), and the latter allows one to compute the one-point terms via Monte Carlo summation (removing the need for a direct sum which has a prohibitive $\mathcal{O}(N_{\rm pix}^{2})$ scaling).

Exploiting spherical harmonic factorizations, the two terms in (4.14) can be written in terms of forward and reverse Fourier-transforms $\mathcal{F}$ and $\mathcal{F}^{-1}$:

	$\displaystyle g_{\ell}^{a}[\bm{y}](\bm{r})$	$\displaystyle\equiv$	$\displaystyle\frac{4\pi}{2\ell+1}\sum_{m=-\ell}^{\ell}\mathcal{F}^{-1}\left[{\Theta^{a}(k)Y^{*}_{\ell m}(\hat{\bm{k}})\mathcal{F}\left[{n\mathsf{H}^{-1}\bm{y}\,Y_{\ell m}}\right](\bm{k})}\right](\bm{r})$		(4.15)
	$\displaystyle\tilde{g}_{\ell}^{a}[\bm{y}](\bm{r})$	$\displaystyle\equiv$	$\displaystyle\frac{4\pi}{2\ell+1}\sum_{m=-\ell}^{\ell}\mathcal{F}^{-1}\left[{\Theta^{a}(k)Y^{*}_{\ell m}(\hat{\bm{k}})\mathcal{F}\left[{n\mathsf{A}^{-1}\bm{y}\,Y_{\ell m}}\right](\bm{k})}\right](\bm{r}).$

The second part of the estimator is a data-independent normalization (or Fisher) matrix, $F_{\alpha\beta}$. This acts to remove correlations between bins and multipoles and can be efficiently estimated via Monte Carlo methods. In the limit of ideal weighting ($\mathsf{H}^{-1}\to\mathsf{C}^{-1}$) and vanishing non-Gaussianity, the bispectrum covariance is equal to $F^{-1}$. As in [50], this takes the form

$\displaystyle F_{\alpha\beta}=\frac{1}{12}\left(\left\langle{\phi_{\alpha}^{i}\mathsf{H}^{-1}_{il}\tilde{\phi}_{\beta}^{l}}\right\rangle-\left\langle{\phi_{\alpha}^{i}}\right\rangle\mathsf{H}^{-1}_{il}\left\langle{\tilde{\phi}_{\beta}^{l}}\right\rangle\right),$

(4.16)

with $\phi_{\alpha}^{i}[\bm{a}]=\mathsf{B}^{ijk}_{,\alpha}\mathsf{H}^{-1}_{jj^{\prime}}\mathsf{H}^{-1}_{kk^{\prime}}a^{j^{\prime}}a^{k^{\prime}}$ and analogously for $\tilde{\phi}$ with $\mathsf{H}^{-1}\to\mathsf{A}^{-1}$. (4.16) can be implemented by applying the linear map $\mathsf{H}^{-1}$ to $\tilde{\phi}$ then summing the result (multiplied by $\phi$) in pixel-space. Once again, the expectations can be computed by summation over Monte Carlo realizations $\bm{a}$ with known covariance $\mathsf{A}$ (e.g., Gaussian random fields).

With the above form for the cumulant derivative (4.3), we can write the $\phi$ field explicitly in terms of Fourier transforms:

	$\displaystyle\phi_{\alpha}^{i}[\bm{a}]$	$\displaystyle=$	$\displaystyle\mathsf{B}_{,\alpha}^{ijk}\mathsf{H}^{-1}_{jj^{\prime}}\mathsf{H}^{-1}_{kk^{\prime}}a^{j^{\prime}}a^{k^{\prime}}$
		$\displaystyle=$	$\displaystyle\frac{n(\bm{r}_{i})}{\Delta_{\alpha}}\int d\bm{r}\,\int_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}e^{i\bm{k}_{1}\cdot\bm{r}_{i}}\left[\Theta^{a}(k_{1})\Theta^{b}(k_{2})\Theta^{c}(k_{3})\mathcal{L}_{\ell}(\hat{\bm{k}}_{3}\cdot\hat{\bm{n}})+\text{5\,perms.}\right]$
			$\displaystyle\qquad\qquad\,\times\,e^{-i(\bm{k}_{123})\cdot\bm{r}}[n\mathsf{H}^{-1}\bm{a}](\bm{k}_{2})[n\mathsf{H}^{-1}\bm{a}](\bm{k}_{3})$
		$\displaystyle=$	$\displaystyle\frac{2n(\bm{r}_{i})}{\Delta_{\alpha}}\left\{\mathcal{F}^{-1}\left[{\Theta^{a}(k)\mathcal{F}\left[{g_{0}^{b}[\bm{a}]g_{\ell}^{c}[\bm{a}]}\right](\bm{k})}\right](\bm{r}_{i})+(a\leftrightarrow b)\right.$
			$\displaystyle\qquad\qquad\left.\,+\,\mathcal{F}^{-1}\left[{\Theta^{c}(k)\mathcal{L}_{\ell}(\hat{\bm{k}}\cdot\hat{\bm{r}}_{i})\mathcal{F}\left[{g_{0}^{a}[\bm{a}]g_{0}^{b}[\bm{a}]}\right](\bm{k})}\right](\bm{r}_{i})\right\},$

with an analogous form for $\tilde{\phi}_{\alpha}$ involving $\tilde{g}_{\ell}^{n}$. The final term involves a Legendre polynomial; using spherical harmonic decompositions, this can be simplified to yield the form:

$\displaystyle\left.\phi_{\alpha}^{i}[\bm{a}]\right|_{\rm III}=\frac{2n(\bm{r}_{i})}{\Delta_{\alpha}}\frac{4\pi}{2\ell+1}\sum_{m=-\ell}^{\ell}Y_{\ell m}(\hat{\bm{r}}_{i})\mathcal{F}^{-1}\left[{\Theta^{c}(k)Y_{\ell m}^{*}(\hat{\bm{k}})\mathcal{F}\left[{g_{0}^{a}[\bm{a}]g_{0}^{b}[\bm{a}]}\right](\bm{k})}\right](\bm{r}_{i}).$

(4.18)

Collecting results, the full estimator for the bispectrum is given by

$\displaystyle\widehat{b}_{\alpha}=\sum_{\beta}F_{\alpha\beta}^{-1}\widehat{b}^{\rm num}_{\beta}.$

(4.19)

This is unbiased for any choice of $\mathsf{H}^{-1}$, unwindowed, and, for $\mathsf{H}^{-1}\approx\mathsf{C}^{-1}$, close-to optimal (partly due to the inclusion of a linear term [cf. 108]). These properties are derived formally in [50]. Both the numerator and Fisher matrix can be efficiently computed using $N_{\rm mc}$ Monte Carlo simulations, with the finite number of simulations incurring an error proportional to $\sqrt{1+1/N_{\rm mc}}$. Whilst the latter is computationally expensive (requiring $\mathcal{O}(N_{\rm bins})$ Fourier transforms), it only has to be estimated once for a given survey geometry. We will discuss the specifics of our implementation in §6. A public Python implementation can be found online.⁶⁶6GitHub.com/OliverPhilcox/Spectra-Without-Windows.

To model the galaxy bispectrum models, we will use the tree-level theory introduced in Ref. [15] (see also [2, 9, 109, 110, 111, 46, 47, 112, 113, 114]). At a redshift $z$, the redshift-space bispectrum is the sum of three contributions:

$$B_{\rm ggg}({\bf k}_{1},{\bf k}_{2},{\bf k}_{3})=B_{211}({\bf k}_{1},{\bf k}_{2},{\bf k}_{3})+B_{\rm ctr}({\bf k}_{1},{\bf k}_{2},{\bf k}_{3})+B_{\rm stoch}({\bf k}_{1},{\bf k}_{2},{\bf k}_{3})\,,$$

(5.1)

where $B_{211}$ is the standard deterministic mode-coupling contribution,

$$\begin{split}B_{211}({\bf k}_{1},{\bf k}_{2},{\bf k}_{3})=2Z_{1}({\bf k}_{1})Z_{1}({\bf k}_{2})Z_{2}({\bf k}_{1},{\bf k}_{2})P_{11}(k_{1})P_{11}(k_{2})+\text{2 cyc.}\,.\end{split}$$

(5.2)

Here $P_{11}(k,z)$ is the linear matter power spectrum at redshift $z$, and the redshift-space perturbation theory kernels are given by [cf., 115]

	$\displaystyle Z_{1}({\bf k})=b_{1}+f\mu^{2}\,,$		(5.3a)
	$\displaystyle Z_{2}({\bf k}_{1},{\bf k}_{2})=\frac{b_{2}}{2}+b_{\mathcal{G}_{2}}\left(\frac{({\bf k}_{1}\cdot{\bf k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}-1\right)+b_{1}F_{2}({\bf k}_{1},{\bf k}_{2})+f\mu^{2}G_{2}({\bf k}_{1},{\bf k}_{2})$
	$\displaystyle\qquad\qquad\quad~{}~{}+\frac{f\mu k}{2}\left(\frac{\mu_{1}}{k_{1}}(b_{1}+f\mu_{2}^{2})+\frac{\mu_{2}}{k_{2}}(b_{1}+f\mu_{1}^{2})\right)\,,$		(5.3b)
	$\displaystyle F_{2}({\bf k}_{1},{\bf k}_{2})=\frac{5}{7}+\frac{1}{2}\left(\frac{({\bf k}_{1}\cdot{\bf k}_{2})}{k_{1}^{2}}+\frac{({\bf k}_{1}\cdot{\bf k}_{2})}{k_{2}^{2}}\right)+\frac{2}{7}\frac{({\bf k}_{1}\cdot{\bf k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}\,,$		(5.3c)
	$\displaystyle G_{2}({\bf k}_{1},{\bf k}_{2})=\frac{3}{7}+\frac{1}{2}\left(\frac{({\bf k}_{1}\cdot{\bf k}_{2})}{k_{1}^{2}}+\frac{({\bf k}_{1}\cdot{\bf k}_{2})}{k_{2}^{2}}\right)+\frac{4}{7}\frac{({\bf k}_{1}\cdot{\bf k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}\,,$		(5.3d)

where we introduced the following angles with respect to the line of sight directions: $\mu_{i}\equiv({\bf k}_{i}\cdot\hat{{\bf z}})/k_{i}$ and $\mu\equiv({\bf k}\cdot\hat{{\bf z}})/k$, ${\bf k}\equiv{\bf k}_{1}+{\bf k}_{2}$. Additionally, $f$ is the logarithmic growth factor,

$$f=\frac{d\mathop{\rm ln}\nolimits D_{+}}{d\mathop{\rm ln}\nolimits a}\,,$$

(5.4)

where $D_{+}$ denotes the usual growth rate and $a$ is the scale factor of the Friedmann-Lemaitre-Robertson-Walker metric. The free coefficients $b_{1}$, $b_{2}$, and $b_{\mathcal{G}_{2}}$ capture linear, quadratic, and tidal bias between galaxies and matter [116, 115, 117, 118, 119, 120, 121].