# Multi-tracer power spectra and bispectra: Formalism

Dionysios Karagiannis,1 Roy Maartens,1,2,3 José Fonseca,4,5,1 Stefano
Camera,6,7,8,1 Chris Clarkson9,1

###### Abstract

The power spectrum and bispectrum of dark matter tracers are key and
complementary probes of the Universe. Next-generation surveys will deliver
good measurements of the bispectrum, opening the door to improved cosmological
constraints and the breaking of parameter degeneracies, from the combination
of the power spectrum and bispectrum. Multi-tracer power spectra have been
used to suppress cosmic variance and mitigate the effects of nuisance
parameters and systematics. We present a bispectrum multi-tracer formalism
that can be applied to next-generation survey data. Then we perform a simple
Fisher analysis to illustrate qualitatively the improved precision on
primordial non-Gaussianity that is expected to come from the bispectrum multi-
tracer. In addition, we investigate the parametric dependence of conditional
errors from multi-tracer power spectra and multi-tracer bispectra, on the
differences between the biases and the number densities of two tracers. Our
results suggest that optimal constraints arise from maximising the ratio of
number densities, the difference between the linear biases, the difference
between the quadratic biases, and the difference between the products
$b_{1}\,b_{\Phi}$ for each tracer, where $b_{\Phi}$ is the bias for the
primordial potential.

##  1 Introduction

Spectacular advances in observations and computations have transformed our
understanding of the Universe. At the same time, this has brought new puzzles
to the fore, such as the still unknown nature of dark matter and dark energy
that are core to the current standard model, and the emergence of tensions
between observations that appear to challenge the foundations of the standard
model. The next generation of large-scale structure surveys will deliver major
advances in precision cosmology, cast new light on current tensions and
anomalies, and open new windows on the primordial and late-time Universe.

Up to now, cosmological information has mainly been extracted via the galaxy
power spectrum, but the galaxy bispectrum is set to play an increasingly
important role, in combination with the power spectrum (see e.g. the recent
work [1, 2, 3, 4, 5, 6, 7] and references therein). Recent measurements of the
BOSS bispectrum [2, 3, 5], including the first measurement of the bispectrum
quadrupole [8], prepare the way for next-generation surveys. These upcoming
surveys are also predicted to detect the bispectrum dipole [9, 10, 11] and
other relativistic effects in the bispectrum [12, 13, 14, 15]. Wide-angle
effects in the galaxy bispectrum [16, 17, 15, 18] will also need to be taken
into account.

The bispectrum is a direct probe of non-Gaussianity, and one of its key
strengths lies in constraints on primordial non-Gaussianity (PNG). Unlike the
power spectrum, it can constrain non-local types of PNG, while in combination
with the power spectrum it can improve constraints on local PNG [19, 20, 21,
22, 23, 24].

In order to further improve the precision of power spectrum constraints, one
can use a multi-tracer analysis, which suppresses cosmic variance and also
reduces the effects of nuisance parameters and systematics (see e.g. [25, 26,
27, 28]). A multi-tracer analysis is particularly effective for constraints on
local PNG and other quantities that have a signal on very large scales, where
cosmic variance is strongest (see e.g. [29, 30, 31, 32, 33, 34, 35, 36, 37,
38, 39, 40, 41]).

This suggests the use of a multi-tracer bispectrum analysis, in combination
with the corresponding multi-tracer power spectra. Here we develop for the
first time a Fisher forecast multi-tracer formalism for the tree-level power
spectra + bispectra combination. Multi-tracer power spectra have been combined
with the auto-bispectrum of one of the tracers by [38]. We extend such
approaches to include all the bispectra.

A multi-tracer forecast has been applied to bispectra alone (in real space) by
[42], on the basis of simplifying assumptions and without giving the
covariances or other details. We fill the gap by providing a systematic
analysis in redshift-space, including all clustering bias and noise parameters
at tree level and giving explicit forms for the covariances.

In this first paper, we do not perform detailed forecasts, which will be
presented in future work. However, we do show qualitatively the improvements
that are potentially delivered by the full power spectrum and bispectrum
multi-tracer of 2 tracers, by investigating the parametric behaviour of the
conditional errors on PNG. To begin with, we do not include wide-angle or
relativistic effects.

##  2 Single-tracer power spectrum and bispectrum

The power spectrum of the Bardeen gauge-invariant primordial gravitational
potential is defined in Fourier space by

${\langle\Phi(\bm{k})\,\Phi(\bm{k}^{\prime})\rangle}=(2\,\pi)^{3}\,\delta_{\rm
D}(\bm{k}+\bm{k}^{\prime})\,P_{\Phi}(k)\;,$ (2.1)

where $P_{\Phi}(k)$ is directly related to the power spectrum of the
primordial curvature perturbations $\zeta$ (during the matter-dominated era,
$\Phi=3\,\zeta/5$), which are generated during inflation. In the case of the
standard single-field slow-roll inflationary scenario, they have a nearly
perfect Gaussian distribution, which means that they can be adequately
characterised by their power spectrum. The primordial perturbations $\Phi$ are
in turn related to the linear dark matter over-density field through the
Poisson equation, ${\delta}(\bm{k},z)=M(k,z)\,\Phi(\bm{k})$, where

$\displaystyle{P}(k,z)$ $\displaystyle=M(k,z)^{2}\,P_{\Phi}(k)\;,$
$\displaystyle M(k,z)$ $\displaystyle=\frac{2\,D(z)}{3\,\Omega_{{\rm
m},0}\,H_{0}^{2}\,{g_{\rm dec}}}\,T(k)\,k^{2}\;.$ (2.2)

Here $D(z)$ is the linear growth factor (since the linear fluid equations
generate a linearly evolved matter density field), normalised to unity at
$z=0$, and $T(k)$ is the matter transfer function normalized to unity at large
scales $k\rightarrow 0$. The factor $g_{\rm dec}$ is the Bardeen potential
growth factor at decoupling, which ensures that $f_{\rm NL}$ is in the CMB
convention [43, 44]. The linear power spectrum is computed with the numerical
Boltzmann code CAMB [45].

Deviations from the standard inflationary model will produce a violation of
Gaussian initial conditions [46]. PNG generates nonzero higher-order
correlators in the primordial curvature perturbations, starting with the
bispectrum, i.e. the Fourier transform of the three-point function,

$\langle\Phi(\bm{k}_{1})\,\Phi(\bm{k}_{2})\,\Phi(\bm{k}_{3})\rangle=2\,(2\,\pi)^{3}\,\delta_{\rm
D}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3})\,B_{\Phi}(k_{1},k_{2},k_{3})\;.$ (2.3)

Higher-order correlators in the primordial density field are a product of
nonlinear interactions during the inflationary or reheating stage. The
primordial bispectrum arising from such interactions is characterized by a
dimensionless amplitude parameter $f_{\rm NL}$, giving the strength of the PNG
signal, and a shape function $F$, i.e. $B_{\Phi}(k_{1},k_{2},k_{3})=f_{\rm
NL}\,F(k_{1},k_{2},k_{3})$. $F$ describes the dependence of the bispectrum on
various triangles, where different inflationary models generate signal that
peaks at distinct triangle configurations. By Eq. 2.2, the leading order PNG
contribution to the matter density bispectrum is

$B_{\rm
PNG}(k_{1},k_{2},k_{3},z)=M(k_{1},z)\,M(k_{2},z)\,M(k_{3},z)\,B_{\Phi}(k_{1},k_{2},k_{3})\;.$
(2.4)

The matter bispectrum has additional contributions from nonlinearity induced
by gravity, even at zeroth order. In other words, the dark matter density
field is intrinsically non-Gaussian, even in the absence of PNG. Therefore, in
order to disentangle the PNG information from late-time nonlinearity, an
adequate description for the latter is needed. Here we model nonlinearity in
the framework of Standard Perturbation Theory (SPT) (see e.g. [47] for a
review). In addition to the nonlinear nature of the dark matter density field,
the hierarchical formation of halos and the way they are populated by luminous
matter induces non-Gaussianity in the distribution of observable structures.
Hence, in order to extract the PNG signal from large-scale structure, we need
the bias relationship between the observed tracers and the matter field.

We follow the perturbative approach of [48, 49, 50], in which the halo over-
density field $\delta_{h}$ is described as a function of all possible local
gravitational observables, which are introduced in the expansion in the form
of renormalised operators. A complete Eulerian bias expansion can be built
from the tensor $\partial_{i}\partial_{j}\Phi$, which contains the trace
$\nabla^{2}\Phi\propto\delta$ and the trace-free tidal field
$s_{ij}=\big{(}\partial_{i}\partial_{j}-\delta_{ij}\,\nabla^{2}/3\big{)}\,\nabla^{-2}\delta$.
Note that the term ‘halo’ here can be replaced by a general matter tracer.

We consider up to second-order terms in the expansion, which are sufficient
for the spatial scales considered here, i.e. much larger scales than those
involved in halo formation. For Gaussian initial conditions, the Eulerian halo
density contrast is [48, 49, 50, 44]:

${\delta^{h}}(\bm{x},\tau)=b_{1}(\tau)\,\delta(\bm{x},\tau)+\varepsilon(\bm{x},\tau)+\frac{b_{2}(\tau)}{2}\,\delta(\bm{x},\tau)^{2}+\frac{b_{s}(\tau)}{2}\,s(\bm{x},\tau)^{2}+\varepsilon_{\delta}(\bm{x},\tau)\,\delta(\bm{x},\tau)\;.$
(2.5)

Here $\tau$ is conformal time, $\bm{x}$ are spatial comoving Eulerian
coordinates, $s^{2}=s_{ij}\,s^{ij}$ is the simplest scalar that can be formed
from the tidal field, $\varepsilon$ is the leading stochastic field [51, 52,
53] and $\varepsilon_{\delta}$ is the stochastic field associated with the
linear bias. These fields take into account the stochastic relation between
the tracer density and any large-scale field. The second-order tidal field
bias coefficient, following the convention in [54], is given by
$b_{s}=-4\,(b_{1}-1)/7$.

In the presence of PNG there is a scale-dependent correction to the linear
bias $b_{1}$, especially in the case of local PNG [55, 56, 57, 58, 59, 60],
since the local PNG bispectrum peaks in squeezed triangles. (A similar scale-
dependent bias correction can be derived for any general nonlocal quadratic
non-Gaussianity template [61, 62, 63, 64].) Then the non-Gaussian bias terms,
linear in $f_{\rm NL}$, are [65]:

${\delta^{h}_{\rm
PNG}}(\bm{x},\tau)=b_{\Psi}(\tau)\,\Psi(\bm{q})+b_{\Psi\delta}(\tau)\,\Psi(\bm{q})\,\delta(\bm{x},\tau)+\varepsilon_{\Psi}(\bm{x},\tau)\,\Psi(\bm{q})\;,$
(2.6)

where $\bm{q}$ are the spatial coordinates in the Lagrangian frame and $\Psi$
is a nonlocal transformation of $\Phi$:

$\displaystyle\Psi(\bm{q})=\int\frac{{\rm
d}^{3}\bm{k}}{(2\,\pi)^{3}}\,k^{\alpha}\,\Phi(\bm{k})\,{\rm e}^{{\rm
i}\,\bm{q}\cdot\bm{k}}\;.$ (2.7)

The parameter $\alpha=2,1,0$ for equilateral, orthogonal and local shapes.
Note that in the local case, $\Psi\rightarrow\Phi$. The stochastic counterpart
of $\Psi$ is $\varepsilon_{\Psi}$. The bias coefficients in Eqs. 2.5, LABEL:
and 2.6 can be derived by using the peak-background split argument (see e.g.
[21]).

The single-tracer power spectrum and bispectrum in redshift space are

$\displaystyle{P^{h}}(\bm{k},z)$ $\displaystyle={\cal
D}_{P}(\bm{k},z)\,Z_{1}(\bm{k},z)^{2}\,P(k,z)+P_{\varepsilon}(z),$ (2.8)
$\displaystyle{B^{h}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},z)$
$\displaystyle={\cal
D}_{B}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},z)\,\Big{[}2\,Z_{1}(\bm{k}_{1},z)\,Z_{1}(\bm{k}_{2},z)\,Z_{2}(\bm{k}_{1},\bm{k}_{2},z)\,P(k_{1},z)\,P(k_{2},z)$
$\displaystyle~{}~{}+2\,\text{perm}\Big{]}+2\,P_{\varepsilon\varepsilon_{\delta}}(z)\,\Big{[}Z_{1}(\bm{k}_{1},z)\,P(k_{1},z)+2\,\text{perm}\Big{]}+B_{\varepsilon}(z)\;.$
(2.9)

At tree level, the first- and second-order kernels $Z_{1,2}$ include redshift
space distortions (RSD) [66, 67, 68] and PNG. In order to deal with redshift
errors and nonlinear RSD, we include the phenomenological factors ${\cal
D}_{P,B}$. Explicit expressions are in Appendix A. The stochastic terms
$P_{\varepsilon},\;P_{\varepsilon\varepsilon_{\delta}},\;B_{\varepsilon}$ are
generated by stochastic bias and their fiducial values are taken to be those
predicted by Poisson statistics [69]:

$P_{\varepsilon}=\frac{1}{{\overline{n}_{h}}}\;,~{}~{}P_{\varepsilon\varepsilon_{\delta}}=\frac{b_{1}}{2\,\overline{n}_{h}}\;,~{}~{}B_{\varepsilon}=\frac{1}{\overline{n}^{2}_{h}}\;.$
(2.10)

##  3 Multi-tracer power spectra and bispectra

In the multi-tracer framework (e.g. [25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
35, 36, 37, 38, 39, 40, 41]) we consider the combined information of all
correlations of different tracer samples in a given redshift bin and sky
patch. The multi-tracer power spectra and bispectra are

$\displaystyle\langle\delta^{I}(\bm{k}_{1})\,\delta^{J}(\bm{k}_{2})\rangle$
$\displaystyle=(2\,\pi)^{3}\,\delta_{\rm
D}(\bm{k}_{1}+\bm{k}_{2})\,P^{IJ}(\bm{k}_{1}),$ (3.1)
$\displaystyle\langle\delta^{I}(\bm{k}_{1})\,\delta^{J}(\bm{k}_{2})\,\delta^{K}(\bm{k}_{3})\rangle$
$\displaystyle=2\,(2\,\pi)^{3}\,\delta_{\rm
D}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3})\,B^{IJK}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})\;,$
(3.2)

where here and elsewhere we omit the $z$-dependence for brevity. The tracer
indices $I,J,K$ range over the tracer labels:
$I,J,K=t,t^{\prime},t^{\prime\prime},\cdots$. Note that the cross-power
spectra are symmetric under permutations of the tracers, while the cross-
bispectra require symmetrisation in order to avoid double counting of
information:

$\displaystyle P^{(IJ)}$ $\displaystyle\equiv P^{IJ},$ (3.3)
$\displaystyle{B}^{(IJK)}$
$\displaystyle=\frac{1}{6}\,\left({B}^{IJK}+{B}^{KJI}+{B}^{KIJ}+{B}^{JIK}+{B}^{JKI}+{B}^{IKJ}\right)\;.$
(3.4)

We will also use the short-hand notation

$\displaystyle P^{I}\equiv P^{II},\quad B^{I}\equiv B^{III}\;,$ (3.5)

where convenient.

We define the estimators:

$\displaystyle\hat{P}^{IJ}(\bm{k})$
$\displaystyle=\frac{1}{V_{s}\,V_{k}}\,\sum_{\bm{q}\in
k}\delta^{I}(\bm{q})\,\delta^{J}(-\bm{q})\;,$ (3.6)
$\displaystyle\hat{B}^{(IJK)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle=\frac{1}{V_{s}\,V_{123}}\sum_{\bm{q}_{1}\in
k_{1}}\sum_{\bm{q}_{2}\in k_{2}}\sum_{\bm{q}_{3}\in k_{3}}\delta_{\rm
K}(\bm{q}_{123})\,\delta^{I}(\bm{q}_{1})\,\delta^{J}(\bm{q}_{2})\,\delta^{K}(\bm{q}_{3})\;,$
(3.7)

where $V_{s}$ is the volume of the survey and
$\bm{q}_{123}\equiv\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3}$. For the power spectra
the sum runs over all wavenumbers $\bm{q}$ in the spherical shell of radius
$k$ and width $\Delta k$. For the bispectra, the Kronecker delta symbol
ensures the formation of ‘fundamental triangles’ with sides $\bm{q}_{i}$,
satisfying the condition $\bm{q}_{123}=0$, that fall in the ‘triangular bin’
defined by the triplet of bin centres ($k_{1}$, $k_{2}$, $k_{3}$) and width
$\Delta k$. The volume in Fourier space of the power spectra shells and the
bispectra fundamental triangle bins is, in the thin shell limit [70],

$\displaystyle V_{k}=4\pi\,k^{2}\,\Delta k\,,\quad
V_{123}\equiv\sum_{\bm{q}_{1}\in k_{1}}\sum_{\bm{q}_{2}\in
k_{2}}\sum_{\bm{q}_{3}\in k_{3}}\delta_{\rm
K}(\bm{q}_{123})={4\pi^{2}}\,k_{1}\,k_{2}\,k_{3}\,\Delta k_{1}\,\Delta
k_{2}\,\Delta k_{3}\,.$ (3.8)

The estimators in Eqs. 3.6, LABEL: and 3.7 in the case of one sample ($I=J=K$)
reduce to the well known auto-correlation results.

### The two-tracer case

From now on we consider the case of two tracers, $t$ and $t^{\prime}$, so that
$I,J,K=t,t^{\prime}$. The data vector of the power spectra is

$\bm{D_{P}}=\Big{[}P^{tt},P^{tt^{\prime}},P^{t^{\prime}t^{\prime}}\Big{]}\;.$
(3.9)

The expected values of the spectra are:

$\langle\hat{P}^{IJ}\rangle\equiv P^{IJ}(\bm{k},z)={\cal
D}_{P}^{IJ}(\bm{k},z)\,Z_{1}^{I}(\bm{k},z)\,Z_{1}^{J}(\bm{k},z)\,P(k,z)+\delta^{IJ}\,P^{I}_{\varepsilon}(z)\,,$
(3.10)

where $Z_{1}^{I}$ and ${\cal D}_{P}^{IJ}$ are given in Appendix A.

The data vector of the bispectrum is

$\bm{D_{B}}=\Big{[}B^{ttt},B^{(ttt^{\prime})},B^{(tt^{\prime}t^{\prime})},B^{t^{\prime}t^{\prime}t^{\prime}}\Big{]}\;.$
(3.11)

In a tree-level description, we consider terms up to second order PT theory,
RSD and bias, while keeping only linear $f_{\rm NL}$ terms in the perturbative
expansion. Then the expected values of the estimator for each bispectrum data
vector entry are:

$\displaystyle\langle\hat{B}^{ttt}\rangle\equiv
B^{ttt}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$ $\displaystyle={\cal
D}_{B}^{ttt}\,\bigg{[}Z_{1}^{t}(\bm{k}_{1})\,Z_{1}^{t}(\bm{k}_{2})\,Z_{1}^{t}(\bm{k}_{3})\,B_{\rm
PNG}(k_{1},k_{2},k_{3})$
$\displaystyle+\Big{\\{}2\,Z_{1}^{t}(\bm{k}_{1})\,Z_{1}^{t}(\bm{k}_{2})\,Z_{2}^{t}(\bm{k}_{1},\bm{k}_{2})\,P(k_{1})\,P(k_{2})+2\,\bm{k}\text{-perm}\Big{\\}}\bigg{]}$
$\displaystyle+2\,P^{t}_{\varepsilon\varepsilon_{\delta}}\,\Big{[}Z_{1}^{t}(\bm{k}_{1})\,P(k_{1})+2\,\bm{k}\text{-perm}\Big{]}+B^{t}_{\varepsilon}\;,$
(3.12) $\displaystyle\langle\hat{B}^{(ttt^{\prime})}\rangle\equiv
B^{(ttt^{\prime})}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle=\frac{1}{3}\,\bigg{[}\Big{\\{}{\cal
D}_{B}^{ttt^{\prime}}\,Z_{1}^{t}(\bm{k}_{1})\,Z_{1}^{t}(\bm{k}_{2})\,Z_{1}^{t^{\prime}}(\bm{k}_{3})+{\cal
D}_{B}^{tt^{\prime}t}\,Z_{1}^{t}(\bm{k}_{1})\,Z_{1}^{t^{\prime}}(\bm{k}_{2})\,Z_{1}^{t}(\bm{k}_{3})$
$\displaystyle+{\cal
D}_{B}^{t^{\prime}tt}Z_{1}^{t^{\prime}}(\bm{k}_{1})\,Z_{1}^{t}(\bm{k}_{2})\,Z_{1}^{t}(\bm{k}_{3})\Big{\\}}\,B_{\rm
PNG}(k_{1},k_{2},k_{3})$ $\displaystyle+2\,\Big{\\{}\big{[}{\cal
D}_{B}^{tt^{\prime}t}\,Z_{1}^{t}(\bm{k}_{1})\,Z_{1}^{t^{\prime}}(\bm{k}_{2})+{\cal
D}_{B}^{t^{\prime}tt}\,Z_{1}^{t^{\prime}}(\bm{k}_{1})\,Z_{1}^{t}(\bm{k}_{2})\big{]}\,Z_{2}^{t}(\bm{k}_{1},\bm{k}_{2})$
$\displaystyle+{\cal
D}_{B}^{ttt^{\prime}}\,Z_{1}^{t}(\bm{k}_{1})\,Z_{1}^{t}(\bm{k}_{2})\,Z_{2}^{t^{\prime}}(\bm{k}_{1},\bm{k}_{2})\Big{\\}}\,P(k_{1})\,P(k_{2})+2\,\bm{k}\text{-perm}\bigg{]}$
$\displaystyle+\frac{1}{3}\,P^{t}_{\varepsilon\varepsilon_{\delta}}\,\Big{[}Z^{t^{\prime}}(\bm{k}_{1})\,P(k_{1})+2\,\bm{k}\text{-perm}\Big{]}\;,$
(3.13) $\displaystyle\langle\hat{B}^{(tt^{\prime}t^{\prime})}\rangle\equiv
B^{(tt^{\prime}t^{\prime})}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle=\frac{1}{3}\,\bigg{[}\Big{\\{}{\cal
D}_{B}^{tt^{\prime}t^{\prime}}\,Z_{1}^{t}(\bm{k}_{1})\,Z_{1}^{t^{\prime}}(\bm{k}_{2})\,Z_{1}^{t^{\prime}}(\bm{k}_{3})+{\cal
D}_{B}^{t^{\prime}t^{\prime}t}\,Z_{1}^{t^{\prime}}(\bm{k}_{1})\,Z_{1}^{t^{\prime}}(\bm{k}_{2})\,Z_{1}^{t}(\bm{k}_{3})$
$\displaystyle+{\cal
D}_{B}^{t^{\prime}tt^{\prime}}\,Z_{1}^{t^{\prime}}(\bm{k}_{1})\,Z_{1}^{t}(\bm{k}_{2})\,Z_{1}^{t^{\prime}}(\bm{k}_{3})\Big{\\}}\,B_{\rm
PNG}(k_{1},k_{2},k_{3})$ $\displaystyle+2\,\Big{\\{}\big{[}{\cal
D}_{B}^{tt^{\prime}t^{\prime}}\,Z_{1}^{t}(\bm{k}_{1})\,Z_{1}^{t^{\prime}}(\bm{k}_{2})+{\cal
D}_{B}^{t^{\prime}tt^{\prime}}\,Z_{1}^{t^{\prime}}(\bm{k}_{1})\,Z_{1}^{t}(\bm{k}_{2})\big{]}\,Z_{2}^{t^{\prime}}(\bm{k}_{1},\bm{k}_{2})$
$\displaystyle+{\cal
D}_{B}^{t^{\prime}t^{\prime}t}\,Z_{1}^{t^{\prime}}(\bm{k}_{1})\,Z_{1}^{t^{\prime}}(\bm{k}_{2})\,Z_{2}^{t}(\bm{k}_{1},\bm{k}_{2})\Big{\\}}P(k_{1})\,P(k_{2})+2\,\bm{k}\text{-perm}\bigg{]}$
$\displaystyle+\frac{1}{3}\,{P^{t^{\prime}}_{\varepsilon\varepsilon_{\delta}}}\,\Big{[}Z^{t}(\bm{k}_{1})\,P(k_{1})+2\,\bm{k}\text{-perm}\Big{]}\;,$
(3.14)
$\displaystyle\langle\hat{B}^{t^{\prime}t^{\prime}t^{\prime}}\rangle\equiv
B^{t^{\prime}t^{\prime}t^{\prime}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle={\cal
D}_{B}^{t^{\prime}t^{\prime}t^{\prime}}\,\bigg{[}Z_{1}^{t^{\prime}}(\bm{k}_{1})\,Z_{1}^{t^{\prime}}(\bm{k}_{2})\,Z_{1}^{t^{\prime}}(\bm{k}_{3})\,B_{\rm
PNG}(k_{1},k_{2},k_{3})$
$\displaystyle+\Big{\\{}2\,Z_{1}^{t^{\prime}}(\bm{k}_{1})\,Z_{1}^{t^{\prime}}(\bm{k}_{2})\,Z_{2}^{t^{\prime}}(\bm{k}_{1},\bm{k}_{2})\,P(k_{1})\,P(k_{2})+2\,\bm{k}\text{-perm}\Big{\\}}\bigg{]}$
$\displaystyle+2\,P^{t^{\prime}}_{\varepsilon\varepsilon_{\delta}}\,\Big{[}Z_{1}^{t^{\prime}}(\bm{k}_{1})\,P(k_{1})+2\,\bm{k}\text{-perm}\Big{]}+B^{t^{\prime}}_{\varepsilon}\;,$
(3.15)

where $Z_{2}^{I}$ and ${\cal D}_{B}^{IJK}$ are given in Appendix A. Note that
we omitted the redshift dependence in all terms, and also the $k$-dependence
in ${\cal D}_{B}^{IJK}$. The stochastic terms have fiducial values:

$P_{\varepsilon}^{I}=\frac{1}{\overline{n}_{I}}\;,~{}~{}P_{\varepsilon\varepsilon_{\delta}}^{I}=\frac{\;b_{1}^{I}}{2\,\overline{n}_{I}}\;,~{}~{}B_{\varepsilon}^{I}=\frac{1}{\overline{n}_{I}^{2}}\;.$
(3.16)

We assume that the cross shot-noise of the power spectra and bispectra is
zero:

$P^{IJ}_{\varepsilon}=P_{\varepsilon\varepsilon_{\delta}}^{IJ}=0~{}~{}(I\neq
J)\;,\quad B^{IJK}_{\varepsilon}=0~{}~{}(I,\,J,\,K~{}\mbox{unequal})\;.$
(3.17)

The extension of the power spectrum and bispectrum multi-tracers from 2 to 3
tracers is given in Appendix B.

##  4 Multi-tracer covariance

The multi-tracer combination of power spectra and bispectra for two tracers
has data vector

$\bm{D_{PB}}=\Big{[}P^{tt},P^{tt^{\prime}},P^{t^{\prime}t^{\prime}},B^{ttt},B^{(ttt^{\prime})},B^{(tt^{\prime}t^{\prime})},B^{t^{\prime}t^{\prime}t^{\prime}}\Big{]}\;.$
(4.1)

The covariance of the data vector contains four blocks, coming from the
covariance of the power spectra, the bispectra and the power spectra –
bispectra covariance. In matrix form:

$\bm{\mathrm{Cov}_{PB}}=\begin{bmatrix}\mathrm{Cov}(\bm{P}_{i},\bm{P}_{j})~{}&&\mathrm{Cov}(\bm{P}_{i},\bm{B}_{{b}})\\\
&&\\\ \mathrm{Cov}(\bm{P}_{i},\bm{B}_{{b}})^{\sf
T}&&\mathrm{Cov}(\bm{B}_{{a}},\bm{B}_{{b}})\end{bmatrix}\;,$ (4.2)

where

$\displaystyle\bm{P}_{i}=\bm{P}(\bm{k}_{i})\;,\quad\bm{B}_{{a}}=\bm{B}(\bm{k}_{a_{1}},\bm{k}_{a_{2}},\bm{k}_{a_{3}})\;.$
(4.3)

We assume that the covariance between two power spectra is Gaussian, i.e. the
off-diagonal contributions are neglected. Then the cross-sample covariance of
the power spectra is a symmetric 3$\times$3 diagonal block matrix:

$\displaystyle\mathrm{Cov}(\bm{P}_{i},\bm{P}_{j})$
$\displaystyle=\begin{pmatrix}\mathrm{Cov}\big{(}P_{i}^{tt},P_{j}^{tt}\big{)}&\mathrm{Cov}\big{(}P_{i}^{tt},P_{j}^{tt^{\prime}}\big{)}&\mathrm{Cov}\big{(}P_{i}^{tt},P_{j}^{t^{\prime}t^{\prime}}\big{)}\\\
&&\\\
&\mathrm{Cov}\big{(}P_{i}^{tt^{\prime}},P_{j}^{tt^{\prime}}\big{)}&\mathrm{Cov}\big{(}P_{i}^{tt^{\prime}},P_{j}^{t^{\prime}t^{\prime}}\big{)}\\\
&&\\\
&&\mathrm{Cov}\big{(}P_{i}^{t^{\prime}t^{\prime}},P_{j}^{t^{\prime}t^{\prime}}\big{)}\\\
\end{pmatrix}$ (4.4)
$\displaystyle=2\,\frac{(2\,\pi)^{3}}{V_{s}\,V_{k}}\,\delta_{ij}\begin{pmatrix}\big{[}P^{tt}(\bm{k}_{i})\big{]}^{2}&P^{tt}(\bm{k}_{i})\,P^{tt^{\prime}}(\bm{k}_{i})&\big{[}P^{tt^{\prime}}(\bm{k}_{i})\big{]}^{2}\\\
&&\\\
&\frac{1}{2}\,\big{\\{}P^{tt}(\bm{k}_{i})\,P^{t^{\prime}t^{\prime}}(\bm{k}_{i})+\big{[}P^{tt^{\prime}}(\bm{k}_{i})\big{]}^{2}\big{\\}}&P^{t^{\prime}t^{\prime}}(\bm{k}_{i})\,P^{tt^{\prime}}(\bm{k}_{i})\\\
&&\\\
&&\big{[}P^{t^{\prime}t^{\prime}}(\bm{k}_{i})\big{]}^{2}\end{pmatrix}\;,$

where each sub-block is a diagonal matrix in $\bm{k}_{i}$.

The covariance of the symmetrised bispectra is a symmetric 4$\times$4 block
matrix:

$\displaystyle\mathrm{Cov}(\bm{B}_{{a}},\bm{B}_{{b}})=\frac{(2\,\pi)^{6}}{3\,V_{s}\,V_{123}}\,s_{123}\,\delta_{ab}\times$
(4.5)
$\displaystyle\begin{pmatrix}\mathrm{Cov}\big{(}B_{{a}}^{ttt},B_{{b}}^{ttt}\big{)}&\mathrm{Cov}\big{(}B_{{a}}^{ttt},B_{{b}}^{(ttt^{\prime})}\big{)}&\mathrm{Cov}\big{(}B_{{a}}^{ttt},B_{{b}}^{(tt^{\prime}t^{\prime})}\big{)}&\mathrm{Cov}\big{(}B_{{a}}^{ttt},B_{{b}}^{t^{\prime}t^{\prime}t^{\prime}}\big{)}\\\
&&&\\\
&\mathrm{Cov}\big{(}B_{{a}}^{(ttt^{\prime})},B_{{b}}^{(ttt^{\prime})}\big{)}&\mathrm{Cov}\big{(}B_{{a}}^{(ttt^{\prime})},B_{{b}}^{(tt^{\prime}t^{\prime})}\big{)}&\mathrm{Cov}\big{(}B_{{a}}^{(ttt^{\prime})},B_{{b}}^{t^{\prime}t^{\prime}t^{\prime}}\big{)}\\\
&&&\\\
&&\mathrm{Cov}\big{(}B_{{a}}^{(tt^{\prime}t^{\prime})},B_{{b}}^{(tt^{\prime}t^{\prime})}\big{)}&\mathrm{Cov}\big{(}B_{{a}}^{(tt^{\prime}t^{\prime})},B_{{b}}^{t^{\prime}t^{\prime}t^{\prime}}\big{)}\\\
&&&\\\
&&&\mathrm{Cov}\big{(}B_{{a}}^{t^{\prime}t^{\prime}t^{\prime}},B_{{b}}^{t^{\prime}t^{\prime}t^{\prime}}\big{)}\end{pmatrix}\;,$

where $s_{123}=6,2,1$ for equilateral, isosceles and scalene triangles
respectively. Each block is an $N_{T}\times N_{T}$ matrix, where $N_{T}$ is
the total number of triangles formed. The index ${{a}}$ represents the
triangle triplet $(\bm{k}_{a_{1}},\bm{k}_{a_{2}},\bm{k}_{a_{3}})$. The
Kronecker delta indicates that each sub-block of the bispectrum-bispectrum
covariance is a diagonal matrix, i.e. we only consider the Gaussian
contribution.

These sub-block matrices are:

$\displaystyle\mathrm{Cov}\big{(}B_{a}^{ttt},B_{b}^{ttt}\big{)}$
$\displaystyle=3\,P^{tt}_{a_{1}}\,P^{tt}_{a_{2}}\,P^{tt}_{a_{3}},$
$\displaystyle\mathrm{Cov}\big{(}B_{a}^{ttt},B_{b}^{(ttt^{\prime})}\big{)}$
$\displaystyle=P^{tt}_{a_{1}}\,P^{tt}_{a_{2}}\,P^{tt^{\prime}}_{a_{3}}+\text{2\,perm}\;,$
$\displaystyle\mathrm{Cov}\big{(}B_{a}^{ttt},B_{b}^{(tt^{\prime}t^{\prime})}\big{)}$
$\displaystyle=P^{tt}_{a_{1}}\,P^{tt^{\prime}}_{a_{2}}\,P^{tt^{\prime}}_{a_{3}}+\text{2\,perm}\;,$
$\displaystyle\mathrm{Cov}\big{(}B_{a}^{ttt},B_{b}^{t^{\prime}t^{\prime}t^{\prime}}\big{)}$
$\displaystyle=3\,P^{tt^{\prime}}_{a_{1}}\,P^{tt^{\prime}}_{a_{2}}\,P^{tt^{\prime}}_{a_{3}},$
$\displaystyle\mathrm{Cov}\big{(}B_{a}^{(ttt^{\prime})},B_{b}^{(ttt^{\prime})}\big{)}$
$\displaystyle=\frac{1}{3}\,\big{(}P^{tt}_{a_{1}}\,P^{tt}_{a_{2}}\,P^{t^{\prime}t^{\prime}}_{a_{3}}+2\,P^{tt}_{a_{1}}\,P^{tt^{\prime}}_{a_{2}}\,P^{tt^{\prime}}_{a_{3}}\big{)}+\text{2\,perm}\;,$
$\displaystyle\mathrm{Cov}\big{(}B_{a}^{(ttt^{\prime})},B_{b}^{(tt^{\prime}t^{\prime})}\big{)}$
$\displaystyle=\frac{1}{3}\,\big{(}P^{tt^{\prime}}_{a_{1}}\,P^{tt^{\prime}}_{a_{2}}\,P^{tt^{\prime}}_{a_{3}}+P^{tt}_{a_{1}}\,P^{tt^{\prime}}_{a_{2}}\,P^{t^{\prime}t^{\prime}}_{a_{3}}+P^{t^{\prime}t^{\prime}}_{a_{1}}\,P^{tt^{\prime}}_{a_{2}}\,P^{tt}_{a_{3}}\big{)}+\text{2\,perm}\;,$
$\displaystyle\mathrm{Cov}\big{(}B_{a}^{(ttt^{\prime})},B_{b}^{t^{\prime}t^{\prime}t^{\prime}}\big{)}$
$\displaystyle=P^{tt^{\prime}}_{a_{1}}\,P^{tt^{\prime}}_{a_{2}}\,P^{t^{\prime}t^{\prime}}_{a_{3}}+\text{2\,perm}\;,$
$\displaystyle\mathrm{Cov}\big{(}B_{a}^{(tt^{\prime}t^{\prime})},B_{b}^{(tt^{\prime}t^{\prime})}\big{)}$
$\displaystyle=\frac{1}{3}\,\big{(}P^{tt}_{a_{1}}\,P^{t^{\prime}t^{\prime}}_{a_{2}}\,P^{t^{\prime}t^{\prime}}_{a_{3}}+2\,P^{tt^{\prime}}_{a_{1}}\,P^{tt^{\prime}}_{a_{2}}\,P^{t^{\prime}t^{\prime}}_{a_{3}}\big{)}+\text{2\,perm}\;,$
$\displaystyle\mathrm{Cov}\big{(}B_{a}^{(tt^{\prime}t^{\prime})},B_{b}^{t^{\prime}t^{\prime}t^{\prime}}\big{)}$
$\displaystyle=P^{tt^{\prime}}_{a_{1}}\,P^{t^{\prime}t^{\prime}}_{a_{2}}\,P^{t^{\prime}t^{\prime}}_{a_{3}}+\text{2\,perm}\;,$
$\displaystyle\mathrm{Cov}\big{(}B_{a}^{t^{\prime}t^{\prime}t^{\prime}},B_{b}^{t^{\prime}t^{\prime}t^{\prime}}\big{)}$
$\displaystyle=3\,P^{t^{\prime}t^{\prime}}_{a_{1}}\,P^{t^{\prime}t^{\prime}}_{a_{2}}\,P^{t^{\prime}t^{\prime}}_{a_{3}}\;.$
(4.6)

Finally, the cross-covariance between the power spectra and bispectra is
itself a 3$\times$4 block matrix:

$\displaystyle\mathrm{Cov}(\bm{P}_{i},\bm{B}_{b})=\frac{(2\,\pi)^{3}}{V_{s}\,V_{k}}$
(4.7) $\displaystyle\times$
$\displaystyle\begin{pmatrix}\mathrm{Cov}\big{(}P_{i}^{tt},B_{b}^{ttt}\big{)}&\mathrm{Cov}\big{(}P_{i}^{tt},B_{b}^{(ttt^{\prime})}\big{)}&\mathrm{Cov}\big{(}P_{i}^{tt},B_{b}^{(tt^{\prime}t^{\prime})}\big{)}&\mathrm{Cov}\big{(}P_{i}^{tt},B_{b}^{t^{\prime}t^{\prime}t^{\prime}}\big{)}\\\
&&&\\\
\mathrm{Cov}\big{(}P_{i}^{tt^{\prime}},B_{b}^{ttt}\big{)}&\mathrm{Cov}\big{(}P_{i}^{tt^{\prime}},B_{b}^{(ttt^{\prime})}\big{)}&\mathrm{Cov}\big{(}P_{i}^{tt^{\prime}},B_{b}^{(tt^{\prime}t^{\prime})}\big{)}&\mathrm{Cov}\big{(}P_{i}^{tt^{\prime}},B_{b}^{t^{\prime}t^{\prime}t^{\prime}}\big{)}\\\
&&&\\\
\mathrm{Cov}\big{(}P_{i}^{t^{\prime}t^{\prime}},B_{b}^{ttt}\big{)}&\mathrm{Cov}\big{(}P_{i}^{t^{\prime}t^{\prime}},B_{b}^{(ttt^{\prime})}\big{)}&\mathrm{Cov}\big{(}P_{i}^{t^{\prime}t^{\prime}},B_{b}^{(tt^{\prime}t^{\prime})}\big{)}&\mathrm{Cov}\big{(}P_{i}^{t^{\prime}t^{\prime}},B_{b}^{t^{\prime}t^{\prime}t^{\prime}}\big{)}\end{pmatrix}\;,$

where

$\displaystyle\mathrm{Cov}\big{(}P_{i}^{tt},B_{b}^{ttt}\big{)}$
$\displaystyle={2\,\big{(}\delta_{ib_{1}}+\delta_{ib_{2}}+\delta_{ib_{3}}\big{)}\,P^{tt}({\bm{k}}_{i})\,B^{ttt}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})}\;,$
$\displaystyle\mathrm{Cov}\big{(}P_{i}^{tt},B_{b}^{(ttt^{\prime})}\big{)}$
$\displaystyle=\frac{2}{3}\,\Big{[}\big{(}\delta_{ib_{1}}+\delta_{ib_{2}}+\delta_{ib_{3}}\big{)}\,P^{tt^{\prime}}(\bm{k}_{i})\,B^{ttt}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})$
$\displaystyle~{}~{}~{}+P^{tt}(\bm{k}_{i})B_{2\delta}^{ttt^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\Big{]}\;,$
$\displaystyle\mathrm{Cov}\big{(}P_{i}^{tt},B_{b}^{(tt^{\prime}t^{\prime})}\big{)}$
$\displaystyle=\frac{2}{3}\,\Big{[}P^{tt}(\bm{k}_{i})\,B_{1\delta}^{tt^{\prime}t^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})+P^{tt^{\prime}}(\bm{k}_{i})\,B_{2\delta}^{ttt^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\Big{]}\;,$
$\displaystyle\mathrm{Cov}\big{(}P_{i}^{tt},B_{b}^{t^{\prime}t^{\prime}t^{\prime}}\big{)}$
$\displaystyle=2\,P^{tt^{\prime}}(\bm{k}_{i})\,B_{1\delta}^{tt^{\prime}t^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\;,$
$\displaystyle\mathrm{Cov}\big{(}P_{i}^{tt^{\prime}},B_{b}^{ttt}\big{)}$
$\displaystyle=\big{(}\delta_{ib_{1}}+\delta_{ib_{2}}+\delta_{ib_{3}}\big{)}\,P^{tt^{\prime}}(\bm{k}_{i})\,B^{ttt}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})+P^{tt}(\bm{k}_{i})\,B_{1\delta}^{ttt^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\;,$
$\displaystyle\mathrm{Cov}\big{(}P_{i}^{tt^{\prime}},B_{b}^{(ttt^{\prime})}\big{)}$
$\displaystyle=\frac{1}{3}\,\bigg{\\{}\big{(}\delta_{ib_{1}}+\delta_{ib_{2}}+\delta_{ib_{3}}\big{)}\,\Big{[}P^{tt^{\prime}}(\bm{k}_{i})\,B^{(ttt^{\prime})}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})$
$\displaystyle~{}~{}~{}+P^{t^{\prime}t^{\prime}}(\bm{k}_{i})\,B^{ttt}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\Big{]}+P^{tt}(\bm{k}_{i})\,B_{2\delta}^{tt^{\prime}t^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\bigg{\\}}\;,$
$\displaystyle\mathrm{Cov}\big{(}P_{i}^{tt^{\prime}},B_{b}^{(tt^{\prime}t^{\prime})}\big{)}$
$\displaystyle=\frac{1}{3}\,\bigg{\\{}\big{(}\delta_{ib_{1}}+\delta_{ib_{2}}+\delta_{ib_{3}}\big{)}\,\Big{[}P^{tt^{\prime}}(\bm{k}_{i})\,B^{(tt^{\prime}t^{\prime})}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})$
$\displaystyle~{}~{}~{}+P^{t^{\prime}t^{\prime}}(\bm{k}_{i})\,B^{ttt}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\Big{]}+P^{t^{\prime}t^{\prime}}(\bm{k}_{i})\,B_{2\delta}^{ttt^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\bigg{\\}}\;,$
$\displaystyle\mathrm{Cov}\big{(}P_{i}^{tt^{\prime}},B_{b}^{t^{\prime}t^{\prime}t^{\prime}}\big{)}$
$\displaystyle=\big{(}\delta_{ib_{1}}+\delta_{ib_{2}}+\delta_{ib_{3}}\big{)}\,P^{tt^{\prime}}(\bm{k}_{i})\,B^{t^{\prime}t^{\prime}t^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})$
$\displaystyle~{}~{}~{}+P^{t^{\prime}t^{\prime}}(\bm{k}_{i})\,B_{1\delta}^{tt^{\prime}t^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\;,$
$\displaystyle\mathrm{Cov}\big{(}P_{i}^{t^{\prime}t^{\prime}},B_{b}^{ttt}\big{)}$
$\displaystyle=2P^{tt^{\prime}}(\bm{k}_{i})\,B_{1\delta}^{ttt^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\;,$
$\displaystyle\mathrm{Cov}\big{(}P_{i}^{t^{\prime}t^{\prime}},B_{b}^{(ttt^{\prime})}\big{)}$
$\displaystyle=\frac{2}{3}\,\Big{[}P^{t^{\prime}t^{\prime}}(\bm{k}_{i})\,B_{1\delta}^{ttt^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})+P^{tt^{\prime}}(\bm{k}_{i})\,B_{2\delta}^{tt^{\prime}t^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\Big{]}\;,$
$\displaystyle\mathrm{Cov}\big{(}P_{i}^{t^{\prime}t^{\prime}},B_{b}^{(tt^{\prime}t^{\prime})}\big{)}$
$\displaystyle=\frac{2}{3}\,\Big{[}\big{(}\delta_{ib_{1}}+\delta_{ib_{2}}+\delta_{ib_{3}}\big{)}\,P^{tt^{\prime}}(\bm{k}_{i})\,B^{t^{\prime}t^{\prime}t^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})$
$\displaystyle~{}~{}~{}+P^{t^{\prime}t^{\prime}}(\bm{k}_{i})\,B_{2\delta}^{tt^{\prime}t^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\Big{]}\;,$
$\displaystyle\mathrm{Cov}\big{(}P_{i}^{t^{\prime}t^{\prime}},B_{b}^{t^{\prime}t^{\prime}t^{\prime}}\big{)}$
$\displaystyle=2\,\big{(}\delta_{ib_{1}}+\delta_{ib_{2}}+\delta_{ib_{3}}\big{)}\,P^{t^{\prime}t^{\prime}}(\bm{k}_{i})\,B^{t^{\prime}t^{\prime}t^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\;.$
(4.8)

Here the partially symmetric bispectra are given by:

$\displaystyle
B_{1\delta}^{ttt^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})$
$\displaystyle=\delta_{ib_{1}}\,B^{t^{\prime}tt}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})+\delta_{ib_{2}}\,B^{tt^{\prime}t}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})$
$\displaystyle~{}~{}~{}+\delta_{ib_{3}}\,B^{ttt^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\;,$
(4.9) $\displaystyle
B_{2\delta}^{ttt^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})$
$\displaystyle=[\delta_{ib_{1}}+\delta_{ib_{2}}]\,B^{ttt^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})+[\delta_{ib_{1}}+\delta_{ib_{3}}]\,B^{tt^{\prime}t}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})$
$\displaystyle~{}~{}~{}+[\delta_{ib_{2}}+\delta_{ib_{3}}]\,B^{t^{\prime}tt}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\;.$
(4.10)

These bispectrum tracer indices are used for convenience and can be
generalised for any 2-tracer combination, e.g. $tt^{\prime}t^{\prime}$. Note
that a nonzero power spectrum-bispectrum cross-covariance is generated by
configurations in which the power spectrum $\bm{k}$ mode coincides at least
with one of the sides of a triangle configuration. The results presented above
for the cross-covariance between power spectrum and bispectrum, consider only
the dominant contribution, while the 5-point connected correlation function
part is neglected (see [71] for a discussion). Note that the block matrix of
Section 4 and each of the sub-blocks are not square matrices.

The expressions for 3 tracers are in Appendix B.

##  5 Fisher information on $f_{\rm NL}$

The previous section presented the formalism for the multi-tracer power
spectra, the multi-tracer bispectra and their combination. In this section, we
conduct separate investigations into the precision from the multi-tracer power
spectra and the multi-tracer bispectra, without combining the information
obtained from the two correlators. A follow-up paper to perform more detailed
forecasts will include the combination of multi-tracers.

The Fisher information on $f_{\rm NL}$ in each redshift bin for the multi-
tracer power spectra is given by

$\displaystyle{F_{f_{\rm NL}}^{\bm{P}}}$ $\displaystyle=\int{\rm
d}\bm{k}\,\frac{\partial\bm{D_{P}}(\bm{k})}{\partial f_{\rm
NL}}\,\bm{\mathrm{Cov}_{P}}^{-1}{(\bm{k})}\,\frac{\partial\bm{D_{P}}^{{\sf
T}}(\bm{k})}{\partial f_{\rm NL}}\;.$ (5.1)

For a rough estimate of the improvements on precision that are produced by the
multi-tracer, we consider the error $\sigma^{\bm{P}}(f_{\rm
NL})=(F^{\bm{P}}_{f_{\rm NL}})^{-1/2}$. Then an analytic expression can be
derived for the single-tracer power spectrum:

$\sigma^{I}_{f_{\rm NL}}=N_{P}^{1/2}\,\int{\rm
d}\bm{k}\,\frac{M(k)\,Z_{1}^{I}(\bm{k})}{2\,{b_{\Psi}^{I}}\,k^{\alpha}}\,\frac{{P^{I}}(\bm{k})}{P^{I}(\bm{k})-P_{\varepsilon}^{I}},\qquad
N_{P}=2\,\frac{(2\,\pi)^{3}}{V_{s}\,V_{k}}\;,$ (5.2)

where we omitted the redshift dependence. Here $\alpha=1,0$ for orthogonal and
local shapes, respectively [55, 56, 72, 61, 73, 44, 74]. Note that equilateral
PNG ($\alpha=2$) cannot be constrained using the galaxy power spectrum, as
shown in [65]. Using the linear redshift space kernel Eq. A.1, the analytic
result Eq. 5.2 shows that, for a fixed number density and $\mu$, the tightest
constraints on $f_{\rm NL}$ can be obtained from a sample that has
$b_{\Psi}^{I}\gg b_{1}^{I}$.

For the multi-tracer case the analytic expression of the Fisher matrix on
$f_{\mathrm{NL}}^{\mathrm{loc}}$ is given by:

$\displaystyle F^{\bm{P}}_{f_{\rm NL}}$ $\displaystyle=\sum_{I,J}\int{\rm
d}\bm{k}\,{\cal
N}^{-1}(\bm{k})\,\bigg{\\{}{4\,P(k)^{2}(b_{\Psi}^{J})^{2}\,\big{[}P_{\varepsilon}^{I}+P_{\varepsilon}^{J}\,R(\bm{k})^{2}\big{]}^{2}\,\left[Z_{1}^{J}(\bm{k})\right]^{2}}$
$\displaystyle+{8\,P(k)^{2}\,P_{\varepsilon}^{J}\,b_{\Psi}^{J}\,\left[Z_{1,\rm
G}^{J}(\bm{k})\,b_{\Psi}^{I}-Z_{1,\rm
G}^{I}(\bm{k})\,b_{\Psi}^{J}\,\right]\,\big{[}P_{\varepsilon}^{I}+P_{\varepsilon}^{J}\,R(\bm{k})^{2}\big{]}\,R(\bm{k})\,Z_{1}^{J}(\bm{k})}$
$\displaystyle+{2\,P(k)^{2}\left[Z_{1,\rm
G}^{J}({\bm{k}})\,b_{\Psi}^{I}-Z_{1,\rm
G}^{I}({\bm{k}})\,b_{\Psi}^{J}\right]^{2}\,\big{[}P_{\varepsilon}^{I}\,P^{J}(\bm{k})+P_{\varepsilon}^{J}\,\big{\\{}P_{\varepsilon}^{J}+P^{J}(\bm{k})\big{\\}}\,R^{2}{(\bm{k})}\big{]}}\bigg{\\}}\;,$
(5.3)

where

$\displaystyle{\cal
N}(\bm{k})=N_{P}\,M(k)^{2}\,k^{-2\alpha}\,\big{[}P_{\varepsilon}^{I}\,P^{J}(\bm{k})+P_{\varepsilon}^{J}\,\big{\\{}P^{J}(\bm{k})-P_{\varepsilon}^{J}\big{\\}}\,R(\bm{k})^{2}\big{]}^{2}\quad\mbox{with}~{}R(\bm{k})=\frac{Z_{1}^{I}(\bm{k})}{Z_{1}^{J}(\bm{k})}\;.$
(5.4)

Here $Z_{1,\rm G}^{I}(\bm{k})$ is the Gaussian version of Eq. A.1, i.e. for
$f_{\rm NL}=0$. Note that this expression is valid for a general linear
redshift kernel (e.g. including non-integrated relativistic corrections).

Our Section 5 generalises the result of [40], which only considers local PNG,
takes the cosmic variance limit and neglects RSD effects. At fixed redshift
bin and scales, Section 5 indicates that precision on $f_{\rm NL}$ is
increased by

$\mbox{increasing}~{}~{}\big{|}Z_{1,\rm
G}^{J}({\bm{k}})\,b_{\Psi}^{I}-Z_{1,\rm
G}^{I}({\bm{k}})\,b_{\Psi}^{J}\big{|}~{}~{}\mbox{and}~{}~{}R(\bm{k}).$ (5.5)

One way to maximise the former is by choosing appropriately the line-of-sight
angle $\mu$ to gain a boost from the presence of RSD, which is also evident in
forecasts that utilise single tracers (e.g. [75]). If we neglect RSD (or fix
$\mu$) and consider local PNG, then Eq. 5.5 implies maximising
$|b_{1}^{J}\,b_{\Phi}^{I}-b_{1}^{I}\,b_{\Phi}^{J}|$, as shown also in [40].

For the multi-tracer bispectra, the Fisher matrix is

$\displaystyle{F_{f_{\rm NL}}^{\bm{B}}}$ $\displaystyle={\int{\rm
d}^{3}\bm{k}}\frac{\partial\bm{D_{B}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})}{\partial
f_{\rm
NL}}\,\bm{\mathrm{Cov}_{B}}^{-1}{(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})}\,\frac{\partial\bm{D_{B}}^{{\sf
T}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})}{\partial f_{\rm NL}}\;.$ (5.6)

Analytical calculations for the bispectrum multi-tracer are long and
uninformative, without assuming some limit, either for the shot noise or the
triangle configurations. Consequently, we perform the bispectrum analysis
numerically.

Figure 1:  Cumulative $1\sigma$ error on $f_{\rm NL}$ as a function of
redshift, for the three types of PNG, comparing single- and multi-tracer
cases, for power spectra and bispectra. Here $t$ and $t^{\prime}$ are tracers
of nominal next-generation surveys.

We consider two nominal spectroscopic surveys, with the following features:
$\bullet$ Overlap sky area = $10\,000\,\deg^{2}$, overlap redshift range
$0.9\leq z\leq 1.7$ (bin size $\Delta z=0.1$).
$\bullet$ Tracer $t$:

high density, $\mathcal{O}(10^{-3})\;(\text{$h$/Mpc})^{3}$;
$\overline{n}_{t}(z)$ similar to a Euclid sample [76];

linear bias $b_{1}^{t}(z)=1.3+0.6\,z$.
$\bullet$ Tracer $t^{\prime}$:

lower density, $\mathcal{O}(10^{-4})\;(\text{$h$/Mpc})^{3}$;
$\overline{n}_{t^{\prime}}(z)$ similar to DESI sample [77];

linear bias $b_{1}^{t^{\prime}}(z)=0.84/D(z)$.
In order to derive the higher-order and non-Gaussian bias coefficients needed
for the tree-level bispectrum, we use a simplified halo-occupation
distribution (HOD) model from [78]. The HOD model has two free parameters
which are calibrated from the values of the linear bias and number densities
at each redshift.

This set-up for the nominal surveys exploits the findings of Section 5 and
will function as the benchmark case. The fiducial cosmology is a flat
$\Lambda$CDM model, with $\omega_{c}=0.1199$, $\omega_{b}=0.02237$,
$h=0.6736$, $n_{s}=0.96488$, $A_{s}=2.1\times 10^{-9}$, as measured by [79].
In addition, we assume non-Gaussian initial conditions, with an amplitude of
$f_{\rm NL}=20$, for all three PNG types considered here.

### $\sigma(f_{\rm NL})(z)$: benchmark case

Figure 1 presents the conditional errors on $f_{\rm NL}$ as a function of
redshift in the benchmark case, for the two nominal tracers as well as their
multi-tracer, using power spectrum and bispectrum. For both correlators, the
multi-tracer delivers an overall improvement relative to each single-tracer
survey.

In the case of local PNG, the multi-tracer power spectrum improvement is
significant across the redshift range, while for the bispectrum, this applies
only in the low-redshift bins. However, the multi-tracer precision on
$f_{\mathrm{NL}}^{\mathrm{loc}}$ is effectively the same for $P^{\rm MT}$ and
$B^{\rm MT}$.

For orthogonal PNG, the improvement from the multi-tracer is significantly
more for $P^{\rm MT}$ than for $B^{\rm MT}$. Note that the performances of the
single tracers $t$ and $t^{\prime}$ are opposite to the local case. Equation
5.2 shows that the single-tracer conditional error from the power spectrum on
$f_{\rm NL}$ is parametrically determined by $b_{1},b_{\Psi}$ and
$P_{\varepsilon}$. Numerical tests and the analytic results show the
constraints are more sensitive to the values of $b_{\Psi}$ and in particular
for the case where $b_{\Psi}^{I}\gg b_{1}^{I}$. For the local case
$b_{\Psi}^{t^{\prime}}>b_{\Psi}^{t}$ and hence we see that tracer $t$
outperforms the constraints from $t^{\prime}$, for both correlators. On the
other hand, for orthogonal PNG, due to the nature of this shape (see [44] for
a discussion) $b_{\Psi}^{t}>b_{\Psi}^{t^{\prime}}$, leading to an opposite
behaviour with respect to the local PNG.

For equilateral PNG, the multi-tracer bispectrum gives only a marginal
improvement over the best tracer performance, as in the orthogonal case, while
the power spectrum does not provide any constraints on
$f_{\mathrm{NL}}^{\mathrm{equil}}$ (see e.g. [65]).

The attributes of the nominal surveys were selected to produce a large
$|b_{1}^{J}\,b_{\Psi}^{I}-b_{1}^{I}\,b_{\Psi}^{J}|$. The results in Fig. 1
indicate that in the general benchmark case – i.e. with no additional
assumptions – this is not enough to achieve an optimal gain from the multi-
tracer approach. Apart from the dependence on PNG type shown in Fig. 1,
additional parameters have a significant effect on the level of improvement
achieved by the multi-tracer. The main tracer-dependent variables that affects
$\sigma(f_{\rm NL})$, besides the linear and PNG biases ($b_{1}^{I}$ and
$b_{\Psi}^{I}$), are the number densities $\overline{n}_{I}$, which in turn
affect the stochastic contributions in Eq. 2.10. In addition, for the
bispectrum, higher-order bias parameters, e.g. $b_{2}^{I}$, will have an
effect on the constraining power.

Further examination is needed on the parametric behaviour of the multi-tracer
power spectrum and bispectrum, in order to determine the tracer properties
that maximise gain from combining their information. Here we confine attention
to local PNG, and we fix the redshift at $z=1$. Note that $b_{\Psi}\rightarrow
b_{\Phi}$ for local PNG, as follows from Eq. 2.7 with $\alpha=0$.

### $\sigma(f_{\mathrm{NL}}^{\mathrm{loc}})$: effect of number densities

We start with the effect of number densities on
$f_{\mathrm{NL}}^{\mathrm{loc}}$ forecasts. We vary $\overline{n}_{t}$, while
fixing all other parameters of tracers $t$ and $t^{\prime}$, including
$\overline{n}_{t^{\prime}}$. The forecast errors on
$f_{\mathrm{NL}}^{\mathrm{loc}}$ are shown in the left panel of Fig. 2 as a
function of the ratio $\overline{n}_{t}/\overline{n}_{t^{\prime}}$, at $z=1$.
In addition to the multi-tracer results, we present the forecasts from simply
adding the Fisher matrices of the two tracers, without considering the cross-
terms in the data vectors, Eqs. 3.9, LABEL: and 3.11, and block covariances,
Eqs. 4.4, LABEL: and 4.5.

The left panel of Fig. 2 shows that for both correlators the multi-tracer
gives maximal precision on $f_{\mathrm{NL}}^{\mathrm{loc}}$ for
$\overline{n}_{t}/\overline{n}_{t^{\prime}}\gtrsim 10$. In other words, the
number density of one tracer should be at least an order of magnitude larger
than that of the other. In fact, if the ratio is beyond that point, the
constraints do not improve further, as we remain shot-noise dominated by the
low-density tracer. This behaviour occurs since only one of the two tracers is
allowed to reach the limit where stochastic contributions can be neglected
(i.e. the cosmic variance limit). The other tracer has a constant shot-noise
component, that ultimately saturates the multi-tracer Fisher information. On
the other hand, the constraints from just summing the Fisher matrices of the
two tracers (i.e. neglecting the multi-tracer cross-terms), exhibit a minimal
change as a function of the number density ratio. This indicates that the
complete multi-tracer approach takes advantage of the increasing number
density in an optimal way, by maximising the contribution of the cross-terms
and hence providing a significant improvement on the
$f_{\mathrm{NL}}^{\mathrm{loc}}$ constraints. In this setup, the power
spectrum seems to be more sensitive to the cosmic variance limit of tracer
$t$, than the bispectrum. Note that the scenario where the number densities of
both tracers are equal (i.e. $\overline{n}_{t}/\overline{n}_{t^{\prime}}=1$),
is not an optimal setup for a multi-tracer application, at least when neither
of the tracers approaches the high-density limit.

Figure 2:  Left: $1\sigma$ conditional error on
$f_{\mathrm{NL}}^{\mathrm{loc}}$ as a function of the ratio
$\overline{n}_{t}/\overline{n}_{t^{\prime}}$, at fixed redshift ($z=1$). We
vary only the number density of sample $t$, while fixing all other parameters
of tracer $t$, together with those of sample $t^{\prime}$. Right: Testing the
cosmic variance limit, by enforcing equal number densities, i.e.
$\overline{n}_{t}=\overline{n}_{t^{\prime}}=\overline{n}$, and varying
$\overline{n}$, at $z=1$. Both panels show the multi-tracer forecasts in
comparison to the simple combination of the two Fisher matrices, without
considering the cross-terms in the data vectors and block covariances. + and
$\times$ indicate the benchmark results of Fig. 1 at $z=1$, for the power
spectrum and bispectrum respectively.

In order to investigate further the limit where the shot-noise can be
neglected, we assume equal number densities. Our aim is to explore the low and
high density limits. The forecast error on $f_{\mathrm{NL}}^{\mathrm{loc}}$ as
a function of number density is presented in the right panel of Fig. 2, for
both correlators. This indicates that the multi-tracer constraints on local
PNG exhibit an inverse power-law behaviour as the number densities of the
tracers approach the cosmic variance limit, $\overline{n}\gg
10^{-3}\;(\text{$h$/Mpc})^{3}$. By contrast, the results of simply adding the
Fisher matrices of the two tracers saturate around the same value, without
exhibiting any further improvement as the number densities increase. This
shows the clear advantage of the multi-tracer approach, where the cross-term
contribution is enhanced significantly, due to the decreasing shot-noise
terms, providing competitive constraints on $f_{\mathrm{NL}}^{\mathrm{loc}}$.

Our power spectrum results clearly show the benefit of the multi-tracer, for
the high density samples, as well as for tracers that approach the cosmic
variance limit. This agrees with the findings of the original work on the
multi-tracer approach [25], as well as the more recent results of [40].
Moreover, here we show that the advantage of the multi-tracer towards the high
density limit can also be extended to the bispectrum, which exhibits an order
of magnitude improvement over the power spectrum in the cosmic variance
limited regime.

Figure 3:  $1\sigma$ conditional error on $f_{\mathrm{NL}}^{\mathrm{loc}}$ as
a function of linear bias difference $b_{1}^{t}-b_{1}^{t^{\prime}}$, at fixed
redshift $z=1$, for the multi-tracer power spectrum and bispectrum. Left:
Varying only $b_{1}^{t}$, while fixing all other tracer-dependent parameters
for both $t$ and $t^{\prime}$. Right: Varying $b_{1}^{t}$, while allowing all
the other bias parameters of tracer $t$ to change in accordance with the
dependence of the HOD model on the linear bias. + and $\times$ indicate the
benchmark results of Fig. 1.

### $\sigma(f_{\mathrm{NL}}^{\mathrm{loc}})$: effect of linear biases

Next we investigate the effect of the linear bias on the multi-tracer power
spectrum and bispectrum. We vary the linear bias of tracer $t$, while fixing
all parameters that are dependent on tracer $t^{\prime}$. Due to the
symmetrisation of the two correlators under tracer permutations, the outcome
of the parametric analysis does not depend on which tracer’s parameters we
chose to vary. Since the bias coefficients of each tracer depend on the HOD
model, whose free parameters are calibrated by the value of the linear bias at
each redshift, varying $b_{1}^{t}$ will affect all the bias coefficients of
tracer $t$. Hence the approach is two-fold: (1) varying only $b_{1}^{t}$,
while fixing all other parameters for both tracers, and (2) allowing the
remaining bias parameters of tracer $t$ to vary as well, through the
dependence of the HOD on $b_{1}^{t}$. The first corresponds to the case where
we treat $b^{1}_{\Phi}$ as being independent of $b_{1}^{t}$. The second
corresponds to the case where we consider the universality relation, i.e.
$b_{\Phi}^{t}$ is determined by $b_{1}^{t}$. The forecast multi-tracer
constraints on $f_{\mathrm{NL}}^{\mathrm{loc}}$, for the two correlators and
at $z=1$, are presented in Fig. 3 as a function of the difference between the
linear biases of the two tracers, $b_{1}^{t}-b_{1}^{t^{\prime}}$.

The left panel displays local PNG constraints after varying only $b_{1}^{t}$
and fixing all other tracer-dependent parameters for both tracers. The results
indicate that when the difference between the two biases $\gtrsim 2$, there is
a significant improvement ($\sim 50\%$) for the multi-tracer power spectrum.
For the multi-tracer bispectrum, the gain is more modest. Beyond that point
the improvement saturates for both correlators. The results show that the
constraining power on $f_{\mathrm{NL}}^{\mathrm{loc}}$ is $\propto
b_{1}^{t}-b_{1}^{t^{\prime}}$. Note that the case where the linear biases of
the tracers are similar, $b_{1}^{t}\sim b_{1}^{t^{\prime}}$, corresponds to
the least tight constraints on $f_{\mathrm{NL}}^{\mathrm{loc}}$, for both
correlators.

Figure 4:  Same as the right panel of Fig. 3, but here
$\sigma(f_{\mathrm{NL}}^{\mathrm{loc}})$ is shown as a function of
$b_{1}^{t}\,b_{\Phi}^{t^{\prime}}-b_{1}^{t^{\prime}}\,b_{\Phi}^{t}$.

The behaviour changes once we allow all other bias coefficients of tracer $t$
to vary as well, through the dependence of the HOD free parameters on the
$b_{1}^{t}$ value. The results are shown in the right panel of Fig. 3. The
case where $b_{1}^{t}\sim b_{1}^{t^{\prime}}$ no longer corresponds to poor
constraints for the mulit-tracer power spectrum, but to the tightest ones. On
the other hand, for the multi-tracer bispectrum, the strongest constraints
correspond to values $b_{1}^{t}-b_{1}^{t^{\prime}}>3$. Beyond this point, the
constraints saturate.

### $\sigma(f_{\mathrm{NL}}^{\mathrm{loc}})$: effect of
$b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}$

Another way of interpreting these results, is to present the forecasts on
$\sigma(f_{\mathrm{NL}}^{\mathrm{loc}})$ as a function of
$b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}$. This is
displayed in Fig. 4, where we allow all parameters of tracer $t$ to vary, via
the dependence of the HOD on $b_{1}^{t}$ (i.e. $b_{\Phi}^{t}(b_{1}^{t})$). The
behaviour exhibited indicates that the constraining power on
$f_{\mathrm{NL}}^{\mathrm{loc}}$ increases as a function of
$b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}$. This is
true for both correlators, where the power spectrum gain is more significant
than the bispectrum as a function of
$b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}$, verifying
the findings of Section 5.

Furthermore, the results of Fig. 4 show that the
$f_{\mathrm{NL}}^{\mathrm{loc}}$ errors from the multi-tracer power spectrum
are minimised for large negative values of
$b_{1}^{t}\,b_{\Phi}^{t^{\prime}}-b_{1}^{t^{\prime}}\,b_{\Phi}^{t}$. This
regime corresponds to $b_{1}^{t}-b_{1}^{t^{\prime}}\sim 0$, as shown in the
right panel of Fig. 3, which explains why $b_{1}^{t}\sim b_{1}^{t^{\prime}}$
provides such tight constraints, as discussed before. For the power spectrum,
both plots verify the findings of Section 5. For the multi-tracer bispectrum
the tightest constraint on $f_{\mathrm{NL}}^{\mathrm{loc}}$ is provided once
$b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}<0$,
indicating the different behaviour of the two correlators as a function of
$b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}$, once we
allow for all bias parameters of tracer $t$ to change and $b_{\Phi}^{t}\propto
b_{1}^{t}$. Moreover, the multi-tracer power spectrum has much less pronounced
improvement for
$b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}\approx 0$,
while for the bispectrum the constraints are close to the optimal results.

Figure 5: 1$\sigma$ conditional error on $f_{\mathrm{NL}}^{\mathrm{loc}}$,
with values given by the colour bar, in the $(b_{\Phi}^{t},\,b_{1}^{t})$
plane, for multi-tracer power spectrum (_left_) and bispectrum (_right_). We
vary $b_{1}^{t}$ and $b_{\Phi}^{t}$ independently, without assuming a
universality relation, while fixing all other parameters for both tracers.
Black dotted lines are contours of constant
$b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}$.

The different parametric behaviour of the multi-tracer power spectrum and
bispectrum can also been seen in the contour plots of Fig. 5. In this case we
vary $b_{1}^{t}$ and $b_{\Phi}^{t}$, without assuming any relation between the
two (i.e. no universality relation), while fixing all other tracer-dependent
parameters for both tracers. The black dashed contour lines indicate the lines
of constant
$b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}$. The
results show that indeed for the multi-tracer power spectrum the case
$b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}\approx 0$
correspond to the least pronounced constraints. The tightest forecasts
originate from the case of maximal
$b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}$ values and
in particular for $-10\leq b_{\Phi}\leq 10$ and for $b_{1}\geq 2$ values. The
behavior of the multi-tracer bispectrum is somewhat similar, only now for the
$b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}\approx 0$
regime, which corresponds as well to the looser constraints, the performance
is better than the power spectrum equivalent. As shown as well in Fig. 4, the
optimal constraints from the multi-tracer bispectrum originate for an
increasing
$|b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}|$,
although the gain is more modest than the power spectrum case. In particular,
to get a significant improvement, growing $b_{1}^{t}$ values and very negative
$b_{\Phi}^{t}$ are required (blue shaded region in the right panel of Fig. 5).

Figure 6:  $1\sigma$ conditional error on $f_{\mathrm{NL}}^{\mathrm{loc}}$ as
a function of the difference $b_{2}^{t}-b_{2}^{t^{\prime}}$, varying
$b_{2}^{t}$, while keeping fixed all other parameters for both tracers. Square
points correspond to the summation of the two Fishers, while dots represent
the full multi-tracer results. $\times$ corresponds to the multi-tracer
bispectrum benchmark results of Fig. 1.

The results presented in Figs. 3, LABEL:, 4, LABEL: and 5 indicate the
potential of the multi-tracer to significantly improve constraints on local
PNG, by carefully choosing the two tracers. For the multi-tracer power
spectrum and bispectrum, increasing $|b_{1}^{t}-b_{1}^{t^{\prime}}|$ or
$|b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}|$ gives
the optimal results, where the latter leads to overall better constraints. The
power spectrum exhibits a much larger gain compared to the bispectrum. This
indicates that increasing
$|b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}|$ is not
enough to take advantage of the full potential of the multi-tracer bispectrum.

### $\sigma(f_{\mathrm{NL}}^{\mathrm{loc}})$: effect of second-order biases

In order to investigate further the parametric behaviour of the multi-tracer
bispectrum, we vary the second-order bias coefficient $b_{2}^{t}$, while all
other parameters for both tracers are fixed. The results are shown in Fig. 6
as a function of the difference $b_{2}^{t}-b_{2}^{t^{\prime}}$ for the multi-
tracer bispectrum, and for the sum of the individual Fisher matrices.
Comparing the two shows once again the importance of utilising the full multi-
tracer in order to optimise constraints on $f_{\mathrm{NL}}^{\mathrm{loc}}$
from the bispectrum information of two tracers.

In addition, varying $b_{2}^{t}$ shows a peak around the case where the
nonlinear bias parameters of tracers $t$ and $t^{\prime}$ are similar, i.e.
the difference $b_{2}^{t}-b_{2}^{t^{\prime}}$ is small. This indicates the
regime where the multi-tracer bispectrum offers the least tight constraints.
Therefore, the constraining power on $f_{\mathrm{NL}}^{\mathrm{loc}}$ is
$\propto|b_{2}^{t}-b_{2}^{t^{\prime}}|$. In fact, increasing the value of
$|b_{2}^{t}-b_{2}^{t^{\prime}}|$ improves the constraints on
$f_{\mathrm{NL}}^{\mathrm{loc}}$ significantly more than in the case of an
increasing $|b_{1}^{t}-b_{1}^{t^{\prime}}|$ or
$|b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}|$ (see
Figs. 3, LABEL: and 4). This implies the importance of the nonlinear bias
parameter in utilising the full potential of the multi-tracer bispectrum.

##  6 Summary and conclusion

In this work we investigate the potential of the multi-tracer power spectrum
and multi-tracer bispectrum analysis in constraining the amplitude of PNG,
with a particular focus on local PNG. We consider two nominal next-generation
surveys and use the Fisher matrix to probe the information on
$f_{\mathrm{NL}}^{\mathrm{loc}}$, within the power spectra and bispectra of
the two tracers. The main goal is to quantify the impact of the various
tracer-dependent parameters on the constraining power on
$f_{\mathrm{NL}}^{\mathrm{loc}}$, in order to maximise the gain from the
multi-tracer technique and provide the tracer-selection criteria for future
galaxy surveys. The main findings can be summarised as follows:

  * •

The complete multi-tracer approach takes advantage of increasing number
density in an optimal way, by maximising the contribution of the cross-terms.
This is particularly important for high-density samples, as well as for
tracers that approach the cosmic variance limit. The scenario where the number
densities of the tracers are similar does not provide the tightest
constraints, except for high-density samples.

  * •

The constraining power of the power spectrum and bispectrum multi-tracers is
$\propto|b_{1}^{t}\,b_{\Phi}^{t^{\prime}}-b_{1}^{t^{\prime}}\,b_{\Phi}^{t}|$,
where the power spectrum exhibits a much larger gain compared to the
bispectrum. In particular the optimal regime is for an increasing difference
between the linear biases of the two tracers, as well as for a negative
$b_{\Phi}$ for one of the two tracers.

  * •

Maximising the quantity
$|b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}|$ is not
enough for to fully exploit the potential of the multi-tracer bispectrum.
Increasing as well the second-order bias difference
$|b_{2}^{t}-b_{2}^{t^{\prime}}|$ provides a significant gain in the bispectrum
constraints, compared to only increasing
$|b_{1}^{t^{\prime}}\,b_{\Phi}^{t}-b_{1}^{t}\,b_{\Phi}^{t^{\prime}}|$.

The results presented here do not consider correlations between the various
parameters, within the Fisher matrix formalism, but address only the case of
the conditional $f_{\mathrm{NL}}^{\mathrm{loc}}$ error. Correlations are
important since the power of the multi-tracer approach can be applied to
reducing or breaking the various degeneracies between parameters. In order to
fully investigate the parametric behaviour of the multi-tracer technique we
need to perform a complete Fisher analysis, including the various parameters
of the model. This is the subject of the follow-up paper of this work.

## Acknowledgments

We thank Daisuke Yamauchi for helpful clarifications of their work [30]. DK
and RM are supported by the South African Radio Astronomy Observatory and the
National Research Foundation (Grant No. 75415). JF is supported by Fundação
para a Ciência e a Tecnologia (FCT) through the research grants
UIDB/04434/2020 and UIDP/04434/2020 and through the Investigador FCT Contract
No. 2020.02633.CEECIND/CP1631/CT0002. SC acknowledges support from the
‘Departments of Excellence 2018-2022’ Grant (L. 232/2016) awarded by the
Italian Ministry of University and Research (mur) and from the ‘Ministero
degli Affari Esteri della Cooperazione Internazionale (maeci) – Direzione
Generale per la Promozione del Sistema Paese Progetto di Grande Rilevanza
ZA18GR02. CC is supported by the UK Science & Technology Facilities Council
Consolidated Grant ST/P000592/1. JF, SC and CC thank the University of the
Western Cape for supporting a visit to Cape Town, during which this work was
initiated. This work made use of the South African Centre for High Performance
Computing, under the project Cosmology with Radio Telescopes, ASTRO-0945.

##  Appendix A Redshift-space kernels and damping factors

The general non-Gaussian redshift kernels up to second order for tracer $I$,
neglecting $\mathcal{O}(f_{\rm NL}^{2})$ terms, are given by (see e.g. [80,
81, 75]):

$\displaystyle Z_{1}^{I}(\bm{k}_{i})$
$\displaystyle=b_{1}^{I}+f\,\mu_{i}^{2}+f_{\rm
NL}\,\frac{b^{I}_{\Phi}\,k_{i}^{\alpha}}{M(k_{i})}\;,$ (A.1) $\displaystyle
Z_{2}^{I}(\bm{k}_{i},\bm{k}_{j})$
$\displaystyle=b_{1}^{I}\,F_{2}(\bm{k}_{i},\bm{k}_{j})+f\,\mu_{ij}^{2}\,G_{2}(\bm{k}_{i},\bm{k}_{j})+\frac{b_{2}^{I}}{2}+\frac{b_{s}^{I}}{2}\,S_{2}(\bm{k}_{i},\bm{k}_{j})$
$\displaystyle~{}~{}~{}~{}+\frac{f\,\mu_{ij}\,k_{ij}}{2}\left[\frac{\mu_{i}}{k_{i}}\,Z_{1}^{I}(\bm{k}_{j})+\frac{\mu_{j}}{k_{j}}\,Z_{1}^{I}(\bm{k}_{i})\right]+f_{\rm
NL}\,\frac{b^{I}_{\Psi\delta}}{2}\,\left[\frac{k_{i}^{\alpha}}{M(k_{i})}+\frac{k_{j}^{\alpha}}{M(k_{j})}\right]$
$\displaystyle~{}~{}~{}~{}-f_{\rm
NL}\,b^{I}_{\Psi}\,N_{2}(\bm{k}_{i},\bm{k}_{j})\;.$ (A.2)

Here $f(z)$ is the linear growth rate,
$\mu_{i}=\hat{\bm{k}}_{i}\cdot\hat{\bm{z}}$, with $\hat{\bm{z}}$ the line-of-
sight direction, $\mu_{ij}=(\mu_{i}k_{i}+\mu_{j}k_{j})/k_{ij}$, and
$k_{ij}^{2}=(\bm{k}_{i}+\bm{k}_{j})^{2}$. The kernels
$F_{2}(\bm{k}_{i},\bm{k}_{j})$ and $G_{2}(\bm{k}_{i},\bm{k}_{j})$ are the
second-order symmetric SPT kernels (e.g. [47]), while
$S_{2}(\bm{k}_{1},\bm{k}_{2})=(\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{2})^{2}-1/3$
is the tidal kernel [82, 54]. The symmetrised $N_{2}$ kernel is

$N_{2}(\bm{k}_{i},\bm{k}_{j})=\frac{1}{2}\,\left[\frac{k_{i}^{\alpha}}{k_{j}^{2}\,M(k_{i})}+\frac{k_{j}^{\alpha}}{k_{i}^{2}\,M(k_{j})}\right]\,\bm{k}_{i}\cdot\bm{k}_{j}\;,$
(A.3)

which encodes the coupling of the PNG potential to the Eulerian-to-Lagrangian
displacement field [72, 80]. The parameter $\alpha=2,1$ for equilateral,
orthogonal shapes [61, 73, 44, 74] and $\alpha=0$ for local shapes [55, 56,
72].

The non-perturbative RSD fingers-of-god effect, which models the reduction in
clustering power due to nonlinear velocities, is described by a
phenomenological model, following [83, 84]. Here we use an exponential damping
factor, which in Fourier space is a radial convolution,
$\delta_{\bm{k}}\rightarrow\delta_{\bm{k}}\exp(-k^{2}\mu^{2}\sigma/2)$.
Furthermore, in the case of an optical survey, redshift errors smear the
galaxy density field along the line-of-sight, which in Fourier space can be
described by
$\delta_{\bm{k}}\rightarrow\delta_{\bm{k}}\exp(-k^{2}\mu^{2}\sigma_{r}/2)$
[85]. The redshift errors propagate into comoving distance errors, where
$\sigma_{r}(z)=c/H(z)\sigma_{z}(z)$ and $\sigma_{z}(z)$ is the redshift error
of the survey. Following the most general approach, we assume that the damping
amplitude is tracer-dependent, so that

$\displaystyle{\cal D}_{P}^{IJ}(\bm{k})$
$\displaystyle=\exp\Big{[}-\frac{1}{2}\,k^{2}\,\mu^{2}\,\Big{\\{}\big{(}\sigma_{P}^{I}\big{)}^{2}+\big{(}\sigma_{r}^{I}\big{)}^{2}+\big{(}\sigma_{P}^{J}\big{)}^{2}+\big{(}\sigma_{r}^{J}\big{)}^{2}\Big{\\}}\Big{]}\;,$
(A.4) $\displaystyle{\cal D}_{B}^{IJK}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle=\exp\Big{[}-\frac{1}{2}\,\Big{\\{}k_{1}^{2}\,\mu_{1}^{2}\,\big{[}\big{(}\sigma_{B}^{I}\big{)}^{2}+\big{(}\sigma_{r}^{I}\big{)}^{2}\big{]}+k_{2}^{2}\,\mu_{2}^{2}\,\big{[}\big{(}\sigma_{B}^{J}\big{)}^{2}+\big{(}\sigma_{r}^{J}\big{)}^{2}\big{]}$
(A.5)
$\displaystyle~{}~{}~{}+k_{3}^{2}\,\mu_{3}^{2}\,\big{[}\big{(}\sigma_{B}^{K}\big{)}^{2}+\big{(}\sigma_{r}^{K}\big{)}^{2}\big{]}\Big{\\}}\Big{]}\;,$

where the redshift dependence has been suppressed for brevity. The amplitude
parameters for RSD and redshift errors, for power spectra and bispectra, are
independent free parameters. The fiducial values for RSD are given by the
velocity dispersion, i.e. $\sigma_{P}^{I}(z)=\sigma_{B}^{I}(z)=\sigma_{v}(z)$.

##  Appendix B 3-tracer case

For three tracers, $I,J,K$ take values $t,t^{\prime},t^{\prime\prime}$, and
the power spectrum data vector is

$\bm{D_{P}}=\Big{[}P^{tt},P^{tt^{\prime}},P^{t^{\prime}t^{\prime}},P^{t^{\prime}t^{\prime\prime}},P^{t^{\prime\prime}t^{\prime\prime}},P^{tt^{\prime\prime}}\Big{]}\;,$
(B.1)

where the expected values of the entries can be derived as in the case of two
tracers, by changing accordingly the tracer-dependent quantities in Eq. 3.10.
Likewise for the bispectrum data vector in the case of three tracers:

$\bm{D_{B}}=\Big{[}B^{ttt},B^{(ttt^{\prime})},B^{(tt^{\prime}t^{\prime})},B^{t^{\prime}t^{\prime}t^{\prime}},B^{(t^{\prime}t^{\prime}t^{\prime\prime})},B^{(t^{\prime}t^{\prime\prime}t^{\prime\prime})},B^{t^{\prime\prime}t^{\prime\prime}t^{\prime\prime}},B^{(ttt^{\prime\prime})},B^{(tt^{\prime\prime}t^{\prime\prime})},B^{(tt^{\prime}t^{\prime\prime})}\Big{]}\;,$
(B.2)

where the expected values for the first 9 entries can be derived from Eqs.
3.12 to 3.15 by exchanging accordingly the tracer indexes in the tracer-
dependent quantities. The expected value of the last entry can be defined as

$\displaystyle\langle\hat{B}^{(tt^{\prime}t^{\prime\prime})}\rangle$
$\displaystyle\equiv
B^{(tt^{\prime}t^{\prime\prime})}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=\frac{1}{6}\,\Big{[}\hat{B}^{tt^{\prime}t^{\prime\prime}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})+\hat{B}^{t^{\prime\prime}t^{\prime}t}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})+\hat{B}^{t^{\prime}t^{\prime\prime}t}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle+\hat{B}^{tt^{\prime\prime}t^{\prime}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})+\hat{B}^{t^{\prime}tt^{\prime\prime}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})+\hat{B}^{t^{\prime\prime}tt^{\prime}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})\Big{]}\;,$
(B.3)

where

$\displaystyle\langle\hat{B}^{tt^{\prime}t^{\prime\prime}}\rangle\equiv\hat{B}^{tt^{\prime}t^{\prime\prime}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle={\cal
D}_{B}^{tt^{\prime}t^{\prime\prime}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})\,\bigg{[}Z_{1}^{t}(\bm{k}_{1})\,Z_{1}^{t^{\prime}}(\bm{k}_{2})\,Z_{1}^{t^{\prime\prime}}(\bm{k}_{3})\,B_{\rm
PNG}(k_{1},k_{2},k_{3})$
$\displaystyle+2\,Z_{1}^{t}(\bm{k}_{1})\,Z_{1}^{t^{\prime}}(\bm{k}_{2})\,Z_{2}^{t^{\prime\prime}}(\bm{k}_{1},\bm{k}_{2})\,P(k_{1})\,P(k_{2})$
$\displaystyle+2\,Z_{1}^{t^{\prime}}(\bm{k}_{2})\,Z_{1}^{t^{\prime\prime}}(\bm{k}_{3})\,Z_{2}^{t}(\bm{k}_{2},\bm{k}_{3})\,P(k_{2})\,P(k_{3})$
$\displaystyle+2\,Z_{1}^{t}(\bm{k}_{1})\,Z_{1}^{t^{\prime\prime}}(\bm{k}_{3})\,Z_{2}^{t^{\prime}}(\bm{k}_{3},\bm{k}_{1})\,P(k_{3})\,P(k_{1})\bigg{]}\;,$
(B.4)

and similarly for the other terms, by appropriately permutating the tracer
indexes. Note that the symmetrised bispectrum of Appendix B has no stochastic
terms due to the assumption
$P_{\varepsilon\varepsilon_{\delta}}^{IJ}=B^{IJK}_{\varepsilon}=0$ .

The data vector of the combined power spectra and bispectra is defined by
combining the data vectors of Eqs. B.1, LABEL: and B.2, creating a 16 element
array. The covariance of this vector is again a 2x2 block matrix, as in the
two-tracer case of Eq. 4.2, but now each individual block will have different
dimensions. To be more precise:

  * •

$\mathrm{Cov}(\bm{P}_{i},\bm{P}_{j})$ is a 6$\times$6 matrix,

  * •

$\mathrm{Cov}(\bm{B}_{i},\bm{B}_{b})$ is a 10$\times$10 matrix,

  * •

$\mathrm{Cov}(\bm{P}_{i},\bm{B}_{b})$ is a 6$\times$10 matrix.

As in the case of two tracers, we consider only the Gaussian part of the
covariances, i.e. each sub-block is a diagonal matrix. The expression for each
sub-block in the $\mathrm{Cov}(\bm{P}_{i},\bm{P}_{j})$ block can be easily
derived from Eq. 4.4, by permuting the tracer indexes accordingly. The same
folows for most of the sub-blocks of $\mathrm{Cov}(\bm{B}_{i},\bm{B}_{b})$,
where all the expressions for the diagonal sub-blocks can be derived from
Section 4 by interchanging the tracer indices accordingly. The only sub-blocks
left to define are those that correspond to the correlations with the
symmetrised bispectrum $B^{(tt^{\prime}t^{\prime\prime})}$:

$\displaystyle\mathrm{Cov}\big{(}B_{a}^{ttt},B_{b}^{(tt^{\prime}t^{\prime\prime})}\big{)}$
$\displaystyle=\frac{1}{2}\Big{(}P^{tt}_{a_{1}}P^{tt^{\prime}}_{a_{2}}P^{tt^{\prime\prime}}_{a_{3}}+P^{tt}_{a_{1}}P^{tt^{\prime\prime}}_{a_{2}}P^{tt^{\prime}}_{a_{3}}\Big{)}+\text{2\,perm}\;,$
$\displaystyle\mathrm{Cov}\big{(}B_{a}^{(ttt^{\prime})},B_{b}^{(tt^{\prime}t^{\prime\prime})}\big{)}$
$\displaystyle=\frac{1}{6}\Big{[}P^{tt}_{a_{1}}\Big{(}P^{tt^{\prime\prime}}_{a_{2}}P^{t^{\prime}t^{\prime}}_{a_{3}}+P^{t^{\prime}t^{\prime}}_{a_{2}}P^{tt^{\prime\prime}}_{a_{3}}\Big{)}+P^{tt}_{a_{1}}\Big{(}P^{tt^{\prime}}_{a_{2}}P^{t^{\prime}t^{\prime\prime}}_{a_{3}}+P^{t^{\prime}t^{\prime\prime}}_{a_{2}}P^{tt^{\prime}}_{a_{3}}\Big{)}$
$\displaystyle~{}~{}~{}+2P^{tt^{\prime}}_{a_{1}}P^{tt^{\prime}}_{a_{2}}P^{tt^{\prime\prime}}_{a_{3}}\Big{]}+\text{2\,perm}\;,$
$\displaystyle\mathrm{Cov}\big{(}B_{a}^{(tt^{\prime}t^{\prime\prime})},B_{b}^{(tt^{\prime}t^{\prime\prime})}\big{)}$
$\displaystyle=\frac{1}{6}\Big{[}P^{tt^{\prime}}_{a_{1}}\Big{(}P^{tt^{\prime\prime}}_{a_{2}}P^{t^{\prime}t^{\prime\prime}}_{a_{3}}+P^{t^{\prime}t^{\prime\prime}}_{a_{2}}P^{tt^{\prime\prime}}_{a_{3}}\Big{)}+\frac{1}{2}P^{tt}_{a_{1}}\Big{(}P^{t^{\prime}t^{\prime}}_{a_{2}}P^{t^{\prime\prime}t^{\prime\prime}}_{a_{3}}+P^{t^{\prime\prime}t^{\prime\prime}}_{a_{2}}P^{t^{\prime}t^{\prime}}_{a_{3}}\Big{)}$
$\displaystyle~{}~{}~{}+P^{tt^{\prime\prime}}_{a_{1}}P^{tt^{\prime\prime}}_{a_{2}}P^{t^{\prime}t^{\prime}}_{a_{3}}+P^{tt}_{a_{1}}P^{t^{\prime}t^{\prime\prime}}_{a_{2}}P^{t^{\prime}t^{\prime\prime}}_{a_{3}}+P^{tt^{\prime}}_{a_{1}}P^{tt^{\prime}}_{a_{2}}P^{t^{\prime\prime}t^{\prime\prime}}_{a_{3}}\Big{]}+\text{2\,perm}.$
(B.5)

The remaining sub-block covariances can be retrieved by changing accordingly
the tracer indices in the above expressions.

The additional expressions needed to derive the power spectra-bispectra cross-
covariance in the case of three tracers are:

$\displaystyle\mathrm{Cov}\big{(}P_{i}^{tt},B_{b}^{(tt^{\prime}t^{\prime\prime})}\big{)}$
$\displaystyle=\frac{1}{3}\Big{[}P^{tt^{\prime\prime}}(\bm{k}_{i})B_{2\delta}^{ttt^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})+P^{tt^{\prime}}(\bm{k}_{i})B_{2\delta}^{ttt^{\prime\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})$
$\displaystyle~{}~{}~{}+P^{tt}(\bm{k}_{i})B_{1\delta}^{tt^{\prime}t^{\prime\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\Big{]},$
$\displaystyle\mathrm{Cov}\big{(}P_{i}^{tt^{\prime}},B_{b}^{(tt^{\prime}t^{\prime\prime})}\big{)}$
$\displaystyle=\frac{1}{6}\Big{[}P^{tt^{\prime}}(\bm{k}_{i})B_{2\delta}^{tt^{\prime}t^{\prime\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})+P^{tt}(\bm{k}_{i})B_{1\delta}^{t^{\prime}t^{\prime}t^{\prime\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})$
$\displaystyle~{}~{}~{}+P^{tt^{\prime\prime}}(\bm{k}_{i})B_{1\delta}^{tt^{\prime}t^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})+P^{t^{\prime}t^{\prime}}(\bm{k}_{i})B_{1\delta}^{ttt^{\prime\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})$
$\displaystyle~{}~{}~{}+P^{t^{\prime}t^{\prime\prime}}(\bm{k}_{i})B_{1\delta}^{ttt^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\Big{]},$
(B.6)

where

$\displaystyle
B_{1\delta}^{tt^{\prime}t^{\prime\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})$
$\displaystyle=\delta_{ib_{1}}\Big{[}B^{tt^{\prime}t^{\prime\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})+B^{tt^{\prime\prime}t^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\Big{]}$
$\displaystyle~{}~{}~{}+\delta_{ib_{2}}\Big{[}B^{t^{\prime}tt^{\prime\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})+B^{t^{\prime\prime}tt^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\Big{]}$
$\displaystyle~{}~{}~{}+\delta_{ib_{3}}\Big{[}B^{t^{\prime\prime}t^{\prime}t}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})+B^{t^{\prime}t^{\prime\prime}t}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\Big{]},$
(B.7) $\displaystyle
B_{2\delta}^{tt^{\prime}t^{\prime\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})$
$\displaystyle=\big{(}\delta_{ib_{1}}+\delta_{ib_{2}}\big{)}\Big{[}B^{tt^{\prime}t^{\prime\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})+B^{t^{\prime}tt^{\prime\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\Big{]}$
$\displaystyle~{}~{}~{}+\big{(}\delta_{ib_{1}}+\delta_{ib_{3}}\big{)}\Big{[}B^{tt^{\prime\prime}t^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})+B^{t^{\prime}t^{\prime\prime}t}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\Big{]}$
$\displaystyle~{}~{}~{}+\big{(}\delta_{ib_{2}}+\delta_{ib_{3}}\big{)}\Big{[}B^{t^{\prime\prime}t^{\prime}t}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})+B^{t^{\prime\prime}tt^{\prime}}(\bm{k}_{b_{1}},\bm{k}_{b_{2}},\bm{k}_{b_{3}})\Big{]}.$
(B.8)

## References

  * [1] S. Novell-Masot, H. Gil-Marín, and L. Verde, GEO-FPT: a model of the galaxy bispectrum at mildly non-linear scales, arXiv:2303.15510.
  * [2] N. S. Sugiyama, D. Yamauchi, T. Kobayashi, T. Fujita, S. Arai, S. Hirano, S. Saito, F. Beutler, and H.-J. Seo, New constraints on cosmological modified gravity theories from anisotropic three-point correlation functions of BOSS DR12 galaxies, arXiv:2302.06808.
  * [3] M. M. Ivanov, O. H. E. Philcox, G. Cabass, T. Nishimichi, M. Simonović, and M. Zaldarriaga, Cosmology with the Galaxy Bispectrum Multipoles: Optimal Estimation and Application to BOSS Data, arXiv:2302.04414.
  * [4] D. Karagiannis, R. Maartens, and L. F. Randrianjanahary, Cosmological constraints from the power spectrum and bispectrum of 21cm intensity maps, JCAP 11 (2022) 003, [arXiv:2206.07747].
  * [5] 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 (2022), no. 4 043517, [arXiv:2112.04515].
  * [6] M. M. Ivanov, O. H. E. Philcox, T. Nishimichi, M. Simonović, M. Takada, and M. Zaldarriaga, Precision analysis of the redshift-space galaxy bispectrum, Phys. Rev. D 105 (2022), no. 6 063512, [arXiv:2110.10161].
  * [7] A. Moradinezhad Dizgah, M. Biagetti, E. Sefusatti, V. Desjacques, and J. Noreña, Primordial Non-Gaussianity from Biased Tracers: Likelihood Analysis of Real-Space Power Spectrum and Bispectrum, JCAP 05 (2021) 015, [arXiv:2010.14523].
  * [8] N. S. Sugiyama, S. Saito, F. Beutler, and H.-J. Seo, A complete FFT-based decomposition formalism for the redshift-space bispectrum, Mon. Not. Roy. Astron. Soc. 484 (2019), no. 1 364–384, [arXiv:1803.02132].
  * [9] C. Clarkson, E. M. de Weerd, S. Jolicoeur, R. Maartens, and O. Umeh, The dipole of the galaxy bispectrum, Mon. Not. Roy. Astron. Soc. 486 (2019), no. 1 L101–L104, [arXiv:1812.09512].
  * [10] R. Maartens, S. Jolicoeur, O. Umeh, E. M. De Weerd, C. Clarkson, and S. Camera, Detecting the relativistic galaxy bispectrum, JCAP 03 (2020) 065, [arXiv:1911.02398].
  * [11] S. Jolicoeur, R. Maartens, E. M. De Weerd, O. Umeh, C. Clarkson, and S. Camera, Detecting the relativistic bispectrum in 21cm intensity maps, JCAP 06 (2021) 039, [arXiv:2009.06197].
  * [12] O. Umeh, S. Jolicoeur, R. Maartens, and C. Clarkson, A general relativistic signature in the galaxy bispectrum: the local effects of observing on the lightcone, JCAP 03 (2017) 034, [arXiv:1610.03351].
  * [13] E. Di Dio, H. Perrier, R. Durrer, G. Marozzi, A. Moradinezhad Dizgah, J. Noreña, and A. Riotto, Non-Gaussianities due to Relativistic Corrections to the Observed Galaxy Bispectrum, JCAP 03 (2017) 006, [arXiv:1611.03720].
  * [14] R. Maartens, S. Jolicoeur, O. Umeh, E. M. De Weerd, and C. Clarkson, Local primordial non-Gaussianity in the relativistic galaxy bispectrum, JCAP 04 (2021) 013, [arXiv:2011.13660].
  * [15] M. Noorikuhani and R. Scoccimarro, Wide-angle and relativistic effects in Fourier-space clustering statistics, Phys. Rev. D 107 (2023), no. 8 083528, [arXiv:2207.12383].
  * [16] E. Di Dio, R. Durrer, R. Maartens, F. Montanari, and O. Umeh, The Full-Sky Angular Bispectrum in Redshift Space, JCAP 04 (2019) 053, [arXiv:1812.09297].
  * [17] K. Garcia and Z. Slepian, Improving the line of sight for the anisotropic 3-point correlation function of galaxies: Centroid and Unit-Vector-Average methods scaling as $\mathscr{O}$ (N2), Mon. Not. Roy. Astron. Soc. 515 (2022), no. 1 1199–1217, [arXiv:2011.03503].
  * [18] K. Pardede, E. Di Dio, and E. Castorina, Wide-angle effects in the galaxy bispectrum, arXiv:2302.12789.
  * [19] D. Karagiannis, A. Slosar, and M. Liguori, Forecasts on Primordial non-Gaussianity from 21 cm Intensity Mapping experiments, JCAP 11 (2020) 052, [arXiv:1911.03964].
  * [20] S. Bharadwaj, A. Mazumdar, and D. Sarkar, Quantifying the Redshift Space Distortion of the Bispectrum I: Primordial Non-Gaussianity, Mon. Not. Roy. Astron. Soc. 493 (2020), no. 1 594–602, [arXiv:2001.10243].
  * [21] D. Karagiannis, J. Fonseca, R. Maartens, and S. Camera, Probing primordial non-Gaussianity with the power spectrum and bispectrum of future 21 cm intensity maps, Phys. Dark Univ. 32 (2021) 100821, [arXiv:2010.07034].
  * [22] G. D’Amico, M. Lewandowski, L. Senatore, and P. Zhang, Limits on primordial non-Gaussianities from BOSS galaxy-clustering data, arXiv:2201.11518.
  * [23] W. R. Coulton, F. Villaescusa-Navarro, D. Jamieson, M. Baldi, G. Jung, D. Karagiannis, M. Liguori, L. Verde, and B. D. Wandelt, Quijote-PNG: The Information Content of the Halo Power Spectrum and Bispectrum, Astrophys. J. 943 (2023), no. 2 178, [arXiv:2206.15450].
  * [24] G. Jung, D. Karagiannis, M. Liguori, M. Baldi, W. R. Coulton, D. Jamieson, L. Verde, F. Villaescusa-Navarro, and B. D. Wandelt, Quijote-PNG: Quasi-maximum likelihood estimation of Primordial Non-Gaussianity in the non-linear halo density field, arXiv:2211.07565.
  * [25] U. Seljak, Extracting primordial non-gaussianity without cosmic variance, Phys. Rev. Lett. 102 (2009) 021302, [arXiv:0807.1770].
  * [26] P. McDonald and U. Seljak, How to measure redshift-space distortions without sample variance, JCAP 10 (2009) 007, [arXiv:0810.0323].
  * [27] N. Hamaus, U. Seljak, and V. Desjacques, Optimal Constraints on Local Primordial Non-Gaussianity from the Two-Point Statistics of Large-Scale Structure, Phys. Rev. D 84 (2011) 083509, [arXiv:1104.2321].
  * [28] L. R. Abramo and K. E. Leonard, Why multi-tracer surveys beat cosmic variance, Mon. Not. Roy. Astron. Soc. 432 (2013) 318, [arXiv:1302.5444].
  * [29] L. D. Ferramacho, M. G. Santos, M. J. Jarvis, and S. Camera, Radio galaxy populations and the multitracer technique: pushing the limits on primordial non-Gaussianity, Mon. Not. Roy. Astron. Soc. 442 (2014), no. 3 2511–2518, [arXiv:1402.2290].
  * [30] D. Yamauchi, K. Takahashi, and M. Oguri, Constraining primordial non-Gaussianity via a multitracer technique with surveys by Euclid and the Square Kilometre Array, Phys. Rev. D 90 (2014), no. 8 083520, [arXiv:1407.5453].
  * [31] D. Alonso and P. G. Ferreira, Constraining ultralarge-scale cosmology with multiple tracers in optical and radio surveys, Phys. Rev. D 92 (2015), no. 6 063525, [arXiv:1507.03550].
  * [32] J. Fonseca, S. Camera, M. Santos, and R. Maartens, Hunting down horizon-scale effects with multi-wavelength surveys, Astrophys. J. Lett. 812 (2015), no. 2 L22, [arXiv:1507.04605].
  * [33] J. Fonseca, R. Maartens, and M. G. Santos, Probing the primordial Universe with MeerKAT and DES, Mon. Not. Roy. Astron. Soc. 466 (2017), no. 3 2780–2786, [arXiv:1611.01322].
  * [34] J. Gleyzes, R. de Putter, D. Green, and O. Doré, Biasing and the search for primordial non-Gaussianity beyond the local type, JCAP 04 (2017) 002, [arXiv:1612.06366].
  * [35] M. Ballardini, W. L. Matthewson, and R. Maartens, Constraining primordial non-Gaussianity using two galaxy surveys and CMB lensing, Mon. Not. Roy. Astron. Soc. 489 (2019), no. 2 1950–1956, [arXiv:1906.04730].
  * [36] Z. Gomes, S. Camera, M. J. Jarvis, C. Hale, and J. Fonseca, Non-Gaussianity constraints using future radio continuum surveys and the multitracer technique, Mon. Not. Roy. Astron. Soc. 492 (2020), no. 1 1513–1522, [arXiv:1912.08362].
  * [37] J.-A. Viljoen, J. Fonseca, and R. Maartens, Multi-wavelength spectroscopic probes: prospects for primordial non-Gaussianity and relativistic effects, JCAP 11 (2021) 010, [arXiv:2107.14057].
  * [38] A. Barreira, Predictions for local PNG bias in the galaxy power spectrum and bispectrum and the consequences for $f_{NL}$ constraints, JCAP 01 (2022), no. 01 033, [arXiv:2107.06887].
  * [39] S. Jolicoeur, R. Maartens, and S. Dlamini, Constraining primordial non-Gaussianity by combining next-generation galaxy and 21 cm intensity mapping surveys, Eur. Phys. J. C 83 (2023), no. 4 320, [arXiv:2301.02406].
  * [40] A. Barreira and E. Krause, Towards optimal and robust $f_{\rm NL}$ constraints with multi-tracer analyses, arXiv:2302.09066.
  * [41] J. M. Sullivan, T. Prijon, and U. Seljak, Learning to Concentrate: Multi-tracer Forecasts on Local Primordial Non-Gaussianity with Machine-Learned Bias, arXiv:2303.08901.
  * [42] D. Yamauchi, S. Yokoyama, and K. Takahashi, Multitracer technique for galaxy bispectrum: An application to constraints on nonlocal primordial non-Gaussianities, Phys. Rev. D 95 (March, 2017) 063530, [arXiv:1611.03590].
  * [43] S. Camera, M. G. Santos, and R. Maartens, Probing primordial non-Gaussianity with SKA galaxy redshift surveys: a fully relativistic analysis, Mon. Not. Roy. Astron. Soc. 448 (2015), no. 2 1035–1043, [arXiv:1409.8286]. [Erratum: Mon.Not.Roy.Astron.Soc. 467, 1505–1506 (2017)].
  * [44] V. Desjacques, D. Jeong, and F. Schmidt, Large-scale galaxy bias, Phys. Rep. 733 (Feb, 2018) 1–193, [arXiv:1611.09787].
  * [45] A. Lewis, A. Challinor, and A. Lasenby, Efficient Computation of Cosmic Microwave Background Anisotropies in Closed Friedmann-Robertson-Walker Models, ApJ 538 (August, 2000) 473–476, [astro-ph/9911177].
  * [46] E. Komatsu, N. Afshordi, N. Bartolo, D. Baumann, J. R. Bond, E. I. Buchbinder, C. T. Byrnes, X. Chen, D. J. H. Chung, A. Cooray, P. Creminelli, N. Dalal, O. Dore, R. Easther, A. V. Frolov, J. Khoury, W. H. Kinney, L. Kofman, K. Koyama, L. Leblond, J.-L. Lehners, J. E. Lidsey, M. Liguori, E. A. Lim, A. Linde, D. H. Lyth, J. Maldacena, S. Matarrese, L. McAllister, P. McDonald, S. Mukohyama, B. Ovrut, H. V. Peiris, A. Riotto, Y. Rodrigues, M. Sasaki, R. Scoccimarro, D. Seery, A. Sefusatti, K. M. Smith, A. A. Starobinsky, P. J. Steinhardt, F. Takahashi, M. Tegmark, A. J. Tolley, L. Verde, B. D. Wandelt, D. Wands, S. Weinberg, M. Wyman, A. P. S. Yadav, and M. Zaldarriaga, Non-Gaussianity as a Probe of the Physics of the Primordial Universe and the Astrophysics of the Low Redshift Universe, astro2010: The Astronomy and Astrophysics Decadal Survey 2010 (2009) 158, [arXiv:0902.4759].
  * [47] F. Bernardeau, S. Colombi, E. Gaztañaga, and R. Scoccimarro, Large-scale structure of the Universe and cosmological perturbation theory, Phys. Rep. 367 (September, 2002) 1–248, [astro-ph/0112551].
  * [48] V. Assassi, D. Baumann, D. Green, and M. Zaldarriaga, Renormalized halo bias, J. Cosmology Astropart. Phys 8 (August, 2014) 056, [arXiv:1402.5916].
  * [49] L. Senatore, Bias in the effective field theory of large scale structures, J. Cosmology Astropart. Phys 11 (November, 2015) 007, [arXiv:1406.7843].
  * [50] M. Mirbabayi, F. Schmidt, and M. Zaldarriaga, Biased tracers and time evolution, J. Cosmology Astropart. Phys 7 (July, 2015) 030, [arXiv:1412.5169].
  * [51] A. Dekel and O. Lahav, Stochastic Nonlinear Galaxy Biasing, ApJ 520 (July, 1999) 24–34, [astro-ph/9806193].
  * [52] A. Taruya and J. Soda, Stochastic Biasing and the Galaxy-Mass Density Relation in the Weakly Nonlinear Regime, ApJ 522 (September, 1999) 46–58, [astro-ph/9809204].
  * [53] T. Matsubara, Stochasticity of Bias and Nonlocality of Galaxy Formation: Linear Scales, ApJ 525 (November, 1999) 543–553, [astro-ph/9906029].
  * [54] T. Baldauf, U. Seljak, V. Desjacques, and P. McDonald, Evidence for quadratic tidal tensor bias from the halo bispectrum, Phys. Rev. D 86 (October, 2012) 083540, [arXiv:1201.4827].
  * [55] N. Dalal, O. Doré, D. Huterer, and A. Shirokov, Imprints of primordial non-Gaussianities on large-scale structure: Scale-dependent bias and abundance of virialized objects, Phys. Rev. D 77 (June, 2008) 123514, [arXiv:0710.4560].
  * [56] A. Slosar, C. Hirata, U. Seljak, S. Ho, and N. Padmanabhan, Constraints on local primordial non-Gaussianity from large scale structure, J. Cosmology Astropart. Phys 2008 (Aug, 2008) 031, [arXiv:0805.3580].
  * [57] S. Matarrese and L. Verde, The Effect of Primordial Non-Gaussianity on Halo Bias, ApJ 677 (April, 2008) L77, [arXiv:0801.4826].
  * [58] L. Verde and S. Matarrese, Detectability of the Effect of Inflationary Non-Gaussianity on Halo Bias, ApJ 706 (November, 2009) L91–L95, [arXiv:0909.3224].
  * [59] N. Afshordi and A. J. Tolley, Primordial non-Gaussianity, statistics of collapsed objects, and the integrated Sachs-Wolfe effect, Phys. Rev. D 78 (December, 2008) 123507, [arXiv:0806.1046].
  * [60] V. Desjacques and U. Seljak, Primordial Non-Gaussianity in the Large-Scale Structure of the Universe, Advances in Astronomy 2010 (2010) 908640, [arXiv:1006.4763].
  * [61] F. Schmidt and M. Kamionkowski, Halo clustering with nonlocal non-Gaussianity, Phys. Rev. D 82 (November, 2010) 103002, [arXiv:1008.0638].
  * [62] R. Scoccimarro, L. Hui, M. Manera, and K. C. Chan, Large-scale bias and efficient generation of initial conditions for nonlocal primordial non-Gaussianity, Phys. Rev. D 85 (April, 2012) 083002, [arXiv:1108.5512].
  * [63] V. Desjacques, D. Jeong, and F. Schmidt, Accurate predictions for the scale-dependent galaxy bias from primordial non-Gaussianity, Phys. Rev. D 84 (September, 2011) 061301, [arXiv:1105.3476].
  * [64] F. Schmidt, D. Jeong, and V. Desjacques, Peak-background split, renormalization, and galaxy clustering, Phys. Rev. D 88 (July, 2013) 023515, [arXiv:1212.0868].
  * [65] V. Assassi, D. Baumann, and F. Schmidt, Galaxy bias and primordial non-Gaussianity, J. Cosmology Astropart. Phys 12 (December, 2015) 043, [arXiv:1510.03723].
  * [66] W. L. W. Sargent and E. L. Turner, A statistical method for determining the cosmological density parameter from the redshifts of a complete sample of galaxies, Astrophys. J. 212 (1977) L3–L7.
  * [67] N. Kaiser, Clustering in real space and in redshift space, MNRAS 227 (July, 1987) 1–21.
  * [68] A. J. S. Hamilton, Linear Redshift Distortions: a Review, in The Evolving Universe (D. Hamilton, ed.), vol. 231 of Astrophysics and Space Science Library, p. 185, 1998.  astro-ph/9708102.
  * [69] F. Schmidt, Towards a self-consistent halo model for the nonlinear large-scale structure, Phys. Rev. D 93 (March, 2016) 063512, [arXiv:1511.02231].
  * [70] E. Sefusatti, M. Crocce, S. Pueblas, and R. Scoccimarro, Cosmology and the bispectrum, Phys. Rev. D 74 (July, 2006) 023522, [astro-ph/0604505].
  * [71] M. Biagetti, L. Castiblanco, J. Noreña, and E. Sefusatti, The covariance of squeezed bispectrum configurations, J. Cosmology Astropart. Phys 2022 (Sept., 2022) 009, [arXiv:2111.05887].
  * [72] T. Giannantonio and C. Porciani, Structure formation from non-Gaussian initial conditions: Multivariate biasing, statistics, and comparison with N-body simulations, Phys. Rev. D 81 (March, 2010) 063530, [arXiv:0911.0017].
  * [73] T. Giannantonio, C. Porciani, J. Carron, A. Amara, and A. Pillepich, Constraining primordial non-Gaussianity with future galaxy surveys, MNRAS 422 (June, 2012) 2854–2877, [arXiv:1109.0958].
  * [74] G. Cabass, E. Pajer, and F. Schmidt, Imprints of Oscillatory Bispectra on Galaxy Clustering, JCAP 09 (2018) 003, [arXiv:1804.07295].
  * [75] D. Karagiannis, A. Lazanu, M. Liguori, A. Raccanelli, N. Bartolo, and L. Verde, Constraining primordial non-Gaussianity with bispectrum and power spectrum from upcoming optical and radio surveys, MNRAS 478 (July, 2018) 1341–1376, [arXiv:1801.09280].
  * [76] Euclid Collaboration, A. Blanchard et al., Euclid preparation: VII. Forecast validation for Euclid cosmological probes, Astron. Astrophys. 642 (2020) A191, [arXiv:1910.09273].
  * [77] DESI Collaboration, A. Aghamousa et al., The DESI Experiment Part I: Science,Targeting, and Survey Design, arXiv:1611.00036.
  * [78] V. Yankelevich and C. Porciani, Cosmological information in the redshift-space bispectrum, Mon. Not. Roy. Astron. Soc. 483 (2019), no. 2 2078–2099, [arXiv:1807.07076].
  * [79] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, S. Basak, R. Battye, K. Benabed, J. P. Bernard, M. Bersanelli, P. Bielewicz, J. J. Bock, J. R. Bond, J. Borrill, F. R. Bouchet, F. Boulanger, M. Bucher, C. Burigana, R. C. Butler, E. Calabrese, J. F. Cardoso, J. Carron, A. Challinor, H. C. Chiang, J. Chluba, L. P. L. Colombo, C. Combet, D. Contreras, B. P. Crill, F. Cuttaia, P. de Bernardis, G. de Zotti, J. Delabrouille, J. M. Delouis, E. Di Valentino, J. M. Diego, O. Doré, M. Douspis, A. Ducout, X. Dupac, S. Dusini, G. Efstathiou, F. Elsner, T. A. Enßlin, H. K. Eriksen, Y. Fantaye, M. Farhang, J. Fergusson, R. Fernandez-Cobos, F. Finelli, F. Forastieri, M. Frailis, A. A. Fraisse, E. Franceschi, A. Frolov, S. Galeotta, S. Galli, K. Ganga, R. T. Génova-Santos, M. Gerbino, T. Ghosh, J. González-Nuevo, K. M. Górski, S. Gratton, A. Gruppuso, J. E. Gudmundsson, J. Hamann, W. Handley, F. K. Hansen, D. Herranz, S. R. Hildebrandt, E. Hivon, Z. Huang, A. H. Jaffe, W. C. Jones, A. Karakci, E. Keihänen, R. Keskitalo, K. Kiiveri, J. Kim, T. S. Kisner, L. Knox, N. Krachmalnicoff, M. Kunz, H. Kurki-Suonio, G. Lagache, J. M. Lamarre, A. Lasenby, M. Lattanzi, C. R. Lawrence, M. Le Jeune, P. Lemos, J. Lesgourgues, F. Levrier, A. Lewis, M. Liguori, P. B. Lilje, M. Lilley, V. Lindholm, M. López-Caniego, P. M. Lubin, Y. Z. Ma, J. F. Macías-Pérez, G. Maggio, D. Maino, N. Mandolesi, A. Mangilli, A. Marcos-Caballero, M. Maris, P. G. Martin, M. Martinelli, E. Martínez-González, S. Matarrese, N. Mauri, J. D. McEwen, P. R. Meinhold, A. Melchiorri, A. Mennella, M. Migliaccio, M. Millea, S. Mitra, M. A. Miville-Deschênes, D. Molinari, L. Montier, G. Morgante, A. Moss, P. Natoli, H. U. Nørgaard-Nielsen, L. Pagano, D. Paoletti, B. Partridge, G. Patanchon, H. V. Peiris, F. Perrotta, V. Pettorino, F. Piacentini, L. Polastri, G. Polenta, J. L. Puget, J. P. Rachen, M. Reinecke, M. Remazeilles, A. Renzi, G. Rocha, C. Rosset, G. Roudier, J. A. Rubiño-Martín, B. Ruiz-Granados, L. Salvati, M. Sandri, M. Savelainen, D. Scott, E. P. S. Shellard, C. Sirignano, G. Sirri, L. D. Spencer, R. Sunyaev, A. S. Suur-Uski, J. A. Tauber, D. Tavagnacco, M. Tenti, L. Toffolatti, M. Tomasi, T. Trombetti, L. Valenziano, J. Valiviita, B. Van Tent, L. Vibert, P. Vielva, F. Villa, N. Vittorio, B. D. Wand elt, I. K. Wehus, M. White, S. D. M. White, A. Zacchei, and A. Zonca, Planck 2018 results. VI. Cosmological parameters, A&A 641 (Sept., 2020) A6, [arXiv:1807.06209].
  * [80] T. Baldauf, U. Seljak, and L. Senatore, Primordial non-Gaussianity in the bispectrum of the halo density field, J. Cosmology Astropart. Phys 4 (April, 2011) 006, [arXiv:1011.1513].
  * [81] M. Tellarini, A. J. Ross, G. Tasinato, and D. Wands, Galaxy bispectrum, primordial non-Gaussianity and redshift space distortions, J. Cosmology Astropart. Phys 6 (June, 2016) 014, [arXiv:1603.06814].
  * [82] P. McDonald and A. Roy, Clustering of dark matter tracers: generalizing bias for the coming era of precision LSS, J. Cosmology Astropart. Phys 8 (August, 2009) 020, [arXiv:0902.0991].
  * [83] J. A. Peacock and S. J. Dodds, Reconstructing the Linear Power Spectrum of Cosmological Mass Fluctuations, MNRAS 267 (April, 1994) 1020, [astro-ph/9311057].
  * [84] W. E. Ballinger, J. A. Peacock, and A. F. Heavens, Measuring the cosmological constant with redshift surveys, MNRAS 282 (October, 1996) 877, [astro-ph/9605017].
  * [85] H.-J. Seo and D. J. Eisenstein, Probing Dark Energy with Baryonic Acoustic Oscillations from Future Large Galaxy Redshift Surveys, ApJ 598 (December, 2003) 720–740, [astro-ph/0307460].