# The Halo Bispectrum Multipoles in Redshift Space

Federico Rizzo,11footnotetext: Corresponding author. Chiara Moretti  Kevin
Pardede  Alexander Eggemeier  Andrea Oddo  Emiliano Sefusatti  Cristiano
Porciani  and Pierluigi Monaco

###### Abstract

We present the analysis of the halo bispectrum in redshift-space in terms of
its multipoles, monopole, quadrupole and hexadecapole, measured from a large
set of simulations. We fit such measurements with a tree-level model in
perturbation theory that depends on linear and nonlinear bias parameters as
well as on the growth rate $f$ of density fluctuations. The likelihood
analysis takes advantage of a very large set of mock catalogs, enabling a
robust estimation of the covariance properties for all multipoles. We compare
the numerical estimate of the covariance matrix to its Gaussian prediction
finding discrepancies of 10% or less for all configurations with the sole
exception of the squeezed triangles in the monopole case. We find the range of
validity of the tree-level model, for the total simulation volume of about
1000$\,h^{-3}\,{\rm Gpc}^{3}$, reaches a maximum wavenumber of $0.08\,h\,{\rm
Mpc}^{-1}$ for the monopole, while it is limited to $0.06$ and $0.045\,h\,{\rm
Mpc}^{-1}$ respectively for quadrupole and hexadecapole. Despite this, the
addition of the quadrupole to the analysis allows for significant improvements
on the determination of the model parameters and specifically on $f$,
similarly to the power spectrum case. Finally, we compare our numerical
estimate for the full covariance with its theoretical prediction in the
Gaussian approximation and find the latter to work remarkably well in the
context of simulation boxes with periodic boundary condition.

##  1 Introduction

As the next generation of spectroscopic galaxy surveys [1, 2, 3], will cover
unprecedented cosmological volumes with a significant observational effort,
increasing attention is currently being payed to the full exploitation of the
information they are expected to provide. Recent measurements and analyses of
higher-order statistics such as the 3-point correlation function in
configuration space [4, 5] or its counterpart in Fourier space, the bispectrum
[6, 7, 8, 9, 10, 11, 12], go precisely in this direction, with the goal of
extending and strengthening the results from the standard analyses of the
2-point correlation function and power spectrum.

However, while the analysis of the power spectrum takes full advantage of
redshift-space distortions by means of a multipoles expansion with respect to
the angle between the wavenumber $\mathbf{k}$ and the line-of-sight (see, e.g.
[13, 14]), in the case of the bispectrum past data analyses have always been
limited to the monopole. On the other hand, the potential offered by the
galaxy bispectrum measured in future surveys to further constrain cosmological
parameters has been explored in several papers [15, 16, 17, 18, 19, 20, 21,
22, 23]. A subset of these works specifically considered the relevance of the
anisotropic bispectrum signal [15, 16, 17, 20, 22] remarking that we can
expect additional information in the higher-order multipoles of the
bispectrum, although the exact extent of such improvement on parameters
constraints, typically of the order of tens of percents, highly depends on the
assumptions on the observable, its covariance and the survey specifications.

A first theoretical modelling of the redshift-space bispectrum at tree-level
in Perturbation Theory can be found in [24] (see [25] and references therein
for earlier work on the matter and galaxy bispectrum in real space). Early
comparisons against measurements of the bispectrum monopole in numerical
simulations are presented in [26, 27, 28] with [27] including as well a first
test of the quadrupole. The analysis of the BOSS data-set of [29, 30, 6]
includes the bispectrum monopole and takes advantage of a phenomenological
model [31], based on fits to simulations, to extend the validity of the tree-
level expression to smaller scales, reaching $0.15\,h\,{\rm Mpc}^{-1}$ with a
5% accuracy on the halo bispectrum monopole at redshift $z=0.55$ (to contrast
$0.06\,h\,{\rm Mpc}^{-1}$ in the case of tree-level PT). A similar approach is
adopted as well in [32] where the monopole and quadrupole of the power
spectrum, bispectrum and integrated trispectrum are compared to simulations.

Ref. [33] goes beyond the tree-level expression presenting a one-loop PT model
for the redshift-space matter bispectrum multipoles (but defined differently
from [27]), including additional corrections along the lines of those
introduced by [34] for the power spectrum. The comparison with numerical
simulations shows an agreement up to $k\sim 0.15$ \- $0.2\,h\,{\rm Mpc}^{-1}$
depending on redshift and configuration shape while the corresponding tree-
level approximation typically fails already around $k\sim 0.07$ \-
$0.08\,h\,{\rm Mpc}^{-1}$ for both the monopole and the quadrupole.

More recently, [8] re-analysed the BOSS bispectrum monopole, adopting a tree-
level model up to $0.1\,h\,{\rm Mpc}^{-1}$, although no comparison with
simulations or details on model validation are provided. A further analysis,
extending the model to include one-loop corrections and corrections due to
primordial non-Gaussianity is presented in [10]. A comparison with large-
volume simulations can be found instead, again for the monopole only, in [23]
for measurements obtained from the very large simulation set already adopted
for the challenge paper [35], corresponding to a cumulative volume of
566$\,h^{-3}\,{\rm Gpc}^{3}$: in this case as well the reach of the tree-level
expression is found to be $k_{\rm max}\sim 0.08\,h\,{\rm Mpc}^{-1}$. The same
pipeline for the bispectrum monopole analysis is applied to the BOSS data in
[9, 11, 12].

It appears that, despite the recent attention, tests of the redshift-space
galaxy bispectrum model have been rather limited. In fact, current literature
is for the most part focused on the bispectrum monopole with only partial
assessments of higher-order multipoles predictions in PT.

The main goal of this paper is to provide a rigorous and extensive comparison
of the tree-level predictions for the halo bispectrum monopole, quadrupole and
hexadecapole (as defined in [27]) against measurements in a very large set of
numerical simulations ($\sim 1,000\,h^{-3}\,{\rm Gpc}^{3}$) while taking
advantage of a robust estimate of their covariance properties from an even
larger set of mock catalogs. It constitutes the natural continuation of a
series of papers exploring in details the challenges of a joint analysis of
the galaxy power spectrum and bispectrum, so far focused on real-space
modelling [36, 37, 38]. We will refer, in particular as Paper I and Paper II
respectively to [36] and [38] as this work shares with these references, to a
large extent, both data-sets and methodology. We test the model by means of a
likelihood analysis in terms of bias parameters along with the growth rate of
perturbations $f$, using the simulation input and real-space results as
reference values. For measurement uncertainties corresponding to the total
combined volume of our simulations, we find that the model provides a valid
description up to a maximum wavenumber of $0.08\,h\,{\rm Mpc}^{-1}$ for the
monopole, $0.06\,h\,{\rm Mpc}^{-1}$ for the quadrupole, and $0.045\,h\,{\rm
Mpc}^{-1}$ for the hexadecapole. We show that, as in the power spectrum case,
the inclusion of the bispectrum quadrupole greatly improves the posteriors
from the monopole alone.

The paper is organised as follows. In section 2 we introduce the theoretical
background for the tree-level prediction of the bispectrum multipoles in
Perturbation Theory. Section 3 describes the numerical simulations and mock
catalogs adopted and the bispectrum estimator. In Section 4 we present the
set-up for our likelihood analyses and in Section 5 the relative results. We
will present our conclusions in Section 6.

##  2 Theoretical background

###  2.1 Model

Given the halo number density contrast
$\delta_{h}(\mathbf{x})\equiv[n_{h}(\mathbf{x})-\bar{n}_{h}]/\bar{n}_{h}$
defined in terms of the number density $n_{h}(\mathbf{x})$ and its expectation
value $\bar{n}=\langle n_{h}(\mathbf{x})\rangle$, and its Fourier
Transform222We adopt the convention for the Fourier Transform
$\delta(\mathbf{k})\equiv\int
d^{3}x\,e^{-i\mathbf{k}\cdot\mathbf{x}}\,\delta(\mathbf{x})\,,$ (2.1) with the
inverse given by
$\delta(\mathbf{x})\equiv\int\frac{d^{3}k}{(2\pi)^{3}}\,e^{i\mathbf{k}\cdot\mathbf{x}}\,\delta(\mathbf{k})\,.$
(2.2) $\delta_{h}(\mathbf{k})$ we can define the halo power spectrum $P_{h}$
and bispectrum $B_{h}$ respectively as

$\displaystyle\langle\delta_{h}(\mathbf{k}_{1})\delta_{h}(\mathbf{k}_{2})\rangle$
$\displaystyle\equiv$
$\displaystyle(2\pi)^{3}\delta_{D}(\mathbf{k}_{12})\,P_{h}(k_{1})\,$ (2.3)
$\displaystyle\langle\delta_{h}(\mathbf{k}_{1})\delta_{h}(\mathbf{k}_{2})\delta_{h}(\mathbf{k}_{3})\rangle$
$\displaystyle\equiv$
$\displaystyle(2\pi)^{3}\delta_{D}(\mathbf{k}_{123})\,B_{h}(k_{1},k_{2},k_{3})\,,$
(2.4)

where the Dirac deltas $\delta_{D}$ results from the assumed statistical
homogeneity and isotropy. For the same reason $P_{h}(k_{1})$ is a function of
one variable, $k_{1}=|\mathbf{k}_{1}|$ and $B_{h}(k_{1},k_{2},k_{3})$ is a
function of the three sides of the triangle formed by $\mathbf{k}_{1}$,
$\mathbf{k}_{2}$ and $\mathbf{k}_{3}$ and independent of its orientation.

In redshift-space, peculiar velocities $\mathbf{v}$ induce distortions in the
galaxy distribution along the line-of-sight (LOS) $\hat{n}$. The observed
position $\mathbf{s}$ will then be related to real position $\mathbf{x}$ by

$\mathbf{s}=\mathbf{x}+\frac{\mathbf{v}\cdot\hat{n}}{a\,H(a)}\hat{n}\,.$ (2.5)

As a result clustering properties, and in particular galaxy correlation
functions estimated in a given region of the sky, will depend on the local
LOS. Since our focus is to test the modelling of the bispectrum based on
measurements in simulation boxes with periodic boundary conditions, we will
assume throughout this work the plane-parallel approximation for redshift-
space distortions and therefore a global, constant LOS. The halo bispectrum
will then be a function of the wavenumbers defining the triangular
configuration $\mathbf{k}_{1}$, $\mathbf{k}_{2}$ and $\mathbf{k}_{3}$ plus the
LOS $\hat{n}$, that is $B_{s}=B_{s}(\textbf{k}_{1},\textbf{k}_{2},\hat{n})$.

Our model for the redshift-space halo bispectrum is the sum of a deterministic
and stochastic contribution, as

$B_{s}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})=B_{s}^{\rm(det)}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})+B_{s}^{\rm(stoch)}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})\,,$
(2.6)

corresponding to the tree-level expression in Perturbation Theory (PT)
resulting from the halo density given, in turn, by the sum of a deterministic
and a stochastic component

$\delta_{s}=\delta_{s}^{\rm(det)}+\delta_{s}^{\rm(stoch)}\,.$ (2.7)

In Fourier space and up to the relevant order the deterministic contributions
are given by

$\delta_{s}^{\rm(det)}(\mathbf{k})=Z_{1}(\mathbf{k})\,\delta_{L}(\mathbf{k})+\int
d^{3}q_{1}d^{3}q_{2}\delta_{D}(\mathbf{k}-\mathbf{q}_{12})Z_{2}(\mathbf{q}_{1},\mathbf{q}_{2})\,\delta_{L}(\mathbf{q}_{1})\,\delta_{L}(\mathbf{q}_{1})\,,$
(2.8)

where $\delta_{L}$ is the linear matter density and the relevant $Z_{n}$
redshift-space kernels are given in terms of the local ($b_{1}$, $b_{2}$) and
tidal ($b_{\mathcal{G}_{2}}$) bias parameters and the linear growth rate $f$
by [39, 24, 26, 27, 40, 41, 42]

$\displaystyle Z_{1}(\mathbf{k})$ $\displaystyle=$ $\displaystyle
b_{1}+f\mu^{2}\,,$ (2.9) $\displaystyle Z_{2}(\mathbf{k}_{1},\mathbf{k}_{2})$
$\displaystyle=$
$\displaystyle\frac{b_{2}}{2}+b_{1}F_{2}(\mathbf{k}_{1},\mathbf{k}_{2})+b_{\mathcal{G}_{2}}S(\mathbf{k}_{1},\mathbf{k}_{2})+f\mu_{12}^{2}G_{2}(\mathbf{k}_{1},\mathbf{k}_{2})+$
(2.10)
$\displaystyle+\frac{f\mu_{12}k_{12}}{2}\left[\frac{\mu_{1}}{k_{1}}Z_{1}(\textbf{k}_{2})+\frac{\mu_{2}}{k_{2}}Z_{1}(\textbf{k}_{1})\right]$

with $F_{2}$ and $G_{2}$ representing the usual matter density and velocity
quadratic kernels and

$S(\mathbf{k}_{1},\mathbf{k}_{2})=\left(\hat{k}_{1}\cdot\hat{k}_{2}\right)^{2}-1\,$
(2.11)

while $\mu_{i}\equiv\mathbf{k}_{i}\cdot\hat{n}/k_{i}$ is the cosine of the
angle formed by the wavenumber $\mathbf{k}_{i}$ with the LOS. We also adopt
the short-hand notation for vector sums
$\mathbf{k}_{12}=\mathbf{k}_{1}+\mathbf{k}_{2}$ so that

$\mu_{12}=\frac{\mathbf{k}_{12}\cdot\hat{n}}{k_{12}}=-\frac{\textbf{k}_{3}\cdot\hat{n}}{k_{3}}=-\mu_{3}\,,$
(2.12)

for a closed triangle with $\mathbf{k}_{123}=0$. The expansions of eq. (2.8)
leads to the tree-level prediction for the bispectrum

$B_{s}^{\rm(det)}(\mathbf{k}_{1},\mathbf{k}_{2},\hat{n})=2\,Z_{1}(\mathbf{k}_{1})\,Z_{1}(\mathbf{k}_{2})Z_{2}(\mathbf{k}_{1},\mathbf{k}_{2})P_{L}(k_{1})P_{L}(k_{2})+{\rm
2~{}perm.}\,$ (2.13)

where $P_{L}(k)$ is the linear matter power spectrum.

The stochastic contribution to $\delta_{s}$ is given instead, following [42]
and their notation, by

$\delta_{s}^{\rm(stoch)}(\mathbf{x})=\epsilon(\mathbf{x})+\epsilon_{\delta}(\mathbf{x})\,\delta(\mathbf{x})+\epsilon_{\eta}(\mathbf{x})\eta(\mathbf{x})\,,$
(2.14)

where $\epsilon$, $\epsilon_{\delta}$ and $\epsilon_{\eta}$ are stochastic
fields uncorrelated to the density perturbations. The composite terms are
limited to those linear in the matter density $\delta$ and in the l.o.s.
derivative of the l.o.s velocity component
$\eta\equiv\partial_{\hat{n}}(\mathbf{v}\cdot\hat{n})$, as these are
responsible for the leading order contributions to the bispectrum. We neglect
any higher-derivative operator in the stochastic contribution and we note that
the last term should appear only due to selection effects [42]. In the large
$k$ limit, we expect to recover the Poisson predictions for the power spectrum
and bispectrum of the stochastic fields, that is [43]

$\displaystyle\langle\epsilon(\mathbf{k}_{1})\epsilon(\mathbf{k}_{2})\rangle$
$\displaystyle\rightarrow\delta_{D}(\mathbf{k}_{12})\,\frac{1}{\bar{n}}\,,$
(2.15)
$\displaystyle\langle\epsilon(\mathbf{k}_{1})\epsilon(\mathbf{k}_{2})\epsilon(\mathbf{k}_{3})\rangle$
$\displaystyle\rightarrow\delta_{D}(\mathbf{k}_{123})\,\frac{1}{\bar{n}^{2}}\,,$
(2.16)
$\displaystyle\langle\epsilon(\mathbf{k}_{1})\epsilon_{\delta}(\mathbf{k}_{2})\rangle$
$\displaystyle\rightarrow\delta_{D}(\mathbf{k}_{123})\,\frac{b_{1}}{2\bar{n}}\,,$
(2.17)
$\displaystyle\langle\epsilon(\mathbf{k}_{1})\epsilon_{\eta}(\mathbf{k}_{2})\rangle$
$\displaystyle\rightarrow\delta_{D}(\mathbf{k}_{123})\,\frac{1}{2\bar{n}}\,,$
(2.18)

where the first only appears in the halo power spectrum, while the last three
all contribute to the halo bispectrum. In principle we can expect independent
departures from the Poisson prediction for all three terms, which in the
large-scale limit can be described in terms of three constant parameters333In
[23] the authors follow [44] in the modelling of the stochastic contribution
assuming
$\delta_{s}^{(stoch)}=d_{1}\,\epsilon_{P}+d_{2}\,b_{1}\,\epsilon_{P}\,\delta+d_{1}\,\epsilon_{P}\,\eta\,,$
(2.19) where the coefficients $d_{1}$ and $d_{2}$ parameterize the corrections
to the Poisson prediction represented by field $\epsilon_{P}$ (for which the
limits (2.15) and (2.16) hold as equalities). The Poisson case is recovered
for $d_{1}=2\,d_{2}=1$. This implies that
$\langle\epsilon\epsilon\rangle=\langle\epsilon\epsilon_{\eta}\rangle$ and
their corrections to Poisson are therefore described by a single degree of
freedom. They also relate $\langle\epsilon\epsilon\rangle$ and
$\langle\epsilon\epsilon\epsilon\rangle$ but it does not seem justified. Such
relation also appears inconsistent with the expansion above and it does not
seem to be supported by the halo model description of [45].

The corresponding stochastic contribution to the bispectrum at tree-level will
then read

$\displaystyle
B_{s}^{\rm(stoch)}(\mathbf{k}_{1},\mathbf{k}_{2},\hat{n})=\frac{1}{\bar{n}}\left[(1+\alpha_{1})\,b_{1}+(1+\alpha_{3})\,f\,\mu^{2}\right]\,Z_{1}(\mathbf{k}_{1})\,P_{L}(k_{1})+2~{}{\rm
perm.}+\frac{1+\alpha_{2}}{\bar{n}^{2}}\,,$ (2.20)

where $\alpha_{i}$ are all parameters vanishing in the Poisson limit444The
notation for the $\alpha_{i}$ parameters is chosen in order to be consistent
with Paper I and Paper II, where $\alpha_{2}$ already appeared as correction
to the $1/\bar{n}^{2}$ term, while $\alpha_{3}$ was not present..

In this work we do not consider any modelling of Finger-of-God effects as we
expect them to be negligible at large scales and for a halo distribution.

###  2.2 Bispectrum multipoles

We adopt the definition of the redshift-space multipoles of the bispectrum
introduced by [27] (and assumed as well by [46] and [16]) where the vector
configurations comprising the domain of
$B_{s}(\mathbf{k}_{1}\mathbf{k}_{2},\hat{n})$ are given in terms of the
variables $k_{1}$, $k_{2}$, $k_{3}$, $\mu_{1}\equiv\cos(\theta)$ and $\xi$,
with $\theta_{1}$ being the angle between $\mathbf{k}_{1}$ and the LOS while
$\xi$ is the azimuthal angle describing a rotation of $\mathbf{k}_{2}$ around
$\mathbf{k}_{1}$.

$B_{s}$ is then expanded in spherical harmonics as

$B_{s}(k_{1},k_{2},k_{3},\theta,\xi)=\sum_{\ell}\sum_{m=-\ell}^{\ell}B_{\ell
m}(k_{1},k_{2},k_{3})\,Y_{\ell}^{m}(\theta,\xi)$ (2.21)

with the coefficients of the expansion given by

$B_{\ell
m}(k_{1},k_{2},k_{3})=\int_{-1}^{+1}\\!\\!d\\!\cos{\theta}\int_{0}^{2\pi}d\xi
B_{s}(k_{1},k_{2},k_{3},\theta,\xi)\,Y_{\ell}^{m}(\theta,\xi)\,.$ (2.22)

We only consider $m=0$ multipoles as the loss of information coming from
excluding the $m\neq 0$ terms is negligible [16]. In this case the spherical
harmonics reduce to Legendre polynomials $\mathcal{L}_{\ell}$ and only depend
on $\mu_{1}\equiv\cos\theta$,

$Y_{\ell}^{0}(\theta,\xi)=\sqrt{\frac{2\ell+1}{4\pi}}\mathcal{L}_{\ell}(\mu_{1})\,$
(2.23)

and the expansion of eq. 2.21 is replaced by

$\frac{1}{2\pi}\int d\xi
B_{s}(k_{1},k_{2},k_{3},\theta,\xi)=\sum_{\ell}\,B_{\ell}(k_{1},k_{2},k_{3})\,{\mathcal{L}}_{\ell}(\mu_{1})$
(2.24)

where

$\displaystyle B_{\ell}(k_{1},k_{2},k_{3})$
$\displaystyle=\sqrt{\frac{2\ell+1}{4\pi}}\,B_{\ell 0}(k_{1},k_{2},k_{3})$
$\displaystyle=(2\ell+1)\frac{1}{2}\int_{-1}^{+1}\\!\\!d\\!\cos{\theta}\Big{[}\frac{1}{2\pi}\int_{0}^{2\pi}d\xi
B_{s}(k_{1},k_{2},k_{3},\theta,\xi)\Big{]}\,{\mathcal{L}}_{\ell}(\cos\theta)\,.$
(2.25)

##  3 Data

###  3.1 N-body simulations

The analysis is performed on redshift-space, halo bispectrum measurements from
the set of 298 Minerva N-body simulations [47] whose real-space counterpart
was already studied in Paper I and Paper II. These follow the evolution of
$1000^{3}$ dark matter particles in a cubic box of side $L=1500\,h^{-1}\,{\rm
Mpc}$ and correspond to a total volume of about $1,000\,h^{-3}\,{\rm
Gpc}^{3}$. Each halo catalog is defined by a minimal mass of $M\simeq
1.12\times 10^{13}\;h^{-1}M_{\odot}$. We refer the reader to Paper I for a
more detailed description of the simulations and of the halo catalog
construction.

Paper II provides us with an estimate of the bias parameters characterising
the halo population obtained as posteriors from the joint analysis of the halo
power spectrum and bispectrum in real space. We will use these here as a
reference for our redshift-space analysis, in addition to the value for the
linear growth rate expected from the fiducial cosmology.

###  3.2 Bispectrum multipoles estimator

Our estimator of the bispectrum multipoles follow the definition of [27] for
the $m=0$ case and constitutes an implementation of the one described in [46]
based on Fast-Fourier Transforms. In our case, however, we assume a constant
line-of-sight $\hat{n}$, corresponding to an exact realisation of the plane-
parallel or distant observer approximation. The estimator reduces therefore to
the following expression

$\hat{B}_{\ell}=(2\ell+1)\frac{k_{f}^{3}}{N_{B}}\sum_{\mathbf{q}_{1}\in
k_{1}}\sum_{\mathbf{q}_{2}\in k_{2}}\sum_{\mathbf{q}_{3}\in
k_{3}}\,\delta_{K}(\mathbf{q}_{123})\,\delta_{s}(\mathbf{q}_{1})\,\delta_{s}(\mathbf{q}_{2})\,\delta_{s}(\mathbf{q}_{3})\,\mathcal{L}_{\ell}(\hat{q}_{1}\cdot\hat{n})\,,$
(3.1)

where the sums, accounting for the discrete nature of the Fourier Transform
$\delta_{s}(\mathbf{q})$ of the halo density in a simulation box, are over all
wavenumbers $\mathbf{q}_{i}$ falling into the bin centered at $k_{i}$ of
radial size $\Delta k$, that is, such that $k_{i}-\Delta
k/2\leq|\mathbf{q}_{i}|<k_{i}+\Delta k/2$. Also, $\delta_{K}(\mathbf{q})$ is a
Kronecker symbol equal to unity for $\mathbf{q}={\bf 0}$ and vanishing
otherwise, while the normalisation factor

$N_{B}(k_{1},k_{2},k_{3})=\sum_{\mathbf{q}_{1}\in
k_{1}}\sum_{\mathbf{q}_{2}\in k_{2}}\sum_{\mathbf{q}_{3}\in
k_{3}}\,\delta_{K}(\mathbf{q}_{123})\,,$ (3.2)

provides the number of “fundamental triangles”
$\\{\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3}\\}$ present in the “triangle
bin” $\\{k_{1},k_{2},k_{3}\\}$. The grid-interpolation of the halo density
$\delta_{s}(\mathbf{q})$ is obtained by means of a fourth-order mass
assignment scheme and adopts the interlacing technique for aliasing reduction
[48]. All bispectrum measurements, unless otherwise stated, assume a
wavenumber bin size $\Delta k=k_{f}$, that is corresponding to the fundamental
frequency characterising the simulation box, $k_{f}\equiv 2\pi/L$. This leads
to the measurement of 1475 triangular configurations up to $k_{\rm
max}=0.1\,h\,{\rm Mpc}^{-1}$ for each multipole555We include “open triangle
bins”, that is those where the bin centers cannot form a closed triangle such
as $\\{k_{1},k_{2},k_{3}\\}=\\{6,3,2\\}k_{f}$ but that nevertheless contain
closed fundamental triplets
$\\{\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3}\\}$. See section 2.2 of Paper
I for a detailed description of the binning definition..

###  3.3 Measurements

Figure 1 shows the mean of the bispectrum multipoles measured from the 298
Minerva N-body simulations for all triangular configurations. In these type of
plots, the ordering of the configurations is determined by increasing values
of $k_{1}$, $k_{2}$, $k_{3}$ which obey the requirement, $k_{1}\geq k_{2}\geq
k_{3}$ (see Paper I for a more detailed explanation). Vertical gray lines mark
the triangle where the value of $k_{1}$ changes, so that all configurations on
the left correspond to triangles made up with sides smaller or equal to such
value of $k_{1}$. All measurements include shot-noise. The bottom half of each
panel shows the relative error on the mean, along with the ratio between the
expected Poisson shot-noise contribution and the overall signal (black, dashed
lines).

Figure 1: Mean of the halo bispectrum multipoles measurements (including shot-
noise) from the 298 N-body Minerva simulations, shown for all triangular
configurations. For each multipole moment, the bottom half of each panel shows
the relative error on the mean (in color) along with the relative (Poisson)
shot-noise contribution (in black).

We notice that the relative error on the mean for the bispectrum monopole is
at the ten-percent level and just slightly smaller at smaller scales. The
shot-noise level is comparable to the statistical error at large scales and it
is larger at smaller scales, as it happens in real-space for this halo
population (see Paper I). The relative error on the quadrupole and
hexadecapole mean is instead of the order of tens of percent and hundreds of
percent, respectively. Because of the first term on the r.h.s. of (2.20), all
multipoles receive a shot-noise contribution. This is comparable to the error
on the mean for the quadrupole and lower by an order of magnitude for the
hexadecapole.

##  4 Covariance

###  4.1 Numerical estimate

As for Paper I and Paper II, the covariance properties of our observables are
estimated from a much larger set of 10,000 measurements from mock halo
catalogs obtained with the Pinocchio code [49, 50, 51]. The mocks share the
same cosmology, box size and resolution with the Minerva simulations and 298
realisations also adopt the same initial conditions. The mass threshold for
the mocks is chosen to reproduce (below percent level) the amplitude of the
large-scale halo power spectrum (including shot-noise) of the numerical
simulations (see Paper I for details). This quantity, in fact, accounts for
the Gaussian contribution to the bispectrum covariance, the leading one for
most triangular configurations (see also [52, 53]).

The covariance matrix for each bispectrum multipole and their cross-covariance
is defined as

$C_{\ell_{1}\ell_{2}}(t_{i},t_{j})\equiv\langle\hat{B}_{\ell_{1}}(t_{i})\,\hat{B}_{\ell_{2}}(t_{j})\rangle-\langle\hat{B}_{\ell_{1}}(t_{i})\rangle\langle\hat{B}_{\ell_{2}}(t_{j})\rangle\,,$
(4.1)

where $t_{i}=\\{k_{1i},k_{2i},k_{3i}\\}$ and
$t_{j}=\\{k_{1j},k_{2j},k_{3j}\\}$ represent two triangle configurations. We
will denote as $\hat{C}_{\ell_{1}\ell_{2}}(t_{i},t_{j})$ its estimate from the
10,000 mock catalogs.

Figure 2: Left column: ratio between the mean bispectrum multipoles from the
298 N-body simulations and the mean of the same quantity from the
corresponding Pinocchio mocks with matching initial conditions. Right column:
ratio between the bispectrum multipoles variance estimated from the
simulations and the one estimated from the corresponding mocks.

The left column of figure 2 shows the ratio between the mean of each
bispectrum multipole measured in the numerical simulations and the mean of the
same quantity measured in the Pinocchio mocks, limited to the 298 mocks with
matching initial conditions. We find for the monopole the same discrepancy, as
large as 7-8% depending on the triangle shape, already encountered in real
space (see Paper I). The noise in the measurements for the higher-order
multipoles, on the other hand, does not allow to clearly identify systematic
differences at the level of 10% or below. The right column of figure 2 shows
instead the ratio between the variance $\Delta B_{\ell}(t_{i})\equiv
C_{\ell\ell}(t_{i},t_{i})$, estimated again from the numerical simulations and
the one estimated from the Pinocchio mocks. Again, the Pinocchio set is
limited here to the 298 realisations with matching seeds. Despite possible
systematics on the observables, the variance is recovered by the Pinocchio
mocks with an error below 10% and no apparent systematic difference for all
multipoles. This is expected, given the close match of the power spectra and
the fact that the leading contribution to the bispectrum covariance is fully
determined by the power spectrum, see Eq. (4.2).

Figure 3: Ratio between the theoretical prediction of the multipoles
covariance $C_{\ell_{1}\ell_{2}}(t_{i},t_{i})$ in the Gaussian approximation
and its numerical estimate from the Pinocchio mocks.

###  4.2 Theoretical Gaussian covariance

In addition to the numerical estimate of the covariance we consider as well
the analytical prediction in the Gaussian approximation. The Gaussian
contribution is present for the variance of each multipole, $\Delta
B_{\ell}\equiv C_{\ell\ell}(t_{i},t_{i})$, but also for the correlation
$C_{\ell_{1}\ell_{2}}(t_{i},t_{i})$ between $\hat{B}_{\ell_{1}}(t_{i})$ and
$\hat{B}_{\ell_{2}}(t_{i})$ with $\ell_{1}\neq\ell_{2}$ but measured for the
same triangle $t_{i}$.

From the definition of the multipoles estimators, eq. (3.1), we can write

$\displaystyle C_{\ell_{1}\ell_{2}}(k_{1},k_{2},k_{3})$
$\displaystyle\equiv\langle\hat{B}_{\ell_{1}}\hat{B}_{\ell_{2}}\rangle-\langle\hat{B}_{\ell_{1}}\rangle\langle\hat{B}_{\ell_{2}}\rangle$
$\displaystyle=(2\ell_{1}+1)(2\ell_{2}+1)\frac{k_{f}^{6}}{N^{2}_{B}}\sum_{\mathbf{q}_{1}\in
k_{1}}\sum_{\mathbf{q}_{2}\in k_{2}}\sum_{\mathbf{q}_{3}\in
k_{3}}\,\delta_{K}(\mathbf{q}_{123})\,\sum_{\mathbf{p}_{1}\in
k_{1}}\sum_{\mathbf{p}_{2}\in k_{2}}\sum_{\mathbf{p}_{3}\in
k_{3}}\,\delta_{K}(\mathbf{p}_{123})$
$\displaystyle\times\,\mathcal{L}_{\ell_{1}}(\mu_{\mathbf{q}_{1}})\,\mathcal{L}_{\ell_{2}}(\mu_{\mathbf{p}_{1}})\,\Big{[}\langle\delta_{\mathbf{q}_{1}}\,\delta_{\mathbf{q}_{2}}\,\delta_{\mathbf{q}_{3}}\delta_{\mathbf{p}_{1}}\,\delta_{\mathbf{p}_{2}}\,\delta_{\mathbf{p}_{3}}\rangle-\langle\delta_{\mathbf{q}_{1}}\,\delta_{\mathbf{q}_{2}}\,\delta_{\mathbf{q}_{3}}\rangle\langle\delta_{\mathbf{p}_{1}}\,\delta_{\mathbf{p}_{2}}\,\delta_{\mathbf{p}_{3}}\rangle\Big{]}\,.$
(4.2)

In the Gaussian approximation, from the expectation values on the r.h.s. of
the equation above, we retain only the contributions depending on the field
power spectrum. Assuming (without loss of generality) that $k_{1}\geq
k_{2}\geq k_{3}$, these are given by

$\displaystyle C_{\ell_{1}\ell_{2}}^{G}(k_{1},k_{2},k_{3})$
$\displaystyle=\frac{(2\ell_{1}+1)(2\ell_{2}+1)}{N^{2}_{B}\,k_{f}^{3}}\sum_{\mathbf{q}_{1}\in
k_{1}}\sum_{\mathbf{q}_{2}\in k_{2}}\sum_{\mathbf{q}_{3}\in
k_{3}}\,\delta_{K}(\mathbf{q}_{123})\,P_{tot}(\mathbf{q}_{1})\,P_{tot}(\mathbf{q}_{2})\,P_{tot}(\mathbf{q}_{3})$
$\displaystyle\times\Big{[}(1+\delta^{K}_{k_{2},k_{3}})\,\mathcal{L}_{\ell_{1}}(\mu_{\mathbf{q}_{1}})\mathcal{L}_{\ell_{2}}(-\mu_{\mathbf{q}_{1}})+(\delta^{K}_{k_{1},k_{2}}+\delta^{K}_{k_{2},k_{3}})\,\mathcal{L}_{\ell_{1}}(\mu_{\mathbf{q}_{1}})\mathcal{L}_{\ell_{2}}(-\mu_{\mathbf{q}_{2}})+$
$\displaystyle+2\,\delta^{K}_{k_{1},k_{3}}\,\mathcal{L}_{\ell_{1}}(\mu_{\mathbf{q}_{1}})\mathcal{L}_{\ell_{2}}(-\mu_{\mathbf{q}_{3}})\Big{]}\,,$
(4.3)

where $P_{tot}(\mathbf{q}_{1})=P(\mathbf{q}_{1})+P_{SN}$ is the anisotropic
halo power spectrum including a shot-noise contribution while
$\delta^{K}_{k_{i},k_{j}}$ is the Kronecker symbol equal to one for
$k_{i}=k_{j}$, and vanishing otherwise. Notice that the terms in the square
brackets correspond, in the case $\ell_{1}=\ell_{2}=0$, to the usual factor
equal to 6, 2 and 1 respectively for equilateral, isosceles and scalene
triangles.

Similarly to the power spectrum variance case, see e.g. [47], we can expand
the anisotropic power spectra in multipoles to obtain

$\displaystyle C_{\ell_{1}\ell_{2}}^{G}$
$\displaystyle=\frac{(2\ell_{1}+1)(2\ell_{2}+1)}{N^{2}_{B}\,k_{f}^{3}}\sum_{\ell_{3},\ell_{4},\ell_{5}}\sum_{\mathbf{q}_{1}\in
k_{1}}\sum_{\mathbf{q}_{2}\in k_{2}}\sum_{\mathbf{q}_{3}\in
k_{3}}\,\delta_{K}(\mathbf{q}_{123})\,P_{tot,\,\ell_{3}}(q_{1})\,P_{tot,\,\ell_{4}}(q_{2})\,P_{tot,\,\ell_{5}}(q_{3})$
$\displaystyle\times\Big{[}(1+\delta^{K}_{k_{2},k_{3}})\,\mathcal{L}_{\ell_{1}}(\mu_{\mathbf{q}_{1}})\mathcal{L}_{\ell_{2}}(-\mu_{\mathbf{q}_{1}})+(\delta^{K}_{k_{1},k_{2}}+\delta^{K}_{k_{2},k_{3}})\,\mathcal{L}_{\ell_{1}}(\mu_{\mathbf{q}_{1}})\mathcal{L}_{\ell_{2}}(-\mu_{\mathbf{q}_{2}})+$
$\displaystyle+2\,\delta^{K}_{k_{1},k_{3}}\,\mathcal{L}_{\ell_{1}}(\mu_{\mathbf{q}_{1}})\mathcal{L}_{\ell_{2}}(-\mu_{\mathbf{q}_{3}})\Big{]}\mathcal{L}_{\ell_{3}}(\mu_{1})\mathcal{L}_{\ell_{4}}(\mu_{2})\mathcal{L}_{\ell_{5}}(\mu_{3})\,.$
(4.4)

This is the expression we adopt in our evaluation of the Gaussian variance,
with the sums over the $k$-shells performed exactly over the discrete
wavenumbers $\mathbf{q}$ defining the Fourier-space density grid666It is
possible to simplify further this expression in the thin-shell approximation
so that  $\displaystyle C_{\ell_{1}\ell_{2}}^{G}$
$\displaystyle\simeq\frac{(2\ell_{1}+1)(2\ell_{2}+1)}{N_{B}\,k_{f}^{3}}\sum_{\ell_{3},\ell_{4},\ell_{5}}P_{tot,\,\ell_{3}}(k_{1})\,P_{tot,\,\ell_{4}}(k_{2})\,P_{tot,\,\ell_{5}}(k_{3})R_{\ell_{1},\ell_{2};\ell_{3}\ell_{4},\ell_{5}}(k_{1},k_{2},k_{3})\,.$
(4.5) where we defined  $\displaystyle
R_{\ell_{1},\ell_{2};\ell_{3}\ell_{4},\ell_{5}}(k_{1},k_{2},k_{3})$
$\displaystyle\equiv\frac{1}{N_{B}}\sum_{\mathbf{q}_{1}\in
k_{1}}\sum_{\mathbf{q}_{2}\in k_{2}}\sum_{\mathbf{q}_{3}\in
k_{3}}\,\delta_{K}(\mathbf{q}_{123})\,$
$\displaystyle\times\Big{[}(1+\delta^{K}_{k_{2},k_{3}})\,\mathcal{L}_{\ell_{1}}(\mu_{\mathbf{q}_{1}})\mathcal{L}_{\ell_{2}}(-\mu_{\mathbf{q}_{1}})+(\delta^{K}_{k_{1},k_{2}}+\delta^{K}_{k_{2},k_{3}})\,\mathcal{L}_{\ell_{1}}(\mu_{\mathbf{q}_{1}})\mathcal{L}_{\ell_{2}}(-\mu_{\mathbf{q}_{2}})+$
$\displaystyle+2\,\delta^{K}_{k_{1},k_{3}}\,\mathcal{L}_{\ell_{1}}(\mu_{\mathbf{q}_{1}})\mathcal{L}_{\ell_{2}}(-\mu_{\mathbf{q}_{3}})\Big{]}\mathcal{L}_{\ell_{3}}(\mu_{1})\mathcal{L}_{\ell_{4}}(\mu_{2})\mathcal{L}_{\ell_{5}}(\mu_{3})\,.$
(4.6) In the continuum limit, we can replace the sums over the shells with
integrals and reduced them to a simple average over the orientation of the
triangle $\left\\{\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3}\right\\}$, that
is  $\frac{1}{N_{B}}\sum_{\mathbf{q}_{1}\in k_{1}}\sum_{\mathbf{q}_{2}\in
k_{2}}\sum_{\mathbf{q}_{3}\in
k_{3}}\,\delta_{K}(\mathbf{q}_{123})\simeq\frac{1}{N_{B}\,k_{f}^{6}}\prod_{i=1}^{3}\int_{k_{i}}d^{3}q_{i}\,\delta_{D}(\mathbf{q}_{123})=\frac{1}{4\pi}\int_{-1}^{1}d\mu_{1}\\!\\!\int
d\xi\,,$ (4.7) with an integrand that is now only a function of powers of
$\mu_{1}$, $\mu_{2}$ and $\mu_{3}$. Then, assumming only even values for
$\ell_{3}$, $\ell_{4}$ and $\ell_{5}$ (and clearly for $\ell_{1}$ and
$\ell_{2}$), we can use the expansion for the Legendre polynomials
$\mathcal{L}(\mu)=\frac{1}{2^{\ell}}\sum_{n=1}^{\ell/2}\frac{(-1)^{n}(2\ell-2n)!}{n!(\ell-n)!(\ell-2n)!}\mu^{\ell-2n}\equiv\sum_{n=1}^{\ell/2}C_{\ell,n}\mu^{\ell-2n}$
(4.8) to get an expression that can be automatically evaluated with a software
allowing for symbolic manipulation,  $\displaystyle
R_{\ell_{1},\ell_{2};\ell_{3}\ell_{4},\ell_{5}}(k_{1},k_{2},k_{3})$
$\displaystyle\simeq\prod_{i=1}^{5}\sum_{n_{i}=1}^{\ell_{i}/2}C_{\ell_{i},n_{i}}\Big{[}(1+\delta^{K}_{k_{2},k_{3}})\mathcal{I}_{\ell_{1}+\ell_{2}+\ell_{3}-2(n_{1}+n_{2}+n_{3}),\,\ell_{4}-2n_{4},\,\ell_{5}-2n_{5}}$
$\displaystyle+(\delta^{K}_{k_{1},k_{2}}+\delta^{K}_{k_{2},k_{3}})\mathcal{I}_{\ell_{1}+\ell_{3}-2(n_{1}+n_{3}),\,\ell_{2}+\ell_{4}-2(n_{2}+n_{4}),\,\ell_{5}-2n_{5}}$
$\displaystyle+2\,\delta^{K}_{k_{1},k_{3}}\,\mathcal{I}_{\ell_{1}+\ell_{3}-2(n_{1}+n_{3}),\,\ell_{4}-2n_{4},\,\ell_{2}+\ell_{5}-2(n_{2}+n_{5})}\Big{]}\,,$
(4.9) , where the integrals $\mathcal{I}_{\alpha,\beta,\gamma}$ are defined in
Appendix A. .

The comparison between the numerical estimate and analytic prediction for
these quantities is shown in figure 3. One can see how the Gaussian prediction
for $C_{00}$, $C_{22}$ and $C_{44}$ is able to describe the measured ones at
the level of 5%, with a slight deficit noticeable in the monopole and
quadrupole case for squeezed triangles, due to the missing non-Gaussian
contribution [52, 53]. The agreement in the case of the cross-covariance
$C_{\ell_{1}\ell_{2}}(t_{i},t_{i})$ is also rather good, with the theory
underestimating the measurements by an overall 10%. In these cases the ratio
can take large values when the denominator is close to zero, as it is the case
particularly for $C_{0,4}$.

Figure 4: Subset of the correlation matrix
$r_{\ell_{1}\ell_{2}}(t_{i},t_{j})$, defined in eq. (4.10) for $\ell=0,2$ and
4, restricted to the first 32 triangular configurations $t_{i}$. The top-right
half is estimated from the full set of 10,000 Pinocchio mocks, while the
bottom-left half is the theoretical prediction in the Gaussian approximation.
The bottom panels compare the predicted (blue) and measured (red) coefficients
$r_{\ell_{1}\ell_{2}}(t_{i},t_{i})$ with $\ell_{1}\neq\ell_{2}$ as a function
of the selected triangles $t_{i}$. The two estimates overlap almost exactly.

Finally, figure 4 shows a subset of the correlation matrix, defined as

$r_{\ell_{1}\ell_{2}}(t_{i},t_{j})=\frac{C_{\ell_{1}\ell_{2}}(t_{i},t_{j})}{\sqrt{C_{\ell_{1}\ell_{1}}(t_{i},t_{i})\,C_{\ell_{2}\ell_{2}}(t_{j},t_{j})}}\,.$
(4.10)

Each data-set $\hat{B}_{\ell}(t_{i})$ is restricted, for illustration
purposes, to its first 32 triangular configurations $t_{i}$, denoted in the
figure in terms of the three sides in units of the fundamental frequency, that
is $\\{k_{1},k_{2},k_{3}\\}/k_{f}$. It follows that while the block-diagonal
matrices show the correlation coefficients $r_{\ell\ell}(t_{i},t_{j})$ for
each multipoles, the off-diagonal matrices depict their relative cross-
covariance $r_{\ell_{1}\ell_{2}}(t_{i},t_{j})$. In addition, the top-right
half is estimated from the full set of 10,000 Pinocchio mocks, while the
bottom-left half is the theoretical prediction in the Gaussian approximation,
vanishing for all elements with $t_{i}\neq t_{j}$. The bottom panels compare
in more detail the predicted and measured coefficients
$r_{\ell_{1}\ell_{2}}(t_{i},t_{i})$ with $\ell_{1}\neq\ell_{2}$ as a function
of the selected triangles $t_{i}$. On these quantities the agreement between
theory and numerical estimates is truly remarkable and extends up to $k_{\rm
max}\sim 0.1\,h\,{\rm Mpc}^{-1}$, that is over the full range of scales that
we will consider in the analysis described in the following sections.

Regarding the structure of the correlation coefficient, it is clear that only
the elements corresponding to the expected Gaussian contributions appear to be
relevant at these large scales. These, however, are not limited to the
diagonal for the full data vector $D=\left\\{B_{0},B_{2},B_{4}\right\\}$, but
obviously include all elements corresponding to the correlation between
different multipoles sharing the same triangles.

##  5 Likelihood analysis

###  5.1 Likelihood function

Following Paper I and Paper II, we fit all measurements together assuming
their independence. This means that our total log-likelihood function
corresponds to the sum of the log-likelihood for each individual realisation,

$\ln{\mathcal{L}}_{\rm
tot}=\sum_{\alpha=1}^{N_{R}}\ln{\mathcal{L}}_{\rm\alpha}\,,$ (5.1)

where $N_{R}$ is the total number of realisations considered.

We work under the assumption of Gaussianity for the individual likelihood
${\mathcal{L}}_{\rm\alpha}$. However, we adopt the proposal of [54] to account
for possible uncertainties in the determination of the precision matrix due to
a limited number of mocks. The log-likelihood for a single realisation is
then, modulo an additive, normalization constant,

$\ln{\mathcal{L}}_{\alpha}=-\frac{N_{M}}{2}\ln\Big{[}1+\frac{\chi^{2}_{\alpha}}{N_{M}-1}\Big{]}\,,$
(5.2)

where $N_{M}$ is the number of mock catalogs used for the numerical estimation
of the covariance matrix (we refer the reader to Paper I for further details).
In this expression $\chi_{\alpha}^{2}$ represents the chi-square for the
individual realisation, given by

$\chi^{2}_{\alpha}=\sum_{i,j=1}^{N_{D}}\left[\hat{D}^{(\alpha)}_{i}-D^{(theory)}_{i}\right]C_{ij}^{-1}\left[\hat{D}^{(\alpha)}_{j}-D^{(theory)}_{j}\right]\,,$
(5.3)

where, in the most general case,
$\hat{D}^{(\alpha)}\equiv\left\\{\hat{B}_{0},\hat{B}_{2},\hat{B}_{4}\right\\}$
is the data vector, of size $N_{D}$, encompassing the three bispectrum
multipoles while $D^{(theory)}$ and $C_{ij}$ are, respectively, the
corresponding theoretical model and covariance matrix.

We should notice that given the large number of 10,000 Pinocchio realisations,
even for the largest data-set corresponding to the joint analysis of the three
bispectrum multipoles up to $k_{\rm max,B}=0.1\,h\,{\rm Mpc}^{-1}$, with a
total of $1,475\times 3\simeq 4,425$ data-points, the difference w.r.t. the
Gaussian case is in fact negligible. The alternative approach of re-scaling
the inverse covariance, as suggested in [55, 56], gives rise to error bars up
to 10% larger, although we have checked that these do not lead to any
appreciable differences in the recovered parameter posteriors that we discuss
below.

###  5.2 Model evaluation

Our main goal is assessing the validity and reach of the tree-level bispectrum
model, eqs. (2.13) and (2.20) leaving the exploration of their potential to
constrain cosmological parameters to a future work. For this reason, in our
likelihood analyses we assume galaxy bias, shot-noise and the growth rate $f$
as the only free parameters. The bispectrum multipoles defined in eq. (2.22)
can be written as a linear combination of several contributions where the
dependence on these parameters can be factorised, leading to a quick
exploration of the likelihood function since each term only needs to be
computed once for the fiducial cosmology.

This allows as well for an exact binning of the theoretical model, taking
advantage of the discrete Fourier-space grid characterising numerical
simulations in boxes with periodic boundary conditions. In this case we can
sum $B_{s}(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3},\hat{n})$ over all
discrete modes $\mathbf{q}_{i}$ forming a close triangle $\mathbf{q}_{123}=0$
and belonging to the bin $\\{k_{1},k_{2},k_{3}\\}$. This leads to

$B_{\ell}^{\rm(binned)}(k_{1},k_{2},k_{3})=\frac{2\ell+1}{N_{B}}\sum_{\mathbf{q}_{1}\in
k_{1}}\sum_{\mathbf{q}_{2}\in k_{2}}\sum_{\mathbf{q}_{3}\in
k_{3}}\delta_{K}(\mathbf{q}_{123})B_{s}(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3},\hat{n})\mathcal{L}_{\ell}(\hat{q}_{1}\cdot\hat{n}),$
(5.4)

where the sums account for the angle-average defining the bispectrum
multipoles.

This approach requires the evaluation of the bispectrum model over a very
large number of triangular configurations, making it unfeasible in a
likelihood analysis where cosmological parameters are explored. An approximate
solution would be to evaluate the model $B_{\ell}$ on a single triangle
defined by effective values of the wavenumbers. This approach, that takes
advantage of the analytical evaluation of the angle integrals in (2.2)
described in Appendix A, is presented in Appendix B along with a
quantification of the systematic errors resulting in the parameters
determination.

###  5.3 Goodness of fit and model selection

We will assess the goodness of the fits that we will perform in terms of the
posterior predictive $p$-value (ppp) and the posterior-averaged reduced chi-
square $\langle\chi^{2}_{\nu}\rangle_{\rm post}$. The ${\rm ppp}$ is defined
in the interval $0$ to $1$ with values of ${\rm ppp}\geq 0.95$ signaling a
failure of the model. The $\langle\chi^{2}_{\nu}\rangle_{\rm post}$ is
compared to the 95 percent (upper) confidence limit associated to a number of
degrees of freedom equal to the total number of data points fitted: when
$\langle\chi^{2}_{\nu}\rangle_{\rm post}$ is greater than this value, the
model fails to describe the data. For the comparison between different models
and assumptions on the bias parameters, we use the Deviance Information
Criterion (DIC) computed from the MCMC simulations. For details on the choice
of these diagnostics, we redirect the reader to Paper I.

##  6 Results

###  6.1 Maximal model

We start with a test of the full model with all seven bias and shot-noise
parameters free. We compare the fit to the bispectrum monopole to the joint
analysis of monopole and higher-order multipoles, assuming the full volume of
the 298 Minerva simulations. The first goal is to identify the set of
parameters that can effectively be determined by our data-set, and the
relative importance of the different multipoles in their constraints. We
assume uniform priors on all parameters, with bounds specified in table 1.

Parameter | Prior (uniform)
---|---
$b_{1}$ |  | $[0.9,3.5]$
$b_{2}$ |  | $[-4,4]$
$b_{\mathcal{G}_{2}}$ |  | $[-4,4]$
$\alpha_{1}$ |  | $[-1,2]$
$\alpha_{2}$ |  | $[-1,2]$
$\alpha_{3}$ |  | $[-1,2]$
$f$ |  | $[0.1,1]$
Table 1: Uniform prior intervals of the model parameters.

The main results are shown in figure 5. The left panels show the marginalised,
1-$\sigma$ posteriors for the model parameters as a function of the maximum
wavenumber included in the triangle selection. Two-dimensional marginalised
contours are shown in the bottom-right panel for the $k_{\rm
max}=0.06\,h\,{\rm Mpc}^{-1}$ case (indicated as the vertical, dotted line in
the other two sub-panels). Finally, the top-right panel shows the posterior-
averaged, reduced chi-square, $\langle\chi_{\nu}^{2}\rangle$ and the posterior
predictive $p$-value.

The tree-level model described in section 2.1 provides a good fit to the data
up to a $k_{\rm max}\simeq 0.08\,h\,{\rm Mpc}^{-1}$ for the monopole, while
the reach is restricted to 0.06 and 0.045$\,h\,{\rm Mpc}^{-1}$, respectively,
when $B_{2}$ and $B_{4}$ are also considered.

For $k_{\rm max}=0.06\,h\,{\rm Mpc}^{-1}$, the combination $B_{0}+B_{2}$
properly recovers the best-fit values of the bias parameters obtained from the
joint analysis of the real-space power spectrum and bispectrum along with the
fiducial value of the linear growth rate $f$, which are shown by the gray,
dashed lines in the bottom-right panel.

The addition of the quadrupole $B_{2}$ greatly reduces the large degeneracy
between the linear bias $b_{1}$ and the growth rate $f$, as well as with the
quadratic bias parameters, characterising the monopole-only constraints. This
qualitatively confirms the expectation that the bispectrum monopole does not
fully capture the information potentially present in the anisotropic,
redshift-space bispectrum [17, 57]. On the other hand, including as well the
hexadecapole leads to no significant improvement for any parameter and we will
drop it from all results that will follow.

Figure 5: Results for the analysis of the whole 298 Minerva simulations data-
set in terms of the full, seven-parameters model. Left panels: marginalised,
1-$\sigma$ posteriors for each parameter as a function of $k_{\rm max}$. Top-
right panels: posterior-averaged, reduced chi-square,
$\langle\chi_{\nu}^{2}\rangle$ and the posterior predictive $p$-value (ppp) as
a function of $k_{\rm max}$. The blue, red and green dashed lines in the
$\langle\chi_{\nu}^{2}\rangle$ panel represent the 95% confidence limits for
the three combinations of multipoles considered. Bottom-right panel: two-
dimensional, marginalised 1-$\sigma$ contours for $k_{\rm max}=0.06\,h\,{\rm
Mpc}^{-1}$ case (corresponding to the vertical line in the other panels). In
all panels, the $B_{0}$-only analysis (blue) is compared to the joint
$B_{0}+B_{2}$ (red) and $B_{0}+B_{2}+B_{4}$ (green). All posteriors are
compared with the results from the joint analysis of the real-space power
spectrum and bispectrum derived in Paper II, whose best-fit values are shown
by the gray, dashed lines.

###  6.2 Shot-noise

It is clear that, despite the large total simulation volume, the data set is
not able to provide meaningful constraints on all shot-noise parameters. This
was true as well, for the same halo catalogs, in real space, even including
power spectrum information (see Paper II).

In this section we compare different options to reduce the shot-noise
parameters to a single one. In addition to the maximal model, characterised by
seven parameters in total, we will consider the following models

  * •

$\alpha_{2}=0$ (6 parameters); this is justified by the posteriors obtained
for the maximal model, which show that this parameter is simply not
constrained by the data and can therefore be set to zero without affecting the
overall fit;

  * •

$\alpha_{3}=-1$ (6 parameters); this corresponds to setting
$\epsilon_{\eta}=0$, as expected under the assumption of no velocity bias and
no selection effects.

Another option is to set $\alpha_{1}=\alpha_{3}$. This is implicit in the
Poisson prediction for the shot-noise of a generic distribution in redshift
space, where both corrections vanish. This prediction is also the outcome of a
count-in-cell estimate of the shot-noise contributions to the bispectrum [58,
59] and corresponds to the standard shot-noise correction often implemented in
bispectrum estimators [46, 60] and implicitly assumed in some data analysis
[29, 6]. We will therefore consider the two, additional 5-parameters models
(both assuming $\alpha_{2}=0$):

  * •

$\alpha_{3}=\alpha_{1}$ and $\alpha_{2}=0$ (5 parameters);

  * •

$\alpha_{3}=-1$ and $\alpha_{2}=0$ (5 parameters).

Figure 6: Top left panel: difference in the DIC of the various shot-noise
models with respect to the maximal one with seven free parameters, as a
function of the largest wavenumber $k_{\rm max}$ for the analysis of the
monopole plus quadrupole. Top right panel: effective numbers of parameters as
a function of $k_{\rm max}$ for the same shot-noise models and the same
dataset, $B_{0}+B_{2}$. Bottom panel: contour plots for $B_{0}+B_{2}$ at
$k_{\rm max}=0.06\,h\,{\rm Mpc}^{-1}$, showing the 1- and 2-$\sigma$ 2D
constraints on the bias and shot-noise parameters for the maximal model
compared to the two 5-parameters models defined by setting
$\alpha_{3}=\alpha_{1}$ and $\alpha_{3}=-1$ (both with $\alpha_{2}=0$). Dashed
lines represent the fiducial value of growth rate $f$ and the best-fit values
for the bias parameters obtained in real space in Paper II.

The top left panel in figure 6 shows a general comparison between all the
models described in the bullet points above in terms of the difference in
their DIC w.r.t. the maximal model with seven parameters, as a function of the
maximum wavenumber $k_{\rm max}$, again for the monopole and quadrupole
analysis. Differences larger than 5 are usually considered relevant. The top
right panel of the same figure shows instead the effective number of
parameters we are able to constrain from the data also as a function of
$k_{\rm max}$. For a value of $k_{\rm max}<0.05\,h\,{\rm Mpc}^{-1}$, we do not
have enough information to determine even the least number of parameters and
the $\Delta$DIC simply favours the simplest models. These are still favoured
up to $k_{\rm max}\sim 0.08\,h\,{\rm Mpc}^{-1}$, where the additional degrees
of freedom of more complex models are probably accounting for missing
nonlinear corrections. The test does not clearly indicate a preference for the
5-parameters model with $\alpha_{3}=\alpha_{1}$ model over the one one with
$\alpha_{3}=-1$, except for $k_{\rm max}>0.08\,h\,{\rm Mpc}^{-1}$, where we
know that none of the models provides a good overall fit anymore. A comparison
of the 2D marginalised posteriors from the monopole and quadrupole analysis at
$k_{\rm max}=0.06\,h\,{\rm Mpc}^{-1}$ is shown for the two 5-parameters models
and the maximal one in the bottom panel of figure 6. Both models improve the
constraints on the growth rate, with minimal differences on the posteriors for
the other parameters. The $\alpha_{3}=-1$ case provides a slightly better
agreement with the fiducial value of $f$ and the real-space estimate of
$b_{1}$.

We will assume the $\alpha_{3}=-1$ (i.e. $\epsilon_{\eta}=0$) and
$\alpha_{2}=0$ case as our default model in all following tests. This implies
the expression for the shot-noise contribution

$\displaystyle
B_{s}^{\rm(stoch)}(\mathbf{k}_{1},\mathbf{k}_{2},\hat{n})=\frac{1+\alpha_{1}}{\bar{n}}\,b_{1}\,Z_{1}(\mathbf{k}_{1})\,P_{L}(k_{1})+2~{}{\rm
perm.}+\frac{1}{\bar{n}^{2}}\,,$ (6.1)

only depending on the parameter $\alpha_{1}$. We assume that this model
provides an accurate description of the stochastic contribution to the
bispectrum, consistent with the large-scale expectation, in a relatively
restricted range ($k\lesssim 0.08\,h\,{\rm Mpc}^{-1}$) where the Poisson limit
($\alpha_{1}=\alpha_{3}=0$) does not apply.

###  6.3 Bias relations

The parameter space can be further reduced by introducing relations among the
bias parameters. In Paper I and Paper II we considered a few of such
relations, either theoretically motivated or from results of numerical
simulations [61, 62, 63, 64, 65]. Of those, we select the two that provide the
best improvement to the fit of the power spectrum and bispectrum in real space
and test them again here in redshift space. The first is the fitting function
for $b_{2}(b_{1},b_{\mathcal{G}_{2}})$ obtained in [62] from separate universe
simulations. The second is the fit to the excursion set prediction for the
tidal bias parameter $b_{\mathcal{G}_{2}}(b_{1})$ proposed in [64, 65]. For
convenience we reproduce these two relations here:

$\displaystyle b_{2}(b_{1},b_{\mathcal{G}_{2}})$
$\displaystyle=0.412-2.142\,b_{1}+0.929\,b_{1}^{2}+0.008\,b_{1}^{3}+\frac{4}{3}b_{\mathcal{G}_{2}}\,,$
(6.2) $\displaystyle b_{\mathcal{G}_{2}}(b_{1})$
$\displaystyle=0.524-0.547\,b_{1}+0.046\,b_{1}^{2}\,.$ (6.3)

Figure 7: Top left panel: difference in the DIC of the two model adopting the
bias relations $b_{2}(b_{1},b_{\mathcal{G}_{2}})$ and
$b_{\mathcal{G}_{2}}(b_{1})$ with respect to the reference, 5-parameter model,
as a function of the largest wavenumber $k_{\rm max}$ for the analysis of the
monopole plus quadrupole. Top right panel: effective numbers of parameters as
a function of $k_{\rm max}$ for the same models and datasets, $B_{0}+B_{2}$.
Bottom panel: contour plots for $B_{0}+B_{2}$ at $k_{\rm max}=0.06\,h\,{\rm
Mpc}^{-1}$, showing the 1- and 2-$\sigma$ 2D constraints on the bias and shot-
noise parameters for the three models. Dashed lines represent the fiducial
value of growth rate $f$ and the best-fit values for the bias parameters
obtained in real space in Paper II.

In figure 7 we compare three cases with our reference 5-parameter model:
applying each of the bias relations $b_{2}(b_{1},b_{\mathcal{G}_{2}})$ and
$b_{\mathcal{G}_{2}}(b_{1})$ individually, as well as the two of them
combined. In the top-left panel we show their difference in the DIC with
respect to the reference model, as a function of the largest wavenumber
$k_{\rm max}$ for the analysis of the monopole plus quadrupole. The top-right
panel shows instead the effective numbers of parameters, again as a function
of $k_{\rm max}$ and for $B_{0}+B_{2}$. The DIC shows a marginal preference
for the $b_{\mathcal{G}_{2}}(b_{1})$ relation, for values of $k_{\rm max}$
close to $0.06\,h\,{\rm Mpc}^{-1}$, whereas the combination of the two
relations quickly becomes disfavoured beyond $k_{\rm max}=0.05\,h\,{\rm
Mpc}^{-1}$.

In the bottom panel of figure 7 we show the contour plots from the analysis at
$k_{\rm max}=0.06\,h\,{\rm Mpc}^{-1}$, showing the 1- and 2-$\sigma$ 2D
constraints on the bias and shot-noise parameters for the three models. All
three cases lead to a reduction in the marginalised posterior constraints,
however, we notice how the application of the tidal bias relation leads to
constraints on $b_{1}$ and $f$ that are systematically offset from the
fiducial values, while the other two cases involving the $b_{2}(b_{1})$
relation significantly reduces any potential tension. We caution that this
outcome might in fact be fortuitous since the $b_{2}(b_{1})$ relation crosses
the $b_{1}$ \- $b_{2}$ contour of the reference model close to the fiducial
value of $b_{1}$ recovered from the real-space analysis of Paper II. We should
stress, in any event, that these systematic differences are only evident due
to the very large cumulative simulation volume: we leave for future work an
assessment for a more realistic volume [66] .

###  6.4 Scale cuts

In section 6.1 we have seen that the range of validity of the tree-level model
is more limited for higher-order multipoles with respect to the monopole. This
suggests that it could be convenient to consider different values of $k_{\rm
max}$ for the different multipoles.

Figure 8: Contour plots for the $B_{0}+B_{2}$ analysis assuming $k_{\rm
max}^{B_{0}}=k_{\rm max}^{B_{2}}=0.06\,h\,{\rm Mpc}^{-1}$ (blue) compared to
the cases with $k_{\rm max}^{B_{0}}=0.08\,h\,{\rm Mpc}^{-1}$ and $k_{\rm
max}^{B_{2}}=0.06\,h\,{\rm Mpc}^{-1}$ (yellow) and $k_{\rm max}^{B_{0}}=k_{\rm
max}^{B_{2}}=0.08\,h\,{\rm Mpc}^{-1}$ (magenta). The gray, dashed mark the
best-fit values from the real-space, joint analysis of power spectrum and
bispectrum of Paper II.

We illustrate this point in figure 8 where we compare the contour plots for
the $B_{0}+B_{2}$ analysis under the following assumptions:

  * •

$k_{\rm max}^{B_{0}}=k_{\rm max}^{B_{2}}=0.06\,h\,{\rm Mpc}^{-1}$ (blue);

  * •

$k_{\rm max}^{B_{0}}=0.08\,h\,{\rm Mpc}^{-1}$ and $k_{\rm
max}^{B_{2}}=0.06\,h\,{\rm Mpc}^{-1}$ (yellow);

  * •

$k_{\rm max}^{B_{0}}=k_{\rm max}^{B_{2}}=0.08\,h\,{\rm Mpc}^{-1}$ (magenta).

In all cases we adopt the reference, 5-parameter model.

We find that extending $k_{\rm max}^{B_{0}}$ to $0.08\,h\,{\rm Mpc}^{-1}$ can
significantly reduce the error, in particular on the bias parameters. However,
for the $k_{\rm max}^{B_{0}}=k_{\rm max}^{B_{2}}=0.08\,h\,{\rm Mpc}^{-1}$
case, where the $\chi^{2}$ for the fit is already above the 95% C.L., we
notice that the $f$-$b_{1}$ contour already shows a discrepancy with the
expected values at more than 2-$\sigma$. We will adopt the scale cuts defined
by $k_{\rm max}^{B_{0}}=0.08\,h\,{\rm Mpc}^{-1}$ and $k_{\rm
max}^{B_{2}}=0.06\,h\,{\rm Mpc}^{-1}$ as our reference choice for most of the
tests in the following sections.

###  6.5 Covariance approximations

In section 4 we directly compared the Gaussian prediction for the bispectrum
covariance with the numerical estimate from the Pinocchio mocks, finding a
remarkable agreement both in the variance $C_{\ell_{1}\ell_{2}}(t_{i},t_{i})$
as in the correlation coefficients $r_{\ell_{1}\ell_{2}}(t_{i},t_{i})$.

Figure 9: Contour plots for the $B_{0}+B_{2}$ reference analysis under
different assumptions for the covariance. The case of the full covariance,
$C_{\ell_{1}\ell_{2}}(t_{i},t_{j})$, estimated from the mocks (blue contours)
is compared to the variance $C_{\ell_{1}\ell_{2}}(t_{i},t_{i})$ estimated from
the mocks (orange) and from the Gaussian theory prediction (magenta). The last
two cases overlap almost exactly. The analysis assumes $k_{\rm
max}^{B_{0}}=0.08\,h\,{\rm Mpc}^{-1}$ and $k_{\rm max}^{B_{2}}=0.06\,h\,{\rm
Mpc}^{-1}$.

Here we check if any residual difference could affect the parameters
determination. The comparison is shown in figure 9 in terms of the contour
plots for the reference analysis of $B_{0}$ and $B_{2}$ assuming,
respectively, $k_{\rm max}^{B_{0}}=0.08\,h\,{\rm Mpc}^{-1}$ and $k_{\rm
max}^{B_{2}}=0.06\,h\,{\rm Mpc}^{-1}$. The full numerical covariance,
$\hat{C}_{\ell_{1},\ell_{2}}(t_{i},t_{j})$ (blue contours) is compared to the
numerical variance $\hat{C}_{\ell_{1},\ell_{2}}(t_{i},t_{i})$ (orange) and to
the Gaussian prediction (magenta). In the last two cases, all elements with
$t_{i}\neq t_{j}$ are set to zero. We notice that the inclusion of such
elements is responsible for constraints only slightly larger, while the
Gaussian prediction reproduces the results from the numerical variance case
almost exactly, with no appreciable differences in the 1D or 2D marginalised
posteriors. This is perhaps not too surprising given that the analysis is
restricted to relatively large-scales.

##  7 Conclusions

This work presents a test of the tree-level prediction in Perturbation Theory
for the halo bispectrum in redshift space with particular attention to its
anisotropic signal as described by higher-order multipoles such as the
quadrupole and hexadecapole. It extends previous results in real space (Paper
I and Paper II), taking advantage of a very large set of 298 N-body
simulations corresponding to a cumulative volume of about 1,000$\,h^{-3}\,{\rm
Gpc}^{3}$ and an even larger set of Pinocchio mocks. The latter provides a
robust estimate of the covariance properties for the full data vector given by
the three bispectrum multipoles. We explore different assumptions on the
observables and related covariance models and assess them in terms of
constraints on bias parameters and the fiducial linear growth rate.

We summarise below our main findings.

  * •

The Pinocchio mocks provide a very good description of the variance estimated
from full numerical simulations with a residual scatter below the 10% level
and no apparent shape-dependence (Fig. 2), for all bispectrum multipoles,
extending previous assessments in real-space (Paper I, [67]).

  * •

The comparison of the likelihood analysis based on the monopole alone with the
joint analyses of $B_{0}+B_{2}$ and $B_{0}+B_{2}+B_{4}$ (Fig. 5), using the
full covariance from the mocks, indicates that the addition of the quadrupole
alone greatly improves the determination of bias parameters and, perhaps not
surprisingly, allows to properly constrain the growth rate $f$; the further
addition of the hexadecapole, instead, leads to no appreaciable improvement.

  * •

For our full simulation volume, the tree-level model provides a good fit for
the bispectrum monopole up to $k_{\rm max}=0.08\,h\,{\rm Mpc}^{-1}$ while the
inclusion of the quadrupole and of the hexadecapole reduce significantly this
range to $0.06$ and $0.045\,h\,{\rm Mpc}^{-1}$, respectively. Indeed, an
optimal configuration for the joint $B_{0}+B_{2}$ analysis should assume
distinct values for $k_{\rm max}$. We tested that better constraints on the
model parameters are obtained assuming $k_{\rm max}^{B_{0}}=0.08\,h\,{\rm
Mpc}^{-1}$ and $k_{\rm max}^{B_{2}}=0.06\,h\,{\rm Mpc}^{-1}$ with respect to
the case where a single, lower value of $k_{\rm max}^{B_{0}}=k_{\rm
max}^{B_{2}}=0.06\,h\,{\rm Mpc}^{-1}$ is adopted to keep within the model
validity range (Fig. 8).

  * •

In general our data, despite the large volume, cannot fully determine all
constant shot-noise parameters describing different departures from the
Poisson expectation. It appears, however, that the stochastic velocity
contribution $\epsilon_{\eta}\eta$ to the halo density, that one can expect
when selection effects are present [68], is indeed disfavoured in our ideal
set-up, at least at the large scales we are exploring (Fig. 6).

  * •

Both the fitting function for the quadratic local bias parameter
$b_{2}(b_{1},b_{\mathcal{G}_{2}})$ of [62] as the relation for the tidal bias
parameter $b_{\mathcal{G}_{2}}(b_{1})$ proposed in [64, 65] (and their
combination) appear to significantly improve the posteriors on $b_{1}$ and $f$
(Fig. 7); the first, however, appears to introduce some bias in the
determination of $b_{1}$, although relevant only because of the large
cumulative volume of our simulations set.

  * •

For our ideal measurements in a simulation box with periodic boundary
conditions, the Gaussian model for the bispectrum multipoles covariance
provides a very good approximation to the numerical estimate. A small
underestimate is noticeable (and expected [52, 53]) for the squeezed
configurations of the bispectrum monopole (Fig. 3). On the other hand a quite
remarkable agreement is obtained in the comparison with the cross-correlation
coefficients (Fig. 4). In terms of the posteriors on bias parameters and $f$
we find no noticeable difference between the Gaussian theory variance and the
numerical variance, while a very small difference is present when compared to
the full numerical covariance (Fig. 9).

As mentioned in the introduction, not many works explored the modelling and
the potential information of the anisotropic bispectrum. Our results
qualitatively confirm earlier Fisher-matrix forecasts [15, 16, 17, 22] in
remarking on the importance of going beyond the analysis of the bispectrum
monopole. The same can be said for [20] and [32], both based instead on a full
likelihood analysis and therefore closer, in principle, to our work. For these
last two references, however, many differences in methodology do not allow a
rigorous, quantitative comparison with our results, in addition to the fact,
of course, that we do not include power spectrum measurements in our data-
vector. We will present a joint analysis of the Minerva-set power spectrum and
bispectrum in redshift space elsewhere [66]. For the time being we limit
ourselves to observe that the inclusion of full anisotropic bispectrum
information will likely be an important step toward a complete exploitation of
cosmological information in spectroscopic galaxy surveys.

## Acknowledgments

We are always grateful to Claudio Dalla Vecchia and Ariel Sanchez for running
and making available the Minerva simulations, performed on the Hydra and
Euclid clusters at the Max Planck Computing and Data Facility (MPCDF) in
Garching. We are grateful to Guido d’Amico and Vincent Desjacques for useful
discussions. The Pinocchio mocks were run on the GALILEO cluster at CINECA,
thanks to an agreement with the University of Trieste. K.P., E.S. and P.M. are
partially supported by the INFN INDARK PD51 grant and acknowledge support from
PRIN MIUR 2015 Cosmology and Fundamental Physics: illuminating the Dark
Universe with Euclid. C.M. acknowledges support from a UK Research and
Innovation Future Leaders Fellowship (MR/S016066/1). For the purpose of open
access, the author has applied a Creative Commons Attribution (CC-BY) licence
to any Author Accepted Manuscript version arising from this submission.

##  Appendix A Analytical evaluation of the bispectrum multipoles

The orientation of the
$\left\\{\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3}\right\\}$ triangle
w.r.t. the LOS in the model for the redshift-space bispectrum
$B_{s}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3},\hat{n})$, eq. (2.6), is
expressed entirely in terms of products of powers of the cosines
$\mu_{i}\equiv\hat{k}_{i}\cdot\hat{n}$ that can be factorised in each
contribution.

The integrals defining the bispectrum multipoles in eq. (2.2) are therefore
limited to angle-averages of such combinations, with the additional factors
from the Legendre polynomials corresponding to additional powers of $\mu_{1}$.
We denote these integrals as

$\mathcal{I}_{\alpha\beta\gamma}=\frac{1}{4\pi}\int_{-1}^{+1}d\mu_{1}\int_{0}^{2\pi}d\xi\,\mu_{1}^{\alpha}\,\mu_{2}^{\beta}(\mu_{1},\xi)\,\mu_{3}^{\gamma}(\mu_{1},\xi)\,,$
(A.1)

where

$\mu_{2}(\mu_{1},\xi)=\mu_{1}\mu_{12}-\sqrt{1-\mu_{1}^{2}}\sqrt{1-\mu^{2}_{12}}\cos\xi$
(A.2)

and

$\mu_{3}(\mu_{1},\xi)=-\frac{k_{1}\mu_{1}+k_{2}\mu_{2}(\mu_{1},\xi)}{k_{3}}\,,$
(A.3)

having introduced $\mu_{12}\equiv\hat{k}_{1}\cdot\hat{k_{2}}$. Since the angle
integration of eq. (A.1) is to be intended as a generic integration over all
orientations, it is easy to see that it should satisfy the following property

$\mathcal{I}_{\alpha,\beta,\gamma}(k_{1},k_{2},k_{3})=\mathcal{I}_{\sigma(\alpha,\beta,\gamma)}[\sigma(k_{1},k_{2},k_{3})]\,,$
(A.4)

where $\sigma(...)$ represents a generic permutation applied, at the same
time, to its indices and arguments777Notice that in general is not true that
$\mathcal{I}_{\alpha,\beta,\gamma}(k_{1},k_{2},k_{3})=\mathcal{I}_{\sigma(\alpha,\beta,\gamma)}(k_{1},k_{2},k_{3})$,
so if we write, e.g., $\mathcal{I}_{200}+\mathcal{I}_{020}$ we implicitly mean
$\mathcal{I}_{200}(k_{1},k_{2},k_{3})+\mathcal{I}_{020}(k_{1},k_{2},k_{3})$
corresponding to the sum of two different quantities. We write, for
illustration purposes the first few as

$\displaystyle\mathcal{I}_{\alpha 00}$
$\displaystyle=\frac{1}{1+\alpha}~{}{\rm for}~{}\alpha~{}{\rm
even~{}(vanishing~{}otherwise)}\,,$ (A.5) $\displaystyle\mathcal{I}_{\alpha
01}$
$\displaystyle=\frac{k_{1}^{2}+k_{3}^{2}-k_{2}^{2}}{2(2+\alpha)\,k_{1}\,k_{3}}~{}{\rm
for}~{}\alpha~{}{\rm odd~{}(vanishing~{}otherwise)}\,,$ (A.6)
$\displaystyle\mathcal{I}_{\alpha 11}$
$\displaystyle=\frac{(2+\alpha)\,k_{1}^{4}-\alpha\,(k_{2}^{2}-k_{3}^{2})^{2}-2\,k_{1}^{2}\,(k_{2}^{2}+k_{3}^{2})}{4(1+\alpha)(3+\alpha)k_{1}^{2}\,k_{2}\,k_{3}}~{}{\rm
for}~{}\alpha~{}{\rm even~{}(vanishing~{}otherwise)}\,,$ (A.7)
$\displaystyle\mathcal{I}_{\alpha 02}$
$\displaystyle=\frac{4\,k_{1}^{2}\,k_{3}^{2}+\alpha\,(k_{1}^{2}+k_{3}^{2}-k_{2}^{2})^{2}}{4(1+\alpha)(3+\alpha)k_{1}^{2}\,k_{3}^{2}}~{}{\rm
for}~{}\alpha~{}{\rm even~{}(vanishing~{}otherwise)}\,.$ (A.8)

We can group all contributions to the bispectrum multipoles according to the
source of quadratic nonlinearity, that is

$B_{\ell}(k_{1},k_{2},k_{3})=B_{\ell}^{(F_{2})}+B_{\ell}^{(b_{2})}+B_{\ell}^{(S_{2})}+B_{\ell}^{(G_{2})}+B_{\ell}^{({\rm
mixed})}\,.$ (A.9)

These contributions can be expressed in terms of the
$\mathcal{I}_{\alpha\beta\gamma}$

$\displaystyle B_{0}^{(F_{2})}$
$\displaystyle=2\,b_{1}^{3}\,\mathcal{I}_{000}\,F_{2}(k_{1},k_{2};k_{3})\,P_{L}(k_{1})\,P_{L}(k_{2})+2~{}{\rm
perm.}+$
$\displaystyle+2\,b_{1}^{2}\,f\,(\mathcal{I}_{200}+\mathcal{I}_{020})\,F_{2}(k_{1},k_{2};k_{3})\,P_{L}(k_{1})\,P_{L}(k_{2})+2~{}{\rm
perm.}+$
$\displaystyle+2\,b_{1}\,f^{2}\,\mathcal{I}_{220}\,F_{2}(k_{1},k_{2};k_{3})\,P_{L}(k_{1})\,P_{L}(k_{2})+2~{}{\rm
perm.}\,,$ (A.10) $\displaystyle B_{0}^{(b2)}$
$\displaystyle=b_{1}^{2}\,b_{2}\,\mathcal{I}_{000}\,P_{L}(k_{1})\,P_{L}(k_{2})+2~{}{\rm
perm.}+$
$\displaystyle+b_{1}\,b_{2}\,f\,(\mathcal{I}_{200}+\mathcal{I}_{020})\,P_{L}(k_{1})\,P_{L}(k_{2})+2~{}{\rm
perm.}+$
$\displaystyle+b_{2}\,f^{2}\,\mathcal{I}_{220}P_{L}(k_{1})P_{L}(k_{2})+2~{}{\rm
perm.}\,,$ (A.11) $\displaystyle B_{0}^{(S_{2})}$
$\displaystyle=2\,b_{1}^{2}\,b_{\mathcal{G}_{2}}\,\mathcal{I}_{000}\,S(k_{1},k_{2};k_{3})\,P_{L}(k_{1})\,P_{L}(k_{2})+2~{}{\rm
perm.}+$
$\displaystyle+2\,b_{1}\,b_{\mathcal{G}_{2}}\,f\,(\mathcal{I}_{200}+\mathcal{I}_{020})\,S(k_{1},k_{2};k_{3})\,P_{L}(k_{1})\,P_{L}(k_{2})+2~{}{\rm
perm.}+$
$\displaystyle+2\,b_{\mathcal{G}_{2}}\,f^{2}\,\mathcal{I}_{220}\,S(k_{1},k_{2};k_{3})\,P_{L}(k_{1})\,P_{L}(k_{2})+2~{}{\rm
perm.}\,,$ (A.12) $\displaystyle B_{0}^{(G_{2})}$
$\displaystyle=2\,b_{1}^{2}\,f\,\mathcal{I}_{002}\,G_{2}(k_{1},k_{2};k_{3})\,P_{L}(k_{1})\,P_{L}(k_{2})+2~{}{\rm
perm.}+$
$\displaystyle+2\,b_{1}\,f^{2}\,(\mathcal{I}_{202}+\mathcal{I}_{022})\,G_{2}(k_{1},k_{2};k_{3})\,P_{L}(k_{1})\,P_{L}(k_{2})+2~{}{\rm
perm.}+$
$\displaystyle+2\,f^{3}\,\mathcal{I}_{222}\,G_{2}(k_{1},k_{2};k_{3})\,P_{L}(k_{1})\,P_{L}(k_{2})+2~{}{\rm
perm.}\,,$ (A.13) $\displaystyle B_{0}^{(\text{mixed})}$
$\displaystyle=-b_{1}^{3}\,f\,\Big{(}\frac{k_{3}}{k_{1}}\mathcal{I}_{101}+\frac{k_{3}}{k_{2}}\mathcal{I}_{011}\Big{)}\,P_{L}(k_{1})\,P_{L}(k_{2})+2~{}{\rm
perm.}-$ $\displaystyle-
b_{1}^{2}\,f^{2}\,\Big{[}\frac{k_{3}}{k_{1}}(\mathcal{I}_{301}+2\mathcal{I}_{121})+\frac{k_{3}}{k_{2}}(\mathcal{I}_{031}+2\mathcal{I}_{211})\Big{]}\,P_{L}(k_{1})\,P_{L}(k_{2})+2~{}{\rm
perm.}-$ $\displaystyle-
b_{1}\,f^{3}\,\Big{[}\frac{k_{3}}{k_{1}}(\mathcal{I}_{141}+2\mathcal{I}_{321})+\frac{k_{3}}{k_{2}}(\mathcal{I}_{411}+2\mathcal{I}_{231})\Big{]}\,P_{L}(k_{1})\,P_{L}(k_{2})+2~{}{\rm
perm.}-$
$\displaystyle-f^{4}\,\Big{(}\frac{k_{3}}{k_{1}}\mathcal{I}_{341}+\frac{k_{3}}{k_{2}}\mathcal{I}_{431}\Big{)}\,P_{L}(k_{1})\,P_{L}(k_{2})+2~{}{\rm
perm.}\,.$ (A.14)

Here the permutations are intended to apply as well on the
$\mathcal{I}_{\alpha\beta\gamma}$ integrals. The shot-noise contribution to
the monopole is given by

$\displaystyle B_{0}^{(\text{shot-noise})}$
$\displaystyle=\frac{1}{\tilde{n}}\Big{\\{}b_{1}^{2}\,\mathcal{I}_{000}\,\Big{[}P_{L}(k_{1})+P_{L}(k_{2})+P_{L}(k_{3})\Big{]}+$
$\displaystyle+b_{1}f\Big{[}\mathcal{I}_{200}\,P(k_{1})+\mathcal{I}_{020}\,P(k_{2})+\mathcal{I}_{002}\,P(k_{3})\Big{]}+$
$\displaystyle+f^{2}\Big{[}\mathcal{I}_{400}P(k_{1})+\mathcal{I}_{040}P(k_{2})P+\mathcal{I}_{004}P(k_{3})\Big{]}\Big{\\}}+\frac{\mathcal{I}_{000}}{\bar{n}^{2}}$
(A.15)

Then the quadrupole and hexadecapole, defined as

$\displaystyle B_{2}(k_{1},k_{2},k_{3})$
$\displaystyle=\int_{-1}^{+1}d\mu_{1}\int_{0}^{2\pi}d\xi
B(k_{1},k_{2},k_{2},\mu_{1},\xi)\Big{[}\frac{1}{4}\sqrt{\frac{5}{\pi}}(3\mu_{1}^{2}-1)\Big{]}$
(A.16) $\displaystyle B_{4}(k_{1},k_{2},k_{3})$
$\displaystyle=\int_{-1}^{+1}d\mu_{1}\int_{0}^{2\pi}d\xi
B(k_{1},k_{2},k_{2},\mu_{1},\xi)\Big{[}\frac{3}{16}\sqrt{\frac{1}{\pi}}(35\mu_{1}^{4}-30\mu_{1}^{2}+3)\Big{]}$
(A.17)

can be readily written starting from eq.s (A)-(A), replacing the
$\mathcal{I}_{\alpha\beta\gamma}$ integrals with the quantities

$\mathcal{J}^{(2)}_{\alpha\beta\gamma}\equiv\frac{1}{4}\sqrt{\frac{5}{\pi}}\left(3\,\mathcal{I}_{\alpha+2\,\beta\,\gamma}-\mathcal{I}_{\alpha\,\beta\,\gamma}\right)\,,$
(A.18)

in the quadrupole case and with

$\mathcal{J}^{(4)}_{\alpha\beta\gamma}\equiv\frac{3}{16\,\sqrt{\pi}}\left(35\,\mathcal{I}_{\alpha+4\,\beta\,\gamma}-30\mathcal{I}_{\alpha+2\,\beta\,\gamma}+3\mathcal{I}_{\alpha\,\beta\,\gamma}\right)\,,$
(A.19)

for the hexadecapole.

##  Appendix B Effective binning of the theoretical model

All of our results assumed an evaluation of the theory predictions
implementing the exact scheme of eq. (5.4). Since this approach can be
numerically quite demanding, particularly in likelihood evaluations extended
to several cosmological parameters, it is worth exploring the systematic
errors induced on the parameter posteriors by the more efficient choice of a
single bispectrum evaluation on the effective wavenumbers, after the
analytical integration over the angles described in Appendix A.

In this case, the theoretical prediction is given by

$B_{\ell}^{\rm eff}(k_{1},k_{2},k_{3})\equiv B_{\ell}(k_{\rm eff,l},k_{\rm
eff,m},k_{\rm eff,s}),$ (B.1)

where the definition of the effective triplet, in general not unique, is based
on “sorted” $\left\\{q_{1},q_{2},q_{3}\right\\}$ triplets as (see Paper
I)888See also [65] and [23] for alternative proposals..

$\displaystyle k_{\rm eff,l}(k_{1},k_{2},k_{3})$
$\displaystyle=\frac{1}{N_{B}}\sum_{\textbf{q}_{1}\in
k_{1}}\sum_{\textbf{q}_{2}\in k_{2}}\sum_{\textbf{q}_{3}\in
k_{3}}\delta_{K}(\textbf{q}_{123})\max(q_{1},q_{2},q_{3})\,,$ $\displaystyle
k_{\rm eff,m}(k_{1},k_{2},k_{3})$
$\displaystyle=\frac{1}{N_{B}}\sum_{\textbf{q}_{1}\in
k_{1}}\sum_{\textbf{q}_{2}\in k_{2}}\sum_{\textbf{q}_{3}\in
k_{3}}\delta_{K}(\textbf{q}_{123}){\rm med}(q_{1},q_{2},q_{3})\,,$
$\displaystyle k_{\rm eff,s}(k_{1},k_{2},k_{3})$
$\displaystyle=\frac{1}{N_{B}}\sum_{\textbf{q}_{1}\in
k_{1}}\sum_{\textbf{q}_{2}\in k_{2}}\sum_{\textbf{q}_{3}\in
k_{3}}\delta_{K}(\textbf{q}_{123})\min(q_{1},q_{2},q_{3})\,.$ (B.2)

For the choice of the bin size $\Delta k=k_{f}$ adopted in our result, the
difference between the two approaches, estimated in terms of the posteriors on
the bias and shot-noise parameters, is completely negligible. For a larger
size of the bin $\Delta k$, useful to reduce the overall size of the data
vector, however, we can find some systematic effect on parameters
determination.

## References

  * [1] R. Laureijs, J. Amiaux, S. Arduini, J. . Auguères, J. Brinchmann, R. Cole, M. Cropper, C. Dabin, L. Duvet, A. Ealet, and et al., Euclid Definition Study Report, ArXiv: 1110.3193 (Oct., 2011) [arXiv:1110.3193].
  * [2] M. Levi, C. Bebek, T. Beers, R. Blum, R. Cahn, D. Eisenstein, B. Flaugher, K. Honscheid, R. Kron, O. Lahav, P. McDonald, N. Roe, D. Schlegel, and representing the DESI collaboration, The DESI Experiment, a whitepaper for Snowmass 2013, ArXiv: 1308.0847 (Aug., 2013) [arXiv:1308.0847].
  * [3] 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, ArXiv e-prints (Dec., 2014) [arXiv:1412.4872].
  * [4] Z. Slepian, D. J. Eisenstein, J. R. Brownstein, C.-H. Chuang, H. Gil-Marín, S. Ho, F.-S. Kitaura, W. J. Percival, A. J. Ross, G. Rossi, H.-J. Seo, A. Slosar, and M. Vargas-Magaña, Detection of baryon acoustic oscillation features in the large-scale three-point correlation function of SDSS BOSS DR12 CMASS galaxies, Mon. Not. R. Astron. Soc.  469 (Aug., 2017) 1738–1751, [arXiv:1607.0609].
  * [5] A. Veropalumbo, I. Sáez Casares, E. Branchini, B. R. Granett, L. Guzzo, F. Marulli, M. Moresco, L. Moscardini, A. Pezzotta, and S. de la Torre, A joint 2- and 3-point clustering analysis of the VIPERS PDR2 catalogue at z 1: breaking the degeneracy of cosmological parameters, Mon. Not. R. Astron. Soc.  507 (Oct., 2021) 1184–1201, [arXiv:2106.1258].
  * [6] H. Gil-Marín, W. J. Percival, L. Verde, J. R. Brownstein, C.-H. Chuang, F.-S. Kitaura, S. A. Rodríguez-Torres, and M. D. Olmstead, The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: RSD measurement from the power spectrum and bispectrum of the DR12 BOSS galaxies, Mon. Not. R. Astron. Soc.  465 (Feb., 2017) 1757–1788, [arXiv:1606.0043].
  * [7] D. W. Pearson and L. Samushia, A Detection of the Baryon Acoustic Oscillation features in the SDSS BOSS DR12 Galaxy Bispectrum, Mon. Not. R. Astron. Soc.  478 (Aug., 2018) 4500–4512, [arXiv:1712.0497].
  * [8] G. d’Amico, J. Gleyzes, N. Kokron, K. Markovic, L. Senatore, P. Zhang, F. Beutler, and H. Gil-Marín, The cosmological analysis of the SDSS/BOSS data from the Effective Field Theory of Large-Scale Structure, Journal of Cosmology and Astro-Particle Physics 2020 (May, 2020) 005, [arXiv:1909.0527].
  * [9] O. H. E. Philcox and M. M. Ivanov, BOSS DR12 full-shape cosmology: $\Lambda$ CDM constraints from the large-scale galaxy power spectrum and bispectrum monopole, Phys. Rev. D  105 (Feb., 2022) 043517, [arXiv:2112.0451].
  * [10] G. D’Amico, M. Lewandowski, L. Senatore, and P. Zhang, Limits on primordial non-Gaussianities from BOSS galaxy-clustering data, arXiv e-prints (Jan., 2022) arXiv:2201.11518, [arXiv:2201.1151].
  * [11] G. Cabass, M. M. Ivanov, O. H. E. Philcox, M. Simonović, and M. Zaldarriaga, Constraints on Single-Field Inflation from the BOSS Galaxy Survey, arXiv e-prints (Jan., 2022) arXiv:2201.07238, [arXiv:2201.0723].
  * [12] G. Cabass, M. M. Ivanov, O. H. E. Philcox, M. Simonović, and M. Zaldarriaga, Constraints on Multi-Field Inflation from the BOSS Galaxy Survey, arXiv e-prints (Apr., 2022) arXiv:2204.01781, [arXiv:2204.0178].
  * [13] S. Alam, M. Ata, S. Bailey, F. Beutler, D. Bizyaev, J. A. Blazek, A. S. Bolton, J. R. Brownstein, A. Burden, C.-H. Chuang, J. Comparat, A. J. Cuesta, K. S. Dawson, D. J. Eisenstein, S. Escoffier, H. Gil-Marín, J. N. Grieb, N. Hand, S. Ho, K. Kinemuchi, D. Kirkby, F. Kitaura, E. Malanushenko, V. Malanushenko, C. Maraston, C. K. McBride, R. C. Nichol, M. D. Olmstead, D. Oravetz, N. Padmanabhan, N. Palanque-Delabrouille, K. Pan, M. Pellejero-Ibanez, W. J. Percival, P. Petitjean, F. Prada, A. M. Price-Whelan, B. A. Reid, S. A. Rodríguez-Torres, N. A. Roe, A. J. Ross, N. P. Ross, G. Rossi, J. A. Rubiño-Martín, S. Saito, S. Salazar-Albornoz, L. Samushia, A. G. Sánchez, S. Satpathy, D. J. Schlegel, D. P. Schneider, C. G. Scóccola, H.-J. Seo, E. S. Sheldon, A. Simmons, A. Slosar, M. A. Strauss, M. E. C. Swanson, D. Thomas, J. L. Tinker, R. Tojeiro, M. V. Magaña, J. A. Vazquez, L. Verde, D. A. Wake, Y. Wang, D. H. Weinberg, M. White, W. M. Wood-Vasey, C. Yèche, I. Zehavi, Z. Zhai, and G.-B. Zhao, The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample, Mon. Not. R. Astron. Soc.  470 (Sept., 2017) 2617–2652, [arXiv:1607.0315].
  * [14] R. Ruggeri, W. J. Percival, H. Gil-Marín, F. Beutler, E.-M. Mueller, F. Zhu, N. Padmanabhan, G.-B. Zhao, P. Zarrouk, A. G. Sánchez, J. Bautista, J. Brinkmann, J. R. Brownstein, F. Baumgarten, C.-H. Chuang, K. Dawson, H.-J. Seo, R. Tojeiro, and C. Zhao, The clustering of the SDSS-IV extended Baryon Oscillation Spectroscopic Survey DR14 quasar sample: measuring the evolution of the growth rate using redshift-space distortions between redshift 0.8 and 2.2, Mon. Not. R. Astron. Soc.  483 (Mar., 2019) 3878–3887, [arXiv:1801.0289].
  * [15] Y.-S. Song, A. Taruya, and A. Oka, Cosmology with anisotropic galaxy clustering from the combination of power spectrum and bispectrum, Journal of Cosmology and Astro-Particle Physics 8 (Aug., 2015) 007, [arXiv:1502.0309].
  * [16] P. Gagrani and L. Samushia, Information Content of the Angular Multipoles of Redshift-Space Galaxy Bispectrum, Mon. Not. R. Astron. Soc.  467 (May, 2017) 928–935, [arXiv:1610.0348].
  * [17] V. Yankelevich and C. Porciani, Cosmological information in the redshift-space bispectrum, Mon. Not. R. Astron. Soc.  (Nov., 2018) [arXiv:1807.0707].
  * [18] A. Chudaykin and M. M. Ivanov, Measuring neutrino masses with large-scale structure: Euclid forecast with controlled theoretical error, Journal of Cosmology and Astro-Particle Physics 2019 (Nov., 2019) 034, [arXiv:1907.0666].
  * [19] C. Hahn, F. Villaescusa-Navarro, E. Castorina, and R. Scoccimarro, Constraining Mν with the bispectrum. Part I. Breaking parameter degeneracies, Journal of Cosmology and Astro-Particle Physics 2020 (Mar., 2020) 040, [arXiv:1909.1110].
  * [20] D. Gualdi and L. Verde, Galaxy redshift-space bispectrum: the importance of being anisotropic, Journal of Cosmology and Astro-Particle Physics 2020 (June, 2020) 041, [arXiv:2003.1207].
  * [21] C. Hahn and F. Villaescusa-Navarro, Constraining Mν with the bispectrum. Part II. The information content of the galaxy bispectrum monopole, Journal of Cosmology and Astro-Particle Physics 2021 (Apr., 2021) 029, [arXiv:2012.0220].
  * [22] N. Agarwal, V. Desjacques, D. Jeong, and F. Schmidt, Information content in the redshift-space galaxy power spectrum and bispectrum, Journal of Cosmology and Astro-Particle Physics 2021 (Mar., 2021) 021, [arXiv:2007.0434].
  * [23] M. M. Ivanov, O. H. E. Philcox, T. Nishimichi, M. Simonović, M. Takada, and M. Zaldarriaga, Precision analysis of the redshift-space galaxy bispectrum, arXiv e-prints (Oct., 2021) arXiv:2110.10161, [arXiv:2110.1016].
  * [24] E. Hivon, F. R. Bouchet, S. Colombi, and R. Juszkiewicz, Redshift distortions of clustering: a Lagrangian approach., Astron. Astrophys.  298 (June, 1995) 643, [astro-ph/9407049].
  * [25] F. Bernardeau, S. Colombi, E. Gaztañaga, and R. Scoccimarro, Large-scale structure of the universe and cosmological perturbation theory, Phys. Rep. 367 (Sept., 2002) 1–3, [astro-ph/0112551].
  * [26] L. Verde, A. F. Heavens, S. Matarrese, and L. Moscardini, Large-scale bias in the universe - ii. redshift-space bispectrum, Mon. Not. R. Astron. Soc.  300 (Nov., 1998) 747–756, [astro-ph/].
  * [27] R. Scoccimarro, H. M. P. Couchman, and J. A. Frieman, The bispectrum as a signature of gravitational instability in redshift space, Astrophys. J.  517 (June, 1999) 531–540, [astro-ph/9808305].
  * [28] R. Scoccimarro, The bispectrum: From theory to observations, Astrophys. J.  544 (Dec., 2000) 597–615, [astro-ph/0004086].
  * [29] H. Gil-Marín, J. Noreña, L. Verde, W. J. Percival, C. Wagner, M. Manera, and D. P. Schneider, The power spectrum and bispectrum of SDSS DR11 BOSS galaxies - I. Bias and gravity, Mon. Not. R. Astron. Soc.  451 (July, 2015) 539–580, [arXiv:1407.5668].
  * [30] H. Gil-Marín, L. Verde, J. Noreña, A. J. Cuesta, L. Samushia, W. J. Percival, C. Wagner, M. Manera, and D. P. Schneider, The power spectrum and bispectrum of SDSS DR11 BOSS galaxies - II. Cosmological interpretation, Mon. Not. R. Astron. Soc.  452 (Sept., 2015) 1914–1921, [arXiv:1408.0027].
  * [31] H. Gil-Marín, C. Wagner, J. Noreña, L. Verde, and W. Percival, Dark matter and halo bispectrum in redshift space: theory and applications, Journal of Cosmology and Astro-Particle Physics 12 (Dec., 2014) 029, [arXiv:1407.1836].
  * [32] D. Gualdi, H. Gil-Marín, and L. Verde, Joint analysis of anisotropic power spectrum, bispectrum and trispectrum: application to N-body simulations, Journal of Cosmology and Astro-Particle Physics 2021 (July, 2021) 008, [arXiv:2104.0397].
  * [33] I. Hashimoto, Y. Rasera, and A. Taruya, Precision cosmology with redshift-space bispectrum: A perturbation theory based model at one-loop order, Phys. Rev. D  96 (Aug., 2017) 043526, [arXiv:1705.0257].
  * [34] A. Taruya, T. Nishimichi, and S. Saito, Baryon acoustic oscillations in 2D: Modeling redshift-space power spectrum from perturbation theory, Phys. Rev. D  82 (Sept., 2010) 063522, [arXiv:1006.0699].
  * [35] T. Nishimichi, G. D’Amico, M. M. Ivanov, L. Senatore, M. Simonović, M. Takada, M. Zaldarriaga, and P. Zhang, Blinded challenge for precision cosmology with large-scale structure: Results from effective field theory for the redshift-space galaxy power spectrum, Phys. Rev. D  102 (Dec., 2020) 123541, [arXiv:2003.0827].
  * [36] A. Oddo, E. Sefusatti, C. Porciani, P. Monaco, and A. G. Sánchez, Toward a robust inference method for the galaxy bispectrum: likelihood function and model selection, Journal of Cosmology and Astro-Particle Physics 2020 (Mar., 2020) 056, [arXiv:1908.0177].
  * [37] D. Alkhanishvili, C. Porciani, E. Sefusatti, M. Biagetti, A. Lazanu, A. Oddo, and V. Yankelevich, The reach of next-to-leading-order perturbation theory for the matter bispectrum, arXiv e-prints (July, 2021) arXiv:2107.08054, [arXiv:2107.0805].
  * [38] A. Oddo, F. Rizzo, E. Sefusatti, C. Porciani, and P. Monaco, Cosmological parameters from the likelihood analysis of the galaxy power spectrum and bispectrum in real space, Journal of Cosmology and Astro-Particle Physics 2021 (Nov., 2021) 038, [arXiv:2108.0320].
  * [39] J. N. Fry and E. Gaztañaga, Biasing and hierarchical statistics in large-scale structure, Astrophys. J.  413 (Aug., 1993) 447–452, [astro-ph/9302009].
  * [40] K. C. Chan, R. Scoccimarro, and R. K. Sheth, Gravity and large-scale nonlocal bias, Phys. Rev. D  85 (Apr., 2012) 083509, [arXiv:1201.3614].
  * [41] T. Baldauf, U. Seljak, V. Desjacques, and P. McDonald, Evidence for quadratic tidal tensor bias from the halo bispectrum, Phys. Rev. D  86 (Oct., 2012) 083540, [arXiv:1201.4827].
  * [42] V. Desjacques, D. Jeong, and F. Schmidt, The galaxy power spectrum and bispectrum in redshift space, Journal of Cosmology and Astro-Particle Physics 2018 (Dec., 2018) 035, [arXiv:1806.0401].
  * [43] F. Schmidt, Towards a self-consistent halo model for the nonlinear large-scale structure, Phys. Rev. D  93 (Mar., 2016) 063512, [arXiv:1511.0223].
  * [44] A. Perko, L. Senatore, E. Jennings, and R. H. Wechsler, Biased Tracers in Redshift Space in the EFT of Large-Scale Structure, ArXiv e-prints (Oct., 2016) [arXiv:1610.0932].
  * [45] D. Ginzburg, V. Desjacques, and K. C. Chan, Shot noise and biased tracers: A new look at the halo model, Phys. Rev. D  96 (Oct., 2017) 083528, [arXiv:1706.0873].
  * [46] R. Scoccimarro, Fast estimators for redshift-space clustering, Phys. Rev. D  92 (Oct., 2015) 083532, [arXiv:1506.0272].
  * [47] J. N. Grieb, A. G. Sánchez, S. Salazar-Albornoz, and C. Dalla Vecchia, Gaussian covariance matrices for anisotropic galaxy clustering measurements, Mon. Not. R. Astron. Soc.  457 (Apr., 2016) 1577–1592, [arXiv:1509.0429].
  * [48] E. Sefusatti, M. Crocce, R. Scoccimarro, and H. M. P. Couchman, Accurate estimators of correlation functions in Fourier space, Mon. Not. R. Astron. Soc.  460 (Aug., 2016) 3624–3636, [arXiv:1512.0729].
  * [49] P. Monaco, T. Theuns, and G. Taffoni, The pinocchio algorithm: pinpointing orbit-crossing collapsed hierarchical objects in a linear density field, Mon. Not. R. Astron. Soc.  331 (Apr., 2002) 587–608, [astro-ph/].
  * [50] P. Monaco, E. Sefusatti, S. Borgani, M. Crocce, P. Fosalba, R. K. Sheth, and T. Theuns, An accurate tool for the fast generation of dark matter halo catalogues, Mon. Not. R. Astron. Soc.  433 (Aug., 2013) 2389–2402, [arXiv:1305.1505].
  * [51] E. Munari, P. Monaco, E. Sefusatti, E. Castorina, F. G. Mohammad, S. Anselmi, and S. Borgani, Improving fast generation of halo catalogues with higher order Lagrangian perturbation theory, Mon. Not. R. Astron. Soc.  465 (Mar., 2017) 4658–4677, [arXiv:1605.0478].
  * [52] A. Barreira, The squeezed matter bispectrum covariance with responses, Journal of Cosmology and Astro-Particle Physics 2019 (Mar., 2019) 008, [arXiv:1901.0124].
  * [53] M. Biagetti, L. Castiblanco, J. Noreña, and E. Sefusatti, The Covariance of Squeezed Bispectrum Configurations, arXiv e-prints (Nov., 2021) arXiv:2111.05887, [arXiv:2111.0588].
  * [54] E. Sellentin and A. F. Heavens, Parameter inference with estimated covariance matrices, Mon. Not. R. Astron. Soc.  456 (Feb., 2016) L132–L136, [arXiv:1511.0596].
  * [55] T. W. Anderson, An introduction to multivariate statistical analysis.  Wiley New York, 1958.
  * [56] J. Hartlap, T. Schrabback, P. Simon, and P. Schneider, The non-Gaussianity of the cosmic shear likelihood or how odd is the Chandra Deep Field South?, Astron. Astrophys.  504 (Sept., 2009) 689–703, [arXiv:0901.3269].
  * [57] D. Gualdi, H. Gil-Marín, M. Manera, B. Joachimi, and O. Lahav, GEOMAX: beyond linear compression for three-point galaxy clustering statistics, Mon. Not. R. Astron. Soc.  497 (July, 2020) 776–792, [arXiv:1912.0101].
  * [58] P. J. E. Peebles, The large-scale structure of the universe.  Princeton, N.J., Princeton University Press, 1980. 435 p., 1980.
  * [59] S. Matarrese, L. Verde, and A. F. Heavens, Large-scale bias in the universe: bispectrum method, Mon. Not. R. Astron. Soc.  290 (Oct., 1997) 651–662, [astro-ph/9706059].
  * [60] N. S. Sugiyama, S. Saito, F. Beutler, and H.-J. Seo, A complete FFT-based decomposition formalism for the redshift-space bispectrum, Mon. Not. R. Astron. Soc.  484 (Mar., 2019) 364–384, [arXiv:1803.0213].
  * [61] R. K. Sheth, K. C. Chan, and R. Scoccimarro, Nonlocal Lagrangian bias, Phys. Rev. D  87 (Apr., 2013) 083002, [arXiv:1207.7117].
  * [62] T. Lazeyras, C. Wagner, T. Baldauf, and F. Schmidt, Precision measurement of the local bias of dark matter halos, Journal of Cosmology and Astro-Particle Physics 2 (Feb., 2016) 018, [arXiv:1511.0109].
  * [63] T. Lazeyras and F. Schmidt, Beyond LIMD bias: a measurement of the complete set of third-order halo bias parameters, Journal of Cosmology and Astro-Particle Physics 9 (Sept., 2018) 008, [arXiv:1712.0753].
  * [64] A. Eggemeier, R. Scoccimarro, M. Crocce, A. Pezzotta, and A. G. Sánchez, Testing one-loop galaxy bias: Power spectrum, Phys. Rev. D  102 (Nov., 2020) 103530, [arXiv:2006.0972].
  * [65] A. Eggemeier, R. Scoccimarro, R. E. Smith, M. Crocce, A. Pezzotta, and A. G. Sánchez, Testing one-loop galaxy bias: Joint analysis of power spectrum and bispectrum, Phys. Rev. D  103 (June, 2021) 123550, [arXiv:2102.0690].
  * [66] C. Moretti, F. Rizzo, K. Pardede, A. Oddo, E. Sefusatti, C. Porciani, and P. Monaco in preparation.
  * [67] M. Colavincenzo, E. Sefusatti, P. Monaco, L. Blot, M. Crocce, M. Lippich, A. G. Sánchez, M. A. Alvarez, A. Agrawal, S. Avila, A. Balaguera-Antolínez, R. Bond, S. Codis, C. Dalla Vecchia, A. Dorta, P. Fosalba, A. Izard, F.-S. Kitaura, M. Pellejero-Ibanez, G. Stein, M. Vakili, and G. Yepes, Comparing approximate methods for mock catalogues and covariance matrices - III: bispectrum, Mon. Not. R. Astron. Soc.  482 (Feb., 2019) 4883–4905, [arXiv:1806.0949].
  * [68] V. Desjacques, D. Jeong, and F. Schmidt, Large-scale galaxy bias, Phys. Rep. 733 (Feb., 2018) 1–193, [arXiv:1611.0978].