NUHEP-TH/22-05

  Abstract We analyze the BOSS power spectrum monopole and quadrupole, and the bispectrum monopole and quadrupole data, using the predictions from the Effective Field Theory of Large-Scale Structure (EFTofLSS). Specifically, we use the one loop prediction for the power spectrum and the bispectrum monopole, and the tree level for the bispectrum quadrupole. After validating our pipeline against numerical simulations as well as checking for several internal consistencies, we apply it to the observational data. We find that analyzing the bispectrum monopole to higher wavenumbers thanks to the one-loop prediction, as well as the addition of the tree-level quadrupole, significantly reduces the error bars with respect to our original analysis of the power spectrum at one loop and bispectrum monopole at tree level. After fixing the spectral tilt to Planck preferred value and using a Big Bang Nucleosynthesis prior, we measure $\sigma_{8}=0.794\pm 0.037$, $h=0.692\pm 0.011$, and $\Omega_{m}=0.311\pm 0.010$ to about $4.7\%$, $1.6\%$, and $3.2\%$, at $68\%$ CL, respectively. This represents an error bar reduction with respect to the power spectrum-only analysis of about $30\%$, $18\%$, and $13\%$ respectively. Remarkably, the results are compatible with the ones obtained with a power-spectrum-only analysis, showing the power of the EFTofLSS in simultaneously predicting several observables. We find no tension with Planck.

 

1 Introduction, Main Results and Conclusion

The SDSS-III Baryon Oscillation Spectroscopic Survey (BOSS) has mapped the clustering of galaxies in the nearby Universe in an unprecedented amount and with great accuracy [1]. Although BOSS’ survey volume is modest with respect to upcoming experiments such as DESI [2] or Euclid [3], the BOSS data are remarkable as they have been revealing a wealth of cosmological information from the large-scale structure of the Universe.

In the last couple of years, the Effective Field Theory of Large-Scale Structure (EFTofLSS) prediction at one-loop order has been used to analyze the BOSS Full Shape (FS) of the galaxy Power Spectrum (PS) [4, 5, 6], and Correlation Function (CF) [7, 8]. The BOSS galaxy-clustering bispectrum monopole using the tree-level prediction was first analyzed in [4] (see [9] for a recent slight generalization). See also [10, 11, 12] for other techniques and analysis using linear theory with higher multipoles. All $\Lambda$CDM cosmological parameters have been measured from these data by only imposing a prior from Big Bang Nucleosynthesis (BBN), reaching a remarkable, and perhaps surprising, precision on some of these. For example, the present amount of matter, $\Omega_{m}$, and the Hubble constant (see also [13, 14] for subsequent refinements) have error bars that are not far from the ones obtained from the Cosmic Microwave Background (CMB) [15]. For clustering and smooth quintessence models, limits on the dark energy equation of state $w$ parameter of $\lesssim 5\%$ have been set using only late-time measurements [14, 16]. This is again quite close to the ones obtained with the CMB [15]. These measurements provide a new, CMB-independent, method for determining the Hubble constant [4], resulting in a measurement that is comparable, if not better, to the one based on the cosmic ladder [17, 18] and CMB. Therefore, this tool has been used to shed light on how some models that were proposed to alleviate the tension in the Hubble measurements (see e.g. [19]) between the CMB and cosmic ladder [20, 21] (see also [22, 23]) actually perform.

Very recently, in [24], we used the one-loop EFTofLSS prediction for the bispectrum to set the first and strong limits on primordial inflationary non-Gaussianities from Large-Scale Structure (LSS) (see also [25, 26] for a contemporary and a subsequent paper, where, once put together, the same shapes are constrained but stopping at the tree-level EFTofLSS prediction, and so obtaining much weaker constraints for the same data). We obtained limits on three of the so-called $f_{\rm NL}$ parameters, $f_{\rm NL}^{\rm equil.}=217\pm 297\,,\ f_{\rm NL}^{\rm orth.}=-64\pm 74\,,\ f_{\rm NL}^{\rm loc.}=49\pm 36$, at $68\%$ confidence level, which are predicted to be produced by some single-clock [27, 28] or multiple fields [29, 30, 31, 32, 33] inflationary models. Perhaps quite surprisingly, those constraints were already quite on par with the ones of the powerful CMB experiment WMAP [34], though largely inferior to the more recent CMB experiment Planck [35]. Significant limits from LSS on just $f_{\rm NL}^{\rm loc.}$ were obtained using the power spectrum only, first in [36], using the so-called non-local bias [37, 38, 39], but the analysis of [24] uses for the first time the bispectrum, obtaining much stronger constraints using the data from the same experiments.

It took an intense and years-long line of study to develop the EFTofLSS from the initial formulation to the level that allows it to be applied to data. It appears to us that it often happens that there is no proper acknowledgment of the many works that were needed to reach this point. For instance, several authors cite Refs. [4, 5] for the ‘model’ to analyze the PS FS, but the EFT model that is used in [4, 5] is essentially the same as the one originally proposed in [40]. We therefore find it fair to add the following footnote in every paper where the EFTofLSS is used to analyze observational data. Even though some of the mentioned papers are not strictly required to analyze the data, we, and we believe probably anybody else, would not have applied the EFTofLSS to data without all these intermediate results that allowed us to overcome the widespread skepticism about the usefulness of the EFTofLSS.

In this paper we upgrade our original analysis of the one-loop power spectrum monopole and quadrupole and tree-level bispectrum monopole [4], to include the full one-loop bispectrum monopole and the tree-level bispectrum quadrupole. We scan over all the $\Lambda$CDM parameters with BBN prior on the baryon abundance, $\Omega_{b}h^{2}$, with the exception of the tilt, $n_{s}$, that we fix to the Planck preferred value.

The development of a pipeline that allows us to analyze the one-loop bispectrum predicted by the EFTofLSS has required much theoretical work, and the explanation of such techniques will be presented in two upcoming papers [90, 91]. While the application to constrain primordial non-Gaussianities was already presented in [24], here, instead, we will just give the essential details and focus on constraints of the $\Lambda$CDM parameters.

Our main results are summarized in fig. 1, where we plot the posteriors on the cosmological parameters that are effectively scanned. This analysis improves the error bars on the $\Lambda$CDM parameters $\sigma_{8}$, $h$, and $\Omega_{m}$ with respect to the power spectrum-only analysis by about 30%, 18%, and 13% respectively, achieving a precision of about $4.7\%$, $1.6\%$, and $3.2\%$ at $68\%$ CL, respectively. Notice also that the results improve significantly upon the ones obtained using instead the tree-level prediction for the bispectrum monopole: in particular, $\sigma_{8}$ is better determined by about $30\%$. Naively, a $30\%$ improvement corresponds to doubling the data volume of the survey. As it can be seen in the same figure, the results are compatible with the ones obtained with a power-spectrum-only analysis. We find no tension with Planck: we measure $\sigma_{8}$, $h$, and $\Omega_{m}$ to values consistent at $0.3\sigma,1.4\sigma,0.5\sigma$, respectively, with the ones of Planck $\nu\Lambda$CDM [15].

best-fit $\text{\rm mean}\pm\sigma$	$\Omega_{m}$	$h$	$\sigma_{8}$	$\omega_{cdm}$	$\ln(10^{10}A_{s})$	$S_{8}$
$P_{\ell}$	0.2984 $0.308\pm 0.012$	0.6763 $0.689^{+0.012}_{-0.014}$	0.8305 $0.819^{+0.049}_{-0.055}$	0.1143 $0.1232\pm 0.0075$	3.123 $3.02\pm 0.15$	0.8283 $0.830^{+0.051}_{-0.060}$
$P_{\ell}+B_{0}^{\mathrm{tree}}$	0.3101 $0.309\pm 0.011$	0.6907 $0.691\pm 0.012$	0.8063 $0.804\pm 0.049$	0.1248 $0.1246\pm 0.0058$	2.98 $2.97\pm 0.13$	0.8197 $0.816^{+0.050}_{-0.057}$
$P_{\ell}+B_{0}^{\mathrm{1loop}}$	0.3210 $0.314\pm 0.011$	0.6956 $0.693\pm 0.011$	0.7882 $0.790^{+0.033}_{-0.037}$	0.1331 $0.1278\pm 0.0061$	2.82 $2.90\pm 0.11$	0.8153 $0.807^{+0.037}_{-0.043}$
$P_{\ell}+B_{0}^{\mathrm{1loop}}+B_{2}^{\mathrm{tree}}$	$0.3082$ $0.311\pm 0.010$	$0.6928$ $0.692\pm 0.011$	$0.7856$ $0.794\pm 0.037$	$0.1258$ $0.1255\pm 0.0057$	$2.88$ $2.94\pm 0.11$	$0.7962$ $0.808\pm 0.041$
Planck	$0.3191^{+0.0085}_{-0.016}$	$0.671^{+0.012}_{-0.0067}$	$0.807^{+0.018}_{-0.0079}$	$0.1201\pm 0.0013$	$3.046\pm 0.015$	$0.832\pm 0.013$

Figure 1: Triangle plots, best-fit values, and relative $68\%$-credible intervals of base cosmological parameters measured from the analysis of BOSS power spectrum multipoles $P_{\ell}$, $\ell=0,2$, at one-loop, bispectrum monopole $B_{0}$ at tree or one-loop level, and bispectrum quadrupole $B_{2}$ at tree-level. Planck $\nu\Lambda$CDM results are shown for comparison.

The paper is organized as follows. In sec. 2 we describe the data products and the measurements we use. In sec. 3 we describe the theory model including the observational aspects. In sec. 4, we present the likelihood we use to describe the data. In sec. 5, we provide some tests for our pipeline. Finally, in sec. 6, we provide some additional details about the main results. Technical aspects and additional materials are relegated to the appendices.

A note of warning:

We end this section of the main results with a final note of warning. It should be emphasized that in performing this analysis, as well as the preceding ones using the EFTofLSS by our group [4, 6, 14, 20, 16, 7, 24], we have assumed that the observational data are not affected by any unknown systematic error or undetected foregrounds. In other words, we have simply analyzed the publicly available data: the two- and three-point functions of the galaxy density in redshift space as measured from the public galaxy catalogues. Given the additional cosmological information that the theoretical modeling by the EFTofLSS allows us to exploit in BOSS data, it might be worthwhile to investigate if potential undetected systematic errors might affect our results. We leave an investigation of these issues to future work.

2 Data

Figure 2: Measurements and best fits of bispectrum monopole $B_{0}$ (top), power spectrum multipoles $P_{\ell}$ (bottom left), and bispectrum quadrupole $B_{2}$ (bottom right) from BOSS (points and error bars) and 2048 Patchy (grey regions) CMASS NGC sky. The bispectrum is shown in bins ordered by their central values forming either an equilateral, isoceles, or scalene triangle, shown in blue, orange, or green, respectively. The bin triangle sides (top panel) are shown either by the bin central values (colored lines) or by their effective values (grey points). The best fit (black points) is shown only for the scales analyzed. The relative error bars (turquoise regions) are shown with the best fit residuals for comparison. While only CMASS NGC is shown for clarity, the best fit depicted here is obtained fitting the full combination $P_{\ell}+B_{0}+B_{2}$ on all BOSS 4 skies.

BOSS DR12 LRG sample.

The main data sample analyzed in this work is the SDSS-III BOSS DR12 luminous red galaxies (LRG) sample [1]. We use the BOSS catalogs DR12 (v5) combined CMASS-LOWZ [92]. To each galaxy we assign the standard FKP weights for optimality together with the correction weights described in [92] for BOSS data and in [93] for the patchy mocks. The inverse covariances are corrected by the Hartlap factor to account for the finite number of mocks used in their estimation [94]. In order to test our analysis pipeline, we will analyze the mean over the 2048 Patchy mocks of CMASS NGC (hereafter referred as ‘Patchy’). We will also make use of the Nseries mocks, which are full $N$-body simulations populated with a Halo Occupation Distribution (HOD) model and selection function similar to the one of BOSS CMASS NGC [1]. We will analyze the mean of the 84 Nseries realizations (hereafter referred as ‘Nseries’). All celestial coordinates are converted to comoving distance assuming $\Omega_{m}^{\rm fid}=0.310$.

BOSS P+B full-shape measurements.

In this work, we analyze the full shape of the power spectrum multipoles $\ell=0,2,$ and of the bispectrum monopole and quadrupole (respectively abbreviated ‘$P_{\ell}$’, ‘$B_{0}$’ and ‘$B_{2}$’). Those measurements are shown in fig. 2 (together with the best fit from our theory model that we discuss later). The estimator for the power spectrum is the standard ‘FKP’ estimator [95], generalized to redshift space in [96, 97, 98]. The bispectrum is estimated using the estimator outlined in [99] (see also [100, 101, 67, 102]). The measurements are obtained using the code Rustico [99]. For the power spectrum, we find excellent agreement between the measurements from Rustico and Nbodykit [103]. We use Nbodykit to measure the window functions as described in [104], with consistent normalization in the power spectrum as discussed in [105, 106, 107].

The configurations of the measurements are the following. We use a box of side length $L_{\rm box}=3500\,(2300){\rm Mpc}/h$ for CMASS (LOWZ), with Piecewise Cubic Spline (PCS) particle assignment scheme and grid interlacing as described in [108]. The grid is consisting of $512^{3}$ cells. The power spectrum is binned in $\Delta k\simeq 0.01h\,{\rm Mpc}^{-1}$. Instead, we bin the bispectrum in $\Delta n=12\,(9)$ units of the fundamental frequency of the box $k_{f}$ for CMASS (LOWZ), starting from the bin centered at $n_{\rm min}=6+\Delta n/2$, up to the one centered on $n_{\rm max}=126\,(69)-\Delta n/2$, which correspond in frequencies to bins of size $\Delta k=0.02154\,(0.02459)h\,{\rm Mpc}^{-1}$, with first and last bins centered on $k_{\rm min}=0.0215\,(0.029)h\,{\rm Mpc}^{-1}$ and $k_{\rm max}=0.215\,(0.176)h\,{\rm Mpc}^{-1}$, respectively. This choice of bin size is motivated to keep the Hartlap factor at a value safely close to $1$ to limit the effect from the bias of the inverse covariance estimator. Given that we have $2048$ patchy mocks at our disposal to estimate the covariance, and, in our analysis, we will analyze $42\,(36)$ $k$-bins in $P_{\ell}$ and $150\,(62)$ triangle bins in $B_{0}$, this makes for the Hartlap factor of about $0.91\,(0.95)$ for CMASS (LOWZ). $B_{2}$, as analyzed at tree-level, only adds $9$ bins per quadrupole (for both CMASS and LOWZ), which lead to the Hartlap factor of the same order of about $0.9$. Importantly, we keep all bins whose centers form a closed triangle. Explicitly, we choose the following bins according to their centers ordered as:

$$\begin{split}(n_{1},n_{2},n_{3})\,,\quad&n_{1},n_{2},n_{3}=n_{\rm min},n_{\rm min}+dn,\dots,n_{\rm max}\,,\\ &\text{if }n_{1}\leq n_{2}\leq n_{3}\text{ and }n_{3}\leq n_{1}+n_{2}\,.\end{split}$$ (1)

It follows that there are several bins that contain fundamental triangles that are not closed. How to properly account for them is discussed in sec. 3.5.

3 Theory model

Our model for the power spectrum multipoles $P_{\ell}$, $\ell=0,2$, the bispectrum monopole, $B_{0}$, and the bispectrum quadrupole, $B_{2}$, consists in the prediction of EFTofLSS at one loop for $P_{\ell}$ and $B_{0}$, and at tree-level for $B_{2}$. We also incorporate a number of observational effects in our modeling to make contact with the measurements.

3.1 EFTofLSS at one loop

In this section, we outline the relevant expressions for the one-loop power spectrum and bispectrum for halos in redshift space $P^{r,h}$ and $B^{r,h}$. The power spectrum and bispectrum are defined as the 2- and 3- point functions of the halo overdensity $\delta_{r,h}$, in Fourier space:

$$\langle\delta_{r,h}(\vec{k};\hat{z})\delta_{r,h}(\vec{k}^{\prime};\hat{z})\rangle = (2\pi)^{3}\delta_{D}(\vec{k}+\vec{k}^{\prime})P^{r,h}(k,\hat{k}\cdot\hat{z})$$ (2)

$$\langle\delta_{r,h}(\vec{k}_{1};\hat{z})\delta_{r,h}(\vec{k}_{2};\hat{z})\delta_{r,h}(\vec{k}_{3};\hat{z})\rangle = (2\pi)^{3}\delta_{D}(\vec{k}_{1}+\vec{k}_{2}+\vec{k}_{3})B^{r,h}(k_{1},k_{2},k_{3},\hat{k}_{1}\cdot\hat{z},\hat{k}_{2}\cdot\hat{z})\ ,$$

where $\hat{z}$ is the line-of-sight direction, and $\delta_{D}$ is the Delta dirac function. After perturbatively expanding the halo overdensity, we arrive at the following expressions for the one-loop power spectrum:

$$P^{r,h}_{1\text{-loop tot.}}=P_{11}^{r,h}+(P_{13}^{r,h}+P_{13}^{r,h,ct})+(P_{22}^{r,h}+P_{22}^{r,h,\epsilon})\ ,$$ (3)

and the one-loop bispectrum:

$$\begin{split}B^{r,h}_{1\text{-loop tot.}}&=B_{211}^{r,h}+(B_{321}^{r,h,(II)}+B_{321}^{r,h,(II),ct})+(B_{411}^{r,h}+B_{411}^{r,h,ct})\\ &\hskip 72.26999pt+(B^{r,h}_{222}+B_{222}^{r,h,\epsilon})+(B_{321}^{r,h,(I)}+B_{321}^{r,h,(I),\epsilon})\ ,\end{split}$$ (4)

where we have grouped perturbation theory contributions with the counterterm contributions that renormalize them in parentheses. The loop integrals are evaluated using the techniques described in [91]. We note that the arguments of the above functions are the same as those given in eq. (2), and we will often drop the arguments below for clarity.

Details for all of the above expressions can be found in app. A, but we summarize the dependence on bias parameters and EFT parameters here for convenience. For the perturbation theory contributions, we have

$$\begin{split}&P_{11}^{r,h}[b_{1}]\ ,\quad P_{13}^{r,h}[b_{1},b_{3},b_{8}]\ ,\quad P_{22}^{r,h}[b_{1},b_{2},b_{5}]\ ,\\ &B_{211}^{r,h}[b_{1},b_{2},b_{5}]\ ,\quad B_{321}^{r,h,(II)}[b_{1},b_{2},b_{3},b_{5},b_{8}]\ ,\quad B_{411}^{r,h}[b_{1},\dots,b_{11}]\ ,\\ &B^{r,h}_{222}[b_{1},b_{2},b_{5}]\ ,\quad B_{321}^{r,h,(I)}[b_{1},b_{2},b_{3},b_{5},b_{6},b_{8},b_{10}]\ ,\end{split}$$ (5)

while for the counterterms, we have

$$\begin{split}&P_{13}^{r,h,ct}[b_{1},c_{h,1},c_{\pi,1},c_{\pi v,1},c_{\pi v,3}]\ ,\quad P_{22}^{r,h,\epsilon}[c^{\rm St}_{1},c^{\rm St}_{2},c^{\rm St}_{3}]\ ,\\ &B_{321}^{r,h,(II),ct}[b_{1},b_{2},b_{5},c_{h,1},c_{\pi,1},c_{\pi v,1},c_{\pi v,3}]\ ,\quad B_{321}^{r,h,{(I),\epsilon}}[b_{1},c^{\rm St}_{1},c^{\rm St}_{2},\{c^{\rm St}_{i}\}_{i=4,\dots,13}]\ ,\\ &B_{411}^{r,h,ct}[b_{1},\{c_{h,i}\}_{i=1,\dots,5},c_{\pi,1},c_{\pi,5},\{c_{\pi v,j}\}_{j=1,\dots,7}]\ ,\quad B_{222}^{r,h,\epsilon}[c^{(222)}_{1},c^{(222)}_{2},c^{(222)}_{5}]\ .\end{split}$$ (6)

Notice that the diagrams $P_{13}^{r,h}$, $B_{321}^{r,h,(II)}$, and $B_{411}^{r,h}$ depend on less biases than the kernels in eq. (LABEL:kernelbiasdep) would suggest. This is because, when considering the particular momentum-configuration of the kernels that enter the loop in eq. (LABEL:loopexpressionsrssbias) and eq. (A.1), they are degenerate with EFT parameters.

To make contact with our measurements described in sec 2, what we analyze in the data are various multipoles with respect to the line-of-sight $\hat{z}$. In particular, we analyze the power-spectrum and bispectrum monopole and quadrupoles. The power-spectrum multipoles are given by

$$P_{\ell}^{r,h}(k)=\frac{2\ell+1}{2}\int_{-1}^{1}\mathrm{d}\mu\,\mathcal{P}_{\ell}(\mu)P^{r,h}(k,\mu)\ ,$$ (7)

where $\mathcal{P}_{\ell}$ are the Legendre polynomials, and $\mu=\hat{k}\cdot\hat{z}$. The bispectrum monopole is the average over the line-of-sight angles [109, 102, 110]

$$B_{0}^{r,h}(k_{1},k_{2},k_{3})=\frac{1}{4\pi}\int_{-1}^{1}\mathrm{d}\mu_{1}\int_{0}^{2\pi}\mathrm{d}\phi\,B^{r,h}(k_{1},k_{2},k_{3},\mu_{1},\mu_{2}(\mu_{1},\phi))\ ,$$ (8)

where $\mu_{i}=\hat{k}_{i}\cdot\hat{z}$, and explicitly, from the triangle conditions:

$$\mu_{2}(\mu_{1},\phi) =\mu_{1}\hat{k}_{1}\cdot\hat{k}_{2}+\sqrt{1-\mu_{1}^{2}}\sqrt{1-(\hat{k}_{1}\cdot\hat{k}_{2})^{2}}\,\sin\phi\ ,$$ (9)

$$\mu_{3}(\mu_{1},\phi) =-k_{3}^{-1}(k_{1}\mu_{1}+k_{2}\mu_{2}(\mu_{1},\phi))\ .$$ (10)

The expectation values of the estimator used by Rustico for the quadrupoles are:

$$\begin{split}&B_{{(2,1)}}^{r,h}(k_{1},k_{2},k_{3})\equiv\frac{5}{4\pi}\int_{-1}^{1}\mathrm{d}\mu_{1}\int_{0}^{2\pi}\mathrm{d}\phi\,\mathcal{P}_{2}(\mu_{1})B^{r,h}(k_{1},k_{2},k_{3},\mu_{1},\mu_{2}(\mu_{1},\phi))\ ,\\ &B_{{(2,2)}}^{r,h}(k_{1},k_{2},k_{3})\equiv\frac{5}{4\pi}\int_{-1}^{1}\mathrm{d}\mu_{1}\int_{0}^{2\pi}\mathrm{d}\phi\,\mathcal{P}_{2}(\mu_{2}(\mu_{1},\phi))B^{r,h}(k_{1},k_{2},k_{3},\mu_{1},\mu_{2}(\mu_{1},\phi))\ ,\\ &B_{{(2,3)}}^{r,h}(k_{1},k_{2},k_{3})\equiv\frac{5}{4\pi}\int_{-1}^{1}\mathrm{d}\mu_{1}\int_{0}^{2\pi}\mathrm{d}\phi\,\mathcal{P}_{2}(\mu_{3}(\mu_{1},\phi))B^{r,h}(k_{1},k_{2},k_{3},\mu_{1},\mu_{2}(\mu_{1},\phi))\ .\end{split}$$ (11)

We work directly in this basis of quadrupoles, that are linear combinations of the $B_{2m}$ coefficients of the spherical-harmonics expansion defined in [109]. We note that if only considering the bispectrum monopole, $c_{2}^{(222)}$ and $c_{5}^{(222)}$ become degenerate, so we redefine $c_{2}^{(222)}\rightarrow c_{2}^{(222)}-c_{5}^{(222)}/6$ ().

3.2 IR-resummation

The IR-resummation is a crucial effect to include in our theory model, in order to correctly reproduce the BAO. For the power spectrum, we use the full resummation of [47] as implemented in Pybird [14]. For the bispectrum instead we rely on a wiggle-no wiggle approximation, following [111]. For the linear bispectrum, the formula we implement is:

$$B_{211}^{r,h}=2K_{1}^{r,h}(\vec{k}_{1};\hat{z})K_{1}^{r,h}(\vec{k}_{2};\hat{z})K_{2}^{r,h}(\vec{k}_{1},\vec{k}_{2};\hat{z})P_{\rm LO}(k_{1})P_{\rm LO}(k_{2})+\text{ 2 perms.}\ ,$$ (12)

where

$$P_{\rm LO}(k)=P_{\rm nw}(k)+(1+k^{2}\Sigma_{\rm tot}^{2})e^{-{k^{2}}\Sigma_{\rm tot}^{2}}P_{\rm w}(k)\,.$$ (13)

Here $P_{\rm w}(k)=P_{11}(k)-P_{\rm nw}(k)$, and $P_{\rm nw}(k)$ is the no-wiggle power spectrum, which we obtain using the sine-transform algorithm described in [112] and detailed in [113]. Then $\Sigma_{\rm tot}^{2}$ is defined by

$$\Sigma^{{2}}_{\rm tot}=-\frac{2}{15}f^{2}\,\delta\Sigma^{2}+\left(1+\frac{1}{3}f(2+f)\right)\Sigma^{2}\,,$$ (14)

$$\Sigma^{2}=\frac{4\pi}{3}\int_{0}^{\Lambda}{\frac{\mathrm{d}q}{(2\pi)^{3}}}P_{\rm nw}(q)\left[1-j_{0}(qx_{\rm osc})+2j_{2}(qx_{\rm osc})\right]\,,$$ (15)

$$\delta\Sigma^{2}=4\pi\int_{0}^{\Lambda}{\frac{\mathrm{d}q}{(2\pi)^{3}}}P_{\rm nw}(q)j_{2}(qx_{\rm osc})\,,$$ (16)

where we choose $\Lambda=1\,h\,{\rm Mpc}^{-1}$ and $x_{\rm osc}=110\,{\rm Mpc}/h$, and $j_{l}$ are the spherical Bessel functions. For the loop, our choice is to only substitute the non-integrated $P_{11}(k)$ with

$$P_{\rm NLO}(k)=P_{\rm nw}(k)+e^{-{k^{2}}\Sigma_{\rm tot}^{2}}P_{\rm w}(k)\,,$$ (17)

while for linear power spectra whose argument are being integrated, we use $P_{11}$. We discuss the goodness of this approximation in sec. 5. Another method of BAO damping for the tree-level bispectrum was given and tested on Patchy mocks in [12].

3.3 Window function

The power spectrum and bispectrum need to be convolved with the window function of the survey. For the power spectrum, this is standard and does not present numerical challenges. However, for the bispectrum this becomes more challenging. Therefore, we resort to an approximation used in [114], which amounts to evaluating the linear bispectrum with the windowed power spectrum. In formula, we have:

$$B_{211}^{r,h}=2K_{1}^{r,h}(\vec{k}_{1};\hat{z})K_{1}^{r,h}(\vec{k}_{2};\hat{z})K_{2}^{r,h}(\vec{k}_{1},\vec{k}_{2};\hat{z})[W\ast P_{11}](\vec{k}_{1})[W\ast P_{11}](\vec{k}_{2})+\text{ 2 perms.}\ ,$$ (18)

where $[W\ast P_{11}](\vec{k})=\int\frac{\mathrm{d}^{3}k^{\prime}}{(2\pi)^{3}}W(\vec{k}-\vec{k}^{\prime})P_{11}(\vec{k}^{\prime})$. Rather than projecting (18) into multipoles, we project (18) as if there was no window function, and for $[W\ast P_{11}](\vec{k})$ we use the following: for the monopole, we use $W\rightarrow W_{00}$ and, for the quadrupole we use $W\rightarrow W_{22}$, where $W_{00}$ and $W_{22}$ are defined in [4]. Because the window is a small effect, we do not apply it to the loop bispectrum. We discuss the goodness of this approximation in sec. 5.

3.4 Alcock-Paczynski effect

To estimate the galaxy spectra from data, a reference cosmology is assumed to transform the measured redshifts and celestial coordinates into three-dimensional cartesian coordinates. The difference between the reference cosmology and the true cosmology produces a geometrical distortion known as the Alcock-Paczynski (AP) effect [115]. We introduce the transverse and parallel distortion parameters:

$$q_{\perp}=\frac{D_{A}(z)H_{0}}{D_{A}^{\rm ref}(z)H_{0}^{\rm ref}}\,,\qquad q_{\parallel}=\frac{H^{\rm ref}(z)/H_{0}^{\rm ref}}{H(z)/H_{0}}\,,$$ (19)

where $D_{A}$ is the angular diameter distance, and the factors of $H_{0}$ are there since our wavenumbers are in units $h\,{\rm Mpc}^{-1}$. In terms of these, the true wavenumber and angle with the line of sight are related to the ones in the reference cosmology by:

$$k=\frac{k^{\rm ref}}{q_{\perp}}\left[1+(\mu^{\rm ref})^{2}\left(\frac{1}{F^{2}}-1\right)\right]^{1/2}\,,\quad\mu=\frac{\mu^{\rm ref}}{F}\left[1+(\mu^{\rm ref})^{2}\left(\frac{1}{F^{2}}-1\right)\right]^{-1/2}\,,$$ (20)

where $F=q_{\parallel}/q_{\perp}$. To match the measured power spectrum multipoles, we do the following integral:

$$P_{\ell}(k^{\rm ref})=\frac{2\ell+1}{2q_{\parallel}q_{\perp}^{2}}\int_{-1}^{1}\mathrm{d}\mu^{\rm ref}\mathcal{L}_{\ell}(\mu^{\rm ref})P(k(k^{\rm ref},\mu^{\rm ref}),\mu(\mu^{\rm ref}))\,.$$ (21)

The formula for the bispectrum is:

$$B_{(\ell,i)}(k_{1}^{\rm ref},k_{2}^{\rm ref},k_{3}^{\rm ref})=\frac{{2\ell+1}}{2q_{\parallel}^{2}q_{\perp}^{4}}\int_{-1}^{1}\mathrm{d}\mu_{1}^{\rm ref}\int_{0}^{2\pi}\frac{\mathrm{d}\phi^{\rm ref}}{2\pi}B(k_{1},k_{2},k_{3},\mu_{1},\mu_{2},\mu_{3})\mathcal{P}_{\ell}(\mu_{i})\,.$$ (22)

For the bispectrum, we only apply the Alcock-Paczynski effect on the tree-level part, as it is a small effect: we find that, within BOSS error bars (on $\Omega_{m}$), it is an effect of at most $\sim 1\sigma$, and accordingly, the change in $\chi^{2}$ is at most $1$ if neglecting it completely. Given the size of the loop terms and counterterms, it is thus safe to neglect it there. We find that we can achieve sufficient numerical accuracy using a nested trapezoidal rule with only $13$ points in $\mu$ and $4$ points in $\phi$, after using symmetries to restrict the integration domain to $\mu_{1}\in[0,1]$, and $\phi\in[-\pi/2,\pi/2]$.

3.5 Binning

For the power spectrum, data are an average over spherical shells in momentum space. The theoretical prediction needs therefore to be averaged over the fundamental modes of the chosen grid. Since our bins have many fundamental modes, in practice we do an integral of the power spectrum over a bin, which is numerically very simple. The effect of binning is anyway small for the power spectrum, with respect to the error bars of our data and simulations.

For the bispectrum, we have an average over fundamental (closed) triangles in a bin of width $\Delta k$ around a central triangle with sides $k_{1}$, $k_{2}$, $k_{3}$. Especially for our chosen bins with $\Delta k\simeq 0.02$, it is important to take into account the binning effects when comparing the theory to the data. The average should be done as a sum over fundamental triangles:

$$B^{r,h}_{{(\ell,i)},\rm bin}(k_{1},k_{2},k_{3})=\frac{2\ell+1}{N_{T}}\sum_{\vec{q}_{1}\in k_{1}}\sum_{\vec{q}_{2}\in k_{2}}\sum_{\vec{q}_{2}\in k_{2}}\delta_{K}(\vec{q}_{1}+\vec{q}_{2}+\vec{q}_{3})B^{r,h}(\vec{q}_{1},\vec{q}_{2},\vec{q}_{3})\mathcal{P}_{\ell}(\mu_{i})\,,$$ (23)

and we note that, here and elsewhere, $B^{r,h}(\vec{q}_{1},\vec{q}_{2},\vec{q}_{3})$ is the full redshift-space bispectrum, i.e. we have suppressed the dependence on $\hat{z}$ for notational convenience. Here $N_{T}$ is the number of fundamental triangles in the bin, $\delta_{K}$ is the Kronecker delta function, and the notation $\vec{q}_{i}\in k_{i}$ means a sum over the fundamental modes $\vec{q}_{i}$ for which $k_{i}-\frac{\Delta k}{2}\leq|\vec{q}_{i}|<k_{i}+\frac{\Delta k}{2}$. Calculating such a sum is numerically very challenging. However, since in each bin there are many fundamental triangles, we expect that an integral approximation should work well. The only caveat is that one needs to integrate only over the closed triangles. In particular, this is very important for bins such that $k_{3}+\Delta k/2>k_{1}+k_{2}-\Delta k$ (remember that our ordering is $k_{1}\leq k_{2}\leq k_{3}$), for which there are configurations of modes that do not form a closed triangle in the bin.

Therefore, we implement the following formula:

$$\begin{split}B^{r,h}_{{(\ell,i)},\rm bin}(k_{1},k_{2},k_{3})&=\frac{{2\ell+1}}{V_{T}}\left(\prod_{i=1}^{3}\int_{V_{i}}\frac{\mathrm{d}^{3}q_{i}}{(2\pi)^{3}}\right)(2\pi)^{3}\delta_{D}^{(3)}(\vec{q}_{1}+\vec{q}_{2}+\vec{q}_{3})B^{r,h}(\vec{q}_{1},\vec{q}_{2},\vec{q}_{3})\mathcal{P_{\ell}}(\mu_{i})\,,\end{split}$$ (24)

where

$$V_{T}\equiv\left(\prod_{i=1}^{3}\int_{V_{i}}\frac{\mathrm{d}^{3}q_{i}}{(2\pi)^{3}}\right)(2\pi)^{3}\delta_{D}^{(3)}(\vec{q}_{1}+\vec{q}_{2}+\vec{q}_{3})\,,$$ (25)

and we used the notation

$$\int_{V_{i}}\frac{\mathrm{d}^{3}q_{i}}{(2\pi)^{3}}=\int_{k_{i}}\frac{\mathrm{d}q_{i}}{2\pi^{2}}q_{i}^{2}\int\frac{\mathrm{d}^{2}\hat{q}_{i}}{4\pi}\ ,\quad\text{where}\quad\int_{k_{i}}\equiv\int_{k_{i}-\frac{\Delta k}{2}}^{k_{i}+\frac{\Delta k}{2}}\mathrm{d}q_{i}\ .$$ (26)

As shown in app. B, we can perform the angular integrals and find

$$B^{r,h}_{{(\ell,i)},\rm bin}(k_{1},k_{2},k_{3})=\frac{1}{V_{T}}\int_{k_{1}}\mathrm{d}q_{1}\int_{k_{2}}\mathrm{d}q_{2}\int_{k_{3}}\mathrm{d}q_{3}\,q_{1}q_{2}q_{3}\,\frac{\beta\left(\Delta_{q}\right)}{8\pi^{4}}B^{r,h}_{{(\ell,i)}}(q_{1},q_{2},q_{3})\,\,,$$ (27)

$$V_{T}=\int_{k_{1}}\mathrm{d}q_{1}\int_{k_{2}}\mathrm{d}q_{2}\int_{k_{3}}\mathrm{d}q_{3}\,q_{1}q_{2}q_{3}\,\frac{\beta\left(\Delta_{q}\right)}{8\pi^{4}}\,,$$ (28)

where $\beta(\Delta_{q})=1/2$ if $q_{1}$, $q_{2}$, $q_{3}$ form a folded triangle, $\beta(\Delta_{q})=1$ for all other (closed) triangles, and $\beta(\Delta_{q})=0$ otherwise.

We apply only the binning in this way to the tree-level part. For efficient numerical evaluation of the integrals in eq. (27), we implement the bispectrum binning as follow. For a given bin centered in $(k_{1},k_{2},k_{3})$, we enforce that $(q_{1},q_{2},q_{3})$ forms a triangle by redefining the integration boundaries: $q_{1}\in[k_{1}-\Delta k/2,k_{1}+\Delta k/2]$, $q_{2}\in[k_{2}-\Delta k/2,k_{2}+\Delta k/2]$, and $q_{3}\in[|k_{1}-k_{2}|,k_{1}+k_{2}]$. Whenever $q_{3}$ can not satisfy this triangle inequality, we drop this configuration. As such, we can drop the $\beta(\Delta_{q})$ function inside the integral. We perform a change of variable $q_{3}\rightarrow\cos(\theta_{12})\equiv(q_{3}^{2}-q_{1}^{2}-q_{2}^{2})/(2q_{1}q_{2})$, such that the integral measure becomes $q_{1}q_{2}q_{3}\ \mathrm{d}q_{1}\mathrm{d}q_{2}\mathrm{d}q_{3}\rightarrow q_{1}^{2}q_{2}^{2}\ \mathrm{d}q_{1}\mathrm{d}q_{2}\mathrm{d}\cos(\theta_{12})$, We then discretize the integration domain in 6 evenly-spaced points in $q_{1}$, 6 in $q_{2}$, and 4 in $\cos(\theta_{12})$. On each point of this grid, we evaluate the $14$ pieces of the tree-level part of the bispectrum. The binning integrals are then performed with a nested trapezoidal rule over the grid. Given that the AP integrals, eq. (22), need to be performed for each of those evaluations, we limit the number of evaluations by first looking for (and storing) the common triangles of the discretized domains over all the bins we need to evaluate. For our $150$ bins in CMASS, this reduces the total number of evaluations by about a factor $1.5$, from $6\cdot 6\cdot 4\cdot 150=21600$ to $13782$. After compilation of the integrand expressions going in the AP integrals, we are able to evaluate the binned bispectrum in our Python code with an overall runtime of $\sim 1$ second on $1$ CPU. The numerical precision of such evaluation has been extensively tested, in particular against Monte-Carlo integrations, and is found to be under control for the data and simulations error bars we analyze in this work.

The loop pieces and counterterms, that are small with respect to the linear term, are instead evaluated on effective wavenumbers. They are defined, as described in [116], by the following averages:

$$k_{\rm eff,1}=\frac{1}{V_{T}}\int_{k_{1}}\frac{\mathrm{d}q_{1}}{2\pi}\int_{k_{2}}\frac{\mathrm{d}q_{2}}{2\pi}\int_{k_{3}}\frac{\mathrm{d}q_{3}}{2\pi}\,q_{1}q_{2}q_{3}\,\beta\left(\Delta_{q}\right)\min(q_{1},q_{2},q_{3})\,,$$ (29)

$$k_{\rm eff,2}=\frac{1}{V_{T}}\int_{k_{1}}\frac{\mathrm{d}q_{1}}{2\pi}\int_{k_{2}}\frac{\mathrm{d}q_{2}}{2\pi}\int_{k_{3}}\frac{\mathrm{d}q_{3}}{2\pi}\,q_{1}q_{2}q_{3}\,\beta\left(\Delta_{q}\right)\mathrm{med}(q_{1},q_{2},q_{3})\,,$$ (30)

$$k_{\rm eff,3}=\frac{1}{V_{T}}\int_{k_{1}}\frac{\mathrm{d}q_{1}}{2\pi}\int_{k_{2}}\frac{\mathrm{d}q_{2}}{2\pi}\int_{k_{3}}\frac{\mathrm{d}q_{3}}{2\pi}\,q_{1}q_{2}q_{3}\,\beta\left(\Delta_{q}\right)\max(q_{1},q_{2},q_{3})\,.$$ (31)

As expected from the size of those terms and the size of the binning effect ($\sim 1\sigma$), we have checked that properly binning the loop instead of evaluating them on these effective wavenumbers lead to negligible shift in the $\min\chi^{2}$ and in the posteriors for the analyses presented in this work.

4 Likelihood

To analyze the data, we start from a Gaussian likelihood, which is multiplied by the prior to arrive at the Bayesian posterior $\mathcal{P}$ over cosmological and bias parameters:

$$-2\ln\mathcal{P}=(T_{i}-D_{i})C^{-1}_{ij}(T_{j}-D_{j})-2\ln\mathcal{P}_{\rm pr}\,,$$ (32)

where $T_{i}$ is the full vector of theory predictions in bin $i$, containing power spectrum multipoles and bispectra, $D_{i}$ the corresponding data measurement in bin $i$, $C_{ij}$ is the full covariance between bins $i$ and $j$, and $\mathcal{P}_{\rm pr}$ is a generic prior on the parameters.

Our theory model depends on cosmological and EFT parameters. It is the case that many EFT parameters appear linearly in the theory model. Denoting them by $g_{\alpha}$, we can write

$$T_{i}=g_{\alpha}T_{G,i}^{\alpha}+T_{NG,i}\ ,$$ (33)

where $T_{G,i}^{\alpha}$ and $T_{NG,i}$ depend non-linearly (that is, at least quadratically) on the other cosmological parameters and three biases for each sky cut. Since we are interested in the marginalized posteriors over cosmological parameters, it is very convenient to do the analytical Gaussian integration over the $g_{\alpha}$. We will also choose a Gaussian prior on them, with covariance $\sigma_{\alpha\beta}$ and mean $\hat{g}_{\alpha}$. Collecting the powers of $g_{\alpha}$, we can write the posterior in the following form:

$$-2\ln\mathcal{P}=g_{\alpha}F_{2,\alpha\beta}g_{\beta}-2g_{\alpha}F_{1,\alpha}+F_{0}\,,$$ (34)

where the $F$’s are defined as:

$$F_{2,\alpha\beta}=T_{G,i}^{\alpha}C^{-1}_{ij}T_{G,j}^{\beta}+\sigma^{-1}_{\alpha\beta}\,,$$ (35)

$$F_{1,\alpha}=-T_{G,i}^{\alpha}C^{-1}_{ij}(T_{NG,j}-D_{j})+\sigma^{-1}_{\alpha\beta}\hat{g}_{\beta}\,,$$ (36)

$$F_{0}=(T_{NG,i}-D_{i})C^{-1}_{ij}(T_{NG,j}-D_{j})+\hat{g}_{\alpha}\sigma^{-1}_{\alpha\beta}\hat{g}_{\beta}-2\ln\Pi\,,$$ (37)

where $\Pi$ is a generic prior on the cosmological and bias parameters non analytically marginalized. In other words, we assume that $\mathcal{P}_{\rm pr}$ is a sum of a Gaussian prior over the $g_{\alpha}$ and a remaining prior on the other parameters. After integrating the $g_{\alpha}$, we have the marginalized posterior:

$$-2\ln\mathcal{P}_{\rm marg}=-F_{1,\alpha}F_{2,\alpha\beta}^{-1}F_{1,\beta}+F_{0}+\ln\det\left(\frac{F_{2}}{2\pi}\right)\,.$$ (38)

Prior.

In our analysis, we vary the cosmological parameters $\omega_{cdm}$, $h$, and $\ln\left(10^{10}A_{s}\right)$ with a flat uninformative prior, while we use a Gaussian prior on $\omega_{b}$ of mean $\omega_{b,\textrm{BBN}}=0.02233$ and standard deviation $\sigma_{\rm BBN}=0.00036$, motivated from Big-Bang Nucleosynthesis (BBN) experiments [117]. We instead fix $n_{s}$ to the truth of the simulations or to the Planck preferred value when analyzing the data [15]. When analyzing the BOSS data, we also fix the neutrino to minimal mass following Planck prescription.

The EFT parameters should instead be restricted to be ${\cal{O}}(1)$ numbers, for consistency of the perturbative expansion. The EFTofLSS is an expansion in the size of fluctuations and derivatives. Both of these are suppressed by a nonlinear scale $k_{\rm NL}\simeq k_{\rm M}\simeq 0.7h\,{\rm Mpc}^{-1}$, where $k_{\rm NL}$ is the nonlinear scale for the matter field, and $k_{\rm M}$ is the typical wavenumber associated to galaxy size. However, it was recognized in [77, 118] that terms involving expectation values of velocity fields, coming from the transformation to redshift space, define a new scale, which we denote by $k_{\rm NL,R}\simeq k_{\rm NL}/\sqrt{8}$. We therefore write down each operator in the EFT expansion with either a $k_{\rm M}$ or a $k_{\rm NL,R}$ suppression, depending on its origin. We then use a Gaussian prior of width 2 centered on 0 on all the EFT parameters that we analytically marginalize, with the following exception: on $c_{h,1},c_{\pi,1},c_{\pi v,1}$, and $c_{2}^{\rm St}$, that already appear in the power spectrum (see their definitions in app. A), we put instead a Gaussian prior of width 4 centered on 0, such that the prior is the same as the ones used in our previous series of analyses with the power spectrum only (see e.g. [4, 6, 7]). For the quadratic biases, we define the linear combinations $c_{2}=(b_{2}+b_{5})/\sqrt{2}$, $c_{4}=(b_{2}-b_{5})/\sqrt{2}$, and we assign on them a Gaussian prior of width 2 centered on 0. Finally, for $b_{1}$, which is positive definite, we use a lognormal prior of mean $0.8$ (since $e^{0.8}=2.23$), and variance 0.8, such that $[0,3.4]$ is the $68\%$ bound for this prior on $b_{1}$. For definiteness, in our prior, we take $k_{\rm NL}=k_{\rm M}=0.7h\,{\rm Mpc}^{-1}$ and $\bar{n}=4\cdot 10^{-4}({\rm Mpc}/h)^{3}$.

When analyzing more than one sky, we can use the information that the bias and EFT parameters should be the same at the same redshift, and their time evolution is expected to be comparable to the growth factor to some small power. This allows us to estimate the variation of $b_{1}$ between CMASS or LOWZ effective redshifts, to be about 20%. Therefore, in our multisky analyses the biases are correlated, which, as explained in the following section, helps to mitigate prior volume effects. In practice, let us consider the 4-sky analysis and the $b_{1}$ parameters, which will be a vector $(b_{1}^{(1)},b_{1}^{(2)},b_{1}^{(3)},b_{1}^{(4)})$, with one $b_{1}^{(i)}$ for each sky. The prior on it is a multivariate lognormal with correlation matrix:

$$\begin{pmatrix}1&\rho_{12}&\rho_{13}&\rho_{12}\rho_{13}\\ \rho_{12}&1&\rho_{12}\rho_{13}&\rho_{13}\\ \rho_{13}&\rho_{12}\rho_{13}&1&\rho_{12}\\ \rho_{12}\rho_{13}&\rho_{13}&\rho_{12}&1\end{pmatrix}\,,$$ (39)

where $\rho_{ij}=1-\epsilon_{ij}^{2}/2$, and we choose $\epsilon_{12}=0.1$, $\epsilon_{13}=0.2$. This formula is motivated by the fact that two variables distributed according to a bivariate normal with correlation $\rho$, the standard deviation of the difference is $\epsilon=\sqrt{2(1-\rho)}$. Our choices of $\epsilon_{ij}$ then reflect that we expect the values of $b_{1}$ to be different only by about $10\%$ between NGC and SGC, given slightly different selection function, and only by about $20\%$ between CMASS and LOWZ, given the redshift evolution of $b_{1}$. We use the same correlation matrix for the Gaussian priors on all the quadruplets $c_{2}$, $c_{4}$ and the $g_{\alpha}$ parameters.

Posterior sampling.

Our analyses are performed using the Metropolis-Hastings sampler as implemented in MontePython 3 [119], with the theory model evaluated using CLASS [120] and PyBird. We declare our MCMC converged when the Gelman-Rubin criterion [121] is $\leq 0.02$. The plots and summary statistics are calculated with the GetDist [122] package.

5 Pipeline validation

For our analyses, we use the following scale cut: $k_{\rm min}=0.01h\,{\rm Mpc}^{-1}$ for all observables, $k_{\rm max}=0.23h\,{\rm Mpc}^{-1}$ for $P_{\ell}$ and $B_{0}$, and $k_{\rm max}=0.08h\,{\rm Mpc}^{-1}$ for $B_{2}$, on CMASS. For LOWZ instead, we use $k_{\rm max}=0.20h\,{\rm Mpc}^{-1}$ for $P_{\ell}$ and $B_{0}$ , following [4]. We keep $k_{\rm max}=0.08h\,{\rm Mpc}^{-1}$ for $B_{2}$ on LOWZ. In this section, we perform multiple checks to validate our method at this scale cut.

5.1 Measuring and fixing phase-space effects

$\sigma_{\rm proj}/\sigma_{\rm stat}$	$\Omega_{m}$	$h$	$\sigma_{8}$	$\omega_{cdm}$	$\ln(10^{10}A_{s})$	$S_{8}$
1 sky, $\sim 100\,V_{\rm 1sky}$	-0.1	-0.14	-0.21	-0.2	-0.07	-0.23
1 sky, $V_{\rm 1sky}$, adjust.	0.13	0.06	0.04	0.15	-0.04	0.08
4 skies, $V_{\rm 4skies}$, adjust.	0.1	0.	-0.05	0.07	-0.06	-0.01

Figure 3: Triangle plots of base cosmological parameters obtained fitting synthetic data analyzed using a covariance with BOSS volume $V_{\rm BOSS}$ or rescaled to a large volume $\sim 100V_{\rm BOSS}$, with prior on the EFT parameters centered on their truth, or with phase-space projection adjustment. Here the synthetic data are corresponding exactly to our model $P_{\ell}+B_{0}+B_{2}$ on the best fit of patchy. ‘1 sky’ or ‘4 skies’ correspond to CMASS NGC or all BOSS skycuts, respectively. The grey lines in the triangle plots represent the truth. We also show the relative deviations $\sigma_{\rm proj}/\sigma_{\rm stat}$ on the base cosmological parameters from the truth from those various analyses. In summary, the addition of a phase-space correction prior to our likelihood allows us to recover unbiased mean in the 1D posteriors of the cosmological parameters of interest.

$\sigma_{\rm proj}/\sigma_{\rm stat}^{\rm data}$		$\Omega_{m}$	$h$	$\sigma_{8}$	$\omega_{cdm}$	$\ln(10^{10}A_{s})$	$S_{8}$
Nseries	$P_{\ell}$	-0.02	0.05	0.08	0.02	0.05	0.07
	$P_{\ell}+B_{0}$	-0.06	-0.03	-0.04	-0.08	0.03	-0.06
	$P_{\ell}+B_{0}+B_{2}$	-0.12	-0.	-0.04	-0.11	0.04	-0.08
1 sky	$P_{\ell}$	-0.15	0.07	-0.11	-0.06	-0.08	-0.15
	$P_{\ell}+B_{0}$	0.07	0.06	0.09	0.11	0.02	0.1
	$P_{\ell}+B_{0}+B_{2}$	0.13	0.06	0.04	0.15	-0.04	0.08
4 skies	$P_{\ell}$	-0.01	0.05	-0.03	0.02	-0.04	-0.03
	$P_{\ell}+B_{0}$	0.05	-0.	0.01	0.03	0.01	0.03
	$P_{\ell}+B_{0}+B_{2}$	0.1	0.	-0.05	0.07	-0.06	-0.01

Table 1: Residual deviations $\sigma_{\rm proj}$ after phase-space projection adjustment measured on synthetic data generated and fitted with our model $P_{\ell}+B_{0}+B_{2}$ with truth given by the best fits of Nseries, Patchy 1 sky, or Patchy 4 skies, relative to BOSS error bars $\sigma_{\rm stat}^{\rm data}$.

Our likelihood has several EFT parameters on top of the cosmological parameters. Some of these appear in the likelihood in a Gaussian way, and we analytically marginalize over them. Performing such a Gaussian integral corresponds to putting these parameters to their best fit values, given all the other parameters and observational data. At this point, we are left with a likelihood which has a non-Gaussian dependence on the EFT and the cosmological parameters.

Now, there is an interesting phenomenon that we would like to describe. Let us analyze data that are generated with our theory model: the EFTofLSS plus observational effects as described in sec 3. We refer to these as ‘synthetic’ data. We generate these synthetic data by choosing the best fit EFT parameters that we find by fitting the average of $2048$ Patchy simulations, on the Patchy cosmology, so that the resulting EFT and cosmological parameters are at realistic values. In this case, the best fit has $\chi^{2}=0$ once we put flat priors on the EFT parameters, and we should clearly recover the correct cosmological parameters. However, as it can be seen from fig. 3, in green, the sampled posteriors show biases in all 1D posteriors of the cosmological parameters, and in particular in $\sigma_{8}$ and $\Omega_{m}$. What is going on?

The first hypothesis is that there could be an error in our pipeline. This hypothesis can be discarded by noticing that if we analyze the data with a covariance that is about 100 times smaller, we recover the truth with exquisite precision (see the blue curve in fig. 3). So, we exclude this hypothesis.

Another reason for the offset of the green curve in fig. 3 could be the prior on the EFT parameters. In fact, while on the synthetic data the EFT parameters have some definite values (which are well within the priors), our Gaussian priors are centered at zero, and so the true value of the EFT parameters are slightly disfavored by the priors. We check if this can be the reason to the offset seen in the posteriors of the cosmological parameters by sampling instead with priors centered around the synthetic truth. We find that the resulting posteriors are close to previous results (grey vs. green in fig. 3), suggesting that the central value of the prior of the bias parameters does not play a substantial role. This means that even if the truth is the maximum likelihood point, the posteriors will not recover it.

Having excluded that the bias in the posteriors on synthetic data is due to an error in our pipeline or due to our priors, we conclude that it must be due to phase-space projection effects. In fact, if the posteriors of the EFT parameters are effectively non-Gaussian (i.e. if the error bars are sufficiently large that the Taylor expansion at second order around the maximum of the posterior is not accurate enough to describe the actual posterior), then, upon marginalization, one can get projection effects on the remaining parameters, which, in this case, are the cosmological ones, even if the maximum likelihood point is the truth. Given the large number of EFT parameters, it is not so surprising that this might be the case. We call this effect ‘phase space effect,’ but it is also known as ‘prior volume effect’ or ‘projection effect.’

We decide to fix this issue with the following procedure. As a measurement of the phase-space effect, for all analyses in this work, we take the shift in the 1D posteriors from the truth obtained fitting synthetic data with the same modeling and covariance. We add a prior of the following form to the $\log$-likelihood of $P_{\ell}+B_{0}(+B_{2})$:

$$\ln\mathcal{P}_{\rm pr}^{\rm ph.\;sp.\;1sky}=-18\left(\frac{b_{1}}{2}\right)+8\left(\frac{\Omega_{m}}{0.31}\right)+14\left(\frac{h}{0.68}\right)\ ,$$ (40)

$$\ln\mathcal{P}_{\rm pr}^{\rm ph.\;sp.\;4sky}=-48\left(\frac{b_{1}}{2}\right)+32\left(\frac{\Omega_{m}}{0.31}\right)+48\left(\frac{h}{0.68}\right)\ ,$$

respectively for 1 sky and 4 skies. As such, upon marginalization, we recover unbiased 1D posteriors from the fit to the synthetic data (see the red curve in fig. 3 and also the associated table). More in detail, in tab. 1, we show the residual deviation from phase-space projection on the base cosmological parameters measured from synthetic data. We see that for all data volume (either the one of CMASS NGC or of all BOSS 4 skies) and cosmologies tested here (either the one of Nseries or the one of Patchy), we find that the residual deviation are negligibly small ($\lesssim 0.15$ or the error bars obtained with BOSS-volume covariance). Since the synthetic data are close to the Patchy ones (and so to the data), and since we expect the phase-space projection effects to be a slowly-varying function of the cosmological and EFT parameters, we add the same phase-space-correcting prior to the likelihood of the BOSS data.

$\Delta_{\rm sys}/\sigma_{\rm stat}$	$\Omega_{m}$	$h$	$\sigma_{8}$	$\omega_{cdm}$	$\ln(10^{10}A_{s})$	$S_{8}$
$P_{\ell}+B_{0}$: base - w/ NNLO	-0.03	-0.09	-0.03	-0.1	0.05	-0.04
$P_{\ell}+B_{0}$: base - w/o $B_{0}$ window	0.11	-0.05	0.01	0.05	-0.01	0.05
$P_{\ell}+B_{0}+B_{2}$: base - w/o $B_{0},B_{2}$ window	0.51	0.09	0.02	0.51	-0.25	0.19

Table 2: Relative shifts $\Delta_{\rm sys}/\sigma_{\rm stat}$ on base cosmological parameters measured from various modeling choices compared to our baseline: inclusion of the NNLO or removal of the window function in the bispectrum.

5.2 Scale cut from NNLO

A simulation-independent way to evaluate the theoretical error as a function of $k_{\rm max}$ is to analyze the data by adding to the theory model a part of the next order terms: for our one-loop model, this part consists in the next-to-next-to-leading-order (NNLO) terms. Such a procedure was successfully applied to estimate the scale cut for the CF [7]. Here we use the same technique. We add the following two-loop counterterms to the EFTofLSS prediction at one-loop for the power spectrum:

$$P_{\rm NNLO}(k,\mu)=\frac{1}{4}c_{r,4}b_{1}^{2}\mu^{4}\frac{k^{4}}{k_{\rm NL,R}^{4}}P_{11}(k)+\frac{1}{4}c_{r,6}b_{1}\mu^{6}\frac{k^{4}}{k_{\rm NL,R}^{4}}P_{11}(k)\,,$$ (42)

and for the bispectrum:

$$B_{\rm NNLO}(k_{1},k_{2},k_{3},\mu,\phi)=2c_{\rm NNLO,1}K_{2}^{r,h}(\vec{k}_{1},\vec{k}_{2};\hat{z})K_{1}^{r,h}(\vec{k}_{2};\hat{z})f\mu_{1}^{2}\frac{k_{1}^{4}}{k_{\rm NL,R}^{4}}P_{11}(k_{1})P_{11}(k_{2})$$

$$+c_{\rm NNLO,2}K_{1}^{r,h}(\vec{k}_{1};\hat{z})K_{1}^{r,h}(\vec{k}_{2};\hat{z})P_{11}(k_{1})P_{11}(k_{2})f\mu_{3}k_{3}\frac{(k_{1}^{2}+k_{2}^{2})}{4k_{1}^{2}k_{2}^{2}k_{\rm NL,R}^{4}}\Big{[}-2\vec{k}_{1}\cdot\vec{k}_{2}(k_{1}^{3}\mu_{1}+k_{2}^{3}\mu_{2})$$

$$\quad+2f\mu_{1}\mu_{2}\mu_{3}k_{1}k_{2}k_{3}(k_{1}^{2}+k_{2}^{2})\Big{]}+\textrm{perm.}\,,$$ (43)

where $k_{\rm NL,R}=k_{\rm NL}/\sqrt{8}$, as discussed in sec. 4. The prefactors $c_{r,4},c_{r,6},c_{\rm NNLO,1}$, and $c_{\rm NNLO,2}$ are given a Gaussian prior centered on zero and of width 2. We then analyze the data as a function of $k_{\rm max}$, and determine the maximum wavenumber by taking the largest $k_{\rm max}$ where the shift in all 1D posteriors of the cosmological parameters with respect to the analysis without these terms is equal to $1/3\cdot\sigma$. This would mean that our results would have become sensitive to these terms that we do not fully compute, and so we need to analyze the data only up to this threshold. For simplicity, rather than determining the $k_{\rm max}$ in this way, we check the effect of these NNLO terms close to the $k_{\rm max}$ that we find in simulations, and check that the effect of the NNLO terms is indeed not too large. The results are presented in tab. 2. We see that the effect is negligibly small, confirming what we find in simulations next, i.e. that our scale cut is appropriate.

$\Delta X/X$		$\Omega_{m}$	$h$	$\sigma_{8}$	$\omega_{cdm}$	$\ln(10^{10}A_{s})$	$S_{8}$
Nseries	$P_{\ell}$	$-0.017\pm 0.048$	$0.003^{+0.022}_{-0.024}$	$0.047^{+0.070}_{-0.086}$	$-0.013^{+0.063}_{-0.071}$	$0.035\pm 0.055$	$0.038^{+0.074}_{-0.092}$
	$P_{\ell}+B_{0}$	$-0.005\pm 0.042$	$0.005\pm 0.019$	$-0.012\pm 0.052$	$0.004^{+0.052}_{-0.058}$	$-0.010\pm 0.040$	$-0.015\pm 0.058$
	$P_{\ell}+B_{0}+B_{2}$	$-0.010\pm 0.041$	$0.006^{+0.018}_{-0.021}$	$-0.009\pm 0.053$	$0.001\pm 0.053$	$-0.007\pm 0.041$	$-0.014\pm 0.059$
Patchy	$P_{\ell}$	$-0.021^{+0.052}_{-0.059}$	$0.005\pm 0.028$	$0.034^{+0.076}_{-0.098}$	$-0.014\pm 0.078$	$0.026^{+0.056}_{-0.066}$	$0.022^{+0.083}_{-0.10}$
	$P_{\ell}+B_{0}$	$-0.011^{+0.044}_{-0.051}$	$0.004\pm 0.023$	$-0.011\pm 0.054$	$-0.004^{+0.052}_{-0.062}$	$-0.006\pm 0.044$	$-0.017\pm 0.058$
	$P_{\ell}+B_{0}+B_{2}$	$-0.012\pm 0.046$	$0.004\pm 0.024$	$-0.004\pm 0.053$	$-0.006^{+0.052}_{-0.059}$	$-0.001\pm 0.043$	$-0.011\pm 0.058$

Figure 4: Triangle plots and relative $68\%$-credible intervals of base cosmological parameters measured from the Nseries and Patchy simulations analyzed using a covariance with CMASS NGC volume. The grey lines in the triangle plots represent the simulation truth.

5.3 Tests of additional modeling effects

Our implementation of the IR-resummation and of the window function is approximate, without a control parameter. We therefore check the accuracy of the two implementations in the following way.

For the window function, the correctness of our approximation has been checked in [123] for the monopole. In fact, as shown in the second line of tab. 2, the difference between the bispectrum computed with our approximation, and the one where we apply no window is within $1/4$ of the error bars obtained on all cosmological parameters from the fit to BOSS data. For the quadrupole, the third line of tab. 2 shows that the difference with applying no window is about 0.5$\sigma$ on the posterior of $\Omega_{m}$ (while negligible for the other cosmological parameters). While this might seem too large an effect to tolerate, one should keep in mind the following. Roughly speaking, the correct window function should consist of applying $3/2$ factors of $W$ to the bispectrum (i.e. one for each field). Applying no window therefore is a radical negligence of all these factors, much worse than the approximation we do (which applies two factors of $W$). We therefore believe that a more reliable estimate of the error associated to our implementation of the window function for the quadrupole is obtained by dividing the effect in tab. 2 by a factor of 4. Even if our estimate were to be wrong by a factor 2, this would make the effect safely negligible. It would be interesting to compare our approach to an analysis using another estimator based on tri-polar spherical harmonics (described in [10] and tested on Patchy mocks in [12]) for which the window functions can be estimated on an equal footing, making its application more straightforward.

Let us now discuss the goodness of our approximate implementation of the IR-resummation of the bispectrum. It should be emphasized that the wiggle/no-wiggle procedure is affected by several uncontrolled approximations (i.e. not controlled by a small parameter, but numerically accidentally small) [62]. On top of those, our formulas neglect the angle dependence of the IR-resummation, and, perhaps even more quantitatively importantly, do not damp the oscillations in the power spectra whose momenta are integrated in the loop integrals, as for example proposed in [111]. We checked that applying the damping for those power spectra leads to a negligible ($\lesssim 0.25$) change in the $\chi^{2}$ when keeping all the parameters of the model fixed. We therefore conclude that neglecting the IR-resummation on the ‘wiggly’ parts from inside the loop integrals is accurate enough for BOSS data. We leave to future work more careful inspection of the remaining approximations in our IR-resummation scheme.

5.4 Tests against simulations

$\sigma_{\rm sys}^{\rm sim}/\sigma_{\rm stat}^{\rm data}$	$\Omega_{m}$	$h$	$\sigma_{8}$	$\omega_{cdm}$	$\ln(10^{10}A_{s})$	$S_{8}$
Nseries $P_{\ell}+B_{0}$	0.02	0.17	0.15	-0.03	0.17	0.17
Nseries $P_{\ell}+B_{0}+B_{2}$	0.16	0.25	0.08	-0.09	0.08	0.16
Patchy $P_{\ell}+B_{0}$	0.27	0.21	0.23	0.05	0.14	0.33
Patchy $P_{\ell}+B_{0}+B_{2}$	0.31	0.2	0.07	0.09	0	0.2

Table 3: Report of systematic errors on base cosmological parameters measured from the Nseries and Patchy simulations. The systematic error, reported relative to the BOSS error bars $\sigma_{\rm stat}^{\rm data}$, is defined as $\sigma_{\rm sys}^{\rm sim}\equiv\max(|\text{mean}-\textrm{truth}|-\sigma_{\rm stat}^{\rm sim}/\sqrt{N_{\rm sim}},0)$. Here $\sigma_{\rm stat}^{\rm sim}/\sqrt{N_{\rm sim}}$ represents the uncertainty from the simulation cosmic variance, which corresponds to about $0.15$ or $0.03$ in $\sigma_{\rm stat}^{\rm data}$ for $N_{\rm sim}=84$ Nseries or $N_{\rm sim}=2048$ Patchy realizations, respectively.

We now test the accuracy of the model by comparing against $N$-body simulations described in sec. 2. This does not only test the effect of the theoretical error due to the next order terms not included in our baseline model or to the approximate IR-resummation, but also of the other observational effects that we model imperfectly, such as the window function. In fig. 4, we show the posteriors from the analysis of the average of 84 Nseries boxes, analyzed with the covariance of one box, such that we also account for the phase space effect. Since the actual cosmic variance associated to this average of 84 boxes is about 1/9 of the posteriors in fig. 4, we measure for each cosmological parameter the theoretical error as the distance of the mean of the posterior to the truth of the simulation minus 1/9 of the standard deviation (we take zero if this number is negative). This allows us to detect theoretical errors larger than 1/9 of a standard deviation of the posterior in fig. 4, which corresponds to about 0.15 of the of the error bars obtained on the BOSS data. Our results show that the theoretical error that we can detect is safely below 1/3 of the error bars obtained on BOSS, as summarized in tab. 3.

In fig. 4, we also present the analogous analysis on the average of 2048 Patchy mocks. In this case, the detectable theoretical error is almost unaffected by cosmic variance. Thus, assuming no systematic error in the Patchy simulations, the minimal detectable theoretical error is practically zero. Also in this case, the theoretical error is safely below 1/3 the error bars obtained on BOSS, as summarized in tab. 3.

We conclude that our analysis pipeline is free from significant systematics for BOSS volume at the scale cuts chosen at the beginning of this section, and we now move on to the analysis of the observational data.

6 Results

	$N_{\rm bin}$ // dof	$\min\chi^{2}$	$\min\chi^{2}/$dof	$p$-value
CMASS NGC	$42+150+9=201$	159.5	0.79	0.99
CMASS SGC	$42+150+9=201$	188.7	0.94	0.72
LOWZ NGC	$36+62+9=107$	98.3	0.92	0.71
LOWZ SGC	$36+62+9=107$	106.4	0.99	0.50
Parameter Prior	$3+41(1+0.1+0.2+0.1\cdot 0.2)\simeq 57$	8.9	-	-
Total	$616-57=559$	$561.9$	1.01	0.46

Table 4: Goodness of fit given by the maximal $\log$-likelihood value $\log\mathcal{L}\equiv-\min\chi^{2}/2$ obtained fitting BOSS 4 skies $P_{\ell}+B_{0}+B_{2}$, and associated $p$-value. For each skycut, we detail the number of bins $N_{\rm bin}=N_{\rm bin}^{P_{\ell}}+N_{\rm bin}^{B_{0}}+N_{\rm bin}^{B_{2}}$, while in ‘Parameter Prior’ we give instead the degrees of freedom (dof). The dof are taken as the sum of 3 varied cosmological parameters (that are not prior dominated) plus an effective number of correlated EFT parameters. The $p$-value are calculated assuming there is no correlation within the data.

Figure 5: Triangle plots of base cosmological parameters measured from the analysis of BOSS power spectrum multipoles $P_{\ell}$, $\ell=0,2$, at one-loop, and bispectrum monopole $B_{0}$ at tree or one-loop level.

When analyzing the BOSS data, we find that there is no additional gain by adding all the three independent quadrupoles after one has been included. We therefore present results including only $B_{(2,3)}^{r,h}$.

In fig. 2, we show the best fit residuals and in tab. 4 the best-fit $\chi^{2}$ and associated $p$-value. The $p$-value is very good and we do not find any concerning systematic behavior in the residuals. In fig. 1, we provide the best-fit parameters, which safely lie within our $68\%$-credible intervals.

The posteriors associated to the analysis of the BOSS data are presented in fig. 1 and fig. 5. They are discussed in the Introduction. In App. C we provide the posteriors for the other non-marginalized parameters as well their confidence interval. One can see that the bispectrum improves their measurement by order 100%.

Appendix A EFTofLSS details

A.1 General expressions

Here we give the details necessary to compute the one-loop power spectrum and one-loop bispectrum of biased tracers in redshift space in the EFTofLSS. The bias expansion for the halo overdensity is given by

$$\begin{split}\delta_{h}(\vec{x},t)=&b_{1}\left(\mathbb{C}_{\delta,1}^{(1)}(\vec{x},t)+\mathbb{C}_{\delta,1}^{(2)}(\vec{x},t)+\mathbb{C}_{\delta,1}^{(3)}(\vec{x},t)+\mathbb{C}_{\delta,1}^{(4)}(\vec{x},t)\right)\\ &+b_{2}\left(\mathbb{C}_{\delta,2}^{(2)}(\vec{x},t)+\mathbb{C}_{\delta,2}^{(3)}(\vec{x},t)+\mathbb{C}_{\delta,2}^{(4)}(\vec{x},t)\right)+b_{3}\left(\mathbb{C}_{\delta,3}^{(3)}(\vec{x},t)+\mathbb{C}_{\delta,3}^{(4)}(\vec{x},t)\right)\\ &+b_{4}\;\mathbb{C}_{\delta,4}^{(4)}(\vec{x},t)+b_{5}\left(\mathbb{C}_{\delta^{2},1}^{(2)}(\vec{x},t)+\mathbb{C}_{\delta^{2},1}^{(3)}(\vec{x},t)+\mathbb{C}_{\delta^{2},1}^{(4)}(\vec{x},t)\right)\\ &+b_{6}\left(\mathbb{C}_{\delta^{2},2}^{(3)}(\vec{x},t)+\mathbb{C}_{\delta^{2},2}^{(4)}(\vec{x},t)\right)+b_{7}\,\mathbb{C}_{\delta^{2},3}^{(4)}(\vec{x},t)+b_{8}\left(\mathbb{C}_{r^{2},2}^{(3)}(\vec{x},t)+\mathbb{C}_{r^{2},2}^{(4)}(\vec{x},t)\right)\\ &+b_{9}\,\mathbb{C}_{r^{2},3}^{(4)}(\vec{x},t)+b_{10}\left(\mathbb{C}_{\delta^{3},1}^{(3)}(\vec{x},t)+\mathbb{C}_{\delta^{3},1}^{(4)}(\vec{x},t)\right)+b_{11}\,\mathbb{C}_{r^{3},2}^{(4)}(\vec{x},t)\\ &+b_{12}\,\mathbb{C}_{\delta^{3},2}^{(4)}(\vec{x},t)+b_{13}\,\mathbb{C}_{r^{2}\delta,2}^{(4)}(\vec{x},t)+b_{14}\,\mathbb{C}_{\delta^{4},1}^{(4)}(\vec{x},t)+b_{15}\,\mathbb{C}_{\delta r^{3},1}^{(4)}(\vec{x},t)\ .\end{split}$$ (44)

In the above, the $\mathbb{C}^{(n)}_{\mathcal{O},\alpha}$ functions are defined by Taylor expanding the operator $\mathcal{O}$ in the fluid line element $\vec{x}_{\rm fl}$ [58], which is given recursively by

$$\vec{x}_{\rm fl}(\vec{x},t,t^{\prime})=\vec{x}-\int_{t^{\prime}}^{t}\frac{dt^{\prime\prime}}{a(t^{\prime\prime})}\vec{v}(\vec{x}_{\rm fl}(\vec{x},t,t^{\prime\prime}),t^{\prime\prime})\ .$$ (45)

Explicitly, writing $\mathcal{O}_{m}$ to represent an operator that is the product of $m$ powers of fluctuations (i.e. $m=3$ for $\delta^{3}$, $\partial^{2}\phi\partial_{i}\partial_{j}\phi\partial^{i}\partial^{j}\phi$, etc.), we have

$$[\mathcal{O}_{m}(\vec{x}_{\rm fl}(t,t^{\prime}),t^{\prime})]^{(n)}=\sum_{\alpha=1}^{n-(m-1)}\left(\frac{D(t^{\prime})}{D(t)}\right)^{\alpha+m-1}\mathbb{C}^{(n)}_{\mathcal{O}_{m},\alpha}(\vec{x},t)\ ,$$ (46)

where the notation $[\dots]^{(n)}$ means to take the $n$-th order term in the perturbative expansion. Explicit expressions for the operators up to third order can be found in [73], with the identification $\mathbb{C}^{(3)}_{r^{2},2}=\mathbb{C}^{(3)}_{s^{2},2}+\mathbb{C}^{(3)}_{\delta^{2},2}/3$, and fourth order operators can be found with eq. (46). Operators at fourth order that are not related to previously used cubic terms are

$$\begin{split}&\mathbb{C}^{(4)}_{r^{3},2}(\vec{x},t)=[r^{3}(\vec{x},t)]^{(4)}-\partial_{i}r^{3}(\vec{x},t)\frac{\partial_{i}\theta(\vec{x},t)}{\partial^{2}}\ ,\\ &\mathbb{C}^{(4)}_{r^{2}\delta,2}(\vec{x},t)=[r^{2}(\vec{x},t)\delta(\vec{x},t)]^{(4)}-\partial_{i}(r^{2}(\vec{x},t)\delta(\vec{x},t))\frac{\partial_{i}\theta(\vec{x},t)}{\partial^{2}}\ ,\\ &\mathbb{C}^{(4)}_{\delta^{4},1}(\vec{x},t)=\delta(\vec{x},t)^{4}\ ,\quad\mathbb{C}^{(4)}_{\delta r^{3},1}(\vec{x},t)=\delta(\vec{x},t)r^{3}(\vec{x},t)\ ,\end{split}$$ (47)

where

$$r_{ij}\equiv\partial_{i}\partial_{j}\phi\ ,\quad r^{2}\equiv r_{ij}r_{ji}\ ,\quad r^{3}\equiv r_{ij}r_{jk}r_{ki}\ ,\quad\partial^{2}\phi\equiv\delta\ ,\quad\text{and}\quad\theta\equiv-\frac{\partial_{i}v^{i}}{faH}\ .$$ (48)

Given the position space halo overdensity $\delta_{h}$, the transformation to redshift space, $\delta_{r,h}$, is accomplished by

$$\begin{split}\delta_{r,h}&=\delta_{h}-\frac{\hat{z}^{i}\hat{z}^{j}}{aH}\partial_{i}\left((1+\delta_{h})v^{j}\right)+\frac{\hat{z}^{i}\hat{z}^{j}\hat{z}^{k}\hat{z}^{l}}{2(aH)^{2}}\partial_{i}\partial_{j}((1+\delta_{h})v^{k}v^{l})\\ &-\frac{\prod_{a=1}^{6}\hat{z}^{i_{a}}}{3!(aH)^{3}}\partial_{i_{1}}\partial_{i_{2}}\partial_{i_{3}}((1+\delta_{h})v^{i_{4}}v^{i_{5}}v^{i_{6}})+\frac{\prod_{a=1}^{8}\hat{z}^{i_{a}}}{4!(aH)^{4}}\partial_{i_{1}}\partial_{i_{2}}\partial_{i_{3}}\partial_{i_{4}}(v^{i_{5}}v^{i_{6}}v^{i_{7}}v^{i_{8}})+\dots\ ,\end{split}$$ (49)

where $\hat{z}$ is the line-of-sight direction. We then expand the overdensity perturbatively as

$$\delta_{r,h}(\vec{k};\hat{z})=\sum_{n}\delta^{(n)}_{r,h}(\vec{k};\hat{z})\ .$$ (50)

For the one-loop power spectrum and the one-loop bispectrum, we need the overdensity to fourth order, $n=4$. The solutions can be written as an expansion in powers of the linear dark-matter overdensity $\delta^{(1)}$ in terms of the symmetric $n$-th order halo kernels $K_{n}^{r,h}$ defined as in

$$\begin{split}\delta_{r,h}^{(1)}(\vec{k};\hat{z})&=K_{1}^{r,h}(\vec{k};\hat{z})\delta^{(1)}(\vec{k})\ ,\\ \delta_{r,h}^{(n)}(\vec{k};\hat{z})&=\int_{\vec{k}_{1},\dots,\vec{k}_{n}}^{\vec{k}}K_{n}^{r,h}(\vec{k}_{1},\dots,\vec{k}_{n};\hat{z})\delta^{(1)}(\vec{k}_{1})\cdots\delta^{(1)}(\vec{k}_{n})\ ,\quad\text{for }n\geq 2\ .\end{split}$$ (52)

For example,

$$K_{1}^{r,h}(\vec{k};\hat{z})=b_{1}+f(\hat{k}\cdot\hat{z})^{2}\ ,$$ (53)

is the famous Kaiser result for linear theory. The explicit expressions for $K_{2}^{r,h}$ and $K_{3}^{r,h}$ are given in [40], while the expression for $K_{4}^{r,h}$ is available in [90] and can be straightforwardly derived from the above expressions. We provide the dependence of the kernels $K_{1,\dots,4}^{r,h}$ on the bias parameters $\{b_{i}\}$ for $i=1,\dots,15$ here for convenience

$$\begin{split}&K_{1}^{r,h}[b_{1}]\ ,\quad K_{2}^{r,h}[b_{1},b_{2},b_{5}]\ ,\quad K_{3}^{r,h}[b_{1},b_{2},b_{3},b_{5},b_{6},b_{8},b_{10}]\ ,\quad\text{and}\quad K_{4}^{r,h}[b_{1},\dots,b_{15}]\ .\end{split}$$ (54)

Given the perturbative expansion above, we can write the observables of interest (in our case the one-loop power spectrum and the one-loop bispectrum) in terms of the kernels $K^{r,h}_{n}$. The relevant quantities for the power spectrum that enter eq. (3) are the tree-level power spectrum

$$P_{11}^{r,h}(k,\hat{k}\cdot\hat{z})=(b_{1}+f(\hat{k}\cdot\hat{z})^{2})^{2}P_{11}(k)\ ,$$ (55)

and the one-loop contributions

$$\begin{split}&P_{22}^{r,h}(k,\hat{k}\cdot\hat{z})=2\int_{\vec{q}}K_{2}^{r,h}(\vec{q},\vec{k}-\vec{q};\hat{z})^{2}P_{11}(q)P_{11}(|\vec{k}-\vec{q}|)\ ,\\ &P_{13}^{r,h}(k,\hat{k}\cdot\hat{z})=6P_{11}(k)K_{1}^{r,h}(\vec{k};\hat{z})\int_{\vec{q}}K_{3}^{r,h}(\vec{q},-\vec{q},\vec{k};\hat{z})P_{11}(q)\ .\end{split}$$ (56)

For the bispectrum, the quantities that enter eq. (4) are the tree-level bispectrum

$$B_{211}^{r,h}=2K_{1}^{r,h}(\vec{k}_{1};\hat{z})K_{1}^{r,h}(\vec{k}_{2};\hat{z})K_{2}^{r,h}(-\vec{k}_{1},-\vec{k}_{2};\hat{z})P_{11}(k_{1})P_{11}(k_{2})+\text{ 2 perms.}\ ,$$ (57)

and the one-loop contributions

$$B^{r,h}_{222}=8\int_{\vec{q}}P_{11}(q)P_{11}(|\vec{k}_{2}-\vec{q}|)P_{11}(|\vec{k}_{1}+\vec{q}|)$$

$$\hskip 108.405pt\times K_{2}^{r,h}(-\vec{q},\vec{k}_{1}+\vec{q};\hat{z})K^{r,h}_{2}(\vec{k}_{1}+\vec{q},\vec{k}_{2}-\vec{q};\hat{z})K^{r,h}_{2}(\vec{k}_{2}-\vec{q},\vec{q};\hat{z})\ ,$$

$$B_{321}^{r,h,(I)}=6P_{11}(k_{1})K_{1}^{r,h}(\vec{k}_{1};\hat{z})\int_{\vec{q}}P_{11}(q)P_{11}(|\vec{k}_{2}-\vec{q}|)$$ (58)

$$\hskip 108.405pt\times K_{3}^{r,h}(-\vec{q},-\vec{k}_{2}+\vec{q},-\vec{k}_{1};\hat{z})K^{r,h}_{2}(\vec{q},\vec{k}_{2}-\vec{q};\hat{z})+\text{ 5 perms.}\ ,$$

$$B_{321}^{r,h,(II)}=6P_{11}(k_{1})P_{11}(k_{2})K_{1}^{r,h}(\vec{k}_{1};\hat{z})K_{2}^{r,h}(\vec{k}_{1},\vec{k}_{2};\hat{z})\int_{\vec{q}}P_{11}(q)K_{3}^{r,h}(\vec{k}_{2},\vec{q},-\vec{q};\hat{z})+\text{ 5 perms.}\ ,$$

$$B_{411}^{r,h}=12P_{11}(k_{1})P_{11}(k_{2})K_{1}^{r,h}(\vec{k}_{1};\hat{z})K_{1}^{r,h}(\vec{k}_{2};\hat{z})\int_{\vec{q}}P_{11}(q)K_{4}^{r,h}(\vec{q},-\vec{q},-\vec{k}_{1},-\vec{k}_{2};\hat{z})+\text{ 2 perms.}\ .$$

Note that on the left-hand sides we have dropped the argument $(k_{1},k_{2},k_{3},\hat{k}_{1}\cdot\hat{z},\hat{k}_{2}\cdot\hat{z})$ on the bispectrum expressions to remove clutter.

Next we move on to the EFT counterterm contributions where, for renormalization up to the one-loop bispectrum, we need the EFT counterterms to second order in fields. Additionally, since we work with biased tracers, we can introduce the counterterms directly into eq. (49), i.e. directly at the level of biased tracers in redshift space, as in [40, 77] for example. We divide the counterterm contributions into two sources, response terms which are proportional to powers of the linear field $\delta^{(1)}$, and stochastic terms which contain randomly fluctuating fields, typically denoted with an ‘$\epsilon$.’ We can write the response terms as

$$\begin{split}&\delta^{(1)}_{r,h,ct}(\vec{k};\hat{z})=K_{1}^{r,h,ct}(\vec{k};\hat{z})\delta^{(1)}(\vec{k})\ ,\\ &\delta^{(2)}_{r,h,ct}(\vec{k};\hat{z})=\int_{\vec{q}_{1},\vec{q}_{2}}^{\vec{k}}K_{2}^{r,h,ct}(\vec{q}_{1},\vec{q}_{2};\hat{z})\delta^{(1)}(\vec{q}_{1})\delta^{(1)}(\vec{q}_{2})\ ,\end{split}$$ (59)

where

$$K_{1}^{r,h,ct}(\vec{k};\hat{z})=\frac{k^{2}}{k_{\rm NL}^{2}}\left(-c_{h,1}+f(\hat{k}\cdot\hat{z})^{2}c_{\pi,1}-{\textstyle{\frac{1}{2}}}f^{2}(\hat{k}\cdot\hat{z})^{4}c_{\pi v,1}-{\textstyle{\frac{1}{2}}}f^{2}(\hat{k}\cdot\hat{z})^{2}c_{\pi v,3}\right)\ ,$$ (60)

and

$$\begin{split}K_{2}^{r,h,ct}(\vec{k}_{1},\vec{k}_{2};\hat{z})=\sum_{i=1}^{14}c_{i}\,e^{K_{2}}_{i}(\vec{k}_{1},\vec{k}_{2};\hat{z})\ ,\end{split}$$ (61)

with

$$c_{i}=\{c_{h,1},c_{h,2},c_{h,3},c_{h,4},c_{h,5},c_{\pi,1},c_{\pi,5},c_{\pi v,1},c_{\pi v,2},c_{\pi v,3},c_{\pi v,4},c_{\pi v,5},c_{\pi v,6},c_{\pi v,7}\}\ ,$$ (62)

and the $e^{K_{2}}_{i}$ functions are given below in sec. A.2. Similarly, we denote the first order stochastic term as $\delta_{1}^{r,h,\epsilon}(\vec{k};\hat{z})$, and the second order ones as

$$\begin{split}&\delta^{(2)}_{r,h,\epsilon}(\vec{k};\hat{z})=\int_{\vec{q}_{1},\vec{q}_{2}}^{\vec{k}}\delta_{2}^{r,h,\epsilon}(\vec{q}_{1},\vec{q}_{2};\hat{z})\delta^{(1)}(\vec{q}_{2})\ .\end{split}$$ (63)

Here, $\delta_{1}^{r,h,\epsilon}(\vec{k};\hat{z})$ and $\delta_{2}^{r,h,\epsilon}(\vec{q}_{1},\vec{q}_{2};\hat{z})$ contain all allowed contractions of tensor stochastic fields $\epsilon^{ij\dots}$, which we assume to be Poisson distributed and do not correlate with the matter field $\delta^{(1)}$. Contractions of the stochastic fields are defined as in [77].

In terms of the kernels above, the response counterterms are

$$\begin{split}&P_{13}^{r,h,ct}(k,\hat{k}\cdot\hat{z})=2K_{1}^{r,h}(\vec{k};\hat{z})K_{1}^{r,h,ct}(-\vec{k};\hat{z})P_{11}(k)\ ,\\ &B_{321}^{r,h,(II),ct}=2P_{11}(k_{1})P_{11}(k_{2})K_{1}^{r,h,ct}(\vec{k}_{1};\hat{z})K_{1}^{r,h}(\vec{k}_{2};\hat{z})K^{r,h}_{2}(-\vec{k}_{1},-\vec{k}_{2};\hat{z})+\text{ 5 perms.}\ ,\\ &B_{411}^{r,h,ct}=2P_{11}(k_{1})P_{11}(k_{2})K_{1}^{r,h}(\vec{k}_{1};\hat{z})K_{1}^{r,h}(\vec{k}_{2};\hat{z})K_{2}^{r,h,ct}(-\vec{k}_{1},-\vec{k}_{2};\hat{z})+\text{ 2 perms.}\ .\end{split}$$ (64)

The two stochastic contributions that involve only $\delta_{1}^{r,h,\epsilon}$ and do not contain any long-wavelength fields are

$$P_{22}^{r,h,\epsilon}=\frac{1}{\bar{n}}\left(c_{1}^{\rm St}+c^{\rm St}_{2}\frac{k^{2}}{k_{\rm NL}^{2}}+c^{\rm St}_{3}\frac{k^{2}}{k_{\rm NL}^{2}}f(\hat{k}\cdot\hat{z})^{2}\right)\ ,$$ (65)

and

$$B_{222}^{r,h,\epsilon}=\frac{1}{\bar{n}^{2}}\left(c_{1}^{(222)}+\frac{1}{k_{\rm NL}^{2}}\left(c^{(222)}_{2}(k_{1}^{2}+k_{2}^{2}+k_{3}^{2})+c_{5}^{(222)}\hat{z}^{i}\hat{z}^{j}\left(k_{1}^{i}k_{2}^{j}+k_{1}^{i}k_{3}^{j}+k_{2}^{i}k_{3}^{j}\right)\right)\right)\ .$$ (66)

The final, mixed response-stochastic, contribution is slightly more complicated, so we first define $\tilde{B}_{321}^{r,h,(I),\epsilon}$ from

$$\begin{split}&(2\pi)^{3}\delta_{D}(\vec{k}_{1}+\vec{k}_{2}+\vec{k}_{3})\tilde{B}_{321}^{r,h,(I),\epsilon}(\vec{k}_{1},\vec{k}_{2},\vec{k}_{3};\hat{z})=\\ &\hskip 72.26999pt\langle\delta_{r,h}^{(1)}(\vec{k}_{1};\hat{z})\delta_{1}^{r,h,\epsilon}(\vec{k}_{2};\hat{z})\delta^{(2)}_{r,h,\epsilon}(\vec{k}_{3};\hat{z})\rangle+\langle\delta_{r,h}^{(1)}(\vec{k}_{1};\hat{z})\delta_{1}^{r,h,\epsilon}(\vec{k}_{3};\hat{z})\delta^{(2)}_{r,h,\epsilon}(\vec{k}_{2};\hat{z})\rangle\ ,\end{split}$$ (67)

so that

$$B_{321}^{r,h,(I),\epsilon}=\tilde{B}_{321}^{r,h,(I),\epsilon}(\vec{k}_{1},\vec{k}_{2},\vec{k}_{3};\hat{z})+\tilde{B}_{321}^{r,h,(I),\epsilon}(\vec{k}_{3},\vec{k}_{1},\vec{k}_{2};\hat{z})+\tilde{B}_{321}^{r,h,(I),\epsilon}(\vec{k}_{2},\vec{k}_{3},\vec{k}_{1};\hat{z})\ .$$ (68)

Then we have

$$\tilde{B}_{321}^{r,h,(I),\epsilon}(\vec{k}_{1},\vec{k}_{2},\vec{k}_{3};\hat{z})=\frac{(b_{1}+f(\hat{k}_{1}\cdot\hat{z})^{2})}{\bar{n}}P_{11}(k_{1})\sum_{i=1}^{13}c^{\rm St}_{i}e_{i}^{\rm St}(\vec{k}_{1},\vec{k}_{2},\vec{k}_{3};\hat{z})\ .$$ (69)

The $e_{i}^{\rm St}$ functions are given below in sec. A.2. Notice that $e_{3}^{\rm St}=0$, so that there is no $c_{3}^{\rm St}$ parameter.

A full description of the fourth order bias expansion and renormalization in redshift space is given in [90]. See [124] for the bias expansion for the real space one-loop bispectrum (and measurement of the scalar amplitude $A_{s}$ from simulations).

A.2 Explicit functions

The functions that enter $K_{2}^{r,h,ct}$ in eq. (61) are given by

$$\begin{split}e_{1}^{K_{2}}&=\frac{-2fk_{2}k_{1}^{3}\mu_{1}\mu_{2}-2fk_{2}^{4}\mu_{1}^{2}+k_{1}^{4}+\left(k_{2}^{2}-k_{3}^{2}\right)k_{1}^{2}}{4k_{2}^{2}k_{\text{NL}}^{2}}+(1\leftrightarrow 2)\ ,\\ e_{2}^{K_{2}}&=-\frac{1}{28k_{1}^{2}k_{2}^{2}k_{\text{NL}}^{2}}\left(k_{3}^{6}+5k_{1}^{2}k_{2}^{2}k_{3}^{2}+7k_{1}^{6}+\left(7k_{2}^{2}-12k_{3}^{2}\right)k_{1}^{4}+3k_{2}^{2}k_{3}^{4}+(1\leftrightarrow 2)\right)\ ,\\ e_{3}^{K_{2}}&=-\frac{k_{3}^{2}}{k_{\rm NL}^{2}}\ ,\quad e_{4}^{K_{2}}=-\frac{k_{3}^{2}\left(k_{1}^{2}+k_{2}^{2}-k_{3}^{2}\right){}^{2}}{4k_{1}^{2}k_{2}^{2}k_{\text{NL}}^{2}}\ ,\quad e_{5}^{K_{2}}=\frac{k_{1}^{2}+k_{2}^{2}-k_{3}^{2}}{2k_{\text{NL}}^{2}}\ ,\\ \end{split}$$ (70)

$$e_{6}^{K_{2}} =\frac{f(k_{1}\mu_{1}+k_{2}\mu_{2})}{4k_{1}^{2}k_{2}^{2}k_{\text{NL}}^{2}}\left(k_{1}^{2}\mu_{1}\left(2fk_{2}k_{1}^{2}\mu_{1}\mu_{2}+2fk_{2}^{3}\mu_{1}\mu_{2}-k_{1}^{3}+\left(k_{3}^{2}-k_{2}^{2}\right)k_{1}\right)+(1\leftrightarrow 2)\right)\ ,$$

$$e_{7}^{K_{2}} =\frac{f\left(k_{1}^{2}-k_{2}^{2}-k_{3}^{2}\right)\left(k_{1}^{2}+k_{2}^{2}-k_{3}^{2}\right)\left(k_{1}^{2}-k_{2}^{2}+k_{3}^{2}\right)\left(k_{1}\mu_{1}+k_{2}\mu_{2}\right){}^{2}}{8k_{1}^{2}k_{2}^{2}k_{3}^{2}k_{\text{NL}}^{2}}\ ,$$ (71)

$$e_{8}^{K_{2}} =\frac{f^{2}(k_{1}\mu_{1}+k_{2}\mu_{2})^{2}}{8k_{1}^{2}k_{2}^{2}k_{\text{NL}}^{2}}\left(-2fk_{1}^{2}k_{2}^{2}\mu_{1}^{2}\mu_{2}^{2}+k_{1}^{2}\mu_{1}^{2}\left(-2fk_{2}k_{1}\mu_{1}\mu_{2}+k_{1}^{2}+k_{2}^{2}-k_{3}^{2}\right)+(1\leftrightarrow 2)\right)\ ,$$

$$\begin{split}e_{9}^{K_{2}}&=\frac{f^{2}(k_{1}\mu_{1}+k_{2}\mu_{2})^{2}}{56k_{1}^{2}k_{2}^{2}k_{3}^{2}k_{\text{NL}}^{2}}\Bigg{(}-2k_{1}k_{2}\left(k_{3}^{4}+5k_{1}^{2}k_{2}^{2}\right)\mu_{1}\mu_{2}+k_{1}^{2}\mu_{1}\Big{(}5k_{1}^{4}\mu_{1}+10k_{2}k_{1}^{3}\mu_{2}\\ &\hskip 72.26999pt-10\left(k_{2}^{2}+k_{3}^{2}\right)k_{1}^{2}\mu_{1}-6k_{2}k_{3}^{2}k_{1}\mu_{2}+5\left(k_{2}^{2}-k_{3}^{2}\right){}^{2}\mu_{1}\Big{)}+(1\leftrightarrow 2)\Bigg{)}\ ,\\ e_{10}^{K_{2}}&=\frac{f^{2}(k_{1}\mu_{1}+k_{2}\mu_{2})^{2}}{8k_{1}^{2}k_{2}^{2}k_{\rm NL}^{2}}\left(k_{1}^{2}k_{2}^{2}+k_{1}^{2}\left(-2fk_{2}k_{1}\mu_{1}\mu_{2}-2fk_{2}^{2}\mu_{1}^{2}+k_{1}^{2}-k_{3}^{2}\right)+(1\leftrightarrow 2)\right)\ ,\\ e_{11}^{K_{2}}&=\frac{f^{2}(k_{1}\mu_{1}+k_{2}\mu_{2})^{2}}{28k_{1}^{2}k_{2}^{2}k_{\rm NL}^{2}}\left(-k_{3}^{4}-12k_{1}^{2}k_{2}^{2}+\left(2k_{3}^{2}-k_{1}^{2}\right)k_{1}^{2}+\left(2k_{3}^{2}-k_{2}^{2}\right)k_{2}^{2}\right)\ ,\\ e_{12}^{K_{2}}&=-\frac{f^{2}\left(\mu_{1}^{2}+\mu_{2}^{2}\right)\left(k_{1}\mu_{1}+k_{2}\mu_{2}\right){}^{2}}{4k_{\text{NL}}^{2}}\ ,\quad e_{13}^{K_{2}}=\frac{f^{2}\left(k_{1}^{2}+k_{2}^{2}-k_{3}^{2}\right)\mu_{1}\mu_{2}\left(k_{1}\mu_{1}+k_{2}\mu_{2}\right){}^{2}}{4k_{1}k_{2}k_{\text{NL}}^{2}}\ ,\\ e_{14}^{K_{2}}&=-\frac{f^{2}\left(k_{1}\mu_{1}+k_{2}\mu_{2}\right){}^{2}}{2k_{\text{NL}}^{2}}\ .\end{split}$$ (72)

In the above, we have used the notation $\mu_{i}\equiv\hat{k}_{i}\cdot\hat{z}$.

The functions that enter the stochastic counterterm $B_{321}^{r,h,(I),\epsilon}$ in eq. (69) are:

$$e_{1}^{\rm St}=f\mu_{1}^{2}-1\ ,$$

$$e_{2}^{\rm St}=-\frac{k_{1}^{2}\left(k_{2}^{2}\left(1-2f\mu_{1}^{2}\right)+k_{3}^{2}\right)+2fk_{2}\left(k_{3}^{2}-k_{2}^{2}\right)k_{1}\mu_{1}\mu_{2}+\left(k_{2}^{2}-k_{3}^{2}\right){}^{2}}{2k_{1}^{2}k_{\text{NL}}^{2}}\ ,$$ (73)

$$e_{3}^{\rm St}=0\ ,$$

$$e_{4}^{\rm St}=-\frac{f^{2}\mu_{1}\left(k_{1}^{3}\mu_{1}\left(2f\mu_{1}^{2}-1\right)+4fk_{2}k_{1}^{2}\mu_{1}^{2}\mu_{2}+k_{1}\mu_{1}\left(k_{2}^{2}\left(4f\mu_{2}^{2}-1\right)+k_{3}^{2}\right)+2k_{2}\left(k_{3}^{2}-k_{2}^{2}\right)\mu_{2}\right)}{4k_{1}k_{\text{NL}}^{2}}\ ,$$

$$e_{5}^{\rm St}=\frac{f^{2}\mu_{1}\left(4fk_{2}k_{1}^{2}\mu_{1}^{2}\mu_{2}+k_{1}\mu_{1}\left(k_{2}^{2}\left(4f\mu_{2}^{2}-1\right)+k_{3}^{2}\right)+k_{1}^{3}\mu_{1}+2k_{2}\left(k_{3}^{2}-k_{2}^{2}\right)\mu_{2}\right)}{4k_{1}k_{\text{NL}}^{2}}\ ,$$

$$e_{6}^{\rm St}=2\ ,\quad e_{7}^{\rm St}=-\frac{k_{2}^{2}+k_{3}^{2}}{k_{\text{NL}}^{2}}\ ,\quad e_{8}^{\rm St}=-\frac{k_{1}^{4}+\left(k_{2}^{2}-k_{3}^{2}\right){}^{2}}{2k_{1}^{2}k_{\text{NL}}^{2}}\ ,$$

$$e_{9}^{\rm St}=-\frac{k_{1}^{2}}{k_{\text{NL}}^{2}}\ ,\quad e_{10}^{\rm St}=-\frac{f\left(k_{1}\mu_{1}+2k_{2}\mu_{2}\right)\left(\left(k_{1}^{2}-k_{2}^{2}+k_{3}^{2}\right)\mu_{1}+2k_{1}k_{2}\mu_{2}\right)}{4k_{1}k_{\text{NL}}^{2}}\ ,$$ (74)

$$e_{11}^{\rm St}=\frac{f\mu_{1}\left(k_{1}\left(k_{1}^{2}+k_{2}^{2}-k_{3}^{2}\right)\mu_{1}+2k_{2}\left(k_{2}^{2}-k_{3}^{2}\right)\mu_{2}\right)}{2k_{1}k_{\text{NL}}^{2}}\ ,\quad e_{12}^{\rm St}=-\frac{2fk_{2}\mu_{2}\left(k_{1}\mu_{1}+k_{2}\mu_{2}\right)}{k_{\text{NL}}^{2}}\ ,$$

$$e_{13}^{\rm St}=\frac{1}{4k_{1}^{2}k_{3}^{2}k_{\text{NL}}^{2}}f\Big{(}k_{1}^{2}\left(k_{1}^{2}-k_{2}^{2}+k_{3}^{2}\right){}^{2}\mu_{1}^{2}+2k_{1}k_{2}\left(k_{1}^{2}-k_{2}^{2}+k_{3}^{2}\right){}^{2}\mu_{1}\mu_{2}$$

$$\hskip 72.26999pt+\left(\left(k_{2}^{2}+k_{3}^{2}\right)k_{1}^{4}-2\left(k_{2}^{2}-k_{3}^{2}\right){}^{2}k_{1}^{2}+\left(k_{2}^{2}-k_{3}^{2}\right){}^{2}\left(k_{2}^{2}+k_{3}^{2}\right)\right)\mu_{2}^{2}\Big{)}\ .$$

All of the above $e_{i}^{\rm St}$ are symmetric when swapping $\vec{k}_{2}$ and $\vec{k}_{3}$, as expected from eq. (LABEL:btilde321). To see it, one must swap $k_{2}\leftrightarrow k_{3}$ and $\mu_{2}\leftrightarrow\mu_{3}$, and then replace $\mu_{3}=-k_{3}^{-1}(k_{1}\mu_{1}+k_{2}\mu_{2})$.

Appendix B Binning formula details

In this appendix, we want to show that the binning formula for the bispectrum

$$B_{(l,i),\rm bin}^{r,h}(k_{1},k_{2},k_{3})=\frac{2l+1}{V_{T}}\left(\prod_{i}\int_{V_{i}}\frac{\mathrm{d}^{3}q_{i}}{(2\pi)^{3}}\right)(2\pi)^{3}\delta_{D}^{(3)}(\vec{q}_{1}+\vec{q}_{2}+\vec{q}_{3})\mathcal{P}_{l}(\mu_{i})B^{r,h}(\vec{q}_{1},\vec{q}_{2},\vec{q}_{3})\ ,$$ (75)

is equivalent to

$$B_{(l,i),\rm bin}^{r,h}(k_{1},k_{2},k_{3})=\frac{1}{V_{T}}\left(\prod_{i}\int_{k_{i}}\mathrm{d}q_{i}\,q_{i}\right)\,\frac{\beta\left(\Delta_{q}\right)}{8\pi^{4}}B^{r,h}_{{(l,i)}}(q_{1},q_{2},q_{3})\,\,,$$ (76)

We do the calculation here for general $l$, which is relevant to us for $l=0,2$.

Given that the bispectrum is a polynomial in $\mu_{1}$ and $\mu_{2}$, and switching to a basis of Legendre polynomials, we can write

$$B^{r,h}(\vec{q}_{1},\vec{q}_{2},\vec{q}_{3})=\sum_{n_{1},n_{2}}B^{r,h}_{n_{1},n_{2}}(q_{1},q_{2},q_{3})\mathcal{P}_{n_{1}}(\mu_{1})\mathcal{P}_{n_{2}}(\mu_{2})\,.$$ (77)

We can focus on the case of $\mathcal{P}_{l}(\mu_{1})$, since the other cases just correspond to a permutation of $\mu_{i}$ in eq. (77). Let us then start by writing the delta function as an integral over plane waves:

$$\begin{split}B^{r,h}_{(l,1),\rm bin}&=\frac{2l+1}{V_{T}}\left(\prod_{i}\int_{V_{i}}\frac{\mathrm{d}^{3}q_{i}}{(2\pi)^{3}}\right)\int\mathrm{d}^{3}x\,e^{i\vec{x}\cdot(\vec{q}_{1}+\vec{q}_{2}+\vec{q}_{3})}\mathcal{P}_{l}(\mu_{1})\\ &\hskip 144.54pt\times\sum_{n_{1},n_{2}}B^{r,h}_{n_{1},n_{2}}(q_{1},q_{2},q_{3})\mathcal{P}_{n_{1}}(\mu_{1})\mathcal{P}_{n_{2}}(\mu_{2})\,.\end{split}$$ (78)

The integral over $\mathrm{d}^{2}\hat{q}_{3}$, is just

$$\int\mathrm{d}^{2}\hat{q}_{3}\,e^{i\vec{q}_{3}\cdot\vec{x}}=4\pi j_{0}(q_{3}x)\,.$$ (79)

And the rest of the exponentials we expand the plane wave:

$$e^{i\vec{q}_{1}\cdot\vec{x}}=\sum_{l_{1}}i^{l_{1}}(2l_{1}+1)j_{l_{1}}(q_{1}x)\mathcal{P}_{l_{1}}(\hat{q}_{1}\cdot\hat{x}).$$ (80)

Putting all of this into eq. (78) we get:

$$\begin{split}B^{r,h}_{(l,1),\rm bin}&=\frac{2l+1}{V_{T}}\sum_{n_{1},n_{2}}\sum_{l_{1},l_{2}}i^{l_{1}+l_{2}}(2l_{1}+1)(2l_{2}+1)\left(\prod_{i=1}^{2}\int_{V_{i}}\frac{\mathrm{d}^{3}q_{i}}{(2\pi)^{3}}\right)\int_{k_{3}}\frac{\mathrm{d}q_{3}}{2\pi^{2}}q_{3}^{2}B^{r,h}_{n_{1},n_{2}}(q_{1},q_{2},q_{3})\\ &\times\int\mathrm{d}^{3}xj_{0}(q_{3}x)j_{l_{1}}(q_{1}x)j_{l_{2}}(q_{2}x)\mathcal{P}_{l_{1}}(\hat{q}_{1}\cdot\hat{x})\mathcal{P}_{l_{2}}(\hat{q}_{2}\cdot\hat{x})\mathcal{P}_{l}(\hat{q}_{1}\cdot\hat{z})\mathcal{P}_{n_{1}}(\hat{q}_{1}\cdot\hat{z})\mathcal{P}_{n_{2}}(\hat{q}_{2}\cdot\hat{z})\,.\end{split}$$ (81)

Now we can do all the angular integrals, using the formula,

$$\int\mathrm{d}^{2}\hat{x}\mathcal{P}_{l_{1}}(\hat{q}_{1}\cdot\hat{x})\mathcal{P}_{l_{2}}(\hat{q}_{2}\cdot\hat{x})=\delta_{l_{1},l_{2}}\,\frac{4\pi}{2l_{1}+1}\mathcal{P}_{l_{1}}(\hat{q}_{1}\cdot\hat{q}_{2})\ ,$$ (82)

to evaluate

$$\begin{split}\int\mathrm{d}^{2}\hat{x}\int\mathrm{d}^{2}&\hat{q}_{1}\int\mathrm{d}^{2}\hat{q}_{2}\,\mathcal{P}_{l_{1}}(\hat{q}_{1}\cdot\hat{x})\mathcal{P}_{l_{2}}(\hat{q}_{2}\cdot\hat{x})\mathcal{P}_{l}(\hat{q}_{1}\cdot\hat{z})\mathcal{P}_{n_{1}}(\hat{q}_{1}\cdot\hat{z})\mathcal{P}_{n_{2}}(\hat{q}_{2}\cdot\hat{z})\\ &=\delta_{l_{1},l_{2}}\,\frac{4\pi}{2l_{1}+1}\int\mathrm{d}^{2}\hat{q}_{1}\int\mathrm{d}^{2}\hat{q}_{2}\,\mathcal{P}_{l_{1}}(\hat{q}_{1}\cdot\hat{q}_{2})\mathcal{P}_{l}(\hat{q}_{1}\cdot\hat{z})\mathcal{P}_{n_{1}}(\hat{q}_{1}\cdot\hat{z})\mathcal{P}_{n_{2}}(\hat{q}_{2}\cdot\hat{z})\\ &=\delta_{l_{1},l_{2}}\,\delta_{l_{1},n_{2}}\,\frac{(4\pi)^{2}}{(2l_{1}+1)^{2}}\int\mathrm{d}^{2}\hat{q}_{1}\,\mathcal{P}_{l_{1}}(\hat{q}_{1}\cdot\hat{z})\mathcal{P}_{l}(\hat{q}_{1}\cdot\hat{z})\mathcal{P}_{n_{1}}(\hat{q}_{1}\cdot\hat{z})\\ &=\delta_{l_{1},l_{2}}\,\delta_{l_{1},n_{2}}\,\frac{(4\pi)^{3}}{(2l_{1}+1)^{2}}\begin{pmatrix}n_{2}&l&n_{1}\\ 0&0&0\end{pmatrix}^{2}\ ,\end{split}$$ (83)

and additionally we used the integral of three Legendre polynomials in terms of the Wigner 3-j symbol in the last line. We are now only left with integrals over the magnitudes:

$$\begin{split}B^{r,h}_{(l,1),\rm bin}&=4\pi\frac{(2l+1)}{V_{T}}\sum_{n_{1},n_{2}}(-1)^{n_{2}}\left(\prod_{i=1}^{3}\int_{k_{i}}\frac{\mathrm{d}q_{i}}{2\pi^{2}}q_{i}^{2}\right)\\ &\quad\times\int_{0}^{\infty}\mathrm{d}xx^{2}j_{n_{2}}(q_{1}x)j_{n_{2}}(q_{2}x)j_{0}(q_{3}x)\begin{pmatrix}n_{2}&l&n_{1}\\ 0&0&0\end{pmatrix}^{2}B^{r,h}_{n_{1},n_{2}}(q_{1},q_{2},q_{3})\,.\end{split}$$ (84)

Further using the following integral over three spherical Bessel functions:

$$\int_{0}^{\infty}\mathrm{d}xx^{2}j_{n_{2}}(q_{1}x)j_{n_{2}}(q_{2}x)j_{0}(q_{3}x)=\frac{\pi}{4q_{1}q_{2}q_{3}}\beta(\hat{q}_{1}\cdot\hat{q}_{2})\mathcal{P}_{n_{2}}\left(\frac{q_{1}^{2}+q_{2}^{2}-q_{3}^{2}}{2q_{1}q_{2}}\right)\ ,$$ (85)

where $\beta(\Delta)=1$ for $-1<\Delta<1$, $\beta(\Delta)=1/2$ for $\Delta=\pm 1$, and $\beta(\Delta)=0$ otherwise, and recognizing that the last Legendre is $\mathcal{P}_{n_{2}}(-\hat{q}_{1}\cdot\hat{q}_{2})=(-1)^{n_{2}}\mathcal{P}_{n_{2}}(\hat{q}_{1}\cdot\hat{q}_{2})$, we can put everything together and we get that eq. (75) reduces to

$$B^{r,h}_{(l,1),\rm bin}=\frac{2l+1}{V_{T}}\sum_{n_{1},n_{2}}\left(\prod_{i=1}^{3}\int_{k_{i}}\mathrm{d}q_{i}q_{i}\right)\frac{\beta(\hat{q}_{1}\cdot\hat{q}_{2})}{8\pi^{4}}\mathcal{P}_{n_{2}}(\hat{q}_{1}\cdot\hat{q}_{2})\begin{pmatrix}n_{2}&l&n_{1}\\ 0&0&0\end{pmatrix}^{2}B^{r,h}_{n_{1},n_{2}}(q_{1},q_{2},q_{3})\,.$$ (86)

To get our formula (76), it is now sufficient to show that the unbinned bispectrum satisfies

$$B^{r,h}_{(l,1)}(q_{1},q_{2},q_{3})=(2l+1)\sum_{n_{1},n_{2}}\mathcal{P}_{n_{2}}(\hat{q}_{1}\cdot\hat{q}_{2})\begin{pmatrix}n_{2}&l&n_{1}\\ 0&0&0\end{pmatrix}^{2}B^{r,h}_{n_{1},n_{2}}(q_{1},q_{2},q_{3})\,.$$ (87)

So next, we write the left hand side of the above explicitly, and expand the redshift space bispectrum, plugging eq. (77) into a generalization of eq. (8) and eq. (LABEL:bispmultipolesused):

$$B^{r,h}_{(l,1)}(q_{1},q_{2},q_{3})=(2l+1)\sum_{n_{1},n_{2}}\int_{-1}^{1}\frac{\mathrm{d}\mu_{1}}{2}\int_{0}^{2\pi}\frac{\mathrm{d}\phi}{2\pi}\mathcal{P}_{l}(\mu_{1})\mathcal{P}_{n_{1}}(\mu_{1})\mathcal{P}_{n_{2}}(\mu_{2})B^{r,h}_{n_{1},n_{2}}(q_{1},q_{2},q_{3})\,.$$ (88)

This can be calculated in a coordinate system in which we fix $\hat{q}_{1}$, $\hat{q}_{2}$ and integrate over $\mathrm{d}^{2}\hat{z}$:

$$B^{r,h}_{(l,1)}(q_{1},q_{2},q_{3})=(2l+1)\sum_{n_{1},n_{2}}\int\frac{\mathrm{d}^{2}\hat{z}}{4\pi}\mathcal{P}_{l}(\hat{q}_{1}\cdot\hat{z})\mathcal{P}_{n_{1}}(\hat{q}_{1}\cdot\hat{z})\mathcal{P}_{n_{2}}(\hat{q}_{2}\cdot\hat{z})B^{r,h}_{n_{1},n_{2}}(q_{1},q_{2},q_{3})\,.$$ (89)

Next, we use that the product of two Legendre polynomials is

$$\mathcal{P}_{l}(\hat{q}_{1}\cdot\hat{z})\mathcal{P}_{n_{1}}(\hat{q}_{1}\cdot\hat{z})=\sum_{L=|n_{1}-l|}^{n_{1}+l}(2L+1)\begin{pmatrix}n_{1}&l&L\\ 0&0&0\end{pmatrix}^{2}\mathcal{P}_{L}(\hat{q}_{1}\cdot\hat{z})\ ,$$ (90)

and plug eq. (90) into eq. (89) to get

$$\begin{split}B^{r,h}_{(l,1)}(q_{1},q_{2},q_{3})&=(2l+1)\sum_{n_{1},n_{2}}\sum_{L=|n_{1}-l|}^{n_{1}+l}(2L+1)\begin{pmatrix}n_{1}&l&L\\ 0&0&0\end{pmatrix}^{2}\\ &\hskip 108.405pt\times\int\frac{\mathrm{d}^{2}\hat{z}}{4\pi}\mathcal{P}_{L}(\hat{q}_{1}\cdot\hat{z})\mathcal{P}_{n_{2}}(\hat{q}_{2}\cdot\hat{z})B^{r,h}_{n_{1},n_{2}}(q_{1},q_{2},q_{3})\\ &=(2l+1)\sum_{n_{1},n_{2}}\begin{pmatrix}n_{1}&l&L\\ 0&0&0\end{pmatrix}^{2}\mathcal{P}_{n_{2}}(\hat{q}_{1}\cdot\hat{q}_{2})B^{r,h}_{n_{1},n_{2}}(q_{1},q_{2},q_{3})\ ,\end{split}$$ (91)

as desired.

For completeness we also calculate the volume

$$\begin{split}V_{T}&=\left(\prod_{i}\int_{V_{i}}\frac{\mathrm{d}^{3}q_{i}}{(2\pi)^{3}}\right)(2\pi)^{3}\delta_{D}^{(3)}(\vec{q}_{1}+\vec{q}_{2}+\vec{q}_{3})=\left(\prod_{i}\int_{V_{i}}\frac{\mathrm{d}^{3}q_{i}}{(2\pi)^{3}}\right)\int\mathrm{d}^{3}x\,e^{i\vec{q}_{1}\cdot\vec{x}}e^{i\vec{q}_{2}\cdot\vec{x}}e^{i\vec{q}_{3}\cdot\vec{x}}\,.\end{split}$$ (92)

We then integrate over the plane waves using eq. (79) and the three Bessel functions using eq. (85) to get

$$\begin{split}V_{T}&=\left(\prod_{i}\int_{k_{i}}\mathrm{d}q_{i}\,q_{i}\right)\frac{\beta(\hat{q}_{1}\cdot\hat{q}_{2})}{8\pi^{4}}\ .\end{split}$$ (93)

Figure 6: Full triangle plots from the analysis of BOSS power spectrum multipoles $P_{\ell}$ at one loop, bispectrum monopole $B_{0}$ at tree level or one loop, and bispectrum quadrupole $B_{2}$ at tree level.

Appendix C Additional parameter posteriors

In fig. 6, we show the full triangle plots obtained fitting BOSS 4 skies $P_{\ell}+B_{0}+B_{2}$. In tab. 5, we show the 68%-credible intervals of $b_{1},c_{2}$, and $c_{4}$ obtained on this same fit.

mean $\pm\sigma$		$b_{1}$	$c_{2}$	$c_{4}$
CMASS NGC	$P_{\ell}$	$2.12\pm 0.18$	$0.95^{+0.45}_{-0.68}$	$0.2\pm 1.7$
	$P_{\ell}+B_{0}$	$2.23\pm 0.13$	$1.27\pm 0.26$	$-0.29^{+0.55}_{-0.61}$
	$P_{\ell}+B_{0}+B_{2}$	$2.19\pm 0.13$	$1.18^{+0.22}_{-0.28}$	$-0.25^{+0.54}_{-0.60}$
CMASS SGC	$P_{\ell}$	$2.13\pm 0.18$	$1.01^{+0.45}_{-0.62}$	$0.2\pm 1.7$
	$P_{\ell}+B_{0}$	$2.27\pm 0.13$	$1.23\pm 0.26$	$-0.32^{+0.56}_{-0.63}$
	$P_{\ell}+B_{0}+B_{2}$	$2.22\pm 0.14$	$1.14^{+0.23}_{-0.27}$	$-0.27\pm 0.60$
LOWZ NGC	$P_{\ell}$	$1.93\pm 0.16$	$0.98^{+0.37}_{-0.47}$	$0.2\pm 1.7$
	$P_{\ell}+B_{0}$	$2.04\pm 0.12$	$1.23^{+0.21}_{-0.24}$	$-0.26\pm 0.64$
	$P_{\ell}+B_{0}+B_{2}$	$2.00\pm 0.12$	$1.14^{+0.20}_{-0.24}$	$-0.27^{+0.60}_{-0.67}$
LOWZ SGC	$P_{\ell}$	$1.93\pm 0.15$	$1.04^{+0.34}_{-0.40}$	$0.2\pm 1.7$
	$P_{\ell}+B_{0}$	$2.05\pm 0.12$	$1.21^{+0.21}_{-0.24}$	$-0.30\pm 0.65$
	$P_{\ell}+B_{0}+B_{2}$	$2.02\pm 0.12$	$1.12^{+0.20}_{-0.24}$	$-0.30^{+0.61}_{-0.68}$

Table 5: 68%-credible intervals of $b_{1},c_{2}$, and $c_{4}$ from the analysis of BOSS power spectrum multipoles $P_{\ell}$ at the one-loop, bispectrum monopole $B_{0}$ at the one-loop, and bispectrum quadrupole $B_{2}$ at tree-level.