.pspdf.pdfps2pdf -dEPSCrop -dNOSAFER #1

# The galaxy bispectrum in the Spherical Fourier-Bessel basis

Joshua N. Benabou<EMAIL_ADDRESS>Berkeley Center for Theoretical
Physics, University of California, Berkeley, CA 94720, U.S.A.  Theoretical
Physics Group, Lawrence Berkeley National Laboratory, Berkeley, CA 94720,
U.S.A.  Adriano Testa  California Institute of Technology, 1200 E. California
Boulevard, Pasadena, CA 91125, USA  Chen Heinrich  California Institute of
Technology, 1200 E. California Boulevard, Pasadena, CA 91125, USA  Henry S.
Grasshorn Gebhardt  California Institute of Technology, 1200 E. California
Boulevard, Pasadena, CA 91125, USA  Jet Propulsion Laboratory, California
Institute of Technology, Pasadena, California 91109, USA  Olivier Doré
California Institute of Technology, 1200 E. California Boulevard, Pasadena, CA
91125, USA  Jet Propulsion Laboratory, California Institute of Technology,
Pasadena, California 91109, USA

###### Abstract

The bispectrum, the three-point correlation in Fourier space, is a crucial
statistic for studying many effects targeted by the next-generation galaxy
surveys, such as primordial non-Gaussianity (PNG) and general relativistic
(GR) effects on large scales. In this work we develop a formalism for the
bispectrum in the Spherical Fourier-Bessel (SFB) basis – a natural basis for
computing correlation functions on the curved sky, as it diagonalizes the
Laplacian operator in spherical coordinates. Working in the SFB basis allows
for line-of-sight effects such as redshift space distortions (RSD) and GR to
be accounted for exactly, i.e without having to resort to perturbative
expansions to go beyond the plane-parallel approximation. Only analytic
results for the SFB bispectrum exist in the literature given the intensive
computations needed. We numerically calculate the SFB bispectrum for the first
time, enabled by a few techniques: We implement a template decomposition of
the redshift-space kernel $Z_{2}$ into Legendre polynomials, and separately
treat the PNG and velocity-divergence terms. We derive an identity to
integrate a product of three spherical harmonics connected by a Dirac delta
function as a simple sum, and use it to investigate the limit of a homogeneous
and isotropic Universe. Moreover, we present a formalism for convolving the
signal with separable window functions, and use a toy spherically symmetric
window to demonstrate the computation and give insights into the properties of
the observed bispectrum signal. While our implementation remains
computationally challenging, it is a step toward a feasible full extraction of
information on large scales via a SFB bispectrum analysis.

###### Contents

  1. I Introduction
  2. II Background
    1. II.1 The SFB formalism
    2. II.2 Observed galaxy SFB power spectrum
    3. II.3 Observed Fourier galaxy bispectrum
    4. II.4 SFB bispectrum definition
  3. III SFB bispectrum in a homogeneous and isotropic Universe
    1. III.1 Key identity for fast computation
    2. III.2 Properties of the signal
    3. III.3 Requirement on the sampling frequency for resolving the oscillations
  4. IV The observed SFB bispectrum
    1. IV.1 Template decomposition of the observed bispectrum
    2. IV.2 Signal computation
    3. IV.3 Properties of the signal
      1. IV.3.1 Oscillations in $k$ and mode coupling
      2. IV.3.2 Limber suppression at low $k$ and high $l$
      3. IV.3.3 Violation of the triangle condition
  5. V Discussion
  6. VI  Acknowledgements
  7. A Useful formulae
    1. A.1 Dirac delta distribution
    2. A.2 Spherical Bessel functions
    3. A.3 Spherical harmonics
    4. A.4 Gaunt factor
    5. A.5 Wigner symbols
  8. B Encoding of observational isotropy by the Gaunt factor
  9. C An identity for integrating a product of three spherical harmonics
  10. D Details of the bispectrum computation
    1. D.1 Legendre expansion coefficients $Z_{l_{1}l_{2}l_{3}}(r,k_{1},k_{2})$
    2. D.2 Derivation of the observed SFB bispectrum Eq. 67
    3. D.3 Simplification of the kernel $\mathcal{V}^{L_{1}L_{2}\ell}_{M_{1}M_{2}m}(k,k_{1},k_{2})$
    4. D.4 Contribution from the $f_{\mathrm{NL}}$ and $G_{2}$ terms
      1. D.4.1 Relation between position-dependent bispectrum and SFB bispectrum
      2. D.4.2 $G_{2}$ contribution
      3. D.4.3 Total contribution from $f_{\mathrm{NL}}$ and $G_{2}$ terms

##  I Introduction

Current and next-generation large-scale-structure (LSS) surveys such as DESI
Aghamousa  _et al._ (2016), Euclid Amendola  _et al._ (2018), SPHEREx Doré
_et al._ (2015) and the Nancy Grace Roman Space Telescope Eifler  _et al._
(2021) will measure the galaxy density field over increasingly larger angular
scales, enabling us to constrain interesting physical effects that become
important on those scales, such as primordial non-Gaussianity Dalal  _et al._
(2008) (PNG) and general relativistic effects McDonald (2009).

While many of our current techniques for estimators and modeling are well-
suited for small-area surveys, they are challenged in larger surveys due to
the breaking down of previously used approximations on the full sky. In
particular, the plane-parallel approximation, which assumes that each galaxy
has the same line-of-sight, breaks down when the galaxy separation becomes
large in a full-sky survey. Additionally, the Newtonian modeling of galaxy
density also breaks down as general relativistic effects that grow as $1/k$
become important on large scales (for details see Refs. Yoo  _et al._ (2012);
Yoo (2010a); Bonvin and Durrer (2011); Challinor and Lewis (2011)).

More precisely, redshift space distortions (RSD) induce effects in the
observed galaxy density field that depend on the line-of-sight (LOS) of
individual galaxies. Estimators assuming a fixed LOS for the entire survey
will inevitably lose information at large galaxy separations. Even if one uses
Yamamoto-like estimators Yamamoto  _et al._ (2006); Bianchi  _et al._ (2015);
Scoccimarro (2015) which assume a fixed LOS for each galaxy pair or triplet,
there is still loss of information as the galaxies could have large angular
separations in a given pair or triplet. The signal picked up by the Yamamoto
estimator also includes wide-angle effects that are usually modeled either
perturbatively as an expansion in the angular separation of the galaxy pair
(i.e an expansion whose zero-order term is the plane-parallel approximation)
Reimberg  _et al._ (2016); Beutler  _et al._ (2019); Benabou  _et al._ (prep);
Pardede  _et al._ (2023), or non-perturbatively via an exact calculation in
the correlation function space Di Dio  _et al._ (2016).

This raises the question of whether the Fourier basis is the optimal basis to
use on the full sky. Indeed, the Fourier basis consists of the eigenfunctions
of the Laplacian in Cartesian coordinates; the spherical Fourier-Bessel (SFB)
basis, consisting of the eigenfunctions of the Laplacian in spherical
coordinates, is a more natural basis for data analysis on the curved sky. The
SFB basis was proposed for studying galaxy surveys since the early 90’s Lahav
(1993); Binney and Quinn (1991), and was applied to data in the context of a
power spectrum analysis in Refs. Heavens and Taylor (1995); Tadros  _et al._
(1999); Percival  _et al._ (2004).

Recently, an important limitation of the SFB analysis has been overcome in
Ref. Samushia (2019), rendering computations of the power spectrum much more
feasible: a boundary condition at the lower end of the redshift range was
introduced, avoiding the need to carry many modes to model vanishing power
outside of the survey footprint, which can introduce numerical instabilities.
Later, the authors of Ref. Gebhardt and Doré (2021) developed a SFB power
spectrum estimator with a public code release, which builds on this
improvement as well as those in Ref. Wang  _et al._ (2020) on pixel window
effects and the separation of angular and radial transforms, making a SFB
analysis feasible for surveys measuring the power spectrum such as Nancy Grace
Roman, SPHEREx and Euclid (see also Grasshorn Gebhardt and Doré (2022, 2023)).

An alternative to the SFB basis called tomographic spherical harmonics (TSH)
has also been explored in the literature, where the galaxy density contrast in
a redshift bin is decomposed into spherical harmonics, and many redshift bins
are used. In the limit of thin bins, neighboring bins are highly correlated,
and the covariance matrix could become nearly degenerate. For thick bins, one
loses information about the radial modes that are smaller than the bin size.
SFB modes, in contrast, are more efficient basis functions since the radial
modes are captured by spherical Bessel functions which are orthogonal to each
other, unlike in the case of the redshift bin decomposition. See Ref. Zhang
_et al._ (2021) for a detailed analysis comparing the SFB and TSH power
spectrum (at $\Delta z=0.1$) for current and future surveys, showing better
$f_{\mathrm{NL}}$ constraints in general for the SFB method.

Limited effort, however, has been dedicated to the study of the SFB
bispectrum. The bispectrum is the $3$-point correlation function in Fourier
space and is of great importance to next-generation surveys. It is shown to be
powerful at breaking parameter degeneracies when combined with the power
spectrum for constraining galaxy bias parameters, neutrino masses and
primordial non-Gaussianities (see e.g. Gualdi and Verde (2020); Heinrich  _et
al._ (2023); Coulton  _et al._ (2019); Ruggeri  _et al._ (2018); Doré  _et
al._ (2015)); the odd-parity bispectrum is also a smoking-gun signature for
general relativistic effects that become more important on large scales
Maartens  _et al._ (2020); Jeong and Schmidt (2020); Clarkson  _et al._
(2019).

A comprehensive derivation of the SFB bispectrum including all first and
second-order GR effects, geometric effects and PNG was achieved in Ref.
Bertacca  _et al._ (2018). However, due to the complexity of the computations
involved, there has not yet been any work numerically evaluating the SFB
bispectrum signal. In fact, while most of the integrals involved in
calculating the signal are three-dimensional and are doable, some of the most
important and interesting ones involving RSD and PNG contributions are four-
dimensional and are intractable to compute naively.

In this paper, we derive a mathematical identity to express the six
dimensional angular integral of three spherical harmonics connected by a
Dirac-delta function as a simple sum, and use it to study the bispectrum
signal in a homogeneous and isotropic Universe, to build intuition for later
understanding the observed bispectrum in a realistic Universe. We also use
this identity to accelerate the computation of RSD and PNG terms contributing
to the observed bispectrum signal. Furthermore, we employ a template
decomposition of the second-order coupling kernel in redshift space $Z_{2}$
into products of Legendre polynomials to evaluate all three-dimensional
integrals. These techniques allow us to calculate and visualize for the first
time the SFB bispectrum signal.

We apply a general formalism we develop for convolution with a separable
window function (in the angular and radial direction) to the toy example of a
spherically symmetric window function to obtain numerical results that we
study in detail. We derive key insights into the properties of the observed
SFB bispectrum in a realistic Universe, highlighting those due to geometric
effects.

The structure of this paper is as follows. In Section II we review the SFB
basis and the modeling of the SFB galaxy power spectrum; we also describe the
modeling of the Fourier space tree-level galaxy bispectrum and define the SFB
bispectrum. Then in Section III we explore the calculation of the SFB
bispectrum in the simplest case of a homogeneous and isotropic Universe,
building up key intuition for interpreting the features of the observed
bispectrum in the next section. In Section IV, we incorporate various
observational effects into the bispectrum, including growth of structure,
galaxy bias, RSD, PNG and the survey window function. We present our template
decomposition technique to enable its calculation, deferring the details of
the derivation to the appendices; we then visualize and analyze the behavior
of the observed bispectrum signal. Finally, we conclude and discuss future
work in Section V.

##  II Background

In this section, we begin by reviewing the SFB formalism and the SFB power
spectrum following Ref. Gebhardt and Doré (2021), to which we refer the
readers for more details. We then review the Fourier space bispectrum, and
describe the modeling of the observed galaxy bispectrum in redshift space
including local PNG. Finally, we define the SFB bispectrum, which we later
compute in Sections III and IV.

###  II.1 The SFB formalism

The spherical Fourier-Bessel mode $\delta_{\ell m}(k)$ of the density contrast
field $\delta({\bm{r}})$ is defined by

$\displaystyle\delta_{\ell m}(k)$
$\displaystyle\equiv\int\mathrm{d}^{3}{\bm{r}}\left[\sqrt{\frac{2}{\pi}}\,k\,j_{\ell}(kr)\,Y^{*}_{\ell
m}({\hat{\bm{r}}})\right]\delta({\bm{r}})\,.$ (1)

The inverse transform is then

$\displaystyle\delta({\bm{r}})$ $\displaystyle=\int\mathrm{d}k\,\sum_{\ell
m}\left[\sqrt{\frac{2}{\pi}}\,k\,j_{\ell}(kr)\,Y_{\ell
m}({\hat{\bm{r}}})\right]\delta_{\ell m}(k)\,,$ (2)

where ${\bm{r}}=r{\hat{\bm{r}}}$ is the position vector, $r$ is the comoving
distance from the origin, and ${\hat{\bm{r}}}$ is the line-of-sight direction.

The spherical Fourier-Bessel modes are related to the Fourier modes via

$\displaystyle\delta_{\ell
m}(k)=\frac{k}{(2\pi)^{\frac{3}{2}}}\,i^{\ell}\int\mathrm{d}^{2}{\hat{\bm{k}}}\,Y^{*}_{\ell
m}({\hat{\bm{k}}})\,\delta({\bm{k}})\,,$ (3)

for which the inverse relation is

$\delta({\bm{k}})=\frac{(2\pi)^{\frac{3}{2}}}{k}\sum_{\ell
m}i^{-\ell}\,Y_{\ell m}({\hat{\bm{k}}})\,\delta_{\ell m}(k)\,.$ (4)

Note that in this paper we use the following convention for the Fourier
transform:

$\displaystyle f({\bm{k}})$
$\displaystyle=\int\mathrm{d}^{3}r\,e^{-i{\bm{k}}\cdot{\bm{r}}}\,f({\bm{r}})\,,$
(5) $\displaystyle f({\bm{r}})$
$\displaystyle=\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\,e^{i{\bm{k}}\cdot{\bm{r}}}\,f({\bm{k}})\,.$
(6)

If unambiguous, we use the same symbol in configuration space (e.g,
$f({\bm{r}})$) as in Fourier space (e.g, $f({\bm{k}})$).

###  II.2 Observed galaxy SFB power spectrum

Let us denote the expansion of the matter density in cosmological perturbation
theory by

$\delta({\bm{k}},r)=\sum_{n=1}^{\infty}D^{n}(r)\delta^{(n)}({\bm{k}})\,,$ (7)

with $D(r)$ the growth factor. Then, in the linear regime, the observed galaxy
density contrast to first order can be modeled as

$\displaystyle\delta^{g,\mathrm{obs},(1)}({\bm{r}})$
$\displaystyle=W({\bm{r}})\,D(r)\int\frac{\mathrm{d}^{3}{\bm{q}}}{(2\pi)^{3}}\,e^{i{\bm{q}}\cdot{\bm{r}}}\,\widetilde{A}_{\mathrm{RSD}}(\mu,q\mu,r)b(r,q)\,\delta^{(1)}({\bm{q}})\,,$
(8)

where $W({\bm{r}})$ is the survey window, $b(r,q)$ is the linear galaxy bias,
$\mu={\hat{\bm{r}}}\cdot{\hat{\bm{q}}}$, and $\delta^{(1)}({\bm{q}})$ is the
matter density contrast in Fourier space. In what follows, the matter density
contrast will always be denoted $\delta$ without any superscript, while the
galaxy density contrast is denoted by $\delta^{g}$.

RSD effects are contained in $\widetilde{A}_{\mathrm{RSD}}(\mu,q\mu,r)$, which
can be modeled as

$\displaystyle\widetilde{A}_{\mathrm{RSD}}(\mu,q\mu,r)=\left(1+\beta\mu^{2}\right)\widetilde{A}_{\mathrm{FoG}}(q\mu)\,,$
(9)

where $\beta=f/b$ and $f=\mathrm{d}\ln D/\mathrm{d}\ln a$ (with $a$ the scale
factor) is the linear growth rate. In this paper we ignore the Fingers-of-God
effect and set $\widetilde{A}_{\mathrm{FoG}}(q\mu)=1$.

Transforming to spherical Fourier-Bessel space we have that

$\displaystyle\delta^{g,\mathrm{obs},(1)}_{\ell m}(k)$
$\displaystyle=\int\mathrm{d}q\sum_{LM}\mathcal{W}_{\ell
m}^{LM}(k,q)\,\delta_{LM}^{(1)}(q)\,,$ (10)

where the observed galaxy density $\delta^{g,\mathrm{obs}}_{\ell m}(k)$ is
related to the matter density $\delta_{LM}(q)$ by the mode coupling matrix
$\mathcal{W}_{\ell m}^{LM}(k,q)$. In our convention, this mode coupling matrix
encodes galaxy physics such as galaxy bias and RSD effects, unequal time
effects such as the growth of structure, and the survey window function
$W(\mathbf{r})$:

$\displaystyle\mathcal{W}_{\ell m}^{LM}(k,q)$
$\displaystyle=\int\mathrm{d}^{2}{\hat{\bm{r}}}\,Y_{LM}({\hat{\bm{r}}})\,Y^{*}_{\ell
m}({\hat{\bm{r}}})\,\mathcal{W}_{\ell}^{L}(k,q,{\hat{\bm{r}}})\,,$ (11)

where

$\displaystyle\mathcal{W}_{\ell}^{L}(k,q,{\hat{\bm{r}}})$
$\displaystyle=\frac{2qk}{\pi}\int\mathrm{d}r\,r^{2}\,W({\bm{r}})\,D(r)\,b(r,q)\,j_{\ell}(kr)\widetilde{A}_{\mathrm{RSD}}(-i\partial_{qr},-iq\partial_{qr},r)j_{L}(qr)\,,$
(12)

where we replace the argument $\mu$ of $\widetilde{A}_{\rm RSD}$ by
$-i\partial_{qr}$ which acts on $e^{i{\bm{q}}\cdot{\bm{r}}}=e^{iqr\mu}$.

Noting that in a homogeneous and isotropic Universe, the matter power spectrum
satisfies

$\langle\delta(\bm{k})\delta^{*}\\!\\!\left(\bm{k}^{\prime}\right)\rangle=(2\pi)^{3}\delta^{D}\left(\bm{k}-\bm{k}^{\prime}\right)P(k)\,,$
(13)

it follows that the 2-point function of the SFB modes is

$\displaystyle\left\langle\delta^{g,\mathrm{obs}}_{\ell
m}(k)\,\delta^{g,\mathrm{obs},*}_{\ell^{\prime}m^{\prime}}(k^{\prime})\right\rangle=\int\mathrm{d}q\sum_{LM}\mathcal{W}_{\ell
m}^{LM}(k,q)\,\mathcal{W}_{\ell^{\prime}m^{\prime}}^{LM,*}(k^{\prime},q)\,P(q)\,.$
(14)

In the full-sky limit where $W(\mathbf{r})=W(r)$, we have that
$\mathcal{W}^{L}_{\ell}(k,q,{\hat{\bm{r}}})\rightarrow\mathcal{W}^{\ell}_{\ell}(k,q)$
is independent of ${\hat{\bm{r}}}$. Let us define
$\mathcal{W}_{\ell}(k,q)\equiv\mathcal{W}^{\ell}_{\ell}(k,q)$. Then the SFB
power spectrum $C_{\ell}(k,k^{\prime})$ defined via

$\displaystyle\left\langle\delta^{g,\mathrm{obs}}_{\ell
m}(k)\,\delta^{g,\mathrm{obs},*}_{\ell^{\prime}m^{\prime}}(k^{\prime})\right\rangle$
$\displaystyle=\delta^{K}_{\ell\ell^{\prime}}\delta^{K}_{mm^{\prime}}\,C_{\ell}(k,k^{\prime})\,,$
(15)

can be expressed as

$\displaystyle C_{\ell}(k,k^{\prime})$
$\displaystyle=\int\mathrm{d}q\,\mathcal{W}_{\ell}(k,q)\,\mathcal{W}^{*}_{\ell}(k^{\prime},q)\,P(q)\,,$
(16)

where

$\displaystyle\mathcal{W}_{\ell}(k,q)=\frac{2qk}{\pi}\int\mathrm{d}r\,r^{2}\,W(r)\,D(r)\,j_{\ell}(kr)\bigg{(}b(r,q)j_{\ell}(qr)-f(r)j_{\ell}^{{}^{\prime\prime}}(qr)\bigg{)},$
(17)

where we use Eq. 9 with $\mu\to-i\partial_{qr}$.

Note that in a homogeneous and isotropic Universe, for which
$b(r,q)=D(r)=\widetilde{A}_{\mathrm{RSD}}=W({\bm{r}})=1$,
$\mathcal{W}_{\ell}(k,q)$ becomes a Dirac delta function and we have that
$C_{l}(k,k^{\prime})=\delta_{D}(k-k^{\prime})P(k)$. In reality, the kernels
$\mathcal{W}_{\ell}(k,q)$ are peaked at $k\approx q$. We show examples of
$\mathcal{W}_{\ell}(k,q)$ for the spherical window
$W(r)=\mathbf{1}_{[0,r_{\mathrm{max}}]}(r)$ for various values of $r_{\rm
max}$ in Fig. 1, where we fix $k=4.18\times 10^{-2}$ and vary $q/k$ for
$\ell=20$. Here and in the remainder of this work we use the Planck 2018
cosmology Aghanim  _et al._ (2020) as our fiducial cosmology. The matter power
spectrum at zero redshift and the linear growth factors $f$ and $D$ are
computed with the Boltzmann code CAMB111https://camb.info/ Lewis  _et al._
(2001). All other calculations are performed in Julia Bezanson  _et al._
(2017).

Figure 1: $\mathcal{W}_{\ell}(k,q)$ for fixed $\ell=20$, $k=4.18\times
10^{-2}\,h\mbox{Mpc}^{-1}$ and for various sizes $r_{\mathrm{max}}$ of the
survey window $W(r)=\mathbf{1}_{[0,r_{\mathrm{max}}]}(r)$.

###  II.3 Observed Fourier galaxy bispectrum

We now review the observed galaxy bispectrum in Fourier space including
observational effects such as the RSD and galaxy bias, but without the window
function convolution. For details on the derivation of the various quantities,
we refer the readers to Ref. Tellarini  _et al._ (2016) and Bernardeau  _et
al._ (2002) which we follow closely. As we will also be concerned with
modeling the effects of PNG in the SFB bispectrum, we will include its effects
in the Fourier bispectrum as well.

We consider PNG of the local type, for which the fluctuations of the potential
are parameterized by

$\displaystyle\Phi_{\mathrm{NG}}(\bm{x})=\varphi(\bm{x})+f_{\mathrm{NL}}\left(\varphi^{2}(\bm{x})-\left\langle\varphi^{2}\right\rangle\right)\,,$
(18)

where $\varphi$ is a primordial Gaussian potential. Using the Poisson
equation, we may relate the long-wavelength Gaussian potential to the linearly
evolved primordial matter density perturbation via

$\displaystyle\Phi_{\mathrm{NG}}(\bm{k})=\frac{\delta^{(1)}(\bm{k},z)}{\alpha(k,z)}\,,$
(19)

where

$\displaystyle\alpha(k,z)=\frac{2k^{2}c^{2}D(z)T(k)}{3H_{0}^{2}\Omega_{\mathrm{m}}}\,.$
(20)

Here $\Omega_{m}$ is the matter density, $H_{0}$ is the Hubble constant, and
$T(k)$ is the transfer function of matter perturbations, normalized to $1$ at
low $k$. Eq. 19-20 are valid in the Newtonian limit on subhorizon scales Wands
and Slosar (2009); Yoo  _et al._ (2009); Yoo (2010b).

In perturbation theory, the observed galaxy density contrast field at position
${\bm{r}}$ is given by

$\displaystyle\delta^{g}({\bm{k}},{\bm{r}})$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}D^{n}(r)\int\frac{\mathrm{d}^{3}{\bm{k}}_{1}}{(2\pi)^{3}}\ldots\int\frac{\mathrm{d}^{3}{\bm{k}}_{n}}{(2\pi)^{3}}\,(2\pi)^{3}\delta_{D}({\bm{k}}_{1}+\dots+{\bm{k}}_{n}-{\bm{k}})\,Z_{n}({\bm{k}}_{1},\ldots,{\bm{k}}_{n},{\bm{r}})\,\delta^{(1)}({\bm{k}}_{1})\cdots\delta^{(1)}({\bm{k}}_{n})\,,$
(21)

where $\delta^{(1)}$ is the linear matter density field, and the $n$-th order
redshift space kernels $Z_{n}$ encode the mode coupling effects from
gravitational evolution, PNG and galaxy biasing. We assume the bivariate
galaxy biasing model

$\displaystyle\delta^{\mathrm{g}}(\bm{x})$
$\displaystyle=b^{\mathrm{E}}_{10}\delta(\bm{x})+b^{\mathrm{E}}_{01}\varphi(\bm{x})+b^{\mathrm{E}}_{20}\left(\delta(\bm{x})\right)^{2}+b^{\mathrm{E}}_{11}\delta(\bm{x})\varphi(\bm{x})+b^{\mathrm{E}}_{02}(\varphi(\bm{x}))^{2}+b_{s_{2}}\left(s^{2}-\left\langle
s^{2}\right\rangle\right)-b^{\mathrm{E}}_{01}n^{2}\,.$ (22)

Above we define the tidal term Catelan  _et al._ (2000); Baldauf  _et al._
(2012)

$\displaystyle
s^{2}(\bm{k})=\int\frac{d\bm{q}}{(2\pi)^{3}}S_{2}(\bm{q},\bm{k}-\bm{q})\delta^{(1)}(\bm{q})\delta^{(1)}(\bm{k}-\bm{q}),$
(23)

and the (non-Gaussian) term encoding displacement of galaxies from their
initial Lagrangian coordinate positions $\bm{q}$

$\displaystyle
n^{2}(\bm{k})=2\int\frac{d\bm{q}}{(2\pi)^{3}}N_{2}(\bm{q},\bm{k}-\bm{q})\frac{\delta^{(1)}(\bm{q})\delta^{(1)}(\bm{k}-\bm{q})}{\alpha(|\bm{k}-\bm{q}|)}\,,$
(24)

where above we use the kernels

$\displaystyle N_{2}\left({\bm{k}}_{1},{\bm{k}}_{2}\right)$
$\displaystyle=\frac{{\bm{k}}_{1}\cdot{\bm{k}}_{2}}{2k_{1}^{2}}\,,$ (25)
$\displaystyle S_{2}\left({\bm{k}}_{1},{\bm{k}}_{2}\right)$
$\displaystyle=\frac{\left({\bm{k}}_{1}\cdot{\bm{k}}_{2}\right)^{2}}{k_{1}^{2}k_{2}^{2}}-\frac{1}{3}\,.$
(26)

The redshift space kernels at first and second order, respectively, are given
by

$\displaystyle Z_{1}({\bm{k}},{\bm{r}})$
$\displaystyle=b^{\mathrm{E}}_{10}+f(r)\mu^{2}+\frac{b^{\mathrm{E}}_{01}}{\alpha(k)},$
(27) $\displaystyle Z_{2}({\bm{k}}_{1},{\bm{k}}_{2},{\bm{r}})$
$\displaystyle=b^{\mathrm{E}}_{10}\left[F_{2}\left({\bm{k}}_{1},{\bm{k}}_{2}\right)+f_{\mathrm{NL}}\frac{\alpha(k)}{\alpha\left(k_{1}\right)\alpha\left(k_{2}\right)}\right]+\left[b^{\mathrm{E}}_{20}-\frac{2}{7}b_{10}^{\mathrm{L}}S_{2}\left({\bm{k}}_{1},{\bm{k}}_{2}\right)\right]$
$\displaystyle+\frac{b^{\mathrm{E}}_{11}}{2}\left[\frac{1}{\alpha\left(k_{1}\right)}+\frac{1}{\alpha\left(k_{2}\right)}\right]+\frac{b^{\mathrm{E}}_{02}}{\alpha\left(k_{1}\right)\alpha\left(k_{2}\right)}-b^{\mathrm{E}}_{01}\left[\frac{N_{2}\left({\bm{k}}_{1},{\bm{k}}_{2}\right)}{\alpha\left(k_{2}\right)}+\frac{N_{2}\left({\bm{k}}_{2},{\bm{k}}_{1}\right)}{\alpha\left(k_{1}\right)}\right]$
$\displaystyle+f(r)\mu^{2}\left[G_{2}\left({\bm{k}}_{1},{\bm{k}}_{2}\right)+f_{\mathrm{NL}}\frac{\alpha(k)}{\alpha\left(k_{1}\right)\alpha\left(k_{2}\right)}\right]+\frac{f(r)^{2}k^{2}\mu^{2}}{2}\frac{\mu_{1}\mu_{2}}{k_{1}k_{2}}+b^{\mathrm{E}}_{10}\frac{f(r)\mu
k}{2}\left(\frac{\mu_{1}}{k_{1}}+\frac{\mu_{2}}{k_{2}}\right)$
$\displaystyle+b^{\mathrm{E}}_{01}\frac{f(r)\mu
k}{2}\left[\frac{\mu_{1}}{k_{1}\alpha\left(k_{2}\right)}+\frac{\mu_{2}}{k_{2}\alpha\left(k_{1}\right)}\right],$
(28)

where $\mu\equiv{\hat{\bm{k}}}\cdot{{\hat{\bm{r}}}}$,
${\bm{k}}\equiv{\bm{k}}_{1}+{\bm{k}}_{2}$,
$\mu_{i}\equiv{\hat{\bm{k}}}_{i}\cdot{{\hat{\bm{r}}}}$, and where the coupling
kernels for the real-space density and velocity-divergence fields are

$\displaystyle F_{2}({\bm{q}}_{1},{\bm{q}}_{2})$
$\displaystyle=\frac{5}{7}+\frac{1}{2}\bigg{(}\frac{q_{1}}{q_{2}}+\frac{q_{2}}{q_{1}}\bigg{)}\frac{{\bm{q}}_{1}\cdot{\bm{q}}_{2}}{q_{1}q_{2}}+\frac{2}{7}\frac{({\bm{q}}_{1}\cdot{\bm{q}}_{2})^{2}}{(q_{1}q_{2})^{2}}\,,$
(29) $\displaystyle G_{2}({\bm{q}}_{1},{\bm{q}}_{2})$
$\displaystyle=\frac{3}{7}+\frac{1}{2}\bigg{(}\frac{q_{1}}{q_{2}}+\frac{q_{2}}{q_{1}}\bigg{)}\frac{{\bm{q}}_{1}\cdot{\bm{q}}_{2}}{q_{1}q_{2}}+\frac{4}{7}\frac{({\bm{q}}_{1}\cdot{\bm{q}}_{2})^{2}}{(q_{1}q_{2})^{2}}\,.$
(30)

From Eq. 27 it follows that we may write the linear galaxy bias appearing in
Eq. 17 as $b(r,q)=b^{\mathrm{E}}_{10}+b^{\mathrm{E}}_{01}/\alpha(q,r)$.

The galaxy bispectrum in Fourier space is defined via

$\displaystyle\left\langle\delta^{g}({\bm{k}}_{1},{\bm{r}}_{1})\delta^{g}({\bm{k}}_{2},{\bm{r}}_{2})\delta^{g}({\bm{k}}_{3},{\bm{r}}_{3})\right\rangle=B_{s}({\bm{k}}_{1},{\bm{k}}_{2},{\bm{k}}_{3},{\bm{r}}_{1},{\bm{r}}_{2},{\bm{r}}_{3})(2\pi)^{3}\delta_{D}({\bm{k}}_{1}+{\bm{k}}_{2}+{\bm{k}}_{3})\,.$
(31)

Working up to second order in the galaxy density field expansion, the tree-
level bispectrum is

$\displaystyle
B_{s}({\bm{k}}_{1},{\bm{k}}_{2},{\bm{k}}_{3},{\bm{r}}_{1},{\bm{r}}_{2},{\bm{r}}_{3})\equiv
2D(r_{1})D(r_{2})D^{2}(r_{3})P(k_{1})P(k_{2})Z_{1}({\bm{k}}_{1},{\bm{r}}_{1})Z_{1}({\bm{k}}_{2},{\bm{r}}_{2})Z_{2}({\bm{k}}_{1},{\bm{k}}_{2},{\bm{r}}_{3})+\mathrm{2\
cyc.\ perm.}\,,$ (32)

where we sum over all cyclic permutations of the subscripts of the quantities
in parentheses. Note that in the absence of RSD, linear growth, galaxy bias,
and PNG, Eq. 32 reduces to the matter bispectrum

$\displaystyle B_{m}({\bm{k}}_{1},{\bm{k}}_{2},{\bm{k}}_{3})\equiv
2P(k_{1})P(k_{2})F_{2}({\bm{k}}_{1},{\bm{k}}_{2})+\mathrm{2\ cyc.\
perm.}\,\,.$ (33)

We follow Ref. Tellarini  _et al._ (2016) to model the Eulerian biases

$\displaystyle b_{10}^{\mathrm{E}}=1+b_{10}^{\mathrm{L}}\,,$ (34)
$\displaystyle b_{01}^{\mathrm{E}}=b_{01}^{\mathrm{L}}\,,$ (35) $\displaystyle
b_{20}^{\mathrm{E}}=\frac{8}{21}b_{10}^{\mathrm{L}}+b_{20}^{\mathrm{L}}\,,$
(36) $\displaystyle
b_{11}^{\mathrm{E}}=b_{01}^{\mathrm{L}}+b_{11}^{\mathrm{L}}\,,$ (37)
$\displaystyle b_{02}^{\mathrm{E}}=b_{02}^{\mathrm{L}}\,,$ (38) in terms of
the Lagrangian biases, which are given by $\displaystyle b_{01}^{\mathrm{L}}$
$\displaystyle=2f_{\mathrm{NL}}\delta_{c}b_{10}^{\mathrm{L}},$ (39)
$\displaystyle b_{11}^{\mathrm{L}}$
$\displaystyle=2f_{\mathrm{NL}}\left(\delta_{c}b_{20}^{\mathrm{L}}-b_{10}^{\mathrm{L}}\right),$
(40) $\displaystyle b_{02}^{\mathrm{L}}$
$\displaystyle=4f_{\mathrm{NL}}^{2}\delta_{c}\left(\delta_{c}b_{20}^{\mathrm{L}}-2b_{10}^{\mathrm{L}}\right)\,,$
(41)

if one assumes a Universal Mass Function, and where $\delta_{c}$ is the
critical overdensity for halo collapse, here set to its value for spherical
collapse $\delta_{c}=1.686$.

Note that only $b_{10}^{E}$ and $b_{20}^{E}$ need to be specified in order to
determine all the other bias parameters. Specifically for our SFB bispectrum
computation later, we will set $b_{10}^{E}=1.8$ and $b_{20}^{E}=0.305$. While
there is no technical obstacle to including redshift-dependent biases in the
SFB calculation, we choose flat biases here for simplicity.

###  II.4 SFB bispectrum definition

We now review the formalism for the SFB bispectrum. We seek to compute the
$3$-point correlation function of the observed galaxy over-density field in
SFB space:

$\displaystyle\left\langle\delta^{g,\mathrm{obs}}_{l_{1}m_{1}}(k_{1})\,\delta^{g,\mathrm{obs}}_{l_{2}m_{2}}(k_{2})\delta^{g,\mathrm{obs}}_{l_{3}m_{3}}(k_{3})\right\rangle\,.$
(42)

In the following, it will be useful to distinguish between two notions of
isotropy, which we term _observational_ isotropy and _intrinsic_ isotropy.
Intrinsic isotropy refers to the statistically isotropic distribution of
galaxies on the largest-scales in real-space. Due to RSD, the galaxy
clustering observed in surveys is not intrinsically isotropic since, in
redshift space, it depends on the angle to a given LOS. On the other hand, the
distribution observed by a full-sky survey remains invariant under rotations
about the observer position. This observational isotropy is only broken by a
survey window which is not spherically symmetric. We show in Appendix B that,
assuming observational isotropy, Eq. 42 is real and proportional to the Gaunt
factor

$\displaystyle\mathcal{G}^{l_{1}l_{2}l_{3}}_{m_{1}m_{2}m_{3}}$
$\displaystyle\equiv\int\mathrm{d}^{2}{\hat{\bm{r}}}\,Y_{l_{1}m_{1}}({\hat{\bm{r}}})\,Y_{l_{2}m_{2}}({\hat{\bm{r}}})\,Y_{l_{3}m_{3}}({\hat{\bm{r}}})\,,$
(43)

which can be expressed in terms of Wigner-$3j$ symbols,

$\displaystyle\mathcal{G}^{l_{1}l_{2}l_{3}}_{m_{1}m_{2}m_{3}}$
$\displaystyle=\left(\frac{(2l_{1}+1)(2l_{2}+1)(2l_{3}+1)}{4\pi}\right)^{\frac{1}{2}}\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
0&0&0\end{pmatrix}\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}\,.$ (44)

The Wigner-$3j$’s ensure the SFB $3$-point function vanishes unless the
following conditions are satisfied: (i) $m_{1}+m_{2}+m_{3}=0$, (ii) triangle
inequality on the $l_{i}$: $l_{i}\geq l_{j}-l_{k}$, and (iii)
$l_{1}+l_{2}+l_{3}$ is even.

In order to rid of the purely geometric information contained in the $m_{i}$,
we compute the “angle-averaged” bispectrum

$\displaystyle B^{\text{SFB}}_{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})$
$\displaystyle\equiv\sum_{m_{1},m_{2},m_{3}}\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}\left\langle\delta^{g,\mathrm{obs}}_{l_{1}m_{1}}(k_{1})\delta^{g,\mathrm{obs}}_{l_{2}m_{2}}(k_{2})\delta^{g,\mathrm{obs}}_{l_{3}m_{3}}(k_{3})\right\rangle\,.$
(45)

Using the orthogonality relation in Eq. 94 then gives

$\displaystyle\left\langle\delta^{g,\mathrm{obs}}_{l_{1}m_{1}}(k_{1})\delta^{g,\mathrm{obs}}_{l_{2}m_{2}}(k_{2})\delta^{g,\mathrm{obs}}_{l_{3}m_{3}}(k_{3})\right\rangle=\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}B^{\text{SFB}}_{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})\,.$
(46)

In the following subsections we will always plot the dimensionless reduced
bispectrum

$\displaystyle
Q^{\text{SFB}}_{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})\equiv\frac{B^{\text{SFB}}_{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})}{P(k_{1})P(k_{2})+P(k_{1})P(k_{3})+P(k_{2})P(k_{3})}\,,$
(47)

which partially projects out the dependence of the signal on $k_{i}$ coming
from the matter power spectrum. Finally, we note that the bispectrum is
invariant under simultaneous cyclic permutations of $(l_{1},l_{2},l_{3})$ and
$(k_{1},k_{2},k_{3})$, which allows us to restrict to $l_{1}\leq l_{2}\leq
l_{3}$.

Before delving into the computation of the SFB bispectrum, let us briefly
remark on its relation to the angular bispectrum in spherical shells (i.e, the
TSH bispectrum) $b_{l_{1}l_{2}l_{3}}(r_{1},r_{2},r_{3})$ of Ref. Di Dio  _et
al._ (2016), defined via

$\displaystyle\mathcal{G}^{l_{1}l_{2}l_{3}}_{m_{1}m_{2}m_{3}}b_{l_{1}l_{2}l_{3}}(r_{1},r_{2},r_{3})\equiv\int
d^{2}{\hat{\bm{r}}}_{1}d^{2}{\hat{\bm{r}}}_{2}d^{2}{\hat{\bm{r}}}_{3}Y^{*}_{l_{1},m_{1}}({\hat{\bm{r}}}_{1})Y^{*}_{l_{2},m_{2}}({\hat{\bm{r}}}_{2})Y^{*}_{l_{3},m_{3}}({\hat{\bm{r}}}_{3})\left\langle\delta^{g,\mathrm{obs}}({\bm{r}}_{1})\delta^{g,\mathrm{obs}}({\bm{r}}_{2})\delta^{g,\mathrm{obs}}({\bm{r}}_{3})\right\rangle\,.$
(48)

Using Eq. 46 in combination with Eq. 1, it follows that

$\displaystyle B^{\text{SFB}}_{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})$
$\displaystyle=\bigg{(}\frac{2}{\pi}\bigg{)}^{\frac{3}{2}}k_{1}k_{2}k_{3}\int\left(\prod_{i}\mathrm{d}r_{i}r_{i}^{2}j_{l_{i}}(k_{i}r_{i})\right)b_{l_{1}l_{2}l_{3}}(r_{1},r_{2},r_{3})\,.$
(49)

We see that the SFB bispectrum and the TSH bispectrum are related by an
invertible linear transformation, and that the multipole indices $l_{i}$ are
the same in the SFB and TSH formalisms (the wavenumbers $k_{i}$ are the same
between SFB and Fourier space). In practice, in the TSH formalism many radial
bins are desirable to fully exploit the large scale radial modes (see Zhang
_et al._ (2021), which studied this for the power spectrum), leading to a
covariance matrix which is difficult to invert, whereas this issue does not
arise in the SFB basis.

##  III SFB bispectrum in a homogeneous and isotropic Universe

###  III.1 Key identity for fast computation

We first examine the bispectrum in the limit of a homogeneous and
intrinsically isotropic Universe (by ignoring the growth of structure, galaxy
bias evolution, redshift-space distortions and window function effects) in
order to study its features and build up intuition for understanding the
observed SFB bispectrum later in Section IV.

We can relate the SFB bispectrum to the Fourier bispectrum using the relation
between the SFB and Fourier modes in Eq. 3:

$\displaystyle\left\langle\delta_{l_{1}m_{1}}(k_{1})\,\delta_{l_{2}m_{2}}(k_{2})\delta_{l_{3}m_{3}}(k_{3})\right\rangle=\frac{k_{1}k_{2}k_{3}}{(2\pi)^{\frac{9}{2}}}i^{l_{1}+l_{2}+l_{3}}\int
d^{2}{\hat{\bm{k}}}_{1}d^{2}{\hat{\bm{k}}}_{2}d^{2}{\hat{\bm{k}}}_{3}Y^{*}_{l_{1},m_{1}}({\hat{\bm{k}}}_{1})Y^{*}_{l_{2},m_{2}}({\hat{\bm{k}}}_{2})Y^{*}_{l_{3},m_{3}}({\hat{\bm{k}}}_{3})\left\langle\delta({\bm{k}}_{1})\,\delta({\bm{k}}_{2})\delta({\bm{k}}_{3})\right\rangle\,.$
(50)

Due to homogeneity and isotropy, the Fourier bispectrum
$B_{m}(k_{1},k_{2},k_{3})$ (Eq. 33) depends only on the lengths $k_{1},k_{2}$
and $k_{3}$, so that we may write

$\displaystyle\left\langle\delta_{l_{1}m_{1}}(k_{1})\,\delta_{l_{2}m_{2}}(k_{2})\delta_{l_{3}m_{3}}(k_{3})\right\rangle$
$\displaystyle=\frac{k_{1}k_{2}k_{3}}{(2\pi)^{\frac{9}{2}}}i^{l_{1}+l_{2}+l_{3}}(2\pi)^{3}B_{m}(k_{1},k_{2},k_{3})\;I_{m_{1}m_{2}m_{3}}^{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3}),$
(51)

where

$I_{m_{1}m_{2}m_{3}}^{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})\equiv\int
d^{2}{\hat{\bm{k}}}_{1}d^{2}{\hat{\bm{k}}}_{2}d^{2}{\hat{\bm{k}}}_{3}Y^{*}_{l_{1},m_{1}}({\hat{\bm{k}}}_{1})Y^{*}_{l_{2},m_{2}}({\hat{\bm{k}}}_{2})Y^{*}_{l_{3},m_{3}}({\hat{\bm{k}}}_{3})\delta_{D}({\bm{k}}_{1}+{\bm{k}}_{2}+{\bm{k}}_{3}).$
(52)

Eq. 52 has typically been written in terms of an integral over the spherical
Bessel functions Bertacca  _et al._ (2018) (for details, see appendix C) as

$\displaystyle
I_{m_{1}m_{2}m_{3}}^{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})=8i^{l_{1}+l_{2}+l_{3}}\int
r^{2}\mathrm{d}rj_{l_{1}}(k_{1}r)j_{l_{2}}(k_{2}r)j_{l_{3}}(k_{3}r)\int
d^{2}{\hat{\bm{r}}}Y^{*}_{l_{1},m_{1}}({\hat{\bm{r}}})Y^{*}_{l_{2},m_{2}}({\hat{\bm{r}}})Y^{*}_{l_{3},m_{3}}({\hat{\bm{r}}})\,.$
(53)

where the second integral is the Gaunt factor (Eq. 43). Here we instead derive
the identity (see derivation in appendix C)

$\displaystyle I_{m_{1}m_{2}m_{3}}^{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})$
$\displaystyle=\frac{4\pi^{\frac{3}{2}}}{k_{1}k_{2}k_{3}}\sqrt{(2l_{1}+1)}\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}\sum_{|m|\leq\mathrm{min}(l_{2},l_{3})}\\!\\!\\!\\!\\!{Y_{l_{2},m}}(\theta_{12},0)Y_{l_{3},-m}(\theta_{13},0)(-1)^{m}\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
0&m&-m\end{pmatrix}\,,$ (54)

which allows us to rapidly compute the angle-averaged bispectrum without any
numerical integration222For this purpose we precompute a lookup table of
$Y_{lm}(\theta,0)$ values and interpolate. Also note that we may halve the
number of terms in the sum by using its invariance under $m\to-m$.

$\displaystyle B^{\text{SFB, iso/homo}}_{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})$
$\displaystyle=B_{m}(k_{1},k_{2},k_{3})i^{l_{1}+l_{2}+l_{3}}\sqrt{2(2l_{1}+1)}$
$\displaystyle\times\sum_{|m|\leq\mathrm{min}(l_{2},l_{3})}Y_{l_{2},m}(\theta_{12},0)Y_{l_{3},-m}(\theta_{13},0)(-1)^{m}\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
0&m&-m\end{pmatrix}\,.$ (55)

Above we define $\theta_{ij}$ as the angle between ${\bm{k}}_{i}$ and
${\bm{k}}_{j}$, such that
$\cos(\theta_{12})\equiv{\hat{\bm{k}}}_{1}\cdot{\hat{\bm{k}}}_{2}=\vartheta(k_{1},k_{2},k_{3})$,
where
$\vartheta(k_{1},k_{2},k_{3})\equiv\frac{k_{3}^{2}-k_{1}^{2}-k_{2}^{2}}{2k_{1}k_{2}}$,
and $\cos(\theta_{13})\equiv\vartheta(k_{1},k_{3},k_{2})$. Furthermore,
denoting the spherical coordinates of ${\hat{\bm{r}}}$ by $(\theta,\phi)$, we
define $Y_{\ell,m}({\hat{\bm{r}}})\equiv Y_{\ell,m}(\theta,\phi)$. It is clear
from Eq. 55 that the SFB bispectrum in an isotropic and homogeneous Universe
is proportional to the Fourier bispectrum by a geometric coupling factor
depending on the $l_{i}$ and the angles between the ${\bm{k}}_{i}$. This
factor also imposes the triangle inequality on the wavenumbers, i.e,
$k_{i}\leq|k_{j}-k_{k}|$, which is relaxed as we shall see in the next section
for the observed SFB bispectrum and is only imposed approximately for
spherically symmetric surveys which extend to sufficiently large redshifts.

The identity Eq. 54 is one of our key results. In addition to trivializing the
computation of the signal in a homogeneous and isotropic Universe, it provides
analytic insight into the geometric features of the observed bispectrum.
Crucially, we will employ this identity to render the computation of the
observed bispectrum tractable. We will discuss these points in IV.2.

Figure 2: The reduced SFB bispectrum signal for an isotropic and homogeneous
Universe, as a function of $k_{2}$ and $k_{3}$ in units of $k_{1}$, which is
fixed to $k_{1}=4.18\times 10^{-2}\ h\mbox{Mpc}^{-1}$ here. Each panel
displays a different triplet $l_{i}\equiv(l_{1},l_{2},l_{3})$. The bispectrum
vanishes identically for configurations $(k_{1},k_{2},k_{3})$ which do not
satisfy the triangle condition. The oscillations are a result of the geometric
coupling in Eq. 56; their number is controlled by the values of $l$. Note that
the colorbar limits are saturated in each panel.

###  III.2 Properties of the signal

In Fig. 2, we show two-dimensional cross-sections of the reduced bispectrum in
an isotropic and homogeneous Universe as a function of $k_{2}/k_{1}$ and
$k_{3}/k_{1}$ for fixed $k_{1}$, for three $l$-triplets $(l_{1},l_{2},l_{3})$.
The most striking feature is the rectangular border outside of which the
signal vanishes; this is the enforcement of the triangle inequality.

Another important feature is that the signal oscillates in the space of
$k_{i}$’s, which is not surprising, given that the products of spherical
harmonics in $I^{l_{1}l_{2}l_{3}}_{m_{1}m_{2}m_{3}}(k_{1},k_{2},k_{3})$,

$\displaystyle Y_{l_{2},m}(\theta_{12},0)Y_{l_{3},-m}(\theta_{13},0)\,,$ (56)

oscillate as the angles between the ${\bm{k}}_{i}$’s vary. Further, the number
of oscillations as one moves from the center of the plot corresponding to an
equilateral $k$-triangle, toward the borders of the rectangular region,
corresponding to degenerate triangles, is higher for larger $l_{i}$ values.

Figure 3: The reduced bispectrum signal in an isotropic and homogeneous
Universe for fixed $k_{1}=4.18\times 10^{-2}\ h\mbox{Mpc}^{-1}$ as in Fig. 2,
but now taking a cross-section along the diagonal $k_{2}=k_{3}$. We show two
equilateral $l$-triplets $l_{1}=l_{2}=l_{3}=l$ with $l=10$ and 30. The plot is
cut at $k_{2}/k_{1}=0.5$ on the left since there is no signal below it where
the triangle condition is violated (this property does not hold for the
observed bispectrum in a realistic Universe, however). At high $k_{2}$, the
SFB bispectrum is proportional to the Fourier-space bispectrum.

We also show a one-dimensional cross section in Fig. 3, taking the diagonal
$k_{2}=k_{3}$, for two different equilateral $l$ shapes
($l_{1}=l_{2}=l_{3}=l$). Perhaps the most important feature to note in this
plot is that the SFB bispectrum effectively reduces to the matter bispectrum
at the limit $k_{2}=k_{3}\gg k_{1}$ in the homogeneous and isotropic Universe.
In this limit, we have $\cos(\theta_{12})=\cos(\theta_{13})\to 0$ and
$B^{\mathrm{iso/homo}}_{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})$ reduces to
$C_{l_{1}l_{2}l_{3}}B_{m}(k_{1},k_{2},k_{3})$ for a prefactor
$C_{l_{1}l_{2}l_{3}}$, which quickly tends towards a constant as $l$ grows.

Similarly, in the limit of degenerate triangles e.g. $k_{3}=k_{2}+k_{1}$, we
have $\cos(\theta_{12})=\cos(\theta_{13})=1$, such that only the $m=0$ term in
Eq. 55 is nonzero, and the bispectrum reduces to

$\displaystyle B^{\text{SFB,
iso/homo}}_{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})=B_{m}(k_{1},k_{2},k_{3})i^{l_{1}+l_{2}+l_{3}}\sqrt{2(2l_{1}+1)}\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
0&0&0\end{pmatrix}\,,$ (57)

where again we have that the cosmological signal and the geometric coupling
separate into a $k$-dependent and an $l$-dependent piece.

###  III.3 Requirement on the sampling frequency for resolving the
oscillations

As we will see in the next section, the observed SFB bispectrum signal has a
similar oscillation pattern in the space of $k$’s as in the isotropic and
homogeneous limit. With the analytic formula Eq. 55 at hand, we can easily
estimate the local frequency of oscillations in $k$-space, and thus the
minimum sampling of $k_{i}$’s required to resolve these oscillations, assuming
that the computation is performed on a uniform cubic grid of
$(k_{1},k_{2},k_{3})$ with spacing $\Delta k$.

For example, consider for a given $k_{1}$, the oscillations along the diagonal
$k_{2}=k_{3}$, as are visible in Fig. 3. Estimating the frequency of the
oscillations amounts to estimating the spacing between the roots of the
associated Legendre polynomials in the products of $Y_{lm}$’s in Eq. 56. On
the diagonal, the lowest point for which the signal is nonzero corresponds to
the degenerate isosceles triangle $k_{2}=k_{3}=k_{1}/2$, and the diagonal
extends to the top right into a squeezed triangle where $k_{2},k_{3}\gg
k_{1}$. On this trajectory, $\theta_{12}$ and $\theta_{13}$ vary from 0 to
$\pi/2$. The associated Legendre polynomial $P^{l}_{m}(\cos(\theta))$ has
$l-|m|$ roots on the range $0<\theta<\pi$, which are symmetric about $\pi/2$,
so there are $(l-|m|)/2$ roots on the ranges we consider.

Consequently, the product of $Y_{l_{2}m}Y_{l_{3}m}$ crosses zero at most
$(l_{2}+l_{3}-2|m|)/2$ times333If $l_{2}=l_{3}$, it crosses zero $\sim
l_{2}/2$ times.. The $m=0$ term has the highest number of roots and may be
used to estimate an upper bound on the spacing of roots. Towards that purpose,
note that $P^{0}_{l}(\cos(\theta))$ is simply the Legendre polynomial
$\mathcal{L}_{l}(\cos(\theta))$. Let $\theta_{1},\dots,\theta_{l}$ be the
sequence of roots of $\mathcal{L}_{l}(\cos(\theta))$ in the interval
$(0,\pi)$, listed in increasing order. Then we have the inequalities on the
location of the roots Szego (1936)

$\displaystyle\frac{\nu-\frac{1}{2}}{l}\pi<\theta_{\nu}<\frac{\nu}{l+1}\pi\quad(\nu=1,2,\cdots,\left\lceil
l/2\right\rceil)\,.$ (58)

Hence, the first value of $k_{2}$ above $k_{2,\rm{min}}=k_{1}/2$ for which the
spherical harmonic product of index $m=0$ vanishes along the diagonal
satisfies

$k_{2}^{\rm root\,\nu=1}<\frac{k_{1}}{2\cos(\pi/(l_{3}+1))},$ (59)

Thus, to have at least $N$ sampling points per oscillation, the sampling
$\Delta k$ must satisfy

$\displaystyle N\Delta k<k_{2}^{\rm
root\,\nu=1}-k_{2,\rm{min}}<k_{1}\left[\frac{1}{2\cos(\frac{\pi}{l_{3}+1})}-\frac{1}{2}\right]\;\xrightarrow{l_{3}\gg\pi}\;\frac{k_{1}}{4}\left(\frac{\pi}{l_{3}}\right)^{2}.$
(60)

We have verified numerically that the estimated $l_{3}^{-2}$ scaling holds for
bispectrum signal. Given the cost of computing the observed bispectrum, this
means that resolving the oscillations of the signal within the triangle
inequality region is challenging at high $l$.

In Section IV.3, we shall see that one property of the observed bispectrum
signal is that for large enough $l$ it is “Limber-suppressed” at low-$k$, in
particular inside the region where the triangle inequality is satisfied. This
means that we actually do not need to resolve the signal close to the borders
of this region, where the local frequency of oscillations is higher, as the
signal contains comparatively little information there.

##  IV The observed SFB bispectrum

In this section, we begin by describing the template decomposition we use in
order to render the computation of the observed SFB bispectrum feasible. We
then give details of the signal computation before studying its properties.

Note that we now incorporate redshift evolution, RSD effects, PNG and survey
window effects. Statistical homogeneity is now broken by the growth of
structure and RSD444Note that in Fourier space, one can still assume
statistical homogeneity by restricting to a given redshift bin and choosing an
effective redshift for the entire bin., and the latter also breaks intrinsic
isotropy. We restrict ourselves to the linear regime $k\lesssim
0.1\,h\mbox{Mpc}^{-1}$ for which the tree-level bispectrum remains valid down
to $z=0$, and therefore do not include Fingers-of-God effects. We also choose
to not model the monopole and dipole here, as they are affected by observer
terms in GR such as the observer potential and peculiar velocity Bertacca  _et
al._ (2018), i.e we restrict ourselves to multipoles $l_{i}\geq 2$.

###  IV.1 Template decomposition of the observed bispectrum

We begin by expressing the observed galaxy density field by applying the
window function to the galaxy density field in redshift space,

$\displaystyle\delta^{g,\mathrm{obs}}({\bm{r}})$
$\displaystyle=W({\bm{r}})\int\frac{\mathrm{d}^{3}{\bm{q}}}{(2\pi)^{3}}\,e^{i{\bm{q}}\cdot{\bm{r}}}\delta^{g}({\bm{q}},{\bm{r}})\,.$
(61)

To second order in the linear matter density field $\delta^{(1)}$, we have
from Eq. 21 that
$\delta^{g,\mathrm{obs}}({\bm{r}})=\delta^{g,\mathrm{obs},(1)}({\bm{r}})+\delta^{g,\mathrm{obs},(2)}({\bm{r}})$,
where

$\displaystyle\delta^{g,\mathrm{obs},(1)}({\bm{r}})$
$\displaystyle=W({\bm{r}})\int\frac{\mathrm{d}^{3}{\bm{q}}}{(2\pi)^{3}}\,e^{i{\bm{q}}\cdot{\bm{r}}}\,D(r)\,Z_{1}(\bm{q},{\bm{r}})\,\delta^{(1)}(\bm{q})\,,$
(62) $\displaystyle\delta^{g,\mathrm{obs},(2)}({\bm{r}})$
$\displaystyle=W({\bm{r}})\int\frac{\mathrm{d}^{3}{\bm{q}}}{(2\pi)^{3}}\,e^{i{\bm{q}}\cdot{\bm{r}}}\,D^{2}(r)\int\frac{1}{(2\pi)^{3}}\mathrm{d}^{3}{\bm{k}}_{1}\,\mathrm{d}^{3}{\bm{k}}_{2}\,Z_{2}({\bm{k}}_{1},{\bm{k}}_{2},{\bm{r}})\,\delta^{(1)}({\bm{k}}_{1})\delta^{(1)}({\bm{k}}_{2})\,\delta_{D}({\bm{k}}_{1}+{\bm{k}}_{2}-{\bm{q}})\,.$
(63)

Transforming the linear density contrast Eq. (62) into SFB space, we retrieve
Eq. (10) where the kernel $\mathcal{W}_{\ell m}^{LM}(k,q)$ encodes (linear)
galaxy physics and RSD. We aim to derive a similar relation for the second-
order density contrast. We now transform Eq. 63 into SFB space using Eq. 1.
Expressing the linear matter density contrast in the SFB basis from Eq. 4, and
writing the Dirac-delta as an integral over complex exponentials, we obtain

$\displaystyle\delta_{\ell m}^{g,(2)}(k)$
$\displaystyle=\sqrt{\frac{2}{\pi}}\,k\int\mathrm{d}^{3}{\bm{r}}\,j_{\ell}(kr)\,Y^{*}_{\ell
m}({\hat{\bm{r}}})\,W({\bm{r}})\,D^{2}(r)\int\mathrm{d}^{3}{\bm{q}}\,e^{i{\bm{q}}\cdot{\bm{r}}}\int\frac{\mathrm{d}^{3}{\bm{k}}_{1}}{(2\pi)^{3}}\int\frac{\mathrm{d}^{3}{\bm{k}}_{2}}{(2\pi)^{3}}\,Z_{2}({\bm{k}}_{1},{\bm{k}}_{2},{\bm{r}})\int\mathrm{d}^{3}{\bm{x}}\,e^{i{\bm{k}}_{1}\cdot{\bm{x}}}\,e^{i{\bm{k}}_{2}\cdot{\bm{x}}}\,e^{-i{\bm{q}}\cdot{\bm{x}}}$
$\displaystyle\quad\times\frac{1}{k_{1}}\sum_{L_{1}M_{1}}i^{-L_{1}}\,Y_{L_{1}M_{1}}({\hat{\bm{k}}}_{1})\,\delta^{(1)}_{L_{1}M_{1}}(k_{1})\frac{1}{k_{2}}\sum_{L_{2}M_{2}}i^{-L_{2}}\,Y_{L_{2}M_{2}}({\hat{\bm{k}}}_{2})\,\delta^{(1)}_{L_{2}M_{2}}(k_{2})\,.$
(64)

Naively inserting the expression for $Z_{2}$ (Eq. 28) into Eq. 64 would
require evaluating high-dimensional angular integrals, which is intractable.
To simplify the calculation, we remark that
$Z_{2}({\bm{k}}_{1},{\bm{k}}_{2},{\bm{r}})$ is nearly a polynomial in
${\hat{\bm{k}}}_{1}\cdot{\hat{\bm{k}}}_{2}$,
${\hat{\bm{k}}}_{1}\cdot{\hat{\bm{r}}}$, and
${\hat{\bm{k}}}_{2}\cdot{\hat{\bm{r}}}$. Indeed, defining $\tilde{Z}_{2}$ such
that

$\displaystyle Z_{2}({\bm{k}}_{1},{\bm{k}}_{2},{\bm{r}})$
$\displaystyle=\tilde{Z}_{2}({\bm{k}}_{1},{\bm{k}}_{2},{\bm{r}})+f_{\mathrm{NL}}\frac{\alpha(k)}{\alpha\left(k_{1}\right)\alpha\left(k_{2}\right)}\left(b^{\mathrm{E}}_{10}+f(r)\mu^{2}\right)+f(r)\mu^{2}G_{2}\left({\bm{k}}_{1},{\bm{k}}_{2}\right)\,,$
(65)

we can decompose $\tilde{Z}_{2}$ into Legendre polynomials in those three
variables and thereby factorize the dependence on the ${\hat{\bm{k}}}_{i}$ and
${\hat{\bm{r}}}$:

$\displaystyle\tilde{Z}_{2}({\bm{k}}_{1},{\bm{k}}_{2},{\bm{r}})$
$\displaystyle=\sum_{l_{1}l_{2}l_{3}}Z_{l_{1}l_{2}l_{3}}(k_{1},k_{2},r)\,\mathcal{L}_{l_{1}}({\hat{\bm{k}}}_{1}\cdot{\hat{\bm{k}}}_{2})\,\mathcal{L}_{l_{2}}({\hat{\bm{k}}}_{1}\cdot{\hat{\bm{r}}})\,\mathcal{L}_{l_{3}}({\hat{\bm{k}}}_{2}\cdot{\hat{\bm{r}}})\,.$
(66)

Importantly, the sum over $l_{1},l_{2},l_{3}$ is finite, and indexed by $9$
triplets $(l_{1},l_{2},l_{3})$ whose corresponding coefficients
$Z_{l_{1}l_{2}l_{3}}$ are listed in Section D.1. As we discuss below, the
above Legendre decomposition permits to reduce the bispectrum to a triple
integral.

Two other terms remain. The term proportional to $f_{\mathrm{NL}}$ in Eq. 65
depends on $k=|{\bm{k}}_{1}+{\bm{k}}_{2}|$ through $\alpha(k)$, so it cannot
be decomposed it in a similar fashion. In principle, the term proportional to
$G_{2}({\bm{k}}_{1},{\bm{k}}_{2})$ in Eq. 65 can be decomposed in this manner;
however, since
$\mu^{2}=({\bm{k}}_{1}\cdot{\hat{\bm{r}}}+{\bm{k}}_{2}\cdot{\hat{\bm{r}}})/(k_{1}^{2}+k_{2}^{2}+2{\bm{k}}_{1}\cdot{\bm{k}}_{2})$
it would render the sum Eq. 66 infinite and slowly-converging. Hence, we
choose to treat the $G_{2}$ and $f_{\mathrm{NL}}$ terms separately (and
exactly), as described in Section D.4.

We summarize the remainder of the derivation here, leaving details to Appendix
D: First, we insert the decomposition Eq. 66 of
$\tilde{Z}_{2}({\bm{k}}_{1},{\bm{k}}_{2},{\bm{r}})$ into Eq. 64, and use the
plane-wave expansion Eq. 85 to decompose the complex exponentials into
spherical harmonics and spherical Bessel functions. We rid of the angular
integrals over ${\hat{\bm{q}}}$ with the orthogonality relation for spherical
harmonics Eq. 80.

Then we apply Wick’s theorem to compute $\left\langle\delta_{\ell
m}^{g,(2)}(k)\,\delta^{g,(1)}_{\ell^{\prime}m^{\prime}}(k^{\prime})\,\delta^{g,(1)}_{\ell^{\prime\prime}m^{\prime\prime}}(k^{\prime\prime})\right\rangle$
in terms of the two-point functions. Finally, proceeding under the assumption
of a spherically symmetric window $W({\bm{r}})=W(r)$, we compute the angle
averaged bispectrum using Eq. 45 and obtain:

$\displaystyle B^{\text{SFB}}_{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})$
$\displaystyle=2\int\mathrm{d}q_{2}\mathcal{W}_{l_{2}}(k_{2},q_{2})\,P(q_{2})\int\mathrm{d}q_{3}\,\mathcal{W}_{l_{3}}(k_{3},q_{3})\,P(q_{3})\
\mathcal{V}_{\mathrm{tot.}}^{l_{1}l_{2}l_{3}}(k_{1},q_{2},q_{3})+\mathrm{2\
cyc.\ perm.}\,,$ (67)

where

$\displaystyle\mathcal{V}_{\mathrm{tot.}}^{l_{1}l_{2}l_{3}}(k_{1},q_{2},q_{3})$
$\displaystyle\equiv\mathcal{V}_{\mathrm{}}^{l_{1}l_{2}l_{3}}(k_{1},q_{2},q_{3})+\mathcal{V}_{f_{\mathrm{NL},G_{2}}}^{l_{1}l_{2}l_{3}}(k_{1},q_{2},q_{3})\,,$
(68)

where the specific forms of $\mathcal{V}_{\mathrm{}}^{l_{1}l_{2}l_{3}}$ and
$\mathcal{V}_{f_{\mathrm{NL},G_{2}}}^{l_{1}l_{2}l_{3}}$ are given by Eqs. 143
and 155 respectively. Note that the kernel $\mathcal{W}_{l}(k,q)$ is already
given by Eq. 17.

Let us briefly comment on the form of the dimensionless kernel Eq. 68. The
first term of Eq. 68 is given by

$\displaystyle\mathcal{V}^{\ell
L_{1}L_{2}}(k,k_{1},k_{2})\equiv(32\pi)^{\frac{3}{2}}kk_{1}k_{2}\sum_{l_{1}l_{2}l_{3}L_{3}L_{4}}g^{L_{1}L_{2}\ell}_{l_{1}l_{2}l_{3}L_{3}L_{4}}J^{\ell
L_{3}L_{4}}_{l_{1}l_{2}l_{3}}(k,k_{1},k_{2})\,,$ (69)

where

$\displaystyle J^{\ell
L_{3}L_{4}}_{l_{1}l_{2}l_{3}}(k,k_{1},k_{2})\equiv\int\mathrm{d}r\,r^{2}\,j_{\ell}(kr)\,j_{L_{3}}(k_{1}r)\,j_{L_{4}}(k_{2}r)\,W(r)\,D^{2}(r)\,Z_{l_{1}l_{2}l_{3}}(k_{1},k_{2},r)\,,$
(70)

is the contribution to SFB mode coupling by cosmological sources (e.g.
redshift evolution, RSD and PNG) and survey window, and where
$g^{L_{1}L_{2}\ell}_{l_{1}l_{2}l_{3}L_{3}L_{4}}$ is a purely geometric mode
coupling coefficient given by Eq. 142.

In analogy to the SFB power spectrum, in which the matter power spectrum is
convolved with kernels $\mathcal{W}_{\ell}(k,q)$ that describe the mode
coupling through the product of two spherical Bessel functions, the kernel
$\mathcal{V}^{l_{1}l_{2}l_{3}}(k_{1},q_{2},q_{3})$ contains a product of three
spherical Bessel functions. This endows the SFB bispectrum with key geometric
features which we discuss shortly. The second term of Eq. 68, given by Eq.
156, contains the contribution from the $f_{\mathrm{NL}}$ and $G_{2}$ terms
and is of a similar form to Eq. 69. Our result matches the general form of the
SFB bispectrum derived in Ref. Bertacca  _et al._ (2018).

###  IV.2 Signal computation

We compute the bispectrum for a uniform grid of $(k_{1},k_{2},k_{3})$ of size
$200^{3}$ with $k_{i}$ between $k_{\mathrm{min}}=4\times 10^{-3}\
h\mbox{Mpc}^{-1}$ (note that future surveys like SPHEREx will be able to probe
down to $10^{-3}\ h\mbox{Mpc}^{-1}$) and $k_{\mathrm{max}}=8\times 10^{-2}\
h\mbox{Mpc}^{-1}$ (to stay within linear regime) with uniform spacing $\Delta
k=3.8\times 10^{-4}\ h\mbox{Mpc}^{-1}$. For the toy window function, we assume
a sphere $W({\bm{r}})=\mathbf{1}_{[0,r_{\mathrm{max}}]}(r)$ with
$r_{\mathrm{max}}=5000\,\mbox{Mpc}h^{-1}$, corresponding to a maximum redshift
$z\sim 4.1$. As a result of the large redshift range chosen here, the kernels
$\mathcal{W}_{l}(k,q)$ in Eq. 67 are highly peaked around $k\approx q$. For
surveys with a smaller redshift extent $r_{\mathrm{max}}$, the kernel
$\mathcal{W}_{l}(k,q)$ would have a smoother peak and lower frequency
oscillations (as in Fig. 1), which would make the computation less
computationally demanding.

Let us now examine more closely the form of the integrals to be computed. Note
first that Eq. 67 is a two-dimensional integral (over $q_{2}$ and $q_{3}$) of
the kernels $\mathcal{V}_{\rm tot}$. The first term in this kernel,
$\mathcal{V}^{l_{1}l_{2}l_{3}}$ (Eq. 140) is a sum of the one-dimensional
Bessel integrals given in Eq. 70. The second term
$\mathcal{V}^{l_{1}l_{2}l_{3}}_{f_{\mathrm{NL}},G_{2}}$ is also effectively a
one-dimensional integral since the computationally-intensive parts of the
integrand, the integrals $\mathcal{W}^{G_{2}}_{l}$,
$\mathcal{W}^{f_{\mathrm{NL}}}_{l}$ and $I_{l_{1}l_{2}l_{3}}$, can be
precomputed on a grid. The precomputation for $I_{l_{1}l_{2}l_{3}}$ requires a
few seconds using the identity Eq. 54, which expresses it as a finite sum with
$l_{3}$ terms.

Finally, the kernels $\mathcal{V}^{l_{1}l_{2}l_{3}}(k_{1},q_{2},q_{3})$
themselves are precomputed on a $q_{2}$-$q_{3}$ grid to be reused for the
various triplets $(k_{1},k_{2},k_{3})$. This, along with the final integration
in Eq. 67, is the bottleneck of the calculation, and limits the number of
$l$-triplets we may feasibly calculate. However, the signal is sufficiently
smooth in $l$ that this might not pose a problem for e.g, a Fisher forecast
exercise. All grid computations are parallelized using Julia’s multithreading
functionality; computing the kernels
$\mathcal{V}^{l_{1}l_{2}l_{3}}(k_{1},q_{2},q_{3})$ requires a few hours per
triplet $(l_{1},l_{2},l_{3})$ with multithreading across $256$ AMD EPYC 7763
CPUs.

Beyond dimensionality, a second numerical concern is the oscillatory nature of
the integrands. We perform all integrals via Gauss-Legendre quadrature. To
accurately integrate
$\mathcal{V}_{\mathrm{tot}}^{l_{1}l_{2}l_{3}}(k_{1},q_{2},q_{3})$, we observe
(empirically) that if the averaging $r$-spacing is $\Delta r$, then we must
impose $k_{1}+q_{2}+q_{3}\lesssim\pi/\Delta r$. We use $\Delta r=5\
\mbox{Mpc}h^{-1}$. Further, to evaluate the bispectrum we must convolve the
kernels $\mathcal{V}_{\mathrm{tot}}^{l_{1}l_{2}l_{3}}(k_{1},q_{2},q_{3})$ with
the $\mathcal{W}_{l_{2}}(k_{2},q_{2})$, both of which oscillate quickly as
$q_{2}$ is varied. This requires a sufficient sampling in the $q_{i}$ space.
As the computation time is quadratic in the number of sampling points for each
$q_{i}$, only modest samplings are feasible; we use $300$ points for each
$q_{i}$, chosen as Gauss-Legendre nodes.

###  IV.3 Properties of the signal

The observed bispectrum signal displays a number of salient features which we
now discuss.

####  IV.3.1 Oscillations in $k$ and mode coupling

A cross-section of the reduced bispectrum for fixed
$k_{1}=0.0418\,h\mbox{Mpc}^{-1}$, for the same set of $l$-triplets as in Fig.
2, is shown in Fig. 4. Perhaps the most striking feature here is that the
patterns of oscillations in $k$-space are similar to those visible in Fig. 2
for the isotropic and homogeneous case. We may understand this from the mode
coupling coefficients $g^{L_{1}L_{2}\ell}_{l_{1}l_{2}l_{3}L_{3}L_{4}}$ in Eq.
69. Numerically we find that they are generally suppressed unless
$(L_{3},L_{4})=(L_{1},L_{2})$, such that the dominant contribution to the
bispectrum signal $B^{\mathrm{SFB}}_{l_{1}l_{2}l_{3}}$ is from integrals of
type Eq. 70 where the spherical Bessel functions have the same indices
$l_{1},l_{2},l_{3}$, as in the isotropic and homogeneous case (Eq. 51 – 53).

Figure 4: The observed reduced SFB bispectrum for a realistic Universe,
assuming a spherically symmetric survey window, for the same set of
$l$-triplets as in Fig. 2. The signal is sampled at $200^{2}$ pairs
$(k_{2},k_{3})$ in each panel, whereas in Fig. 2 there are $400^{2}$ sampled
pairs.

####  IV.3.2 Limber suppression at low $k$ and high $l$

Figure 5: A illustration of the Limber suppression at low $k$, which becomes
visible for the reduced observed SFB bispectrum signal within the triangle
inequality region for sufficiently large $l$. Here we show the signal for
$l_{1}=l_{2}=l_{3}=90$, while still fixing $k_{1}=4.18\times
10^{-2}\,h\mbox{Mpc}^{-1}$. The onset of Limber suppression is indicated by
the gray dashed lines, where we expect the signal to be suppressed according
to the Limber approximation $k_{i}=(l_{i}+\frac{1}{2})/r_{\mathrm{max}}$. The
border of the region where the triangle inequality on $(k_{1},k_{2},k_{3})$
holds is shown by the solid lines.

In Fig. 5, we increase the $l$ values to $l_{1}=l_{2}=l_{3}=90$, and see that
the overall amplitude of the oscillations decreases by roughly an order of
magnitude relative to the leftmost panel of Fig. 4. To understand the
suppression with increasing $l$, which is generic, we note that for fixed $r$,
the spherical Bessel function $j_{l}(kr)$ is proportional to $(kr)^{l}$ for
small $k$; for large $k$, it oscillates with an amplitude proportional to
$(kr)^{-1}$. Further, in the Limber approximation555Calling this the _Limber
approximation_ is standard in cosmology. However, the term is slightly
misleading, since the original approximation by Limber (1954) was in
configuration space, and only applied to Fourier-space by Kaiser (1998). The
resulting approximation effectively is that of Eq. 75 (LoVerde and Afshordi,
2008). (Eq. 75), for a fixed $l$ and $k$, the Bessel function $j_{l}(kr)$ is
peaked around $r\sim(l+\frac{1}{2})/k$ with a peak value equal to
$\sim\sqrt{\pi/(2l)}$ (LoVerde and Afshordi, 2008).

Physically, we may understand the suppression as follows. For fixed
$(k_{1},k_{2},k_{3})$, the bispectrum at higher $l_{i}$ probes higher
redshifts. If the survey window has finite radial extent, these higher
redshift contributions to the signal are necessarily smaller, and vanish once
the redshift exceeds the extent of the survey. By contrast, when the survey
window has infinite size as in the homogeneous and isotropic Universe, for
every $l_{i}$ there is a corresponding redshift which contributes non-
negligibly to the signal, such that there is no such suppression at high
$l_{i}$.

Given that the spherical Bessel functions go as $(kr)^{l}$ when $k\lesssim
l/r$, we should also expect a sharp suppression of the bispectrum at low $k$.
This suppression is not visible in the first panel of Fig. 4 because this
effect is only relevant for $l\gtrsim k_{\mathrm{min}}r_{\mathrm{max}}=20$ in
our fiducial setup. On the other hand it visible in Fig. 5, for
$l_{1}=l_{2}=l_{3}=90$. The onset of the low $k$ suppression is inside of the
region where the triangle inequality holds, and is indicated by the dotted
gray lines.

Recall from the previous sections that the frequency of the oscillations
increases as we approach the triangle inequality boundary, making the
computation increasingly difficult close to the boundary with higher sampling
needed to resolve these oscillations. For small $l$, we have large spacings
$\Delta k\propto l^{-2}$ which are manageable. For large $l$, the Limber
suppression is helpful in the sense that it is not necessary to compute the
signal at the boundary since it is suppressed there by several orders of
magnitude. This is true both in a Fisher analysis or in a real data analysis
where we can simply choose to ignore this part of the data vector as it
contains almost no information.

####  IV.3.3 Violation of the triangle condition

We saw that in an isotropic and homogeneous Universe, the bispectrum signal
vanishes identically when $(k_{1},k_{2},k_{3})$ does not satisfy the triangle
inequality. This is a consequence of (intrinsic) isotropy. In the case of the
observed bispectrum, where this isotropy is broken by e.g, redshift space
distortions, this is no longer true, though it holds approximately for wider
survey windows. For lower values of $r_{\mathrm{max}}$, the signal strength is
non-negligible even when $(k_{1},k_{2},k_{3})$ violate the triangle
inequality.

This is due to two reasons. First, the kernels $\mathcal{W}_{l}(k,q)$ (Eq. 17)
are less peaked for smaller $r_{\mathrm{max}}$, such that for
$(k_{1},k_{2},k_{3})$ which violate the triangle inequality, the integral Eq.
67 can pick up a non-negligible contribution from the kernel
$\mathcal{V}^{l_{1}l_{2}l_{3}}_{\rm tot}(k_{1},q_{2},q_{3})$ for
$(k_{1},q_{2},q_{3})$ which do form a triangle. Secondly, for smaller
$r_{\mathrm{max}}$, $\mathcal{V}^{l_{1}l_{2}l_{3}}_{\rm
tot}(k_{1},q_{2},q_{3})$ itself can take non-negligible values for
$(k_{1},q_{2},q_{3})$ which do not form a triangle (indeed, the Bessel
integrals Eq. 70 do not come with any triangle condition). As a result, the
triangle inequality is broken in the observed SFB bispectrum, with the
violation more severe at lower $r_{\rm max}$. See Fig. 6 which illustrates
this effect using $r_{\rm max}=500,1000$ and $2500$ $\,\mbox{Mpc}h^{-1}$.

Figure 6: A demonstration of the triangle condition violation in a realistic
Universe – the observed reduced bispectrum signal for different values of the
redshift extent $r_{\mathrm{max}}=500,1000$ and $2500\,\mbox{Mpc}h^{-1}$
(corresponding to $z\sim 0.18$, $0.37$ and $1.2$ respectively) at fixed
$k_{1}=4.18\times 10^{-2}\ h\mbox{Mpc}^{-1}$ and
$(l_{1},l_{2},l_{3})=(4,6,8)$. For smaller $r_{\rm max}$, the violation of the
triangle condition is more severe due to less peaked $\mathcal{W}_{l}$ kernels
(see Section IV.3.3 for more explanations).

##  V Discussion

In this paper, we computed the SFB bispectrum signal for the first time,
discussing how to account for redshift space distortions and primordial non-
Gaussianity. Starting with a toy example of the homogeneous and isotropic
Universe, we built up intuition for later understanding some key features of
the observed bispectrum convolved with a toy spherically symmetric window.

To render the computation tractable, we leveraged a decomposition of the
second-order redshift space kernel into products of three Legendre
polynomials, which allowed us to express the bispectrum, modulo RSD and PNG
terms, as a triple integral. Furthermore, we derived an identity to express as
a simple sum the $6$-dimensional angular integral of three spherical harmonics
(or equivalently, the one-dimensional integral of a product of three spherical
Bessel functions on an infinite interval). This enabled us to rapidly compute
and study the signal in the case of a homogeneous and isotropic Universe, and
to accelerate the calculation of the RSD and PNG contributions to the observed
bispectrum.

Even with these techniques and the various numerical optimizations we
employed, computing the SFB bispectrum is clearly expensive: for each triplet
of multipole indices $(l_{1},l_{2},l_{3})$ and of wavenumbers
$(k_{1},k_{2},k_{3})$, we need to evaluate triple integrals with oscillatory
integrands. Our method requires $\mathcal{O}(100)$ CPU hours for each
$l$-triplet to compute the signal on a grid of $200^{3}$ $(k_{1},k_{2},k_{3})$
triplets. In a realistic data analysis, one would need to calculate the signal
for different cosmologies in a Monte-Carlo Markov Chain (MCMC) on the order of
seconds. We note however that we have chosen a very large redshift extent
$z\lesssim 4$, corresponding to $r_{\rm max}=5000\,\mbox{Mpc}h^{-1}$, for
which the integration is the most challenging. Surveys with smaller redshift
extent would require less time for computation.

There are several possibilities to accelerate and improve the accuracy of the
computation that we leave for future work. For example, as the local frequency
of the oscillations in the signal can be estimated as a function of
$(k_{i},l_{i})$ from the isotropic and homogeneous bispectrum as in Section
III.3, the signal can instead be sampled on a suitable non-uniform grid of
$(k_{1},k_{2},k_{3})$. To improve upon the Gauss-Legendre quadrature method
for integrating the spherical Bessel product Eq. 70, the $3$-dimensional
generalization of the FFTLog method of Grasshorn Gebhardt and Jeong (2018)
could prove superior. Leveraging cache-friendly memory layouts may also help
to speed up some of the linear algebra operations involved in computing these
Bessel integrals on a grid in $k$-space.

For the purpose of calculations in a MCMC analysis, one could also explore
decomposing the bispectrum dependence on cosmological parameters into
precomputed templates, and varying the coefficients of the template
corresponding to the varying cosmology. If it becomes impossible to directly
compute the signal for each point in the cosmology parameter space sampled
during the MCMC, the use of emulators could also aid in minimizing the
evaluation time.

Another natural extension to this work is to incorporate more physical and
observational effects in the computation. The most important physical effects
on large-scales that we have not included here are general-relativistic (GR)
effects. The authors of Ref. Bertacca  _et al._ (2018) detailed how to
incorporate them in the SFB bispectrum; in principle they may be evaluated
numerically within the same framework described in our work, i.e through
(many) additional terms in the first order kernels $\mathcal{W}_{l}(k,q)$ and
in the second-order kernels $\mathcal{V}^{l_{1}l_{2}l_{3}}$, after a template
decomposition into Legendre polynomials as in Section IV.1. On smaller scales,
more detailed modeling of RSD (e.g. Fingers-of-God and Alcock-Paczynski
effects) and of the nonlinear regime would be needed. In addition, we have
used a spherically symmetric window function to demonstrate the calculation,
while realistic window convolution is still to be explored.

To feasibly use the SFB bispectrum to analyze survey data, a number of missing
pieces would still need to be filled in. In particular, it would be necessary
to develop an efficient SFB bispectrum estimator, e.g, by building off of
techniques developed by Gebhardt and Doré (2021); Grasshorn Gebhardt and Doré
(2023) for an SFB power spectrum estimator. As allowing for a survey window of
arbitrary geometry in the modeling of the signal would greatly increase the
computational cost, one may explore accounting for it in the estimator, e.g,
by using a windowless estimator which directly returns the window-deconvolved
bispectrum as pioneered in Ref. Philcox (2021) for the bispectrum multipoles.

Moreover, a realistic covariance matrix for the SFB bispectrum beyond the
Gaussian approximation also needs to be developed, including complexities due
to window function convolution as well as non-Gaussian covariance. If the
window effects can be reliably removed at the estimator level, then the
covariance would be significantly simplified. For the non-Gaussian covariance,
an approximation similar in form to that proposed in Ref. Biagetti  _et al._
(2022) may be applicable for the SFB bispectrum, where the non-Gaussian part
of the covariance is dominated by the product of two bispectra sharing the
same large scale – a good approximation for squeezed configurations and also
tested to be good enough for other configurations in Ref. Biagetti  _et al._
(2022) in the context of Fourier-space bispectrum. Alternatively, to
incorporate all complexities at once, one may also develop mocks to compute
the mock-based covariance by averaging over many realizations, once a fast SFB
bispectrum estimator exists. This method would include wide-angle effects
directly for mocks with large enough angular area, while it could be
challenging to incorporate all GR effects into the mocks.

Since an advantage of the SFB formalism is that it avoids the loss of
information due to assuming inexact lines-of-sights for individual galaxy
triplets, it would also be interesting to evaluate more quantitatively now
this information gain, for example by comparing to the standard bispectrum
multipole formalism in the local plane-parallel approximation and to
perturbative corrections thereof as in Pardede  _et al._ (2023). Our work to
enable the computation of the signal will allow for such a study to be
conducted. With a suitable scheme to interpolate the signal in the space of
multipoles $l_{i}$, it could be feasible to conduct a Fisher forecast for
various cosmological parameters of interest, such as $f_{\mathrm{NL}}$ or RSD
parameters.

Note that this loss of information may be small for surveys with small angular
extent, but more important for full-sky surveys like SPHEREx. Currently, with
the exception of the TSH formalism, only perturbative approaches to modeling
wide-angle effects in the bispectrum have been proposed Pardede  _et al._
(2023); Noorikuhani and Scoccimarro (2023), expanding from the global plane-
parallel approximation. Thus, the SFB bispectrum remains the only method to
fully account for all large scale effects non-perturbatively on the largest
angular scales while preserving the potential of retaining all information
contained in the radial modes.

Another advantage of the SFB formalism is that some of the GR terms, which are
mostly radial effects along the line-of-sight, become easily disentangled from
other effects. In particular, the monopole and the dipole terms in the SFB
formalism contain all the observer terms in GR arising from the potential and
velocity at the observer position. Some of these terms may be quantified via
other means before being subtracted (e.g. the velocity term in the dipole),
while others are intrinsically undetectable (e.g. observer potential) and may
need to be modeled through constrained realization if they affect observables
of interest.

Modeling these terms in a Cartesian framework amounts to propagating these
effects to every mode (and every order if a perturbative expansion from the
plane-parallel approximation is used), which would propagate potential
systematics into every measured mode. In contrast, in a spherical framework
such as SFB, there is a clear radial and angular separation that allows for
the isolation of such terms into just the monopole and dipole, which may then
be discarded or tested separately for systematics.

While the TSH formalism provides a similar advantage, it requires many radial
bins required to resolve the large scale radial modes (which is important to
do for measuring $f_{\mathrm{NL}}$), which introduces highly correlated
neighboring radial bins, and leads to numerically instabilities during
covariance inversion. The SFB method is therefore a trade-off between
extracting the maximal amount of information and the cost of computing the
signal. In this regard, we have made a step forward by rendering the SFB
signal computable and studying its various features.

This is merely the beginning of more efforts to follow to make the calculation
of the SFB bispectrum feasible for next-generation surveys. With increasing
computational power in the future, along with more sophisticated numerical and
mathematical techniques, the SFB bispectrum may become a key formalism that
will allow us to extract all of the possible information from a full-sky
galaxy survey.

##  VI  Acknowledgements

This research used resources of the National Energy Research Scientific
Computing Center (NERSC). We thank Johannes Blaschke for invaluable technical
advice on running and optimizing our Julia code on NERSC, Vincent Desjacques
for enlightening discussions, and Robin Wen for helpful feedback on the
manuscript. J.B. is supported in part by the DOE Early Career Grant
DESC0019225. A.T. is also grateful for the support provided by the Walter
Burke Institute for Theoretical Physics. C.H. thanks her Maker for providing
the amazing people to work together, the opportunities to collaborate
fruitfully, and all the resources and any insights that made this work
possible. This material is based upon work supported by the U.S. Department of
Energy, Office of Science, Office of High Energy Physics, under Award Number
DE-SC0011632.

## References

  * Aghamousa  _et al._ (2016) A. Aghamousa _et al._ (DESI), (2016), arXiv:1611.00036 [astro-ph.IM] .
  * Amendola  _et al._ (2018) L. Amendola _et al._ , Living Rev. Rel. 21, 2 (2018), arXiv:1606.00180 [astro-ph.CO] .
  * Doré  _et al._ (2015) O. Doré, J. Bock, M. Ashby, P. Capak, A. Cooray, R. de Putter, T. Eifler, N. Flagey, Y. Gong, S. Habib, K. Heitmann, C. Hirata, W.-S. Jeong, R. Katti, P. Korngut, E. Krause, D.-H. Lee, D. Masters, P. Mauskopf, G. Melnick, B. Mennesson, H. Nguyen, K. Öberg, A. Pullen, A. Raccanelli, R. Smith, Y.-S. Song, V. Tolls, S. Unwin, T. Venumadhav, M. Viero, M. Werner, and M. Zemcov, “Cosmology with the spherex all-sky spectral survey,” (2015), arXiv:1412.4872 [astro-ph.CO] .
  * Eifler  _et al._ (2021) T. Eifler _et al._ , Mon. Not. Roy. Astron. Soc. 507, 1514 (2021), arXiv:2004.04702 [astro-ph.CO] .
  * Dalal  _et al._ (2008) N. Dalal, O. Dore, D. Huterer, and A. Shirokov, Phys. Rev. D 77, 123514 (2008), arXiv:0710.4560 [astro-ph] .
  * McDonald (2009) P. McDonald, Journal of Cosmology and Astroparticle Physics 2009, 026–026 (2009).
  * Yoo  _et al._ (2012) J. Yoo, N. Hamaus, U. Seljak, and M. Zaldarriaga, Physical Review D 86 (2012), 10.1103/physrevd.86.063514.
  * Yoo (2010a) J. Yoo, Physical Review D 82 (2010a), 10.1103/physrevd.82.083508.
  * Bonvin and Durrer (2011) C. Bonvin and R. Durrer, Physical Review D 84 (2011), 10.1103/physrevd.84.063505.
  * Challinor and Lewis (2011) A. Challinor and A. Lewis, Physical Review D 84 (2011), 10.1103/physrevd.84.043516.
  * Yamamoto  _et al._ (2006) K. Yamamoto, M. Nakamichi, A. Kamino, B. A. Bassett, and H. Nishioka, Publications of the Astronomical Society of Japan 58, 93–102 (2006).
  * Bianchi  _et al._ (2015) D. Bianchi, H. Gil-Marín, R. Ruggeri, and W. J. Percival, Monthly Notices of the Royal Astronomical Society: Letters 453, L11–L15 (2015).
  * Scoccimarro (2015) R. Scoccimarro, Physical Review D 92 (2015), 10.1103/physrevd.92.083532.
  * Reimberg  _et al._ (2016) P. Reimberg, F. Bernardeau, and C. Pitrou, Journal of Cosmology and Astroparticle Physics 2016, 048–048 (2016).
  * Beutler  _et al._ (2019) F. Beutler, E. Castorina, and P. Zhang, JCAP 2019, 040 (2019), arXiv:1810.05051 [astro-ph.CO] .
  * Benabou  _et al._ (prep) J. N. Benabou, I. Sands, C. Heinrich, H. S. Grasshorn Gebhardt, and O. Doré, (in prep).
  * Pardede  _et al._ (2023) K. Pardede, E. D. Dio, and E. Castorina, “Wide-angle effects in the galaxy bispectrum,” (2023), arXiv:2302.12789 [astro-ph.CO] .
  * Di Dio  _et al._ (2016) E. Di Dio, R. Durrer, G. Marozzi, and F. Montanari, Journal of Cosmology and Astroparticle Physics 2016, 016–016 (2016).
  * Lahav (1993) O. Lahav, “Spherical harmonic reconstruction of cosmic density and velocity fields,” (1993), arXiv:astro-ph/9309030 [astro-ph] .
  * Binney and Quinn (1991) J. Binney and T. Quinn, Monthly Notices of the Royal Astronomical Society 249, 678 (1991).
  * Heavens and Taylor (1995) A. F. Heavens and A. N. Taylor, Monthly Notices of the Royal Astronomical Society 275, 483–497 (1995).
  * Tadros  _et al._ (1999) H. Tadros, W. E. Ballinger, A. N. Taylor, A. F. Heavens, G. Efstathiou, W. Saunders, C. S. Frenk, O. Keeble, R. McMahon, S. J. Maddox, S. Oliver, M. Rowan-Robinson, W. J. Sutherland, and S. D. M. White, Monthly Notices of the Royal Astronomical Society 305, 527–546 (1999).
  * Percival  _et al._ (2004) W. J. Percival, D. Burkey, A. Heavens, A. Taylor, S. Cole, J. A. Peacock, C. M. Baugh, J. Bland-Hawthorn, T. Bridges, R. Cannon, M. Colless, C. Collins, W. Couch, G. Dalton, R. De Propris, S. P. Driver, G. Efstathiou, R. S. Ellis, C. S. Frenk, K. Glazebrook, C. Jackson, O. Lahav, I. Lewis, S. Lumsden, S. Maddox, P. Norberg, B. A. Peterson, W. Sutherland, and K. Taylor, Monthly Notices of the Royal Astronomical Society 353, 1201–1218 (2004).
  * Samushia (2019) L. Samushia, “Proper fourier decomposition formalism for cosmological fields in spherical shells,” (2019), arXiv:1906.05866 [astro-ph.CO] .
  * Gebhardt and Doré (2021) H. S. G. Gebhardt and O. Doré, “Superfab: a fabulous code for spherical fourier-bessel decomposition,” (2021), arXiv:2102.10079 [astro-ph.CO] .
  * Wang  _et al._ (2020) M. Wang, S. Avila, D. Bianchi, R. Crittenden, and W. Percival, Journal of Cosmology and Astroparticle Physics 2020, 022–022 (2020).
  * Grasshorn Gebhardt and Doré (2022) H. S. Grasshorn Gebhardt and O. Doré, JCAP 01, 038 (2022), arXiv:2109.13352 [astro-ph.CO] .
  * Grasshorn Gebhardt and Doré (2023) H. Grasshorn Gebhardt and O. Doré, (2023), arXiv:2310.17677 [astro-ph.IM] .
  * Zhang  _et al._ (2021) Y. Zhang, A. R. Pullen, and A. S. Maniyar, Physical Review D 104 (2021), 10.1103/physrevd.104.103523.
  * Gualdi and Verde (2020) D. Gualdi and L. Verde, Journal of Cosmology and Astroparticle Physics 2020, 041–041 (2020).
  * Heinrich  _et al._ (2023) C. Heinrich, O. Dore, and E. Krause, (2023), arXiv:2311.13082 [astro-ph.CO] .
  * Coulton  _et al._ (2019) W. R. Coulton, J. Liu, M. S. Madhavacheril, V. Böhm, and D. N. Spergel, Journal of Cosmology and Astroparticle Physics 2019, 043–043 (2019).
  * Ruggeri  _et al._ (2018) R. Ruggeri, E. Castorina, C. Carbone, and E. Sefusatti, Journal of Cosmology and Astroparticle Physics 2018, 003–003 (2018).
  * Maartens  _et al._ (2020) R. Maartens, S. Jolicoeur, O. Umeh, E. M. D. Weerd, C. Clarkson, and S. Camera, Journal of Cosmology and Astroparticle Physics 2020, 065–065 (2020).
  * Jeong and Schmidt (2020) D. Jeong and F. Schmidt, Phys. Rev. D 102, 023530 (2020), arXiv:1906.05198 [astro-ph.CO] .
  * Clarkson  _et al._ (2019) C. Clarkson, E. M. de Weerd, S. Jolicoeur, R. Maartens, and O. Umeh, Mon. Not. Roy. Astron. Soc. 486, L101 (2019), arXiv:1812.09512 [astro-ph.CO] .
  * Bertacca  _et al._ (2018) D. Bertacca, A. Raccanelli, N. Bartolo, M. Liguori, S. Matarrese, and L. Verde, Physical Review D 97 (2018), 10.1103/physrevd.97.023531.
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Lewis  _et al._ (2001) A. Lewis, A. Challinor, and A. Lasenby, The Astrophysical Journal 538, 473 (2001).
  * Bezanson  _et al._ (2017) J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, SIAM Review 59, 65 (2017).
  * Tellarini  _et al._ (2016) M. Tellarini, A. J. Ross, G. Tasinato, and D. Wands, Journal of Cosmology and Astroparticle Physics 2016, 014–014 (2016).
  * Bernardeau  _et al._ (2002) F. Bernardeau, S. Colombi, E. Gaztañaga, and R. Scoccimarro, Physics Reports 367, 1–248 (2002).
  * Wands and Slosar (2009) D. Wands and A. Slosar, Phys. Rev. D 79, 123507 (2009), arXiv:0902.1084 [astro-ph.CO] .
  * Yoo  _et al._ (2009) J. Yoo, A. L. Fitzpatrick, and M. Zaldarriaga, Phys. Rev. D 80, 083514 (2009), arXiv:0907.0707 [astro-ph.CO] .
  * Yoo (2010b) J. Yoo, Phys. Rev. D 82, 083508 (2010b), arXiv:1009.3021 [astro-ph.CO] .
  * Catelan  _et al._ (2000) P. Catelan, C. Porciani, and M. Kamionkowski, Mon. Not. Roy. Astron. Soc. 318, 39 (2000), arXiv:astro-ph/0005544 .
  * Baldauf  _et al._ (2012) T. Baldauf, U. Seljak, V. Desjacques, and P. McDonald, Phys. Rev. D 86, 083540 (2012), arXiv:1201.4827 [astro-ph.CO] .
  * Szego (1936) G. Szego, Transactions of the American Mathematical Society 39, 1 (1936).
  * Limber (1954) D. N. Limber, Astrophys. J. 119, 655 (1954).
  * Kaiser (1998) N. Kaiser, Astrophys. J. 498, 26 (1998), arXiv:astro-ph/9610120 [astro-ph] .
  * LoVerde and Afshordi (2008) M. LoVerde and N. Afshordi, Phys. Rev. D 78, 123506 (2008), arXiv:0809.5112 [astro-ph] .
  * Grasshorn Gebhardt and Jeong (2018) H. S. Grasshorn Gebhardt and D. Jeong, Physical Review D 97 (2018), 10.1103/physrevd.97.023504.
  * Philcox (2021) O. H. Philcox, Physical Review D 104 (2021), 10.1103/physrevd.104.123529.
  * Biagetti  _et al._ (2022) M. Biagetti, L. Castiblanco, J. Noreña, and E. Sefusatti, Journal of Cosmology and Astroparticle Physics 2022, 009 (2022).
  * Noorikuhani and Scoccimarro (2023) M. Noorikuhani and R. Scoccimarro, Phys. Rev. D 107, 083528 (2023), arXiv:2207.12383 [astro-ph.CO] .
  * Slepian (2018) Z. Slepian, “On decoupling the integrals of cosmological perturbation theory,” (2018), arXiv:1812.02728 [astro-ph.CO] .
  * Mehrem  _et al._ (1991) R. Mehrem, J. T. Londergan, and M. H. Macfarlane, Journal of Physics A: Mathematical and General 24, 1435 (1991).
  * Fabrikant (2013) V. I. Fabrikant, Quarterly of Applied Mathematics 71, 573 (2013).
  * Dong  _et al._ (2015) H. Dong, J. S. Kim, and G. S. Chirikjian, Journal of Computational Biology 22, 787–805 (2015).
  * Chellino and Slepian (2023) J. Chellino and Z. Slepian, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 479 (2023), 10.1098/rspa.2023.0138.
  * Jackson and Maximon (1972) A. Jackson and L. Maximon, Siam Journal on Mathematical Analysis - SIAM J MATH ANAL 3 (1972), 10.1137/0503043.

##  Appendix A Useful formulae

###  A.1 Dirac delta distribution

For a continuously differentiable function $g$ with simple roots $\\{x_{i}\\}$
we have

$\displaystyle\delta_{D}(g(x))=\sum_{i}\frac{\delta_{D}(x-x_{i})}{|g^{\prime}(x_{i})|}\,.$
(71)

For any function $f({\bm{k}})$,

$\displaystyle\int
k^{2}\mathrm{d}{k}\,\mathrm{d}^{2}{\hat{\bm{k}}}\,\delta^{D}({\bm{k}}-{\bm{k}}^{\prime})\,f({\bm{k}})=\int\mathrm{d}{k}\,\delta^{D}(k-k^{\prime})\int\mathrm{d}^{2}{\hat{\bm{k}}}\,\delta^{D}({\hat{\bm{k}}}-{\hat{\bm{k}}}^{\prime})\,f({\bm{k}})\,.$
(72)

Therefore,

$\displaystyle\delta^{D}({\bm{k}}-{\bm{k}}^{\prime})$
$\displaystyle=k^{-2}\,\delta^{D}(k-k^{\prime})\,\delta^{D}({\hat{\bm{k}}}-{\hat{\bm{k}}}^{\prime})\,.$
(73)

###  A.2 Spherical Bessel functions

To first order the Bessel function $J_{\nu}(x)$ may be approximated by a Dirac
delta as LoVerde and Afshordi (2008)

$\displaystyle J_{\nu}(kr)$
$\displaystyle\simeq\delta^{D}\\!\left(kr-\nu\right)\,.$ (74)

Therefore, for a spherical Bessel function
$j_{\ell}(x)=\sqrt{\pi/2x}\,J_{\ell+\frac{1}{2}}(x)$ we have

$\displaystyle j_{\ell}(kr)$
$\displaystyle\simeq\sqrt{\frac{\pi}{2rk}}\,\delta^{D}\\!\left(r-\frac{\ell+\frac{1}{2}}{k}\right)\,,$
(75)

to first order. In the cosmology literature, a version of this is often called
_Limber’s approximation_ (Limber, 1954). Spherical Bessel functions satisfy
the orthogonality relation

$\displaystyle\delta^{D}(k-k^{\prime})$
$\displaystyle=\frac{2kk^{\prime}}{\pi}\int_{0}^{\infty}\mathrm{d}{r}\,r^{2}\,j_{\ell}(kr)\,j_{\ell}(k^{\prime}r)\,.$
(76)

###  A.3 Spherical harmonics

Spherical harmonics can be expressed in terms of a complex exponential and
real associated Legendre functions $\mathrm{P}_{\ell}^{m}(x)$ as

$\displaystyle Y_{\ell m}({\hat{\bm{r}}})$
$\displaystyle=e^{im\phi}\left(\frac{(\ell-m)!(2\ell+1)}{4\pi\left(\ell+m\right)!}\right)^{\frac{1}{2}}\mathrm{P}_{\ell}^{m}\\!\left(\cos\theta\right).$
(77)

The associated Legendre functions are even or odd according to the index,

$\displaystyle P_{\ell}^{m}(-x)=(-1)^{\ell+m}P_{\ell}^{m}(x)\,.$ (78)

The completeness relation is

$\displaystyle\sum_{\ell m}Y_{\ell m}({\hat{\bm{r}}})\,Y^{*}_{\ell
m}({\hat{\bm{r}}}^{\prime})$
$\displaystyle=\delta^{D}({\hat{\bm{r}}}-{\hat{\bm{r}}}^{\prime})\,.$ (79)

The spherical harmonics satisfy the orthogonality relation

$\displaystyle\int\mathrm{d}{\Omega}_{{\hat{\bm{r}}}}\,Y_{\ell
m}({\hat{\bm{r}}})\,Y^{*}_{\ell^{\prime}m^{\prime}}({\hat{\bm{r}}})=\delta^{K}_{\ell\ell^{\prime}}\delta^{K}_{mm^{\prime}}\,.$
(80)

For a rotation $\mathcal{R}$ about the origin that sends the unit vector
$\mathbf{r}$ to $\mathbf{r}^{\prime}$, we have

$\displaystyle
Y_{\ell,m}(\mathbf{r}^{\prime})=\sum_{m^{\prime}=-\ell}^{\ell}[D^{(\ell)}_{mm^{\prime}}({\mathcal{R}})]^{*}Y_{\ell,m^{\prime}}({\mathbf{r}}),$
(81)

where $[D^{(\ell)}_{mm^{\prime}}({\mathcal{R}})]^{*}$ is the complex conjugate
of an entry of the Wigner $D$-matrix. The Wigner $D$-matrix is a unitary
square matrix of dimension $2j+1$. If $\mathcal{R}$ is defined by proper Euler
angles $\alpha,\beta,\gamma$ in the $z$-$y$-$z$ convention, we have the
property

$\displaystyle
D^{\ell}_{m0}(\mathcal{R})=\sqrt{\frac{4\pi}{2\ell+1}}Y_{\ell,m}^{*}(\beta,\alpha)\,,$
(82)

and also the relation

$\displaystyle
D^{j}_{mk}(\mathcal{R})D^{j^{\prime}}_{m^{\prime}k^{\prime}}(\mathcal{R})=\sum_{J=|j-j^{\prime}|}^{j+j^{\prime}}\langle
jmj^{\prime}m^{\prime}|J\left(m+m^{\prime}\right)\rangle\langle
jkj^{\prime}k^{\prime}|J\left(k+k^{\prime}\right)\rangle
D^{J}_{\left(m+m^{\prime}\right)\left(k+k^{\prime}\right)}(\mathcal{R})\,,$
(83)

where $\langle j_{1}m_{1}j_{2}m_{2}|j_{3}m_{3}\rangle$ is a Clebsch-Gordan
coefficient. The latter is related to the Wigner $3j$ symbols by:

$\displaystyle\langle j_{1}\,m_{1}\,j_{2}\,m_{2}|J\,M\rangle$
$\displaystyle=(-1)^{-j_{1}+j_{2}-M}\sqrt{2J+1}\begin{pmatrix}j_{1}&j_{2}&J\\\
m_{1}&m_{2}&-M\end{pmatrix}.$ (84)

We may also expand plane waves in terms of spherical Bessels and spherical
harmonics,

$\displaystyle e^{i{\bm{q}}\cdot{\bm{r}}}$
$\displaystyle=4\pi\sum_{\ell^{\prime},m^{\prime}}i^{\ell^{\prime}}j_{\ell^{\prime}}(qr)\,Y^{*}_{\ell^{\prime}m^{\prime}}({\hat{\bm{q}}})\,Y_{\ell^{\prime}m^{\prime}}({\hat{\bm{r}}})\,,$
(85)

from which it follows, using Eq. 80, that

$\displaystyle\int\mathrm{d}^{2}{\hat{\bm{r}}}Y^{*}_{l,m}({\hat{\bm{r}}})e^{i{\bm{q}}\cdot{\bm{r}}}({\hat{\bm{q}}}\cdot{\hat{\bm{r}}})^{\alpha}=4\pi
i^{l}Y^{*}_{l,m}({\hat{\bm{q}}})(-i\partial_{qr})^{\alpha}j_{l}(qr)\,.$ (86)

The Legendre polynomials can be expressed as a sum over spherical harmonics as

$\displaystyle\mathcal{L}_{\ell}({\hat{\bm{k}}}\cdot{\hat{\bm{r}}})$
$\displaystyle=\frac{4\pi}{2\ell+1}\sum_{m}Y_{\ell m}({\hat{\bm{k}}})\,Y_{\ell
m}^{*}({\hat{\bm{r}}})\,.$ (87)

We also have the identities

$\displaystyle Y^{*}_{\ell m}({\hat{\bm{r}}})$
$\displaystyle=(-1)^{m}Y_{\ell,-m}({\hat{\bm{r}}})\,,$ (88) $\displaystyle
Y_{\ell m}(-{\hat{\bm{r}}})$
$\displaystyle=(-1)^{\ell}Y_{\ell,m}({\hat{\bm{r}}})\,.$ (89)

###  A.4 Gaunt factor

The Gaunt factor is

$\displaystyle\mathcal{G}^{\ell LL_{1}}_{mMM_{1}}$
$\displaystyle\equiv\int\mathrm{d}^{2}{\hat{\bm{r}}}\,Y_{\ell
m}({\hat{\bm{r}}})\,Y_{LM}({\hat{\bm{r}}})\,Y_{L_{1}M_{1}}({\hat{\bm{r}}})\,,$
(90)

and it can be expressed in terms of Wigner-$3j$ symbols,

$\displaystyle\mathcal{G}^{\ell LL_{1}}_{mMM_{1}}$
$\displaystyle=\left(\frac{(2\ell+1)(2L+1)(2L_{1}+1)}{4\pi}\right)^{\frac{1}{2}}\begin{pmatrix}\ell&L&L_{1}\\\
0&0&0\end{pmatrix}\begin{pmatrix}\ell&L&L_{1}\\\ m&M&M_{1}\end{pmatrix}\,.$
(91)

Hence, a product of two spherical harmonics can be reduced to a linear
combination of spherical harmonics by

$\displaystyle
Y_{l_{1}m_{1}}({\hat{\bm{r}}})Y_{l_{2}m_{2}}({\hat{\bm{r}}})=\sum_{L}(-1)^{M}\mathcal{G}_{m_{1}m_{2}-M}^{l_{1}l_{2}L}Y_{LM}({\hat{\bm{r}}})\,,$
(92)

where $M=m_{1}+m_{2}$. Using this identity, one can derive by recursion the
analogous integral to Eq. 90 for any number of spherical harmonics. For four
spherical harmonics we have (as in Appendix A of Ref. Slepian (2018)):

$\displaystyle\int\mathrm{d}^{2}{\hat{\bm{r}}}\
Y_{l_{1}m_{1}}({\hat{\bm{r}}})Y_{l_{2}m_{2}}({\hat{\bm{r}}})Y_{l_{3}m_{3}}({\hat{\bm{r}}})Y_{l_{4}m_{4}}({\hat{\bm{r}}})=\sum_{L}(-1)^{M}\mathcal{G}^{l_{1}l_{2}L}_{m_{1}m_{2}-M}\mathcal{G}^{Ll_{3}l_{4}}_{Mm_{3}m_{4}}\,.$
(93)

###  A.5 Wigner symbols

The Wigner $3j$ symbols obey an orthogonality relation

$\displaystyle\sum_{mM}\begin{pmatrix}\ell&L&L_{1}\\\
m&M&M_{1}\end{pmatrix}\begin{pmatrix}\ell&L&L_{2}\\\ m&M&M_{2}\end{pmatrix}$
$\displaystyle=\frac{\delta^{K}_{L_{1}L_{2}}\delta^{K}_{M_{1}M_{2}}\delta^{T}(\ell,L,L_{1})}{2L_{1}+1}\,,$
(94)

where $\delta^{T}(\ell,L,L_{1})$ enforces the triangle relation that is obeyed
by the Wigner 3$j$-symbols, i.e they vanish unless $|\ell-L|\leq
L_{1}\leq\ell+L$ and $m+M+M_{1}=0$.

The Wigner $3j$’s acquire a phase for $m_{i}\rightarrow-m_{i}$:

$\displaystyle\begin{pmatrix}j_{1}&j_{2}&j_{3}\\\
-m_{1}&-m_{2}&-m_{3}\end{pmatrix}=(-1)^{j_{1}+j_{2}+j_{3}}\begin{pmatrix}j_{1}&j_{2}&j_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}\,.$ (95)

We also have the identity

$\displaystyle\begin{pmatrix}j_{1}&j_{2}&j_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}\begin{Bmatrix}j_{1}&j_{2}&j_{3}\\\
l_{1}&l_{2}&l_{3}\end{Bmatrix}=\sum_{m^{\prime}_{1}m^{\prime}_{2}m^{\prime}_{3}}(-1)^{l_{1}+l_{2}+l_{3}+m^{\prime}_{1}+m^{\prime}_{2}+m^{\prime}_{3}}\begin{pmatrix}j_{1}&l_{2}&l_{3}\\\
m_{1}&m^{\prime}_{2}&-m^{\prime}_{3}\end{pmatrix}\begin{pmatrix}l_{1}&j_{2}&l_{3}\\\
-m^{\prime}_{1}&m_{2}&m^{\prime}_{3}\end{pmatrix}\begin{pmatrix}l_{1}&l_{2}&j_{3}\\\
m^{\prime}_{1}&-m^{\prime}_{2}&m_{3}\end{pmatrix}\,,$ (96)

where a Wigner $6j$-symbol appears on the LHS. Lastly, we also have

$\displaystyle\begin{pmatrix}j_{13}&j_{23}&j_{33}\\\
m_{13}&m_{23}&m_{33}\end{pmatrix}\begin{Bmatrix}j_{11}&j_{12}&j_{13}\\\
j_{21}&j_{22}&j_{23}\\\
j_{31}&j_{32}&j_{33}\end{Bmatrix}=\sum_{m_{r1},m_{r2},r=1,2,3}\begin{pmatrix}j_{11}&j_{12}&j_{13}\\\
m_{11}&m_{12}&m_{13}\end{pmatrix}\begin{pmatrix}j_{21}&j_{22}&j_{23}\\\
m_{21}&m_{22}&m_{23}\end{pmatrix}\begin{pmatrix}j_{31}&j_{32}&j_{33}\\\
m_{31}&m_{32}&m_{33}\end{pmatrix}$
$\displaystyle\times\begin{pmatrix}j_{11}&j_{21}&j_{31}\\\
m_{11}&m_{21}&m_{31}\end{pmatrix}\begin{pmatrix}j_{12}&j_{22}&j_{32}\\\
m_{12}&m_{22}&m_{32}\end{pmatrix}\,,$ (97)

where a Wigner $9j$-symbol appears on the LHS.

##  Appendix B Encoding of observational isotropy by the Gaunt factor

Here we show that the $3$-point function of the SFB modes of an
observationally-isotropic real-valued field $\delta(\mathbf{r})$ is real and
proportional to the Gaunt factor. To see this, note that in real space, the
3-point function of $\delta$ can only depend on the distances to each point
and the angles on the sky. Therefore, we may expand it in Legendre polynomials
as

$\displaystyle\langle\delta(\bm{r}_{1})\delta(\bm{r}_{2})\delta(\bm{r}_{3})\rangle$
$\displaystyle=\sum_{L_{1}L_{2}L_{3}}f_{L_{1}L_{2}L_{3}}(r_{1},r_{2},r_{3})\,\mathcal{L}_{L_{1}}(\hat{\bm{r}}_{1}\cdot\hat{\bm{r}}_{2})\,\mathcal{L}_{L_{2}}(\hat{\bm{r}}_{2}\cdot\hat{\bm{r}}_{3})\,\mathcal{L}_{L_{3}}(\hat{\bm{r}}_{3}\cdot\hat{\bm{r}}_{1})$
(98)

The Legendre polynomials may be further decomposed into sums over spherical
harmonics via Eq. 87. We may then transform Eq. 98 to spherical harmonic space
to obtain

$\displaystyle\langle\delta_{\ell_{1}m_{1}}(r_{1})\delta_{\ell_{2}m_{2}}(r_{2})\delta_{\ell_{3}m_{3}}(r_{3})\rangle$
$\displaystyle=\int
d^{2}\hat{r}_{1}\,d^{2}\hat{r}_{2}\,d^{2}\hat{r}_{3}\,Y^{*}_{\ell_{1}m_{1}}(\hat{\bm{r}}_{1})\,Y^{*}_{\ell_{2}m_{2}}(\hat{\bm{r}}_{2})\,Y^{*}_{\ell_{3}m_{3}}(\hat{\bm{r}}_{3})$
$\displaystyle\quad\times\sum_{L_{1}L_{2}L_{3}}A_{L_{1}L_{2}L_{3}}\,f_{L_{1}L_{2}L_{3}}(r_{1},r_{2},r_{3})\sum_{M_{1}M_{2}M_{3}}(-1)^{M_{1}+M_{2}+M_{3}}\,Y^{*}_{L_{1}-M_{1}}(\hat{\bm{r}}_{1})Y^{*}_{L_{1}M_{1}}(\hat{\bm{r}}_{2})$
$\displaystyle\quad\times\,Y^{*}_{L_{2}-M_{2}}(\hat{\bm{r}}_{2})Y^{*}_{L_{2}M_{2}}(\hat{\bm{r}}_{3})\,Y^{*}_{L_{3}-M_{3}}(\hat{\bm{r}}_{3})Y^{*}_{L_{3}M_{3}}(\hat{\bm{r}}_{1})$
$\displaystyle=\sum_{L_{1}L_{2}L_{3}}A_{L_{1}L_{2}L_{3}}\,f_{L_{1}L_{2}L_{3}}(r_{1},r_{2},r_{3})\sum_{M_{1}M_{2}M_{3}}(-1)^{M_{1}+M_{2}+M_{3}}\,\mathcal{G}^{\ell_{1}L_{1}L_{3}}_{m_{1}-M_{1}M_{3}}\,\mathcal{G}^{\ell_{2}L_{2}L_{1}}_{m_{2}-M_{2}M_{1}}\,\mathcal{G}^{\ell_{3}L_{3}L_{2}}_{m_{3}-M_{3}M_{2}}\,.$
(99)

where we define $A_{L_{1}L_{2}L_{3}}\equiv\prod_{i}4\pi/(2L_{i}+1)$. Using the
identity Eq. 96 to evaluate the sum of Gaunt factors, we may write

$\displaystyle\langle\delta_{\ell_{1}m_{1}}(r_{1})\delta_{\ell_{2}m_{2}}(r_{2})\delta_{\ell_{3}m_{3}}(r_{3})\rangle$
$\displaystyle=\sqrt{(4\pi)^{3}(2l_{1}+1)(2l_{2}+1)(2l_{3}+1)}\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}$
$\displaystyle\times\sum_{L_{1}L_{2}L_{3}}f_{L_{1}L_{2}L_{3}}(r_{1},r_{2},r_{3})$
$\displaystyle(-1)^{L_{1}+L_{2}+L_{3}}\begin{Bmatrix}l_{1}&l_{2}&l_{3}\\\
L_{1}&L_{2}&L_{3}\end{Bmatrix}\begin{pmatrix}l_{1}&L_{1}&L_{3}\\\
0&0&0\end{pmatrix}\begin{pmatrix}l_{2}&L_{2}&L_{1}\\\
0&0&0\end{pmatrix}\begin{pmatrix}l_{3}&L_{3}&L_{2}\\\ 0&0&0\end{pmatrix}.$
(100)

The Wigner $3j$ symbols inside the sum over the $L_{i}$ impose that
$l_{1}+l_{2}+l_{3}$ be even, hence the sum is proportional to the symbol
$\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\ 0&0&0\end{pmatrix}$ , and
$\langle\delta_{\ell_{1}m_{1}}(r_{1})\delta_{\ell_{2}m_{2}}(r_{2})\delta_{\ell_{3}m_{3}}(r_{3})\rangle$,
to the Gaunt factor $\mathcal{G}^{\l_{1}l_{2}l_{3}}_{m_{1}m_{2}m_{3}}$, which
encodes the isotropy. The $3$-point function of the SFB modes is then obtained
by applying the basis transformation Eq. 1, hence it is real.

##  Appendix C An identity for integrating a product of three spherical
harmonics

Here we derive the identity

$\displaystyle I_{m_{1}m_{2}m_{3}}^{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})$
$\displaystyle\equiv\int
d^{2}{\hat{\bm{k}}}_{1}d^{2}{\hat{\bm{k}}}_{2}d^{2}{\hat{\bm{k}}}_{3}Y^{*}_{l_{1},m_{1}}({\hat{\bm{k}}}_{1})Y^{*}_{l_{2},m_{2}}({\hat{\bm{k}}}_{2})Y^{*}_{l_{3},m_{3}}({\hat{\bm{k}}}_{3})\delta_{D}({\bm{k}}_{1}+{\bm{k}}_{2}+{\bm{k}}_{3})$
$\displaystyle=\frac{4\pi^{\frac{3}{2}}}{k_{1}k_{2}k_{3}}\sqrt{(2l_{1}+1)}\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}\sum_{|m|\leq\mathrm{min}(l_{2},l_{3})}\\!\\!\\!\\!\\!{Y_{l_{2},m}}(\theta_{12},0)Y_{l_{3},-m}(\theta_{13},0)(-1)^{m}\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
0&m&-m\end{pmatrix}\,,$ (101)

which we use to calculate the SFB bispectrum in the homogeneous and isotropic
Universe, as well as to accelerate part of the calculation of the observed
bispectrum.

We need only to show this for $k_{1},k_{2},k_{3}$ which satisfy the triangle
inequality, as otherwise the integral clearly vanishes. Above, we define
$\theta_{ij}$ as the angle between ${\bm{k}}_{i}$ and ${\bm{k}}_{j}$, such
that
$\cos(\theta_{12})=\vartheta(k_{1},k_{2},k_{3})\equiv\frac{k_{3}^{2}-k_{1}^{2}-k_{2}^{2}}{2k_{1}k_{2}}$
and $\cos({\theta_{13}})=\vartheta(k_{1},k_{3},k_{2})$, and we denote for the
unit vector ${\hat{\bm{r}}}$ of spherical angles $(\theta,\phi)$,
$Y_{\ell,m}(\theta,\phi)\equiv Y_{\ell,m}({\hat{\bm{r}}})$. The angle-averaged
form of Eq. 101 is denoted $I_{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})$ such that

$\displaystyle
I_{m_{1}m_{2}m_{3}}^{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})=\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}I_{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})\,.$ (102)

We begin by using Eq. 71 and Eq. 73 to write

$\displaystyle\delta_{D}({\bm{k}}_{1}+{\bm{k}}_{2}+{\bm{k}}_{3})$
$\displaystyle=k_{3}^{-2}\delta_{D}(|{\bm{k}}_{1}+{\bm{k}}_{2}|-k_{3})\delta_{D}(\widehat{{\bm{k}}_{1}+{\bm{k}}_{2}}+{\hat{\bm{k}}}_{3})$
$\displaystyle=(k_{1}k_{2}k_{3})^{-1}\delta_{D}({\hat{\bm{k}}}_{1}\cdot{\hat{\bm{k}}}_{2}-\vartheta(k_{1},k_{2},k_{3}))\delta_{D}(\widehat{{\bm{k}}_{1}+{\bm{k}}_{2}}+{\hat{\bm{k}}}_{3})\,,$
(103)

such that, after integration over $d^{2}{\hat{\bm{k}}}_{3}$, Eq. 101 becomes

$\displaystyle
I_{m_{1}m_{2}m_{3}}^{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})=(k_{1}k_{2}k_{3})^{-1}(-1)^{l_{3}}\int
d^{2}{\hat{\bm{k}}}_{1}d^{2}{\hat{\bm{k}}}_{2}Y^{*}_{l_{1},m_{1}}({\hat{\bm{k}}}_{1})Y^{*}_{l_{2},m_{2}}({\hat{\bm{k}}}_{2})Y^{*}_{l_{3},m_{3}}((\widehat{{\bm{k}}_{1}+{\bm{k}}_{2}}))\delta_{D}({\hat{\bm{k}}}_{1}\cdot{\hat{\bm{k}}}_{2}-\vartheta(k_{1},k_{2},k_{3}))\,,$
(104)

where we used the parity property Eq. 89. We then integrate over
${\hat{\bm{k}}}_{2}$ by rotating it through an angle $\varphi_{2}$ around
${\hat{\bm{k}}}_{1}$, as the Dirac delta fixes
$\cos(\theta_{12})=\vartheta(k_{1},k_{2},k_{3})$. Then $\cos(\pi-\theta_{13})$
also remains fixed and $\widehat{{\bm{k}}_{1}+{\bm{k}}_{2}}$ rotates about
${\hat{\bm{k}}}_{1}$ by the same angle $\varphi_{2}$.

We denote by $\mathcal{R}({\hat{\bm{k}}}_{1})$ the rotation sending the axis
$\mathbf{\hat{z}}$ to ${\hat{\bm{k}}}_{1}$. Using the rotation formula for
spherical harmonics Eq. 81 and integrating over
${\hat{\bm{k}}}_{1}\cdot{\hat{\bm{k}}}_{2}$, Eq. 104 becomes

$\displaystyle
I_{m_{1}m_{2}m_{3}}^{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})=(k_{1}k_{2}k_{3})^{-1}(-1)^{l_{3}}\int
d^{2}{\hat{\bm{k}}}_{1}d\varphi_{2}Y^{*}_{l_{1},m_{1}}$
$\displaystyle({\hat{\bm{k}}}_{1})\bigg{[}\sum_{m_{2}^{\prime}=-l_{2}}^{l_{2}}D^{(l_{2})}_{m_{2},m_{2}^{\prime}}({\mathcal{R}({\hat{\bm{k}}}_{1})})Y^{*}_{l_{2},m_{2}^{\prime}}(\theta_{12},\varphi_{2})\bigg{]}$
$\displaystyle\times\bigg{[}\sum_{m_{3}^{\prime}=-l_{3}}^{l_{3}}D^{(l_{3})}_{m_{3},m_{3}^{\prime}}({\mathcal{R}({\hat{\bm{k}}}_{1})})Y^{*}_{l_{3},m_{3}^{\prime}}(\pi-\theta_{13},\varphi_{2})\bigg{]}\,,$
(105)

where $D^{(\ell)}_{m,m^{\prime}}(\mathcal{R}({\hat{\bm{k}}}_{1}))$ are Wigner
$D$-matrix elements. As $Y^{*}_{\ell,m}(\theta,\phi)$ is proportional to
$e^{-im\phi}$, integrating the pairwise products of spherical harmonics over
$d\varphi_{2}$ gives factors
$2\pi\delta^{K}_{m_{2}^{\prime},-m_{3}^{\prime}}$. Hence

$\displaystyle
I_{m_{1}m_{2}m_{3}}^{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})=2\pi(k_{1}k_{2}k_{3})^{-1}(-1)^{l_{3}}\sum_{m_{2}^{\prime}=-\mathrm{min}(l_{2},l_{3})}^{\mathrm{min}(l_{2},l_{3})}$
$\displaystyle
Y^{*}_{l_{2},m_{2}^{\prime}}(\theta_{12},0)Y^{*}_{l_{3},-m_{2}^{\prime}}(\pi-\theta_{13},0)$
$\displaystyle\times\int
d^{2}{\hat{\bm{k}}}_{1}Y^{*}_{l_{1},m_{1}}({\hat{\bm{k}}}_{1})D^{(l_{2})}_{m_{2},m_{2}^{\prime}}({\mathcal{R}({\hat{\bm{k}}}_{1})})D^{(l_{3})}_{m_{3},-m_{2}^{\prime}}({\mathcal{R}({\hat{\bm{k}}}_{1})})\,,$
(106)

where $D^{(\ell)}_{mm^{\prime}}({\mathcal{R}})$ is the Wigner $D$-matrix.
Using Eq. 83, the integral in Eq. 106 becomes

$\displaystyle\int
d^{2}{\hat{\bm{k}}}_{1}Y^{*}_{l_{1},m_{1}}({\hat{\bm{k}}}_{1})D^{(l_{2})}_{m_{2},m_{2}^{\prime}}({\mathcal{R}({\hat{\bm{k}}}_{1})})D^{(l_{3})}_{m_{3},-m_{2}^{\prime}}({\mathcal{R}({\hat{\bm{k}}}_{1})})$
$\displaystyle=\sum_{J=|l_{2}-l_{3}|}^{l_{2}+l_{3}}\langle
l_{2}m_{2}l_{3}m_{3}|J\left(m_{2}+m_{3}\right)\rangle\langle
l_{2}m_{2}^{\prime}l_{3}(-m_{2}^{\prime})|J0\rangle$ $\displaystyle\times\int
d^{2}{\hat{\bm{k}}}_{1}Y^{*}_{l_{1},m_{1}}({\hat{\bm{k}}}_{1})D^{J}_{m_{2}+m_{3},0}({\mathcal{R}({\hat{\bm{k}}}_{1})})$
$\displaystyle=(-1)^{m_{1}}\sqrt{\frac{4\pi}{2l_{1}+1}}\langle
l_{2}m_{2}l_{3}m_{3}|l_{1}\left(-m_{1}\right)\rangle\langle
l_{2}m_{2}^{\prime}l_{3}(-m_{2}^{\prime})|l_{1}0\rangle\,,$ (107)

where the second equation above follows from the first by applying the
identity Eq. 82 to write

$\displaystyle
D^{J}_{m_{2}+m_{3},0}({\mathcal{R}({\hat{\bm{k}}}_{1})})=\sqrt{\frac{4\pi}{2J+1}}Y^{*}_{J,m_{2}+m_{3}}({\hat{\bm{k}}}_{1})\,,$
(108)

and then using the orthogonality of spherical harmonics Eq. 80. We obtain the
final expression Eq. 101 by expressing the Clebsch-Gordan coefficients in
terms of Wigner $3j$’s via Eq. 84 and inserting Eq. 107 into Eq. 106. In
writing Eq. 101 we have also removed the complex conjugations from the
spherical harmonics as they are real, and used the parity of associated
Legendre polynomials Eq. 78.

It is also instructive to rewrite the integral Eq. 101 to make explicit the
consequence of isotropy. Using that

$\displaystyle\delta_{D}({\bm{k}}_{1}+{\bm{k}}_{2}+{\bm{k}}_{3})=\frac{1}{(2\pi)^{3}}\int\mathrm{d}^{3}{\bm{r}}e^{i({\bm{k}}_{1}+{\bm{k}}_{2}+{\bm{k}}_{3})\cdot{\bm{r}}}\,,$
(109)

and Eq. 86 (for $\alpha=0$), Eq. 101 becomes

$\displaystyle
I_{m_{1}m_{2}m_{3}}^{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})=8i^{l_{1}+l_{2}+l_{3}}\int
r^{2}\mathrm{d}rj_{l_{1}}(k_{1}r)j_{l_{2}}(k_{2}r)j_{l_{3}}(k_{3}r)\int
d^{2}{\hat{\bm{r}}}Y^{*}_{l_{1},m_{1}}({\hat{\bm{r}}})Y^{*}_{l_{2},m_{2}}({\hat{\bm{r}}})Y^{*}_{l_{3},m_{3}}({\hat{\bm{r}}})\,,$
(110)

where the angular integral is the Gaunt factor
$\mathcal{G}^{l_{1}l_{2}l_{3}}_{m_{1}m_{2}m_{3}}$ encoding the isotropy. The
above radial integral of the product of three spherical Bessel functions has
been evaluated analytically and by recursion in Mehrem  _et al._ (1991);
Fabrikant (2013); Dong  _et al._ (2015); Chellino and Slepian (2023), though
typically with methods requiring more computation than Eq. 101. During the
writing of this manuscript the authors discovered Ref. Jackson and Maximon
(1972), which also evaluated the radial integral, leading to a result
equivalent to Eq. 101 via an alternate derivation.

##  Appendix D Details of the bispectrum computation

###  D.1 Legendre expansion coefficients $Z_{l_{1}l_{2}l_{3}}(r,k_{1},k_{2})$

In this subsection we derive the coefficients
$Z_{l_{1}l_{2}l_{3}}(r,k_{1},k_{2})$ introduced in Eq. 66, which we reproduce
here:

$\displaystyle\tilde{Z}_{2}({\bm{k}}_{1},{\bm{k}}_{2},{\bm{r}})$
$\displaystyle=\sum_{l_{1}l_{2}l_{3}}Z_{l_{1}l_{2}l_{3}}(k_{1},k_{2},r)\,\mathcal{L}_{l_{1}}({\hat{\bm{k}}}_{1}\cdot{\hat{\bm{k}}}_{2})\,\mathcal{L}_{l_{2}}({\hat{\bm{k}}}_{1}\cdot{\hat{\bm{r}}})\,\mathcal{L}_{l_{3}}({\hat{\bm{k}}}_{2}\cdot{\hat{\bm{r}}})$
$\displaystyle=\sum_{\\{l_{i}m_{i}\\}}\frac{(4\pi)^{3}\,Z_{l_{1}l_{2}l_{3}}(k_{1},k_{2},r)}{\left(2l_{1}+1\right)\left(2l_{2}+1\right)\left(2l_{3}+1\right)}\,Y_{l_{1}m_{1}}({\hat{\bm{k}}}_{1})\,Y^{*}_{l_{1}m_{1}}({\hat{\bm{k}}}_{2})\,Y_{l_{2}m_{2}}({\hat{\bm{k}}}_{1})\,Y^{*}_{l_{2}m_{2}}({\hat{\bm{r}}})\,Y_{l_{3}m_{3}}({\hat{\bm{k}}}_{2})\,Y^{*}_{l_{3}m_{3}}({\hat{\bm{r}}})\,.$
(111)

We first note that the form of $Z_{2}({\bm{k}}_{1},{\bm{k}}_{2},{\bm{r}})$
(see Eq. 28) implies that it can be written as a polynomial in $\mu_{1}$,
$\mu_{2}$, ${\hat{\bm{k}}}_{1}\cdot{\hat{\bm{k}}}_{2}$, and
$\mu=(k_{1}\mu_{1}+k_{2}\mu_{2})/k$, except for the terms proportional to
$f_{\mathrm{NL}}$:

$\displaystyle
f_{\mathrm{NL}}\frac{\alpha(k)}{\alpha\left(k_{1}\right)\alpha\left(k_{2}\right)}\left(b^{\mathrm{E}}_{10}+f(r)\mu^{2}\right)\,.$
(112)

Furthermore, except for the term

$\displaystyle f(r)\mu^{2}G_{2}\left({\bm{k}}_{1},{\bm{k}}_{2}\right)\,,$
(113)

$Z_{2}({\bm{k}}_{1},{\bm{k}}_{2},{\bm{r}})$ is a polynomial of in $\mu_{1}$,
$\mu_{2}$ and ${\hat{\bm{k}}}_{1}\cdot{\hat{\bm{k}}}_{2}$ alone, which allows
us to write the decomposition Eq. 66. In fact, Eq. 113 could also be included
in this decomposition, by writing

$\displaystyle\mu^{2}=\frac{(\mu
k)^{2}}{k_{1}^{2}+k_{2}^{2}+2k_{1}k_{2}{\hat{\bm{k}}}_{1}\cdot{\hat{\bm{k}}}_{2}}=\frac{(k_{1}\mu_{1}+k_{2}\mu_{2})^{2}}{k_{1}^{2}+k_{2}^{2}}\sum_{n\geq
0}\bigg{(}-\frac{2k_{1}k_{2}}{k_{1}^{2}+k_{2}^{2}}{\hat{\bm{k}}}_{1}\cdot{\hat{\bm{k}}}_{2}\bigg{)}^{n}\,.$
(114)

However, then the sum over $l_{1},l_{2},l_{3}$ in Eq. 66 becomes infinite, and
truncating this slowly-converging series after even a small number of terms
greatly increases the number of non-vanishing coefficients
$Z_{l_{1}l_{2}l_{3}}$. Hence, we opt to treat Eq. 113 along with Eq. 112
separately (and exactly) as described in Appendix D.4.

Leaving out these terms, we have only $9$ coefficients
$Z_{l_{1}l_{2}l_{3}}(r,k_{1},k_{2})$, which are listed below by triplet
$(l_{1},l_{2},l_{3})$. We have dropped the implicit $r$-dependence for
brevity.

$\displaystyle(0,0,0)$
$\displaystyle:\quad\frac{1}{3}b^{\mathrm{E}}_{10}f+\frac{17}{21}b^{\mathrm{E}}_{10}+\frac{1}{2}b^{\mathrm{E}}_{20}+\frac{1}{9}f^{2}+\frac{b^{\mathrm{E}}_{01}f}{6\alpha(k_{1})}+\frac{b^{\mathrm{E}}_{01}f}{6\alpha(k_{2})}+\frac{b^{\mathrm{E}}_{11}}{2\alpha(k_{1})}+\frac{b^{\mathrm{E}}_{11}}{2\alpha(k_{2})}+\frac{b^{\mathrm{E}}_{02}}{\alpha(k_{1})\alpha(k_{2})}$
(115) $\displaystyle(0,0,2)$
$\displaystyle:\quad\frac{1}{9}f(3b^{\mathrm{E}}_{10}+2f)+f\frac{b^{\mathrm{E}}_{01}}{3\alpha(k_{1})}$
(116) $\displaystyle(0,1,1)$
$\displaystyle:\quad\bigg{(}\frac{k_{1}}{2k_{2}}+\frac{k_{2}}{2k_{1}}\bigg{)}f(b^{\mathrm{E}}_{10}+\frac{3}{5}f)+b^{\mathrm{E}}_{01}f\bigg{(}\frac{k_{2}}{2k_{1}\alpha(k_{2})}+\frac{k_{1}}{2k_{2}\alpha(k_{1})}\bigg{)}$
(117) $\displaystyle(0,1,3)$ $\displaystyle:\quad f^{2}\frac{k_{2}}{5k_{1}}$
(118) $\displaystyle(0,2,0)$
$\displaystyle:\quad\frac{1}{9}f(3b^{\mathrm{E}}_{10}+2f)+f\frac{b^{\mathrm{E}}_{01}}{3\alpha(k_{2})}$
(119) $\displaystyle(0,2,2)$ $\displaystyle:\quad\frac{4}{9}f^{2}$ (120)
$\displaystyle(0,3,1)$ $\displaystyle:\quad f^{2}\frac{k_{1}}{5k_{2}}$ (121)
$\displaystyle(1,0,0)$ $\displaystyle:\quad
b^{\mathrm{E}}_{10}\bigg{(}\frac{k_{1}}{2k_{2}}+\frac{k_{2}}{2k_{1}}\bigg{)}+b^{\mathrm{E}}_{01}\bigg{(}\frac{k_{2}}{2k_{1}\alpha(k_{2})}+\frac{k_{1}}{2k_{2}\alpha(k_{1})}\bigg{)}$
(122) $\displaystyle(2,0,0)$ $\displaystyle:\quad\frac{4}{21}$ (123)

###  D.2 Derivation of the observed SFB bispectrum Eq. 67

Here we detail the remainder of the derivation of the observed bispectrum Eq.
67 after the Legendre expansion of the second-order redshift-space kernel has
been performed as in Section D.1.

We first insert the decomposition of
$\tilde{Z}_{2}({\bm{k}}_{1},{\bm{k}}_{2},{\bm{r}})$ (Eq. 111) into
$\delta_{lm}^{g,(2)}$ (Eq. 64), and use the identity in Eq. 85 to simplify the
complex exponentials. We rid of the angular integrals over ${\hat{\bm{q}}}$
with the orthogonality relation for spherical harmonics Eq. 80. After
rearranging for the angular and radial integrals and assuming a separable
window $W({\bm{r}})=W(r)W({\hat{\bm{r}}})$, we obtain

$\displaystyle\delta_{\ell m}^{g,(2)}(k)$
$\displaystyle=\sqrt{\frac{2}{\pi}}\,k\int\mathrm{d}r\,r^{2}\,j_{\ell}(kr)\,W(r)\,D^{2}(r)\,4\pi\int\frac{\mathrm{d}k_{1}\,k_{1}}{(2\pi)^{3}}\int\frac{\mathrm{d}k_{2}\,k_{2}}{(2\pi)^{3}}\sum_{\\{l_{i}m_{i}\\}}\frac{(4\pi)^{3}\,Z_{l_{1}l_{2}l_{3}}(k_{1},k_{2},r)}{\left(2l_{1}+1\right)\left(2l_{2}+1\right)\left(2l_{3}+1\right)}$
$\displaystyle\quad\times
r^{2}\,(4\pi)^{3}\sum_{\\{L_{i}M_{i}\\}}\sum_{LM}i^{-L_{1}-L_{2}+L_{3}+L_{4}}\,j_{L_{3}}(k_{1}r)\,j_{L_{4}}(k_{2}r)\,\delta^{(1)}_{L_{1}M_{1}}(k_{1})\,\delta^{(1)}_{L_{2}M_{2}}(k_{2})$
$\displaystyle\quad\times\frac{\pi}{2r^{2}}\,\mathcal{G}^{L_{3}}_{M_{3}}{}^{L_{4}}_{M_{4}}{}^{L}_{M}\int\mathrm{d}^{2}{\hat{\bm{r}}}\,W({\hat{\bm{r}}})\,Y_{LM}({\hat{\bm{r}}})\,Y^{*}_{\ell
m}({\hat{\bm{r}}})\,Y^{*}_{l_{2}m_{2}}({\hat{\bm{r}}})\,Y^{*}_{l_{3}m_{3}}({\hat{\bm{r}}})$
$\displaystyle\quad\times\int\mathrm{d}^{2}{\hat{\bm{k}}}_{1}\,Y_{l_{1}m_{1}}({\hat{\bm{k}}}_{1})\,Y_{l_{2}m_{2}}({\hat{\bm{k}}}_{1})\,Y_{L_{3}M_{3}}({\hat{\bm{k}}}_{1})\,Y_{L_{1}M_{1}}({\hat{\bm{k}}}_{1})\int\mathrm{d}^{2}{\hat{\bm{k}}}_{2}\,Y^{*}_{l_{1}m_{1}}({\hat{\bm{k}}}_{2})\,Y_{l_{3}m_{3}}({\hat{\bm{k}}}_{2})\,Y_{L_{4}M_{4}}({\hat{\bm{k}}}_{2})\,Y_{L_{2}M_{2}}({\hat{\bm{k}}}_{2})\,.$
(124)

We may rewrite this more compactly as

$\displaystyle\delta_{\ell m}^{g,(2)}(k)$
$\displaystyle=\int\mathrm{d}k_{1}\int\mathrm{d}k_{2}\sum_{L_{1}M_{1}}\sum_{L_{2}M_{2}}\mathcal{V}^{\ell
L_{1}L_{2}}_{mM_{1}M_{2}}(k,k_{1},k_{2})\,\delta^{(1)}_{L_{1}M_{1}}(k_{1})\,\delta^{(1)}_{L_{2}M_{2}}(k_{2})\,,$
(125)

where we use $\mathcal{V}$ to denote the second-order coupling kernel in the
SFB bispectrum

$\displaystyle\mathcal{V}^{\ell L_{1}L_{2}}_{mM_{1}M_{2}}(k,k_{1},k_{2})$
$\displaystyle\equiv\sqrt{\frac{\pi}{2}}\,2^{8}\pi\,kk_{1}k_{2}\int\mathrm{d}r\,r^{2}\,j_{\ell}(kr)\,W(r)\,D^{2}(r)\sum_{l_{1}l_{2}l_{3}}\frac{Z_{l_{1}l_{2}l_{3}}(k_{1},k_{2},r)}{\left(2l_{1}+1\right)\left(2l_{2}+1\right)\left(2l_{3}+1\right)}$
$\displaystyle\quad\times\sum_{L_{3}L_{4}}\,j_{L_{3}}(k_{1}r)\,j_{L_{4}}(k_{2}r)\,\mathcal{C}^{L_{1}}_{M_{1}}{}^{L_{2}}_{M_{2}}{}^{\ell}_{m}{}^{,L_{3}L_{4}}_{,l_{1}l_{2}l_{3}}\,.$
(126)

Here
${C}^{L_{1}}_{M_{1}}{}^{L_{2}}_{M_{2}}{}^{\ell}_{m}{}^{,L_{3}L_{4}}_{,l_{1}l_{2}l_{3}}$
is a mode coupling coefficient

$\displaystyle\mathcal{C}^{L_{1}}_{M_{1}}{}^{L_{2}}_{M_{2}}{}^{\ell}_{m}{}^{,L_{3}L_{4}}_{,l_{1}l_{2}l_{3}}$
$\displaystyle=i^{-L_{1}-L_{2}+L_{3}+L_{4}}\sum_{m_{1}m_{2}m_{3}}\sum_{M_{3}M_{4}}\sum_{LM}(-1)^{M+m_{1}}\mathcal{G}^{L_{3}L_{4}L}_{M_{3}M_{4}M}\mathcal{H}^{L\ell
l_{2}l_{3}}_{-Mmm_{2}m_{3}}\mathcal{H}^{l_{1}l_{2}L_{3}L_{1}}_{m_{1}m_{2}M_{3}M_{1}}\mathcal{H}^{l_{1}l_{3}L_{4}L_{2}}_{-m_{1}m_{3}M_{4}M_{2}}\,.$
(127)

and

$\displaystyle\mathcal{H}^{l_{1}l_{2}l_{3}l_{4}}_{m_{1}m_{2}m_{3}m_{4}}\equiv\int\mathrm{d}^{2}{\hat{\bm{r}}}\
W({\hat{\bm{r}}})Y_{l_{1}m_{1}}({\hat{\bm{r}}})Y_{l_{2}m_{2}}({\hat{\bm{r}}})Y_{l_{3}m_{3}}({\hat{\bm{r}}})Y_{l_{4}m_{4}}({\hat{\bm{r}}})\,.$
(128)

is the integral over four spherical harmonics and the angular part of the
window function.

Using the form of $\delta^{g,(1)}_{\ell m}(k)$ in Eq. 10, the terms
contributing to the tree-level $3$-point correlation function of the SFB modes
are

$\displaystyle\left\langle\delta_{\ell
m}^{g,(2)}(k)\,\delta^{g,(1)}_{\ell^{\prime}m^{\prime}}(k^{\prime})\,\delta^{g,(1)}_{\ell^{\prime\prime}m^{\prime\prime}}(k^{\prime\prime})\right\rangle$
$\displaystyle=\int\mathrm{d}k_{1}\int\mathrm{d}k_{2}\sum_{L_{1}M_{1}}\sum_{L_{2}M_{2}}\mathcal{V}^{\ell
L_{1}L_{2}}_{mM_{1}M_{2}}(k,k_{1},k_{2})\int\mathrm{d}q^{\prime}\sum_{L^{\prime}M^{\prime}}\mathcal{W}_{\ell^{\prime}m^{\prime}}^{L^{\prime}M^{\prime}}(k^{\prime},q^{\prime})\int\mathrm{d}q^{\prime\prime}\sum_{L^{\prime\prime}M^{\prime\prime}}\mathcal{W}_{\ell^{\prime\prime}m^{\prime\prime}}^{L^{\prime\prime}M^{\prime\prime}}(k^{\prime\prime},q^{\prime\prime})$
$\displaystyle\quad\times\left\langle\delta^{(1)}_{L_{1}M_{1}}(k_{1})\,\delta^{(1)}_{L_{2}M_{2}}(k_{2})\,\delta^{(1)}_{L^{\prime}M^{\prime}}(q^{\prime})\,\delta^{(1)}_{L^{\prime\prime}M^{\prime\prime}}(q^{\prime\prime})\right\rangle\,,$
(129)

along with the two other terms with cyclically permuted superscript indices.
Noting that the SFB power spectrum for the constant-time slice matter density
contrast is homogeneous and isotropic, i.e.

$\displaystyle\left\langle\delta^{(1)}_{L_{1}M_{1}}(k_{1})\,\delta^{(1)}_{L_{2}M_{2}}(k_{2})\right\rangle$
$\displaystyle=\delta^{K}_{L_{1}L_{2}}\delta^{K}_{M_{1}-M_{2}}(-1)^{M_{2}}\delta^{D}(k_{1}-k_{2})\,P(k_{1})\,,$
(130)

we may apply Wick’s theorem,

$\displaystyle\left\langle\delta_{\ell
m}^{g,(2)}(k)\,\delta^{g,(1)}_{\ell^{\prime}m^{\prime}}(k^{\prime})\,\delta^{g,(1)}_{\ell^{\prime\prime}m^{\prime\prime}}(k^{\prime\prime})\right\rangle$
$\displaystyle=\int\mathrm{d}k_{1}\sum_{L_{1}M_{1}}(-1)^{M_{1}}\mathcal{V}^{\ell
L_{1}L_{1}}_{mM_{1}-M_{1}}(k,k_{1},k_{1})\,P(k_{1})\int\mathrm{d}q^{\prime}\sum_{L^{\prime}M^{\prime}}(-1)^{M^{\prime}}\mathcal{W}_{\ell^{\prime}m^{\prime}}^{L^{\prime}M^{\prime}}(k^{\prime},q^{\prime})\,\mathcal{W}_{\ell^{\prime\prime}m^{\prime\prime}}^{L^{\prime}-M^{\prime}}(k^{\prime\prime},q^{\prime})\,P(q^{\prime})$
$\displaystyle\quad+\int\mathrm{d}q^{\prime}\int\mathrm{d}q^{\prime\prime}\sum_{L^{\prime}M^{\prime}}\sum_{L^{\prime\prime}M^{\prime\prime}}(-1)^{M^{\prime}+M^{\prime\prime}}\mathcal{V}^{\ell
L^{\prime}L^{\prime\prime}}_{m-M^{\prime}-M^{\prime\prime}}(k,q^{\prime},q^{\prime\prime})\,\mathcal{W}_{\ell^{\prime}m^{\prime}}^{L^{\prime}M^{\prime}}(k^{\prime},q^{\prime})\,\mathcal{W}_{\ell^{\prime\prime}m^{\prime\prime}}^{L^{\prime\prime}M^{\prime\prime}}(k^{\prime\prime},q^{\prime\prime})\,P(q^{\prime})\,P(q^{\prime\prime})$
$\displaystyle\quad+\int\mathrm{d}q^{\prime}\int\mathrm{d}q^{\prime\prime}\sum_{L^{\prime\prime}M^{\prime\prime}}\sum_{L^{\prime}M^{\prime}}(-1)^{M^{\prime}+M^{\prime\prime}}\mathcal{V}^{\ell
L^{\prime\prime}L^{\prime}}_{m-M^{\prime\prime}-M^{\prime}}(k,q^{\prime\prime},q^{\prime})\,\mathcal{W}_{\ell^{\prime}m^{\prime}}^{L^{\prime}M^{\prime}}(k^{\prime},q^{\prime})\,\mathcal{W}_{\ell^{\prime\prime}m^{\prime\prime}}^{L^{\prime\prime}M^{\prime\prime}}(k^{\prime\prime},q^{\prime\prime})\,P(q^{\prime})\,P(q^{\prime\prime})\,.$
(131)

From here on we assume a spherically symmetric window $W({\bm{r}})=W(r)$. Then
we have $\mathcal{W}_{\ell m}^{LM}(k,q)=\delta^{K}_{\ell
L}\delta^{K}_{mM}\,\mathcal{W}_{\ell}(k,q)$. Hence,

$\displaystyle\left\langle\delta_{\ell
m}^{g,(2)}(k)\,\delta^{g,(1)}_{\ell^{\prime}m^{\prime}}(k^{\prime})\,\delta^{g,(1)}_{\ell^{\prime\prime}m^{\prime\prime}}(k^{\prime\prime})\right\rangle$
$\displaystyle=\delta^{K}_{\ell^{\prime}\ell^{\prime\prime}}\delta^{K}_{m^{\prime}-m^{\prime\prime}}(-1)^{m^{\prime}}\,C_{\ell^{\prime}}(k^{\prime},k^{\prime\prime})\int\mathrm{d}k_{1}\sum_{L_{1}M_{1}}(-1)^{M_{1}}\mathcal{V}^{\ell
L_{1}L_{1}}_{mM_{1}-M_{1}}(k,k_{1},k_{1})\,P(k_{1})$
$\displaystyle\quad+\int\mathrm{d}q^{\prime}\,\mathcal{W}_{\ell^{\prime}}(k^{\prime},q^{\prime})\,P(q^{\prime})\int\mathrm{d}q^{\prime\prime}\,\mathcal{W}_{\ell^{\prime\prime}}(k^{\prime\prime},q^{\prime\prime})\,P(q^{\prime\prime})(-1)^{m^{\prime}+m^{\prime\prime}}\left[\mathcal{V}^{\ell\ell^{\prime}\ell^{\prime\prime}}_{m-m^{\prime}-m^{\prime\prime}}(k,q^{\prime},q^{\prime\prime})+\mathcal{V}^{\ell\ell^{\prime\prime}\ell^{\prime}}_{m-m^{\prime\prime}-m^{\prime}}(k,q^{\prime\prime},q^{\prime})\right].$
(132)

In Eq. 125 swapping $L_{1}\leftrightarrow L_{2},k_{1}\leftrightarrow
k_{2},M_{1}\leftrightarrow M_{2}$ gives $\mathcal{V}^{\ell
L_{1}L_{2}}_{mM_{1}M_{2}}(k,k_{1},k_{2})=\mathcal{V}^{\ell
L_{2}L_{1}}_{mM_{2}M_{1}}(k,k_{2},k_{1})$, thus the two terms in brackets in
Eq. 132 are equal. In Appendix D.3 we simplify
$\mathcal{V}^{\ell\ell^{\prime}\ell^{\prime\prime}}_{m-m^{\prime}-m^{\prime\prime}}$
and show that it is proportional to the Gaunt factor
$\mathcal{G}^{\ell\ell^{\prime}\ell^{\prime\prime}}_{mm^{\prime}m^{\prime\prime}}$.
Hence we define the angle-averaged quantity

$\displaystyle\mathcal{V}^{\ell\ell^{\prime}\ell^{\prime\prime}}(k,k^{\prime},k^{\prime\prime})\equiv\sum_{mm^{\prime}m^{\prime\prime}}\begin{pmatrix}\ell&&\ell^{\prime}&&\ell^{\prime\prime}\\\
m&&m^{\prime}&&m^{\prime\prime}\end{pmatrix}(-1)^{m^{\prime}+m^{\prime\prime}}\mathcal{V}^{\ell\ell^{\prime}\ell^{\prime\prime}}_{m-m^{\prime}-m^{\prime\prime}}(k,k^{\prime},k^{\prime\prime})\,,$
(133)

whose simplified expression is given by Eq. 143.

In fact, the first line of Eq. 132 always vanishes when $\ell>0$, which we may
see by the following quick argument. By the remarks at the beginning of
Section II.3 on observational isotropy, Eq. 132 must be proportional to the
Gaunt factor
$\mathcal{G}^{\ell\ell^{\prime}\ell^{\prime\prime}}_{mm^{\prime}m^{\prime\prime}}$,
such that Eq. 132 is nonzero only if $\ell,\ell^{\prime},\ell^{\prime\prime}$
form a triangle. This condition is already imposed by the second line of Eq.
132. Let $\ell>0$ and set $\ell^{\prime}=\ell^{\prime\prime}=0$ such that the
triangle condition is violated. As $C_{0}(k^{\prime},k^{\prime\prime})\neq 0$,
the integral in the first line of Eq. 132 must vanish, and consequently the
first line must vanish for any $(\ell,\ell^{\prime},\ell^{\prime\prime})$ with
$\ell>0$. The authors of Ref. Bertacca  _et al._ (2018) demonstrate this
directly with a lengthy derivation.

Finally, after angle-averaging with Eq. 45 and re-indexing for clarity, the
SFB bispectrum, ignoring contributions from the $f_{\mathrm{NL}}$ and $G_{2}$
terms, is given by

$\displaystyle
B^{\text{SFB,}\,\tilde{Z}_{2}}_{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})$
$\displaystyle=2\int\mathrm{d}q_{2}\mathcal{W}_{l_{2}}(k_{2},q_{2})\,P(q_{2})\int\mathrm{d}q_{3}\,\mathcal{W}_{l_{3}}(k_{3},q_{3})\,P(q_{3})\
\mathcal{V}^{l_{1}l_{2}l_{3}}(k_{1},q_{2},q_{3})+\mathrm{2\ cyc.\ perm.\,.}$
(134)

The contribution from the $f_{\mathrm{NL}}$ and $G_{2}$ terms is given in
Appendix D.4 (Eq. 156). It is of a similar form to the contribution from the
terms in $\tilde{Z}_{2}$.

In principle, it is possible to repeat the above derivation while relaxing the
assumption that $W({\bm{r}})$ is spherically symmetric (but keeping the
assumption that the radial and angular dependencies are separable) by
decomposing $W({\hat{\bm{r}}})$ into spherical harmonics. However, as in this
case we can no longer leverage observational isotropy of the signal, the
calculation is significantly more complicated (and expensive), so we leave the
details to a future work.

Lastly, it is useful to verify that we retrieve Eq. 55 in the limit of an
isotropic and homogeneous Universe, i.e by setting $D=b_{1}=W=1$ and $f=0$. In
this case, the second term in Eq. 67 vanishes, and we may take
$W_{\ell}(k,q)=\delta_{D}(k-q)$, such that

$\displaystyle
B^{\text{SFB}}_{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})=2P(k_{2})P(k_{3})\
\mathcal{V}^{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})+\mathrm{2\ cyc.\ perm.}\,.$
(135)

As it is somewhat tedious to demonstrate equivalence with Eq. 55 analytically,
we omit the details here666For a lengthy proof along these lines see Ref.
Bertacca  _et al._ (2018).; as a test of our code, we verify numerically that
Eq. 135 and Eq. 55 are identical in the considered limit.

###  D.3 Simplification of the kernel
$\mathcal{V}^{L_{1}L_{2}\ell}_{M_{1}M_{2}m}(k,k_{1},k_{2})$

We now show that, under the assumption of a spherically symmetric window
$W({\bm{r}})=W(r)$, the kernel $\mathcal{V}^{\ell L_{1}L_{2}}_{mM_{1}M_{2}}$
is proportional to the Gaunt factor
$\mathcal{G}^{L_{1}L_{2}\ell}_{M_{1}M_{2}-m}$, and compute its angle-averaged
expression, defined by Eq. 133. For compactness, we will denote the Wigner
coefficient $\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}$ by
$K^{\ell_{1}\ell_{2}\ell_{3}}_{m_{1}m_{2}m_{3}}$ and the coefficient of
proportionality between Gaunt factors and Wigner coefficients by

$\displaystyle f_{\ell_{1}\ell_{2}\ell_{3}}$
$\displaystyle\equiv\left(\frac{(2\ell_{1}+1)(2\ell_{2}+1)(2\ell_{3}+1)}{4\pi}\right)^{\frac{1}{2}}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
0&0&0\end{pmatrix}\,,$ (136)

such that
$\mathcal{G}^{\ell_{1}\ell_{2}\ell_{3}}_{m_{1}m_{2}m_{3}}=f_{\ell_{1}\ell_{2}\ell_{3}}K^{\ell_{1}\ell_{2}\ell_{3}}_{m_{1}m_{2}m_{3}}$.

First, we use the identity Eq. 93 to evaluate the integral over four spherical
harmonics Eq. 128, such that we may express Eq. 127 in terms of Gaunt
factors777Note that if we allow for a generic angular dependence of the window
$W({\hat{\bm{r}}})$, it is still possible to evaluate the integral Eq. 93 by
decomposing $W({\hat{\bm{r}}})$ into spherical harmonics and using the
generalization of the integral Eq. 127 to a product of five spherical
harmonics. However this leads to an explosion in the number of terms needed to
compute the kernel $\mathcal{V}^{\ell L_{1}L_{2}}_{mM_{1}M_{2}}$, so in
practice only spherically symmetric windows are currently computationally
feasible.. We also use that the sum of the lower indices of the Gaunt factor
must vanish, to write

$\displaystyle\mathcal{C}^{L_{1}}_{M_{1}}{}^{L_{2}}_{M_{2}}{}^{\ell}_{m}{}^{,L_{3}L_{4}}_{,l_{1}l_{2}l_{3}}$
$\displaystyle=i^{-L_{1}-L_{2}+L_{3}+L_{4}}\sum_{L,H_{1},H_{2},H_{3}}\sum_{m_{1},m_{2},m_{3},M_{3},M_{4},M}(-1)^{m+m_{1}+m_{2}+m_{3}}$
$\displaystyle\quad\times\mathcal{G}^{L_{3}L_{4}L}_{M_{3}M_{4}M}\mathcal{G}^{L\ell
H_{1}}_{-Mm-
N_{1}}\mathcal{G}^{H_{1}l_{2}l_{3}}_{N_{1}m_{2}m_{3}}\mathcal{G}^{l_{1}l_{2}H_{2}}_{m_{1}m_{2}-N_{2}}\mathcal{G}^{H_{2}L_{3}L_{1}}_{N_{2}M_{3}M_{1}}\mathcal{G}^{l_{1}l_{3}H_{3}}_{-m_{1}m_{3}-N_{3}}\mathcal{G}^{H_{3}L_{4}L_{2}}_{N_{3}M_{4}M_{2}}\,,$
(137)

where the sums over $N_{1}$, $N_{2}$, and $N_{3}$ have only one term each, and
are thus not explicitly written. With Eq. 95, and changing signs on the
summation variables, we can use Eq. 96 to simplify the inner sum

$\displaystyle\sum_{m_{1}m_{2}m_{3}}(-1)^{m_{1}+m_{2}+m_{3}}\mathcal{G}^{H_{1}l_{2}l_{3}}_{N_{1}m_{2}m_{3}}\mathcal{G}^{l_{1}l_{2}H_{2}}_{m_{1}m_{2}-N_{2}}\mathcal{G}^{l_{1}l_{3}H_{3}}_{-m_{1}m_{3}-N_{3}}$
$\displaystyle=f_{H_{1}l_{2}l_{3}}f_{H_{2}l_{1}l_{2}}f_{H_{3}l_{1}l_{3}}(-1)^{H_{1}+l_{2}+l_{3}+H_{3}+l_{1}+l_{3}}\sum_{m_{1}m_{2}m_{3}}(-1)^{m_{1}+m_{2}+m_{3}}\,K^{H_{1}l_{2}l_{3}}_{-N_{1}m_{2}-m_{3}}K^{l_{1}H_{3}l_{3}}_{-m_{1}-N_{3}m_{3}}K^{l_{1}l_{2}H_{2}}_{m_{1}-m_{2}-N_{2}}$
$\displaystyle=f_{H_{1}l_{2}l_{3}}f_{H_{2}l_{1}l_{2}}f_{H_{3}l_{1}l_{3}}(-1)^{l_{1}+l_{2}+l_{3}}K^{H_{1}H_{2}H_{3}}_{-N_{1}-N_{2}-N_{3}}\begin{Bmatrix}H_{1}&H_{2}&H_{3}\\\
l_{1}&l_{3}&l_{2}\end{Bmatrix}\,,$ (138)

where we also used that $l_{1}+l_{2}+H_{2}$ must be even. Then, after
expressing all Gaunt factors in terms of Wigner coefficients, Eq. 137 becomes

$\displaystyle\mathcal{C}^{L_{1}}_{M_{1}}{}^{L_{2}}_{M_{2}}{}^{\ell}_{m}{}^{,L_{3}L_{4}}_{,l_{1}l_{2}l_{3}}$
$\displaystyle=i^{-L_{1}-L_{2}+L_{3}+L_{4}}\sum_{L,H_{1},H_{2},H_{3}}f_{H_{1}l_{2}l_{3}}f_{H_{2}l_{1}l_{2}}f_{H_{3}l_{1}l_{3}}(-1)^{l_{1}+l_{2}+l_{3}}\begin{Bmatrix}H_{1}&H_{2}&H_{3}\\\
l_{1}&l_{3}&l_{2}\end{Bmatrix}(-1)^{m}$ $\displaystyle\times
f_{L_{3}L_{4}L}f_{L\ell
H_{1}}f_{H_{2}L_{3}L_{1}}f_{H_{3}L_{4}L_{2}}\sum_{M_{3},M_{4},M}{K}^{L_{3}L_{4}L}_{M_{3}M_{4}M}{K}^{L\ell
H_{1}}_{-Mm-
N_{1}}{K}^{H_{2}L_{3}L_{1}}_{N_{2}M_{3}M_{1}}{K}^{H_{3}L_{4}L_{2}}_{N_{3}M_{4}M_{2}}K^{H_{1}H_{2}H_{3}}_{-N_{1}-N_{2}-N_{3}}$
$\displaystyle=i^{-L_{1}-L_{2}+L_{3}+L_{4}}{K}^{L_{1}L_{2}\ell}_{M_{1}M_{2}-m}(-1)^{m}(-1)^{l_{1}+l_{2}+l_{3}}$
$\displaystyle\times\sum_{L,H_{1},H_{2},H_{3}}f_{H_{1}l_{2}l_{3}}f_{H_{2}l_{1}l_{2}}f_{H_{3}l_{1}l_{3}}f_{L_{3}L_{4}L}f_{L\ell
H_{1}}f_{H_{2}L_{3}L_{1}}f_{H_{3}L_{4}L_{2}}(-1)^{H_{1}+H_{2}+H_{3}}$
$\displaystyle\times\begin{Bmatrix}H_{1}&H_{2}&H_{3}\\\
l_{1}&l_{3}&l_{2}\end{Bmatrix}\begin{Bmatrix}H_{2}&L_{3}&L_{1}\\\
H_{3}&L_{4}&L_{2}\\\ H_{1}&L&\ell\end{Bmatrix}\,,$ (139)

where we used Eq. 97 to obtain the last line.

Finally, substituting Eq. 139 in Eq. 126 and collecting constant factors, we
obtain:

$\displaystyle\mathcal{V}^{\ell
L_{1}L_{2}}_{mM_{1}M_{2}}(k,k_{1},k_{2})\equiv(32\pi)^{\frac{3}{2}}kk_{1}k_{2}\begin{pmatrix}L_{1}&L_{2}&\ell\\\
M_{1}&M_{2}&-m\end{pmatrix}(-1)^{m}\sum_{l_{1}l_{2}l_{3}L_{3}L_{4}}g^{L_{1}L_{2}\ell}_{l_{1}l_{2}l_{3}L_{3}L_{4}}J^{\ell
L_{3}L_{4}}_{l_{1}l_{2}l_{3}}(k,k_{1},k_{2})\,,$ (140)

where we have defined

$\displaystyle J^{\ell
L_{3}L_{4}}_{l_{1}l_{2}l_{3}}(k,k_{1},k_{2})\equiv\int\mathrm{d}r\,r^{2}\,j_{\ell}(kr)\,j_{L_{3}}(k_{1}r)\,j_{L_{4}}(k_{2}r)\,W(r)\,D^{2}(r)\,Z_{l_{1}l_{2}l_{3}}(k_{1},k_{2},r)\,,$
(141)

and

$\displaystyle
g^{L_{1}L_{2}\ell}_{l_{1}l_{2}l_{3}L_{3}L_{4}}\equiv\frac{(-1)^{l_{1}+l_{2}+l_{3}}}{\left(2l_{1}+1\right)\left(2l_{2}+1\right)\left(2l_{3}+1\right)}i^{-L_{1}-L_{2}+L_{3}+L_{4}}\sum_{L,H_{1},H_{2},H_{3}}f_{H_{1}l_{2}l_{3}}f_{H_{2}l_{1}l_{2}}f_{H_{3}l_{1}l_{3}}f_{L_{3}L_{4}L}f_{L\ell
H_{1}}f_{H_{2}L_{3}L_{1}}f_{H_{3}L_{4}L_{2}}$
$\displaystyle\times\begin{Bmatrix}H_{1}&H_{2}&H_{3}\\\
l_{1}&l_{3}&l_{2}\end{Bmatrix}\begin{Bmatrix}H_{2}&L_{3}&L_{1}\\\
H_{3}&L_{4}&L_{2}\\\ H_{1}&L&\ell\end{Bmatrix}\,,$ (142)

which is real. There are $9$ triplets $(l_{1},l_{2},l_{3})$ in the Legendre
decomposition of the kernel $Z_{2}$, excluding contributions from the kernel
$G_{2}$ and from terms proportional to $f_{\mathrm{NL}}$. As a result, for
fixed $L_{1},L_{2},\ell$ there are at most $49$ terms in the sum of Eq. 140.
Finally, angle-averaging with Eq. 133, we obtain:

$\displaystyle\mathcal{V}^{\ell
L_{1}L_{2}}(k,k_{1},k_{2})\equiv(32\pi)^{\frac{3}{2}}kk_{1}k_{2}\sum_{l_{1}l_{2}l_{3}L_{3}L_{4}}g^{L_{1}L_{2}\ell}_{l_{1}l_{2}l_{3}L_{3}L_{4}}J^{\ell
L_{3}L_{4}}_{l_{1}l_{2}l_{3}}(k,k_{1},k_{2})\,.$ (143)

###  D.4 Contribution from the $f_{\mathrm{NL}}$ and $G_{2}$ terms

Here we address the terms in the kernel $Z_{2}$ (Eq. 28) that were left out
when decomposing $Z_{2}$ as polynomial in
${\hat{\bm{k}}}_{1}\cdot{\hat{\bm{k}}}_{2}$,
${\hat{\bm{k}}}_{1}\cdot{\hat{\bm{r}}}$, and
${\hat{\bm{k}}}_{2}\cdot{\hat{\bm{r}}}$ (Eq. 66). These terms are

$\displaystyle
f_{\mathrm{NL}}\frac{\alpha(k,r)}{\alpha\left(k_{1},r\right)\alpha\left(k_{2},r\right)}\left(b^{\mathrm{E}}_{10}+f(r)\mu^{2}\right)\,,$
(144)

and

$\displaystyle f(r)\mu^{2}G_{2}\left({\bm{k}}_{1},{\bm{k}}_{2}\right)\,.$
(145)

To account for them, it will be advantageous to first write the SFB bispectrum
in an alternate form (Eq. 147), which makes clear the relation of the observed
bispectrum to the bispectrum of an isotropic and homogeneous Universe.

####  D.4.1 Relation between position-dependent bispectrum and SFB bispectrum

The SFB bispectrum is obtained from the position-dependent Fourier-space
bispectrum by first transforming Eq. 31 to configuration space using Eq. 61,
and then transforming into SFB space using Eq. 1. We get

$\displaystyle\left\langle\delta^{\mathrm{obs}}_{g,l_{1}m_{1}}(k_{1})\delta^{\mathrm{obs}}_{g,l_{2}m_{2}}(k_{2})\delta^{\mathrm{obs}}_{g,l_{3}m_{3}}(k_{3})\right\rangle=\bigg{(}\frac{2}{\pi}\bigg{)}^{\frac{3}{2}}k_{1}k_{2}k_{3}\int\bigg{(}\prod_{i}\frac{1}{(2\pi)^{3}}r_{i}^{2}\mathrm{d}r_{i}q_{i}^{2}\mathrm{d}q_{i}\mathrm{d}^{2}{\hat{\bm{r}}}_{i}\mathrm{d}^{2}{\hat{\bm{q}}}_{i}W({\bm{r}}_{i})e^{i{\bm{q}}_{i}\cdot{\bm{r}}_{i}}j_{l_{i}}(k_{i}r_{i})Y^{*}_{l_{i},m_{i}}({\hat{\bm{r}}}_{i})\bigg{)}$
$\displaystyle\times
B_{s}({\bm{q}}_{1},{\bm{q}}_{2},{\bm{q}}_{3},{\bm{r}}_{1},{\bm{r}}_{2},{\bm{r}}_{3})(2\pi)^{3}\delta_{D}({\bm{q}}_{1}+{\bm{q}}_{2}+{\bm{q}}_{3})\,.$
(146)

In the absence of RSD, linear growth, galaxy bias, and window, the observed
SFB bispectrum (Eq. 146) reduces to the SFB bispectrum in an isotropic and
homogeneous Universe (Eq. 50), using Eq. 76. Unlike in the isotropic and
homogeneous case, however, fixing the lengths $q_{i}$ and imposing
${\bm{q}}_{1}+{\bm{q}}_{2}+{\bm{q}}_{3}=0$ fixes the angles
${\hat{\bm{q}}}_{i}\cdot{\hat{\bm{q}}}_{j}$ but does not determine
$B_{s}({\bm{q}}_{1},{\bm{q}}_{2},{\bm{q}}_{3},{\bm{r}}_{1},{\bm{r}}_{2},{\bm{r}}_{3})$,
which depends on the $9$ angles
$\mu_{ij}\equiv{\hat{\bm{q}}}_{i}\cdot{\hat{\bm{r}}}_{j}$ to the three lines
of sight ${\hat{\bm{r}}}_{j}$.

Assuming a radial window $W({\bm{r}})=W(r)$, the angle-averaged bispectrum is
given by

$\displaystyle B^{\text{SFB,
obs}}_{l_{1}l_{2}l_{3}}(k_{1},k_{2},k_{3})=\frac{1}{(2\pi)^{6}}\bigg{(}\frac{2}{\pi}\bigg{)}^{\frac{3}{2}}k_{1}k_{2}k_{3}\int\bigg{(}\prod_{i}r_{i}^{2}\mathrm{d}r_{i}q_{i}^{2}\mathrm{d}q_{i}W(r_{i})j_{l_{i}}(k_{i}r_{i})\bigg{)}\mathcal{I}_{l_{1}l_{2}l_{3}}^{\text{ang.}}(q_{1},q_{2},q_{3},r_{1},r_{2},r_{3})\,,$
(147)

where we have defined the angle-averaged angular integral

$\displaystyle\mathcal{I}_{l_{1}l_{2}l_{3}}^{\text{ang.}}(q_{1},q_{2},q_{3},r_{1},r_{2},r_{3})$
$\displaystyle\equiv\sum_{m_{1}m_{2}m_{3}}\begin{pmatrix}l_{1}&l_{2}&l_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}\int\bigg{(}\prod_{i}\mathrm{d}^{2}{\hat{\bm{r}}}_{i}\mathrm{d}^{2}{\hat{\bm{q}}}_{i}e^{i{\bm{q}}_{i}\cdot{\bm{r}}_{i}}Y^{*}_{l_{i},m_{i}}({\hat{\bm{r}}}_{i})\bigg{)}$
$\displaystyle\quad\times