Measuring $f_{\mathrm{NL}}$ with the SPHEREx Multi-tracer Redshift Space Bispectrum

We consider here the local type primordial non-Gaussianity parameterized by $f_{\mathrm{NL}}$

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

(1)

where $\varphi$ is an auxiliary primordial Gaussian potential. Using the Poisson equation, we can relate the primordial potential to the linearly evolved primordial matter density perturbation $\delta_{\rm m,p}$ as (valid on subhorizon scales in the Newtonian limit)

$$\Phi(\bm{k})=\frac{\delta_{\text{m,p}}(\bm{k},z)}{\alpha(k,z)},$$

(2)

where

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

(3)

and where $D(z)$ is the linear growth factor normalized at $z=0$, $T(k)$ the transfer function, $\Omega_{m}$ the matter density and $H_{0}$ the Hubble constant. On large scales where $T(k)=1$, $\alpha$ scales as $k^{2}/H^{2}$.

Working in perturbation theory, we expand the linear matter density contrast $\delta_{\text{m}}(\bm{k})$ as

$$\delta_{\text{m}}(\bm{k})=\delta_{\text{m}}^{(1)}(\bm{k})+\delta_{\text{m}}^{(2)}(\bm{k})+\delta_{\text{m}}^{(3)}(\bm{k})+\ldots\;.$$

(4)

Given Eq. 2 the linearly evolved primordial matter density field up to second order is then given by

$\displaystyle\delta_{\text{m,p}}(\bm{k},z)$

$\displaystyle=$

$\displaystyle\alpha(k,z)\,\Phi(\bm{k})=\delta_{\text{m,p}}^{(1)}(\bm{k},z)+f_{\mathrm{NL}}\,\delta_{\text{m,p}}^{(2)}(\bm{k},z),$

where

$$\delta_{\text{m,p}}^{(1)}(\bm{k},z)=\alpha(k,z)\varphi(\bm{k}),$$

(6)

and

$$\delta_{\text{m,p}}^{(2)}(\bm{k},z)=\alpha(k,z)\,\int\frac{\,\mathrm{d}^{3}q}{(2\pi)^{3}}\varphi(\bm{q})\varphi(\bm{k}-\bm{q}).$$

(7)

Gravitational evolution also contributes to nonlinear coupling starting at second order, so that the non-linearly evolved matter density field receives an additional term proportional to $F_{2}$:

	$\displaystyle\delta^{(1)}_{\text{m}}(\bm{k},z)$	$\displaystyle=$	$\displaystyle\delta_{\text{m,p}}^{(1)}(\bm{k},z),$		(8)
	$\displaystyle\delta^{(2)}_{\text{m}}(\bm{k},z)$	$\displaystyle=$	$\displaystyle\int\frac{\,\mathrm{d}^{3}q}{(2\pi)^{3}}\delta^{(1)}_{\text{m,p}}(\bm{q})\delta^{(1)}_{\text{m,p}}(\bm{k}-\bm{q})F_{2}(\bm{q},\bm{k}-\bm{q})$		(9)
		$\displaystyle+$	$\displaystyle f_{\mathrm{NL}}\,\delta_{\text{m,p}}^{(2)}(\bm{k},z),$

where $F_{2}$ is the second-order mode coupling kernel

$$F_{2}(\bm{k}_{1},\bm{k}_{2})=\frac{5}{7}+\frac{1}{2}\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}k_{2}}\left(\frac{k_{1}}{k_{2}}+\frac{k_{2}}{k_{1}}\right)+\frac{2}{7}\frac{(\bm{k}_{1}\cdot\bm{k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}.$$

(10)

Galaxy surveys observe the density of galaxies, which are biased tracers of the underlying dark matter density. We follow Ref. Tellarini et al. (2016) in our modeling for the galaxy density and the galaxy bispectrum in the section to follow. We consider the following bivariate bias model for the Eulerian galaxy density constrast, where the dependence on $\varphi$ arises in the presence of non-Gaussianities:

	$\displaystyle\delta_{\mathrm{g}}(\bm{x})$	$\displaystyle=$	$\displaystyle b_{10}\delta_{\rm m}(\bm{x})+b_{01}\varphi(\bm{x})$		(11)
		$\displaystyle+$	$\displaystyle\frac{1}{2}b_{20}\left(\delta_{\rm m}(\bm{x})\right)^{2}+b_{11}\delta_{\rm m}(\bm{x})\varphi(\bm{x})$
		$\displaystyle+$	$\displaystyle\frac{1}{2}b_{02}\left(\varphi(\bm{x})\right)^{2}+\frac{1}{2}b_{s_{2}}(s^{2}-\langle s^{2}\rangle)-b_{01}n^{2}.$

Here we have the tidal term Catelan et al. (2000); Baldauf et al. (2012)

$$s^{2}(\bm{k})=\int\frac{d\bm{q}}{(2\pi)^{3}}\mathcal{S}_{2}(\bm{q},\bm{k}-\bm{q})\delta_{\rm m}^{(1)}(\bm{q})\delta_{\rm m}^{(1)}(\bm{k}-\bm{q}),$$

(12)

and the non-Gaussian shift term due to the displacement of galaxies with respect to their initial positions $\bm{q}$ in Lagrangian coordinates

$$n^{2}(\bm{k})=2\int\frac{d\bm{q}}{(2\pi)^{3}}\mathcal{N}_{2}(\bm{q},\bm{k}-\bm{q})\frac{\delta_{\rm m}^{(1)}(\bm{q})\delta_{\rm m}^{(1)}(\bm{k}-\bm{q})}{\alpha(|\bm{k}-\bm{q}|)},$$

(13)

where

$$\mathcal{S}(\bm{k}_{1},\bm{k}_{2})=\frac{(\bm{k}_{1}\cdot\bm{k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}-\frac{1}{3},$$

(14)

and

$$\mathcal{N}(\bm{k}_{1},\bm{k}_{2})=\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{2k_{1}^{2}}.$$

(15)

We note that we have ignored the stochastic contributions to $\delta_{\mathrm{g}}$ in Eq. 11, and also do not marginalize over any potential deviations of the shot noise from Poisson predictions (see Ref. Rizzo et al. (2023) for example for the full stochastic contributions to the tree-level bispectrum).

Taking only the first order terms in Eq. 11 we have

$$\delta_{\mathrm{g}}(\bm{k})=\left(b_{10}+\frac{b_{01}}{\alpha(k)}\right)\delta_{\text{m,p}}(\bm{k}).$$

(16)

We will expound on the specific modeling and values we take for each bias parameter in Section IV, where we will assume an universal mass function, for which the relationship between $b_{10}$ and $b_{01}$ becomes

$$b_{01}=2f_{\mathrm{NL}}\delta_{c}(b_{10}-1),$$

(17)

and we recover the well-known linear order result Dalal et al. (2008)

$$\delta_{\mathrm{g}}(\bm{k})=\left(b_{10}+\frac{2f_{\mathrm{NL}}\delta_{\text{c}}(b_{10}-1)}{\alpha(k)}\right)\delta_{\text{m,p}}(\bm{k}),$$

(18)

where $\delta_{c}=1.686$ is the threshold for spherical collapse.

We note that recent studies have shown how the universal mass function assumption for obtaining the relation Eq. 17 for $b_{01}$ (also known as $b_{\phi}f_{\mathrm{NL}}$) could be inaccurate and could bias constraints on $f_{\mathrm{NL}}$ (e.g. Barreira (2020)). So in a realistic analysis, one could marginalize over theoretically-informed priors for $b_{\phi}$ or choose to constrain the combination $b_{\phi}f_{\mathrm{NL}}$ instead Barreira (2020, 2022). Multi-tracer analysis can also help to improve the $f_{\mathrm{NL}}$ constraints for suitably chosen galaxy samples (e.g. by maximizing the combination $|b_{10}^{B}b_{01}^{A}-b_{10}^{A}b_{01}^{B}|$ in a two-tracer analysis) Barreira and Krause (2023).

Let us define the multi-tracer bispectrum $B_{\rm ggg}^{ABC}$ as

$$\langle\delta_{g}^{A}(\bm{k}_{1})\delta_{g}^{B}(\bm{k}_{2})\delta_{g}^{C}(\bm{k}_{3})\rangle=(2\pi)^{3}\,\delta_{D}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3})B_{ggg}^{ABC}(\bm{k}_{1},\bm{k}_{2}),$$

(19)

where $A,B,C$ denote the galaxy samples. Now the wavevectors $\bm{k}_{1},\bm{k}_{2},\bm{k}_{3}$ will be associated with samples $A,B,C$ respectively. We only use the galaxy bispectrum with samples from the same redshift bin $i$ centered at $z_{i}$. Note that $A,B,C=1..n_{b}$, so that there are $n_{b}=n_{\rm tracer}^{3}$ multitracer combinations for each redshift bin.

To model the bispectrum in redshift space, we include the linear Kaiser effects on large scales, and the damping of small scales due to redshift errors, ignoring for now the Alcock-Pazynski effects.

The redshift errors $\sigma_{z,i}^{A}=\tilde{\sigma}_{z}^{A}(1+z_{i})$ decreases our ability to measure modes parallel to the line-of-sight. We choose to model this effect with a Gaussian suppression using

$$F^{ABC}_{i}(\bm{k}_{1},\bm{k}_{2})=e^{-\frac{1}{2}\left[k_{1}^{2}\mu_{1}^{2}(\sigma_{p,i}^{A})^{2}+k_{2}^{2}\mu_{2}^{2}(\sigma_{p,i}^{B})^{2}+k_{3}^{2}\mu_{3}^{2}(\sigma_{p,i}^{C})^{2}\right]},$$

(20)

and

$$\sigma_{p,i}^{A}(z)=\tilde{\sigma}_{z}^{A}(1+z_{i})\frac{2\pi c}{H(z)}.$$

(21)

Note that here $\sigma_{p,i}(z)$ can no longer be factored out for the multi-tracer bispectrum for which each galaxy sample may have a different redshift error.

The multi-tracer bispectrum in redshift space is then modeled as

$$B^{ABC}_{ggg}(\bm{k}_{1},\bm{k}_{2}|z_{i})=F^{ABC}_{i}(\bm{k}_{1},\bm{k}_{2})\,\left[2Z_{1}^{A}(\bm{k}_{1})Z_{1}^{B}(\bm{k}_{2})Z_{2}^{C}(\bm{k}_{1},\bm{k}_{2})\,P(k_{1})P(k_{2})+2\mathrm{\ cycl.\ perm.}\right],$$

(22)

where

	$\displaystyle Z_{1}^{X}(\bm{k}_{1})=\;$	$\displaystyle b_{10}^{X}$	$\displaystyle(1+\beta\mu_{1}^{2})+\frac{b_{01}^{X}}{\alpha(k_{1})}$		(23)
	$\displaystyle Z_{2}^{X}(\bm{k}_{1},\bm{k}_{2})=\;$	$\displaystyle b_{10}^{X}$	$\displaystyle\left[F_{2}(\bm{k}_{1},\bm{k}_{2})+f_{\mathrm{NL}}\frac{\alpha(k)}{\alpha(k_{1})\alpha(k_{2})}\right]+\frac{b_{20}^{X}}{2}+\frac{1}{2}b_{s_{2}}^{X}S_{2}(\bm{k}_{1},\bm{k}_{2})$		(24)
		$\displaystyle+$	$\displaystyle\frac{b_{11}^{X}}{2}\left[\frac{1}{\alpha(k_{1})}+\frac{1}{\alpha(k_{2})}\right]+\frac{b_{02}^{X}}{2}\frac{1}{\alpha(k_{1})\alpha(k_{2})}-b_{01}^{X}\left[\frac{N_{2}(\bm{k}_{1},\bm{k}_{2})}{\alpha(k_{2})}+\frac{N_{2}(\bm{k}_{2},\bm{k}_{1})}{\alpha(k_{1})}\right]$
		$\displaystyle+$	$\displaystyle f\mu^{2}\left[G_{2}(\bm{k}_{1},\bm{k}_{2})+f_{\mathrm{NL}}\frac{\alpha(k)}{\alpha(k_{1})\alpha(k_{2})}\right]+\frac{1}{2}f\mu k\left(\frac{\mu_{1}}{k_{1}}Z_{1}^{X}(\bm{k}_{2})+\frac{\mu_{2}}{k_{2}}Z_{1}^{X}(\bm{k}_{1})\right),$

where $\mu_{i}=\hat{\bm{k}}_{i}\cdot\hat{\bm{n}}$, where $\hat{\bm{n}}$ is the line-of-sight, $k=k_{3}$ and $\mu=-\mu_{3}$, and $G_{2}$ is the second-order velocity kernel

$$G_{2}(\bm{k}_{1},\bm{k}_{2})=\frac{3}{7}+\frac{1}{2}\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}k_{2}}\left(\frac{k_{1}}{k_{2}}+\frac{k_{2}}{k_{1}}\right)+\frac{4}{7}\frac{(\bm{k}_{1}\cdot\bm{k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}.$$

(25)

The last term in Eq. 24 is often also written as

$\displaystyle\frac{1}{2}f\mu k\left(\frac{\mu_{1}}{k_{1}}Z_{1}^{X}(\bm{k}_{2})+\frac{\mu_{2}}{k_{2}}Z_{1}^{X}(\bm{k}_{1})\right)=f^{2}\mu^{2}k^{2}\frac{\mu_{1}\mu_{2}}{2k_{1}k_{2}}+b_{10}^{X}\frac{f\mu k}{2}\left(\frac{\mu_{1}}{k_{1}}+\frac{\mu_{2}}{k_{2}}\right)+b_{01}^{X}\frac{f\mu k}{2}\left[\frac{\mu_{1}}{k_{1}\alpha(k_{2})}+\frac{\mu_{2}}{k_{2}\alpha(k_{1})}\right].$

(26)

Let $\tilde{B}^{ABC}(k_{1},k_{2},k_{3})$ represent the binned bispectrum over triangle shapes and orientations with bin centers denoted by $(k_{1},k_{2},k_{3})$ and $(\theta,\phi)$ respectively. If we approximate the binned bispectrum by the value at the bin center, then the Fisher matrix for the multi-tracer bispectrum in a single redshift bin can be written as

$\displaystyle F_{ij}$

$\displaystyle=$

$\displaystyle\sum_{(k_{1},k_{2},k_{3})}\sum_{(\theta,\phi)}\sum_{(ABC)}\sum_{(A^{\prime}B^{\prime}C^{\prime})}\frac{\partial\tilde{B}^{ABC}(k_{1},k_{2},k_{3},\theta,\phi)}{\partial p_{i}}\left[\mathrm{\tilde{Cov}}\right]^{-1}\frac{\partial\tilde{B}^{A^{\prime}B^{\prime}C^{\prime}}(k_{1},k_{2},k_{3},\theta,\phi)}{\partial p_{j}},$

(32)

where the sum is over $n_{\rm shape}$ allowed triangle shape bins with centers denoted by $(k_{1},k_{2},k_{3})$ and $n_{\rm ori}$ orientation bins with centers denoted by $(\theta,\phi)$, and there is no correlation between different triangle shapes and orientations in the Gaussian approximation for the covariance matrix. There is however a correlation between the different multi-tracer combinations, so the covariance $\tilde{\mathrm{Cov}}$ is a $n_{b}\times n_{b}$ matrix where $n_{b}=n_{\rm tracers}^{3}$ is the number of multi-tracer combinations. It can be written as

$$\mathrm{\tilde{Cov}}=\mathrm{Cov}\,\frac{V}{N_{\rm modes}},$$

(33)

where in the Gaussian approximation we have for a single mode

			$\displaystyle\langle\delta_{g}^{A}(\bm{k}_{1})\,\delta_{g}^{B}(\bm{k}_{2})\,\delta_{g}^{C}(\bm{k}_{3})\,\delta_{g}^{A^{\prime}}(\bm{k}^{\prime}_{1})\,\delta_{g}^{B^{\prime}}(\bm{k}^{\prime}_{2})\,\delta_{g}^{C^{\prime}}(\bm{k}^{\prime}_{3})\,\rangle$		(34)
		$\displaystyle\approx$	$\displaystyle\left(P_{gg}^{AA^{\prime}}(\bm{k}_{1})+\frac{\delta_{AA^{\prime}}}{\bar{n}_{g}^{A}}\right)\left(P_{gg}^{BB^{\prime}}(\bm{k}_{2})+\frac{\delta_{BB^{\prime}}}{\bar{n}_{g}^{B}}\right)\left(P_{gg}^{CC^{\prime}}(\bm{k}_{3})+\frac{\delta_{CC^{\prime}}}{\bar{n}_{g}^{C}}\right)\delta_{D}(\bm{k}_{1}+\bm{k}^{\prime}_{1})\delta_{D}(\bm{k}_{2}+\bm{k}^{\prime}_{2})\delta_{D}(\bm{k}_{3}+\bm{k}^{\prime}_{3})$		(35)
		$\displaystyle\equiv$	$\displaystyle\mathrm{Cov}\left[B_{ggg}^{ABC}(\bm{k}_{1},\bm{k}_{2}),B_{ggg}^{A^{\prime}B^{\prime}C^{\prime}}(\bm{k}^{\prime}_{1},\bm{k}^{\prime}_{2})\right]\;\delta_{D}(\bm{k}_{1}+\bm{k}^{\prime}_{1})\delta_{D}(\bm{k}_{2}+\bm{k}^{\prime}_{2})\delta_{D}(\bm{k}_{3}+\bm{k}^{\prime}_{3}),$		(36)

and where

$$\frac{V}{N_{\rm modes}}\,=\frac{(2\pi)^{5}}{V(\mathrm{d}k_{1}\,\mathrm{d}k_{2}\,\mathrm{d}k_{3}\,k_{1}k_{2}k_{3})\,(d\mu d\phi)\,\beta}.$$

(37)

Here the number of modes is the number of closed triangles $N_{\rm modes}$ within a triangle shape and orientation bin given a survey volume which sets the fundamental frequency $k_{F}\equiv(2\pi)V^{-1/3}$, where $V^{-1/3}$ is the volume of the survey in the given redshift bin. In the limit that the bin width $\Delta k_{i}\gg k_{F}$, the following expression is a good approximation

			$\displaystyle N_{\rm modes}=\frac{K_{\Delta}}{k_{F}^{6}}$
		$\displaystyle=$	$\displaystyle\frac{V^{2}}{(2\pi)^{6}}\left[8\pi^{2}k_{1}k_{2}k_{3}\Delta k_{1}\Delta k_{2}\Delta k_{3}\beta\left(\frac{1}{4\pi}\Delta\mu_{1}\Delta\phi_{12}\right)\right],$

where $8\pi^{2}k_{1}k_{2}k_{3}\Delta k_{1}\Delta k_{2}\Delta k_{3}\beta$ is the number of closed triangles in a triangle shape bin denoted by the bin centers ($k_{1},k_{2},k_{3}$) with bin widths $\Delta k_{1},\Delta k_{2}$ and $\Delta k_{3}$, and $\beta=0.5$ for degenerate triangles and 1 otherwise. The factor $\frac{1}{4\pi}\Delta\mu_{1}\Delta\phi_{12}$ corresponds to the fraction of triangles in a fixed triangle shape bin that falls into the orientation bin with centers $(\mu,\phi_{12})$.

Note that we do not include the factor $s_{B}$ = 6, 2, 1 for equilateral, isoceles and scalene triangles respectively, often used in the orientation-averaged case. Since here we have triangle orientation bins, we instead have that triangles of the same equilateral or isoceles shape but different orientations can be correlated. These would appear as off-diagonal elements of the $n_{\rm ori}\times n_{\rm ori}$ covariance matrix block for a fixed triangle shape. Furthermore, the $s_{B}$ factor is not valid for multi-tracers, since the other possible contractions in Eq. 34 are not necessarily equal to the first contraction we wrote down in Eq. 36 in the multi-tracer case. For our calculations however, we ignore the effects from those off-diagonal terms because a bin with an equilateral or isoceles triangle as its representative center contains many triangles that are not equilateral or isoceles.

Finally, the length of our multi-tracer bispectrum data vector for each redshift bin is $n=n_{b}\times n_{\rm shape}\times n_{\rm ori}=125\times 10^{4}\times 25\approx 10^{7}$ for $n_{k}=50$, $n_{\theta}=5$, $n_{\phi}=5$. The exact number of triangle shapes $n_{\rm shape}$ varies for different redshift bins, since the values of $k_{\rm min}=2\pi V^{-1/3}$ and $k_{\rm max}=0.2\ h\text{Mpc}^{-1}(1+z)$ are dependent on the redshift bin. We note also that we have used an approximate $k_{\rm min}$ value corresponding to a cubic volume, whereas in reality, the redshift bins are shaped as spherical shells in the full sky limit, or part of a spherical shell when the survey window function is applied. A precise mode counting treatment with the exact window function will be evaluated in a future work.

For the bispectrum multipoles, we have the following Fisher matrix

	$\displaystyle F_{ij}=\int dk_{1}dk_{2}dk_{3}k_{1}k_{2}k_{3}\sum_{(l,m)}\sum_{(l^{\prime},m^{\prime})}\sum_{(ABC)}\sum_{(A^{\prime}B^{\prime}C^{\prime})}$
	$\displaystyle\frac{\partial B_{lm}^{ABC*}(k_{1},k_{2},k_{3})}{\partial p_{i}}\left[\frac{(2\pi)^{5}}{V}\mathrm{Cov}\right]^{-1}\frac{\partial B_{l^{\prime}m^{\prime}}^{A^{\prime}B^{\prime}C^{\prime}}(k_{1},k_{2},k_{3})}{\partial p_{j}}.$

Let $\tilde{B}_{lm}^{ABC}(k_{1},k_{2},k_{3})$ represent the bispectrum binned over triangle shapes that fall into the bin specified by the centers $(k_{1},k_{2},k_{3})$. Let this binned value be approximated by the bispectrum value at this center configuration, and we have that the Fisher matrix can be approximated as

	$\displaystyle F_{ij}=$		$\displaystyle\sum_{(k_{1},k_{2},k_{3})}\sum_{(l,m)}\sum_{(l^{\prime},m^{\prime})}\sum_{(ABC)}\sum_{(A^{\prime}B^{\prime}C^{\prime})}$
			$\displaystyle\frac{\partial\tilde{B}_{lm}^{ABC*}(k_{1},k_{2},k_{3})}{\partial p_{i}}\left[\mathrm{\tilde{Cov}}\right]^{-1}\frac{\partial\tilde{B}_{l^{\prime}m^{\prime}}^{A^{\prime}B^{\prime}C^{\prime}}(k_{1},k_{2},k_{3})}{\partial p_{j}}.$

The outer sum is taken over $n_{\rm shape}$ unique triangle shapes $(k_{1},k_{2},k_{3})$, where $k_{1},k_{2}$ and $k_{3}$ are values of the $k$-bin centers that satisfy the triangle inequality ($0<k_{1}\leq k_{2}\leq k_{3}\leq k_{1}+k_{2}$). There is also a sum over $n_{lm}$ pairs of $(l,m)$ values, where $l$ and $m$ are integers satisfying $0\leq l\leq l_{\rm max}$ and $-l\leq m\leq l$ (and similarly for the $(l^{\prime},m^{\prime})$ pairs), as well as a sum over $n_{b}=n_{\rm tracers}^{3}$ multi-tracer combinations ($ABC$) where $A,B,C=1..n_{\rm tracers}$ (and similarly for ($A^{\prime}B^{\prime}C^{\prime}$)).

Because there is no correlation between different triangle shapes (in the Gaussian covariance), whereas there is a correlation between different $(l,m)$ pairs and multi-tracer combinations, the covariance $\tilde{\mathrm{Cov}}$ of the binned bispectrum is a $N\times N$ matrix, where $N=n_{b}n_{lm}$. It can be calculated using

$$\mathrm{\tilde{Cov}}=\mathrm{Cov}\,\frac{(2\pi)^{5}s_{B}}{V\mathrm{d}k_{1}\mathrm{d}k_{2}\mathrm{d}k_{3}k_{1}k_{2}k_{3}\beta},$$

(41)

where

	$\displaystyle\mathrm{Cov}[B_{lm}^{ABC*}$	,	$\displaystyle B_{l^{\prime}m^{\prime}}^{A^{\prime}B^{\prime}C^{\prime}}]=\frac{1}{(4\pi)^{2}}\int\mathrm{d(cos}\theta)\mathrm{d}\phi\,$
		$\displaystyle\times$	$\displaystyle Y_{lm}^{*}(\theta,\phi)Y_{l^{\prime}m^{\prime}}(\theta,\phi)\left(P^{AA^{\prime}}(\bm{k}_{1})+\frac{1}{\bar{n}^{A}_{g}}\right)$
		$\displaystyle\times$	$\displaystyle\left(P^{BB^{\prime}}(\bm{k}_{2})+\frac{1}{\bar{n}^{B}_{g}}\right)\left(P^{CC^{\prime}}(\bm{k}_{3})+\frac{1}{\bar{n}^{C}_{g}}\right),$

and $\beta=0.5$ for degenerate triangles and 1 otherwise. The factor of $1/(4\pi)^{2}$ comes from the normalization of the bispectrum multipoles in Eq. LABEL:eq:blm_definition.

The marginalized error on parameter $p_{i}$ is obtained using

$$\sigma_{p_{i}}=\left[F^{-1/2}\right]_{ii}.$$

(43)

The length of the data vector is slightly reduced here, with $n=n_{b}\times n_{\rm shape}\times n_{\rm lm}$ for a single redshift bin, where $n_{\rm lm}=25$ when all $(l,m)$ pairs are used up to $l_{\rm max}=4$, which can be reduced $n_{\rm lm}=2$ without much loss of information if only the even $l$ and $m=0$ modes are used up to $l_{\rm max}=2$, as we will show later in Section V.

We briefly summarize the galaxy bias modeling we chose and refer the readers to Ref. Tellarini et al. (2016) for more details. Recall that the galaxy density field up to second order is described by

	$\displaystyle\delta_{\mathrm{g}}(\bm{x})$	$\displaystyle=$	$\displaystyle b_{10}\delta_{\rm m}(\bm{x})+b_{01}\varphi(\bm{x})$		(44)
		$\displaystyle+$	$\displaystyle b_{20}\left(\delta_{\rm m}(\bm{x})\right)^{2}+b_{11}\delta_{\rm m}(\bm{x})\varphi(\bm{x})$
		$\displaystyle+$	$\displaystyle b_{02}\left(\varphi(\bm{x})\right)^{2}+b_{s_{2}}(s^{2}-\langle s^{2}\rangle)-b_{01}n^{2},$

where we have now absorbed the factors of 1/2 into $b_{20}$, $b_{02}$ and $b_{s_{2}}$. There are a total of 6 different kind of bias parameters to model: $b_{10}$, $b_{01}$, $b_{20}$, $b_{11}$, $b_{02}$ and $b_{s_{2}}$, each having 55 distinct values for the 5 galaxy samples and 11 redshift bins.

As previously noted, we assume an universal mass function and use the following relation for $b_{01}$:

$$b_{01}=2f_{\mathrm{NL}}\delta_{c}(b_{10}-1).$$

(45)

For $b_{20}$, we treat it as a function of $b_{10}$ using a fit from simulations in Ref. Lazeyras et al. (2016),

$$b_{20}^{\rm Lezeyras}(b_{10})=\frac{1}{2}(0.412-2.143\,b_{10}+0.929\,b_{10}^{2}+0.008\,b_{10}^{3}).$$

(46)

Similarly for $b_{s_{2}}$, we use a relation to $b_{10}$ obtained from a fit to simulation results from Ref. Baldauf et al. (2012),

$$b_{s_{2}}^{\rm Saito}(b_{10})=-\frac{2}{7}\,(b_{10}-1).$$

(47)

For $b_{11}$ and $b_{02}$, we can relate them to any given values of $b_{10},b_{20}$ and $f_{\mathrm{NL}}$ using,

$$b_{11}=b_{11}^{L}+b_{01}^{L},$$

(48)

and

$$b_{02}=f_{\mathrm{NL}}^{2}4\delta_{c}(\delta_{c}b_{20}^{L}-2b_{10}^{L}),$$

(49)

where the $L$ superscript denotes Lagrangian bias, and we suppressed the $E$ superscript for the Eulerian biases, and the Lagrangian quantities may be related to the Eulerian ones as follows:

$$b_{11}^{L}=2f_{\mathrm{NL}}(\delta_{c}b_{20}^{L}-b_{10}^{L}),$$

(50)

$$b_{20}^{L}=b_{20}-\frac{8}{21}(b_{10}-1),$$

(51)

$$b_{10}^{L}=b_{10}-1,$$

(52)

and

$$b_{01}^{L}=b_{01},$$

(53)

The above results were obtained under the assumption of the universal mass function, conservation of the galaxy number density in a given volume, spherical collapse and no velocity bias between galaxies and matter.

In the Fisher forecast, we vary the parameters $b_{10}$ for each sample and redshift bin when taking the derivative with respect to $b_{10}$, and let all other bias values vary according to the relations described above.

The survey parameters have not changed significantly since the original SPHEREx forecast. There are eleven redshift bins ranging from $z=0$ to $4.6$ (see definitions in Table 2 in Appendix A), with the galaxies in each redshift bin divided by their redshift uncertainties falling in the bins $\tilde{\sigma}_{z}=\sigma_{z}/(1+z)=0-0.003-0.01-0.03-0.1-0.2$. We use the maximum value of the redshift bin in our forecast for more conservative results: $\tilde{\sigma}_{z}^{A}=0.003,0.01,0.03,0.1$ and $0.2$ for $A=1$ to 5 respectively. So a lower sample number means better redshift uncertainties. Note that in the original bispectrum forecast, it was the mean value of $\tilde{\sigma}_{z}^{A}$ that was used, leading to slightly more optimistic but also more realistic results.

The details of the procedure for obtaining the fiducial galaxy number densities and biases are found in Ref. Doré et al. (2014), which we briefly summarize here. First, the simulated SPHEREx galaxy catalog (based on COSMOS Scoville et al. (2007); Ilbert et al. (2009)) was piped through a template-fitting based photometric redshift measurement pipeline in order to produce a redshift and a redshift error estimate for each galaxy. These estimates were used to derive the galaxy number density for the galaxy sample in each redshift bin and $\tilde{\sigma}_{z}$ bin.

The method of abundance matching was then used to obtain an estimate of the linear galaxy bias for each galaxy sample: The galaxies with the best redshift errors are matched to the host halos with the largest total mass. More specifically, the mass function in Ref. Tinker et al. (2008) was used to find the minimum halo mass, for which the halo bias is found using the fitting formula in Ref. Tinker et al. (2010) and set to the linear galaxy bias. The fiducial values for the linear biases and number densities are available at the SPHEREx public Github repository¹¹1SPHEREx public github repository: https://github.com/SPHEREx/Public-products/blob/master/galaxy_density_v28_base_cbe.txt, and also reproduced in Appendix A for convenience in Tables 3 and 4 respectively.

To illustrate how much of the total constraint is captured by a set of bispectrum multipoles, we look at the improvement ratio defined as

$$\mathrm{Improvement\ ratio}=\frac{\sigma_{f_{\mathrm{NL}}}(B_{\rm full})}{\sigma_{f_{\mathrm{NL}}}(B_{lm})},$$

(54)

where $B_{\rm full}$ stands for the full constraint obtained by using the Fourier space bispectrum with all the triangle shape and orientation bins, whereas $B_{lm}$ stands for bispectrum multipoles with a specific set of $l,m$ values. In this notation, the improvement ratio will always be less than one for the bispectrum multipoles, and higher means better constraints.