The BOSS bispectrum analysis at one loop from the Effective Field Theory of Large-Scale Structure

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.


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 Λ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, Ω 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 ≲ 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 NL parameters, f equil.NL = 217 ± 297 , f orth.NL = −64 ± 74 , f loc.NL = 49 ± 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 loc.
NL 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. 1n 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 ΛCDM parameters with BBN prior on the baryon abundance, Ω 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 Λ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 ΛCDM parameters σ 8 , h, and Ω 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. 2Notice also that the results improve significantly upon the ones obtained using instead the tree-level prediction for the bispectrum monopole: in particular, σ 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 σ 8 , h, and Ω m to values consistent at 0.3σ, 1.4σ, 0.5σ, respectively, with the ones of Planck νΛCDM [15].
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 some theoretical developments of the EFTofLSS, such as a careful understanding of renormalization [42,54,55] (including rather-subtle aspects such as lattice-running [42] and a better understanding of the velocity field [44,56]), of several ways for extracting the value of the counterterms from simulations [42,57], and of the non-locality in time of the EFTofLSS [44,46,58].These theoretical explorations also include an enlightening study in 1+1 dimensions [57].An IR-resummation of the long displacement fields had to be performed in order to reproduce the Baryon Acoustic Oscillation (BAO) peak, giving rise to the so-called IR-Resummed EFTofLSS [59,60,61,62,63].Accounts of baryonic effects were presented in [64,65].The dark-matter bispectrum has been computed at one-loop in [66,67], the one-loop trispectrum in [68], and the displacement field in [69].The lensing power spectrum has been computed at two loops in [70].Biased tracers, such as halos and galaxies, have been studied in the context of the EFTofLSS in [58,71,72,73,40,74,75] (see also [76]), the halo and matter power spectra and bispectra (including all cross correlations) in [58,72].Redshift space distortions have been developed in [59,77,40].Neutrinos have been included in the EFTofLSS in [78,79], clustering dark energy in [80,52,81,82], and primordial non-Gaussianities in [72,83,84,85,77,86].Faster evaluation schemes for the calculation of some of the loop integrals have been developed in [87].Comparison with high-quality N -body simulations to show that the EFTofLSS can accurately recover the cosmological parameters have been performed in [4,6,88,89]. 2Here and in the rest of this work, we quote parameter constraints as the Bayesian 68% credible interval from the one-dimensional marginalized posterior.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.

Data
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]. 3 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]. 4 We will analyze the mean of the 84 Nseries realizations (hereafter referred as 'Nseries').All celestial coordinates are converted to comoving distance assuming Ω fid m = 0.310.
BOSS P+B full-shape measurements.In this work, we analyze the full shape of the power spectrum multipoles ℓ = 0, 2, and of the bispectrum monopole and quadrupole (respectively abbreviated 'P ℓ ', '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]. 5For the power spectrum, we find excellent agreement between the measurements from Rustico and Nbodykit [103]. 6We 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 box = 3500 (2300)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 ∆k ≃ 0.01h Mpc −1 .Instead, we bin the bispectrum in ∆n = 12 (9) units of the fundamental frequency of the box k f for CMASS (LOWZ), starting from the bin centered at n min = 6 + ∆n/2, up to the one centered on n max = 126 which correspond in frequencies to bins of size ∆k = 0.02154 (0.02459)h Mpc −1 , with first and last bins centered on k min = 0.0215 (0.029)h Mpc −1 and k max = 0.215 (0.176)h 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 ℓ 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: (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 δ r,h , in Fourier space: where ẑ is the line-of-sight direction, and δ 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-loop tot.= P r,h 11 + (P r,h 13 + P r,h,ct

13
) + (P r,h 22 + P r,h,ϵ 22 ) , and the one-loop bispectrum: ) , 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 while for the counterterms, we have , c ] .
Notice that the diagrams P r,h 13 , B r,h,(II) 321 , and B r,h 411 depend on less biases than the kernels in eq. ( 54) would suggest.This is because, when considering the particular momentumconfiguration of the kernels that enter the loop in eq. ( 56) and eq.( 58), 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 ẑ.In particular, we analyze the powerspectrum and bispectrum monopole and quadrupoles.The power-spectrum multipoles are given by where P ℓ are the Legendre polynomials, and µ = k • ẑ.The bispectrum monopole is the average over the line-of-sight angles [109,102,110] where µ i = ki • ẑ, and explicitly, from the triangle conditions: The expectation values of the estimator used by Rustico for the quadrupoles are: 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 and c (222) 5 become degenerate, so we redefine c

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: where Here P w (k) = P 11 (k) − P nw (k), and P nw (k) is the no-wiggle power spectrum, which we obtain using the sine-transform algorithm described in [112] and detailed in [113].Then Σ 2 tot is defined by where we choose Λ = 1 h Mpc −1 and x osc = 110 Mpc/h, and j l are the spherical Bessel functions. 9For the loop, our choice is to only substitute the non-integrated P 11 (k) with 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].

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: where Rather than projecting (18) into multipoles, we project (18) as if there was no window function, and for [W * P 11 ]( ⃗ k) we use the following: for the monopole, we use W → W 00 and, for the quadrupole we use W → 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.

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: where D A is the angular diameter distance, and the factors of H 0 are there since our wavenumbers are in units h 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: where F = q ∥ /q ⊥ .To match the measured power spectrum multipoles, we do the following integral: The formula for the bispectrum is: 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 Ω m ), it is an effect of at most ∼ 1σ, and accordingly, the change in χ 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 µ and 4 points in ϕ, after using symmetries to restrict the integration domain to µ 1 ∈ [0, 1], and ϕ ∈ [−π/2, π/2].

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 ∆k around a central triangle with sides k 1 , k 2 , k 3 .Especially for our chosen bins with ∆k ≃ 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: and we note that, here and elsewhere, B r,h (⃗ q 1 , ⃗ q 2 , ⃗ q 3 ) is the full redshift-space bispectrum, i.e.
we have suppressed the dependence on ẑ for notational convenience.Here N T is the number of fundamental triangles in the bin, δ K is the Kronecker delta function, and the notation ⃗ q i ∈ k i means a sum over the fundamental modes ⃗ q i for which 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 + ∆k/2 > k 1 + k 2 − ∆k (remember that our ordering is k 1 ≤ k 2 ≤ 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: where and we used the notation , where As shown in app.B, we can perform the angular integrals and find where β(∆ q ) = 1/2 if q 1 , q 2 , q 3 form a folded triangle, β(∆ q ) = 1 for all other (closed) triangles, and β(∆ 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: and . Whenever q 3 can not satisfy this triangle inequality, we drop this configuration.As such, we can drop the β(∆ q ) function inside the integral.We perform a change of variable q 3 → cos(θ 12 ) ≡ (q 2 3 − q 2 1 − q 2 2 )/(2q 1 q 2 ), such that the integral measure becomes q 1 q 2 q 3 dq 1 dq 2 dq 3 → q 2 1 q 2 2 dq 1 dq 2 d cos(θ 12 ), We then discretize the integration domain in 6 evenly-spaced points in q 1 , 6 in q 2 , and 4 in cos(θ 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 • 6 • 4 • 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 ∼ 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. 10They are defined, as described in [116], by the following averages: As expected from the size of those terms and the size of the binning effect (∼ 1σ), we have checked that properly binning the loop instead of evaluating them on these effective wavenumbers lead to negligible shift in the min χ 2 and in the posteriors for the analyses presented in this work.

Likelihood
To analyze the data, we start from a Gaussian likelihood, which is multiplied by the prior to arrive at the Bayesian posterior P over cosmological and bias parameters: 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 P 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 α , we can write where T α G,i and T N G,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 α .We will also choose a Gaussian prior on them, with covariance σ αβ and mean ĝα . 11Collecting the powers of g α , we can write the posterior in the following form: where the F 's are defined as: where Π is a generic prior on the cosmological and bias parameters non analytically marginalized.In other words, we assume that P pr is a sum of a Gaussian prior over the g α and a remaining prior on the other parameters.After integrating the g α , we have the marginalized posterior: Prior.In our analysis, we vary the cosmological parameters ω cdm , h, and ln (10 10 A s ) with a flat uninformative prior, while we use a Gaussian prior on ω b of mean ω b,BBN = 0.02233 and standard deviation σ 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. 12he EFT parameters should instead be restricted to be 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 NL ≃ k M ≃ 0.7h Mpc −1 , where k NL is the nonlinear scale for the matter field, and k 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 NL,R ≃ k NL / √ 8.We therefore write down each operator in the EFT expansion with either a k M or a k 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 π,1 , c πv,1 , and c St 2 , 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 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 NL = k M = 0.7h Mpc −1 and n = 4 • 10 −4 (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 ), with one b (i) 1 for each sky.The prior on it is a multivariate lognormal with correlation matrix: where standard deviation of the difference is ϵ = 2(1 − ρ).Our choices of ϵ 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 α 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 ≤ 0.02.The plots and summary statistics are calculated with the GetDist [122] package.

Pipeline validation
For our analyses, we use the following scale cut: k min = 0.01h Mpc −1 for all observables, k max = 0.23h Mpc −1 for P ℓ and B 0 , and k max = 0.08h Mpc −1 for B 2 , on CMASS.For LOWZ instead, we use k max = 0.20h Mpc −1 for P ℓ and B 0 , following [4].We keep k max = 0.08h Mpc −1 for B 2 on LOWZ. 13In this section, we perform multiple checks to validate our method at this scale cut.

Measuring and fixing phase-space effects
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 13 For comparison purpose, we sometimes fit B 0 using the tree-level prediction instead.When doing so, B 0 is denoted B tree 0 (to distinguish from B 1loop 0 ) and is fitted up to k max = 0.08h Mpc −1 for both CMASS and LOWZ.a covariance with BOSS volume V BOSS or rescaled to a large volume ∼ 100V 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 ℓ + 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 σ proj /σ 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.
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 χ 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 σ 8 and Ω 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 phasespace 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 ℓ + B 0 (+B , respectively for 1 sky and 4 skies. 14As 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 14 When analyzing the power spectrum multipoles P ℓ alone, we put the following prior instead:  are negligibly small (≲ 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.

Scale cut from NNLO
A simulation-independent way to evaluate the theoretical error as a function of k 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: and for the bispectrum: where k NL,R = k NL / √ 8, as discussed in sec.4. The prefactors c r,4 , c r,6 , c NNLO,1 , and c NNLO,2 are given a Gaussian prior centered on zero and of width 2. We then analyze the data as a function of k max , and determine the maximum wavenumber by taking the largest k 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 • σ.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 max in this way, we check the effect of these NNLO terms close to the k 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.

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σ on the posterior of Ω 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 (≲ 0.25) change in the χ 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.

Tests against simulations
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  For each skycut, we detail the number of bins N bin = N P ℓ bin + N B0 bin + N B2 bin , 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.
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.

Results
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 r,h (2,3) .In fig.2, we show the best fit residuals and in tab. 4 the best-fit χ 2 and associated pvalue.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  Given the position space halo overdensity δ h , the transformation to redshift space, δ r,h , is accomplished by where ẑ is the line-of-sight direction.We then expand the overdensity perturbatively as 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 δ (1) in terms of the symmetric n-th order halo kernels K r,h n defined as in 15   δ For example, is the famous Kaiser result for linear theory.The explicit expressions for K r,h 2 and K r,h 3 are given in [40], while the expression for K r,h 4 is available in [90] and can be straightforwardly derived from the above expressions.We provide the dependence of the kernels K r,h 1,...,4 on the bias parameters {b i } for i = 1, . . ., 15 here for convenience .(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 and the one-loop contributions 15 We have introduced the following notation ⃗ k1,..., ⃗ kn Additionally, P 11 is the dark-matter linear power spectrum.
For the bispectrum, the quantities that enter eq. ( 4) are the tree-level bispectrum and the one-loop contributions Note that on the left-hand sides we have dropped the argument (k 1 , k 2 , k 3 , k1 • ẑ, k2 • ẑ) 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 δ (1) , and stochastic terms which contain randomly fluctuating fields, typically denoted with an 'ϵ.' We can write the response terms as where and with The e St i functions are given below in sec.A.2. Notice that e St 3 = 0, so that there is no c St 3 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 r,h,ct 2 in eq. ( 61) are given by e In the above, we have used the notation µ i ≡ ki • ẑ.

C Additional parameter posteriors
In fig.6, we show the full triangle plots obtained fitting BOSS 4 skies P ℓ + B 0 + B 2 .In tab. 5, we show the 68%-credible intervals of b 1 , 2 , and c 4 obtained on this same fit.

( 69 )Figure 2 :
Figure 2: Measurements and best fits of bispectrum monopole B 0 (top), power spectrum multipoles P ℓ (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 ℓ + B 0 + B 2 on all BOSS 4 skies.

Figure 3 :
Figure 3: Triangle plots of base cosmological parameters obtained fitting synthetic data analyzed using

Figure 4 :
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.

and fig. 5 .Figure 5 :
Figure 5: Triangle plots of base cosmological parameters measured from the analysis of BOSS power spectrum multipoles P ℓ , ℓ = 0, 2, at one-loop, and bispectrum monopole B 0 at tree or one-loop level.

2 Figure 6 :
Figure6: Full triangle plots from the analysis of BOSS power spectrum multipoles P ℓ at one loop, bispectrum monopole B 0 at tree level or one loop, and bispectrum quadrupole B 2 at tree level.

Table 1 :
and we choose ϵ 12 = 0.1, ϵ 13 = 0.2.This formula is motivated by the fact that two variables distributed according to a bivariate normal with correlation ρ, the Residual deviations σ proj after phase-space projection adjustment measured on synthetic data generated and fitted with our model P ℓ + 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 σ data stat .

Table 2 :
Relative shifts ∆ sys /σ 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.

Table 3 :
Report of systematic errors on base cosmological parameters measured from the Nseries and sim represents the uncertainty from the simulation cosmic variance, which corresponds to about 0.15 or 0.03 in σ data stat for N sim = 84 Nseries or N sim = 2048 Patchy realizations, respectively.

Table 4 :
Goodness of fit given by the maximal log-likelihood value log L ≡ − min χ 2 /2 obtained fitting BOSS 4 skies P ℓ + B 0 + B 2 , and associated p-value.