Cosmology with the Redshift-Space Galaxy Bispectrum Monopole at One-Loop Order

We study the cosmological information content of the redshift-space galaxy bispectrum monopole at one-loop order in perturbation theory. We incorporate all effects necessary for comparison to data: fourth-order galaxy bias, infrared resummation (accounting for the non-linear evolution of baryon acoustic oscillations), ultraviolet counterterms, non-linear redshift-space distortions, stochastic contributions, projection, and binning effects. The model is implemented using FFTLog, and validated with the PT Challenge suite of $N$-body simulations, whose large volume allows for high-precision tests. Focusing on the mass fluctuation amplitude, $\sigma_8$, and galaxy bias parameters, we find that including one-loop corrections allow us to significantly extend the range of scales over which the bispectrum can be modeled, and greatly tightens constraints on bias parameters. However, this does not lead to noticeable improvements in the $\sigma_8$ errorbar due to the necessary marginalization over a large number of nuisance parameters with conservative priors. Analyzing a BOSS-volume likelihood, we find that the addition of the one-loop bispectrum may lead to improvements on primordial non-Gaussianity constraints by $\lesssim 30\%$ and on $\sigma_8$ by $\approx 10\%$, though we caution that this requires pushing the analysis to short scales where the galaxy bias parameters may not be correctly recovered; this may lead to biases in the recovered parameter values. We conclude that restrictive priors from simulations or higher-order statistics such as the bispectrum multipoles will be needed in order to realize the full information content of the galaxy bispectrum.

Whilst the above references have been pivotal to the development of a bispectrum model, few contain all the necessary ingredients to allow for robust comparison of theory and observation. In particular, one must account for the backreaction of short-scale physics on the large-scale bispectrum [66,67,85], long-wavelength displacements [86][87][88][89][90][91][92][93], and survey geometry [94,95], all of which can lead to biases in derived parameters if not properly accounted for. In [80] a complete model for the tree-level (leading-order) bispectrum of galaxies in redshift-space was presented and validated, including all the above effects (see also [19]). This allows for precise modelling of the angle-averaged bispectrum monopole, and has facilitated a number of analyses constraining ΛCDM parameters [15] and primordial non-Gaussianity [96,97]. However, this model was restricted to relatively large scales (k < 0.08 h −1 Mpc at z = 0.61). If we wish to further exploit the constraining power of the bispectrum, we must push to smaller scales, by extending the perturbation theory to next order. Whilst [98] has recently demonstrated some work in this direction, a full model for the one-loop bispectrum (including all relevant phenomena such as projection effects) has not yet been presented and validated with simulations.
In this work, we present a complete and systematic computation of the redshift-space galaxy bispectrum monopole at one-loop order. This includes all effects necessary to compare with observational data: deterministic contributions, counterterms, bias renormalization, stochasticity, bin-averaging, and coordinate distortions. This involves the galaxy density at fourth-order: we systematically account for all bias operators (following [79]), and include full treatment of all necessary redshift-space counterterms, ensuring a convergent Taylor series. Our model necessarily depends on a number of free parameters: these account for the unknown complexities of ultraviolet physics (such as galaxy formation physics and feedback), and ensure physical robustness. Efficient computation of the one-loop bispectrum is non-trivial; as such, we devote a significant portion of this work to discussing its practical computation with the FFTLog algorithm [99]. We compare the theoretical predictions to real-and redshift-space bispectra obtained from the PT Challenge simulations [100], which serve both to validate the approach and to assess the information content of the one-loop bispectrum model. Though we restrict to the measurement of σ 8 and primordial non-Gaussianity parameters, one can constrain a variety of other phenomena with the bispectrum, and, further still, our methodology can be extended to other correlators including the bispectrum multipoles [65,101] and the recently-detected trispectrum [102,103].
The remainder of this paper is structured as follows. The theoretical model is presented in §2, before its implementation is outlined in §3. In §4 we give details of the data and analysis choices used to validate the model, before presenting the results of likelihood analyses using the real-and redshift-space galaxy bispectrum in §5 & 6 respectively. §7 comments on the method's applicability to current datasets, with a summary and discussion given in §8. Finally, various technical details are presented in the Appendices, including: A the perturbation theory kernels, B details of the bispectrum integration routines, C discussion of the redshift-space counterterms, and D derivation of the stochastic bispectrum components. Appendix E is devoted to prior volume effects. The key plots of this work are Fig. 1, showing the one-loop bispectrum components, and Fig. 4, displaying the utility of the bispectrum for a BOSS-like survey.

Theoretical Model for the One-Loop Bispectrum
In this work, we analyze the power spectrum and bispectrum of biased tracers (i.e. galaxies) in redshift space at one-loop order. Whilst the one-loop power spectrum and tree-level bispectrum have been described in detail before [e.g. 19,80], a complete model for the oneloop bispectrum has not been presented before (though some aspects can be found in [98]) and will be discussed below, with additional technical details found in the appendices. Here, we will restrict to Gaussian initial conditions; extension to primordial non-Gaussianity is discussed in §5. 2. In the EFTofLSS, the bispectrum is comprised of the following terms at one-loop order [e.g., 66,85,104]: where the first and second terms give the tree-level and one-loop bispectrum in Eulerian perturbation theory, B ct is the derivative and counterterm contribution, and B stoch encodes stochasticity. This is strictly a function of five variables: three lengths, {k 1 , k 2 , k 3 } and two angles, {µ 1 , µ 2 }, for µ i ≡k i ·n with line-of-sightn (hereafter LoS), noting that k 1 µ 1 + k 2 µ 2 + k 3 µ 3 = 0. In real-space, this reduces to just three variables: {k 1 , k 2 , k 3 }.

Bias Expansion
To compute the bispectrum within Eulerian perturbation theory, our first step is to express real-space galaxy density field, δ g , in terms of a basis of bias operators, i.e. all combinations of the density and velocity fields (δ and θ) consistent with the relevant symmetries up to a given order in perturbation theory [79,[105][106][107][108][109][110]. For the one-loop bispectrum, we require terms up to fourth-order (δ 4 L ), and here use the basis of Galileon operators proposed in [79]: where curly brackets separate operators of different order and the bias parameters are marked in color. In (2.2), we drop any terms that do not appear in the one-loop bispectrum; these are all composite local evolution operators such as δ 4 and δ 2 G 2 (Φ v ). Here we have ignored both higher-derivative operators (which we return to below) and bias renormalization, which is discussed in Appendix B. The Galileon operators are defined by where Φ v ≡ ∇ −2 θ is the velocity potential, equal to the Newtonian potential Φ ≡ ∇ −2 δ at leading order. These can be simply generalized to functions of multiple potentials, with (2.2) involving the LPT potentials ϕ 1,2 , satisfying Up to third order, this is equivalent to the bias expansion used in [106] and previous works [e.g., 10,80], with the relations 1 Utilizing (2.2), and expanding each operator in terms of the linear density field δ (1) ≡ δ L , we can define the n-th order contributions to the galaxy density field: where the real-space kernels K n are given in Appendix A.1 and depend on the bias parameters given above. Furthermore, this generalizes to the redshift-space density field, δ s (k), using the well-known mapping [e.g., 111] where f is the logarithmic growth rate, u z (q) = (iµ q /q)θ(q) is the Fourier-space LoS velocity field, and µ q ≡q ·n, for LoS vectorn. The associated kernels, analogous to (2.6), are labelled as Z n and defined for n ≤ 4 in Appendix A.2.

Deterministic Contributions
Utilizing the redshift-space kernels of Appendix A.2, the tree-level bispectrum, B 211 ≡ δ can be written where P L (k) is the linear power spectrum (though see the below discussion on infrared resummation). This depends on the bias parameters {b 1 , b 2 , γ 2 }, as well as the growth rate, f (z). The one-loop terms can be written as loop integrals over the linear power spectrum, and come in four flavors [e.g., 66,104]: × P L (q)P L (|k 1 + q|)P L (|k 2 − q|), × P L (q)P L (|k 2 − q|) + 5 perm., where the B II 321 spectrum is similar to the P 13 (k 1 ) contribution to the one-loop power spectrum. Computation of the loop integrals can be performed via explicit numerical integration or with the FFTLog method [99]; we discuss the latter in §3, with details presented in Appendix B. As well as the tree-level biases, these spectra involve the higherorder parameters {b 3 , γ × 2 , γ 3 , γ 21 , γ × 21 , γ 211 , γ 22 , γ 31 }, of which only γ 21 appears in the oneloop power spectrum.

Counterterms
To ensure a self-consistent theoretical model, we require a set of counterterms, which account for non-idealities in fluid equations (via the viscous stress tensor), and absorb the unknown ultraviolet (UV, q k) behavior of the loop integrals in (2.9) [e.g., 66,67,85]. For the one-loop bispectrum in real-space, these operators are degenerate with derivative operators in the bias expansion, such as ∇ 2 δ. Furthermore, the redshift-space bispectrum contains additional counterterms that appear after the renormalization of contract operators in the perturbative mapping of (2.7); these are discussed in detail in Appendix C.
The overall bispectrum counterterm contribution can be written as where F ctr 2 (k 1 , k 2 ) is the real-space counterterm kernel [79]: 2), and we choose the non-linear scale k NL = 0.45h Mpc −1 [44,112,113]. (2.10) additionally involves the µ-dependent redshift-space kernel Z ctr 2 defined by (2.12) as derived in Appendix C, with k 3z ≡ k 3 µ 3 . In principle, two combinations of C 1 , C 2 and C 5 are constrained by the power spectrum, so only one parameter out of three is independent here. In practice, however, we did not find any difference between imposing the power spectrum constraints on C 1 , C 2 , C 5 or treating them as free parameters. This is why we proceed with keeping them free in what follows. In total, the one-loop bispectrum counterterm depends on 14 free parameters, {β B,i }, and {C i } in addition to the one-loop power spectrum counterterms. 2 Notably, many of the counterterms appearing in (2.12) are degenerate at the bispectrum monopole level; nevertheless, we prefer to keep all of them in the model, and marginalize over them within physically motivated priors. This is done for two main reasons. First, terms with different powers of µ can, in principle, be distinguished even at the bispectrum monopole level thanks to the Alcock-Paczynski projection effect [114], which is described below. Second, the degeneracy between these terms can be broken with higher order angular multipole moments of the bispectrum [37,41], which we will analyze in the future.

Stochasticity
Contributions to the bispectrum are also sourced by the non-deterministic part of the density field [106][107][108][109][110], i.e. that uncorrelated with δ L . At tree-level, this gives two terms, ∝ 1/n, P (k)/n (arising from Poissonian shot-noise with sample densityn), whilst at oneloop order, we must keep contributions suppressed by (k/k halo ) 2 , where k −1 halo is some characteristic halo size. From [79], we have the following form at next-to-leading order in real-space: depending on another five free parameters {A shot,0 , B shot } and {A shot,1 , S 0 , S 1 }, which cannot be constrained with the one-loop power spectrum. In the Poisson limit, A shot,0 = B shot = 1, with all higher-order terms (arising, for example, from halo exclusion) vanish.
In redshift-space, significantly more dependencies arise. A systematic derivation of these is presented in Appendix D and yields the following expression: (2.14) where This expression shares the parameter P shot with the power spectrum, but includes an additional 12 nuisance coefficients: {{S n }, {A shot,n }, B shot }. P shot is defined as a constant rescaling of the stochastic power spectrum [80], Note that in the absence of projection effects the counterterms A shot,1 and A shot,2 are fully degenerate. Therefore, for the purposes of this study we will set A shot,2 = 0.

Infrared Resummation
An additional complication arises from the effects of long-wavelength displacements, which can be consistently treated using "infrared resummation". A rigorous derivation of this was presented in [86,87] in the context of time-sliced perturbation theory [115], and, at tree-level order, can be implemented by replacing the linear power spectrum P L with its IR-resummed equivalent, i.e. 3 where P w and P nw are the wiggly and smooth parts of the power spectrum respectively. This has the effect of damping the oscillatory component by a k-and µ-dependent factor. The damping scales are given in terms of the broadband power spectrum as where r BAO is the sound-horizon scale and k S ∼ 0.1h Mpc −1 . At one-loop order, the IR-resummed bispectrum can be written schematically as where B[P ] indicates that the bispectrum should be evaluated using the power spectrum P and we have dropped the counterterms and stochasticity [86]. In this case, the loop corrections become more complex, since the damping factor, Σ 2 tot is a function of the redshift-space angles µ. To allow for efficient computation via the FFTLog procedure (Appendix B), we here adopt the isotropic approximation for the one-loop terms, dropping any µ-dependence in Σ 2 tot inside the integral. This is expected to be a good approximation in practice, and is exact for the real-space case. Note that we keep the full redshiftdependent damping function in the tree-level expressions, i.e. the isotropic templates are used only for the computations of the one-loop corrections.

Loop Integrals
We now discuss how to compute the one-loop bispectrum. The most difficult part of this is evaluating the loop integrals appearing in (2.9): in this work, these are computed via the FFTLog procedure [99], the subtleties of which are described in Appendix B. In essence, the real-space computation proceeds by first writing the integration kernels (products of Z n ) as polynomials in k 2 i , q 2 , and |k i ± q| 2 (or their reciprocals). By expanding the linear (or IRresummed) power spectrum as a sum over complex polynomials, i.e. P L (k) ∼ c m k ν+iηm for frequency η m and FFTLog 'bias' ν, the various terms in (2.9) take the form (using B 222 as an example) for some complex ν i . The integral can be evaluated using techniques borrowed from quantum field theory, and reduces the calculation to a tensor multiplication, noting that all cosmological information is encoded within c m . In redshift space, the appearance of angles,q ·n inside the integral make this more challenging; however, it can be evaluated using similar tricks to the one-loop power spectrum [cf. 112]), as discussed in Appendix B. Following the above tricks, the bispectrum takes the following schematic form, again taking B 222 as an example: where the i index runs over all combinations of bias and f (z), denoted θ i . 4 Additionally, we have expanded in terms of the redshift-space angles {µ, χ ≡ 1 − µ 2 cos φ} (of which there are 47 non-trivial combinations); these are related to the µ i angles via: We adopt this basis rather than the more familiar choice of {µ 1 , µ 2 }, since it avoids pathologies for flattened triangles (whence k 1 ≈ k 2 + k 3 , and µ 1 ≈ −µ 2 ). The underlying shapes, B (i,j,k) , appearing in (3.2) are independent of both redshift-space angles and bias parameters, and depend only on the form of the linear power spectrum, k 1 , x = k 2 3 /k 2 1 and y = k 2 2 /k 2 1 , assuming k 1 ≥ k 2 ≥ k 3 . Two options arise for using the bispectrum templates B (i,j,k) in Monte Carlo Markov Chain (MCMC) analyses: (a) they may be computed once for a fixed linear power spectrum, (b) they may be computed as a tensor multiplication (cf. 3.1) at each step in the MCMC chain, feeding in the relevant linear power spectrum (and thus c m coefficients) each iteration. Whilst (b) is the approach usually adopted for the one-loop power spectrum, we will here adopt (a) for the one-loop bispectrum. This has the effect of fixing cosmology in the bispectrum loops (except for σ 8 , which acts as a global rescaling, modulo a small effect concerning the IR resummation amplitude, which we ignore in this work), and is chosen on computational grounds, since the size of the necessary FFTLog matrices becomes very large. 5 Explicitly, we compute the bispectrum templates, B (i,j,k) , for a grid of values of {x, y, k 1 } (treating flattened triangles with √ x + √ y = 1 separately to avoid divergences), then use these to construct a three-dimensional linear interpolator for each shape. The resulting bispectra have been compared to results from explicit (and computationally intensive) numerical integration for a range of values of bias and triangle shapes and found to be in excellent agreement. Full details of the above steps are given in Appendix B.3. We additionally publicly release all our analysis code: this can be found at GitHub.com/OliverPhilcox/OneLoopBispectrum.

Bin Integration
To robustly compare theory and data, we must integrate the model across some set of bins. Following [80], this is achieved via the integral: where B 0 is the bispectrum monopole of (2.21) and N 123 = 8πk 1 k 2 k 3 (∆k) 3 V 2 /(2π) 6 for bin center (k 1 , k 2 , k 3 ) and width ∆k. As in [80], this is strictly exact only in the narrow-bin limit, and can be corrected by "discreteness weights" as in the former work. In practice, we compute the set of bispectrum templates B (i,j,k) (k 1 , x, y) for a range of values of k 1 , x, y (see Appendix B.3) then perform the bin-averaging by linearly interpolating these values, dropping any triangles that do not satisfy the triangle conditions The integration is performed using Gauss-Legendre quadrature, as for the angular integrals. Finally, we note that we can perform bin integration either within the MCMC chains or as a pre-processing step (allowing us to use bin-averaged templates in the later analysis). We use the latter option for the purposes of this paper. 5 To see this, note that the matrix in (3.1) has size N 3 freq , for N freq FFTLog frequencies. Taking N freq = 64, with 47 angular combinations, O(50) bias parameter combinations, and computing the matrix for 10 choices of each of x and y (noting that k scales out), we find ∼ 5 × 10 10 elements, or ∼ 50 GB in (complex) single precision.

Free Parameters
Our full model for the one-loop galaxy power spectrum and bispectrum depends on the following 44 free parameters (i.e. Wilson coefficients): where parameters appearing only in the power spectrum (following the definitions of [80]), only in the bispectrum, and in both spectra, are shown in blue, black and purple respectively. The three lines give bias parameters, UV counterterms, and stochasticity parameters respectively. Note that here we switch to the power spectrum biases b G 2 and b Γ 3 instead of γ 2 and γ 21 to ease the comparison with previous works [80,116]; these are related via (2.5). Whilst performing an MCMC analysis in this high-dimensional space may seem a formidable task, we note that all parameters except {b 1 , b 2 , b G 2 } enter the theory model linearly, and can thus be analytically marginalized, following [117]. This is exact, and will be applied to all analyses presented in this work, significantly reducing computational cost. Since the parameter b Γ 3 is of physical interest in power spectrum analyses, we opt to marginalize over this explicitly, alongside the quadratic biases. For the purposes of the analytic marginalization, we assume the following priors on the bispectrum nuisance parameters: all means are zeros, and the expectation values given by 10 for all bias parameters, 10 for all real-space counterterms and one-loop stochastic contributions, 20 for redshift-space counterterms (in order to account for enhancements caused by short-scale non-linear redshift-space distortions, known as fingers-of-God [118]), and 20 for redshift-space one-loop stochastic contributions. For the tree-level stochastic counterterms, following [80] we assume standard deviations of 5 for the dimensionless B stoch , P shot , and A shot parameters. The power spectrum nuisance priors match [80,116]. Note that our nuisance parameters are normalized in such a way that their physical values are expected to be O(1) numbers from the naturalness arguments. In this sense our physically-motivated choice of nuisance parameter is conservative, as we allow them to be as large as O(10).

Numerical Results
Before proceeding to use the one-loop bispectra to perform parameter inference, we first consider the form of the spectra themselves. Plotting the bispectrum is a challenge itself: the monopole exists in the three-dimensional simplex of {k 1 , k 2 , k 3 }, and we have contributions from a wide variety of nuisance parameter combinations. For the purpose of visualization, we will fix the bias parameters to simple local-in-Lagrangian space predictions, based on [119]: assuming the bias to be described only by linear and quadratic terms In Fig. 1 we plot the deterministic (Eulerian PT) bispectrum contributions assuming the above bias relations with b L 1 = 1, b L 2 = 0.3 and f (z) = 0.7, as well as distortion parameters α = α ⊥ = 1 and the best-fit PT Challenge input power spectrum (cf. §4). For both equilateral and squeezed triangles we observe a similar form: the one-loop corrections are suppressed on large scales (by k/k NL ) but become large as k increases, with the B I 321 piece exceeding tree-level theory by k ∼ 0.1h Mpc −1 . We find significant cancellation between the various one-loop components (which all depend on the same biases), which is expected from the IR cancellation of loop integrals. Note that the high-k behavior is further modified by the counterterms (scaling as k 2 P 2 L (k)) and stochasticity (scaling as k 0 and P L (k) at leading-order). The individual shapes of the bispectrum components are generally non-trivial, with oscillatory signatures seen in B 411 and, to a lesser extent, B I,II 321 . The smooth nature of B 222 (expected since the three power spectra are all inside the q integral) implies that a smaller number of FFTLog frequencies can likely be used in its computation, which may expedite the template computation, and suggests that this has only weak cosmology dependence. From the deterministic contributions alone, it is clear that the one-loop bispectrum is a significant fraction of B tree for all k 0.1h Mpc −1 , and thus its inclusion is necessary if we wish to model the bispectrum beyond the softest modes.

Data and Analysis Details
The dataset used in this paper is the PT Challenge suite [100], comprising high-resolution N -body simulations at z = 0.61 with a total volume of 566 h −3 Gpc 3 . Galaxies are allocated via a BOSS-like halo occupation prescription, and various summary statistics computed using a fiducial cosmology with Ω m = 0.3. In all our analyses, we use the redshift-space power spectrum multipoles, P (k), and the real-space power spectrum proxy Q 0 , both of which were studied in detail in [116]. In this work, we additionally add the bispectra in both real-and redshift-space, with the comparison allowing us to assess the relative importance of redshift-space distortions in the one-loop bispectrum.
The relevant bispectra are computed as described in [80], which studied the tree-level bispectrum likelihood. We bin the bispectrum data in wavenumber bins of width ∆k = 0.01h Mpc −1 , and use only triangles whose bin centers satisfy momentum conservation. 6 For k max = 0.15, 0.2, 0.3, we find a total of 372, 825, and 2600 independent triangle configurations, respectively, and note that, unlike [80], we do not include the very first bin in the analysis, i.e. we fix k min = 0.01h Mpc −1 for the bispectrum. This matches Dashed lines indicate negative contributions, and we show results for two types of triangle: equilateral, with k 1 = k 2 = k 3 , and squeezed, with k 2 = 0.9k 1 and k 3 = 0.2k 1 . For illustration, we assume coevolution biases following [119], with Lagrangian biases b L 1 = 1, b L 2 = 0.3 and a growth factor f (z) = 0.7. We do not include the contributions from stochasticity or counterterms in this plot, but note that all bias operators have been renormalized. the analyses of the actual surveys like BOSS, where the very first bin is often affected by systematics including stellar contamination [2,15].
Our theory model for the power spectrum matches that of [80,116], and we make use of the publicly available code class-pt [112] to compute the power spectrum models. 7 Similarly, our theoretical model for the bispectrum is discussed in detail in Section 2, and implemented using the FFTLog prescription using Mathematica -we refer the reader to Appendix B for technical details.
An important part of the likelihoods are the covariance matrices, encoding both errors and correlations. As in previous works, we here adopt the Gaussian tree-level approximation for the analytic covariance matrices of power spectra and bispectra, neglecting any cross-correlation between the two statistics. For sufficiently large scales these assumptions are well justified [24,63,65,80]; at smaller-scales, and in the presence of non-uniform survey geometry, a mock-based approach will probably be needed, such as in [15], most likely in combination with some compression scheme [e.g., 117].
The mock galaxy clustering data from the PT Challenge simulations are analyzed within the Bayesian framework. Here, we perform a global MCMC analysis using the publicly available sampler Montepython [120,121] varying the clustering amplitude σ 8 , f equil NL (the amplitude of equilateral primordial non-Gaussianity [96]) and the EFT nuisance parameters. Since the true value of σ 8 in the simulations remains blinded, we will show results only for the fractional error on σ 8 . As noted above, we will marginalize over all physical nuisance parameters given in §3. 3. This is in contrast with some bispectrum studies that aim to fix certain nuisance parameters, such as asserting coevolution relations for Lagrangian biases [81]. Indeed, for some particular purposes, i.e. fits of σ 8 , it may be sufficient to keep fewer parameters in the fit. However, such approximations are unwarrantedtheir validity can break down for other types of analyses. Therefore, we prefer to explicitly vary all physical nuisance parameters in the fit. By virtue of analytic marginalization [117], this is done at no computational cost.

Results: Real-space
We now present results from the above analyses, focusing first on the combination of the redshift-space power spectrum and real-space bispectrum. Though not quite matching observational setups (where the power spectrum and bispectrum are both observed in redshift-space), this analysis will allow us to understand the impact of redshift-space distortions.
To obtain the real space model, we set f = 0 in all calculations and retain only EFT operators that do not depend on the LoS angles, giving a one-loop model fully equivalent to that used in [79,81]. Our data vector contains the power spectrum multipoles, the real-space analog and the bispectrum monopole, i.e. [P (k), Q 0 (k), B(k 1 , k 2 , k 3 )] and we restrict to the z = 0.61 snapshot of the PT Challenge simulations. In most analyses we use P (k) up to k max = 0.16h Mpc −1 , and Q 0 in the range 0.16h Mpc −1 ≤ k < 0.40h Mpc −1 , as validated in [80,100,116]. We explore the impact of varying the bispectrum k max below.

Clustering amplitude and bias parameters
We first focus on measuring the mass clustering amplitude σ 8 and leading galaxy bias These appear both in the one-loop power spectrum and bispectrum models, and hence can be tightly constrained by the data. Unlike σ 8 , the true values of bias parameters in the simulations are unknown. As such, we take their best-fit values at a certain k max (where the one-loop model can be trusted) as a proxy for their true values. This k max is measured as in [100] (see also [113]) by determining at what scale cut posteriors for at least one parameter become biased w.r.t. analyses with lower k max . We take best-fit values of bias parameters at the last stable k max as ground truth, p true . Following this, our parameter measurements are quoted as ∆p = p − p true , to avoid unblinding the results.
We fit the real-space bispectrum data for the following choices of scale-cut: Posterior distributions for the clustering amplitude, σ 8 , and certain nuisance parameters extracted from MCMC analyses of the power spectrum multipoles and the oneloop real-space bispectrum. The power spectrum likelihood is the same for all cases, whilst we vary k B max for the bispectrum, as indicated in the caption. Corresponding marginalized parameter contours for k max = 0.21h Mpc −1 are given in Tab. 1. we find that the posterior on σ 8 remains unbiased up to k B max = 0.21h Mpc −1 , with a shift of 1 − 2% observed for k B max = 0.23h Mpc −1 and k B max = 0.25h Mpc −1 , which becomes significant relative to the PT Challenge error-bars.
However, at k B max ≥ 0.23h Mpc −1 , we observe that the nuisance parameters become biased w.r.t. measurements at low scale cuts, for example, we find a visible tension between the b 2 posterior at k B max = 0.23h Mpc −1 and k B max = 0.19h Mpc −1 . In addition, we see that parameters the b G 2 and b Γ 3 become biased. The optimal values of these parameters scale along the degeneracy direction b G 2 + 0.34b Γ 3 ≈ const, which closely matches the degeneracy combination imposed by the power spectrum, b G 2 + 0.4b Γ 3 [11].
In contrast with the k B max ≥ 0.21 h Mpc −1 picture, the results at all choices of k B max ≤ 0.21h Mpc −1 are fully consistent, implying that k B max = 0.21h Mpc −1 should be chosen as a baseline scale cut. This is somewhat larger than the one-loop power spectrum scale-cut of k max = 0.16h Mpc −1 ; whilst this might appear unusual, we note that the power spectrum contains significantly higher signal-to-noise, and is subject to redshift-space complexities, both of which decrease k P max . We use best-fit values of nuisance parameters from the baseline P + Q 0 + B(k B max = 0.21h Mpc −1 ) analysis as ground truth values in the below. P Parameter 68% limits 0.005 ± 0.050 Table 1. One-dimensional marginalized contraints on low-order bias parameters and the clustering amplitude σ 8 extracted from the PT Challenge data-set. We display results obtained using only the power spectrum multipoles P (left panel, cf. [100]), and those including the power spectrum, Q 0 and the one-loop real-space bispectrum likelihood with k B max = 0.21h Mpc −1 (right panel). The one-loop bispectrum is the main new feature of this work. Most parameters are normalized to their true values, to avoid unblinding the simulation. In real-space, the addition of the bispectrum significantly tightens posteriors on bias parameters (by at least an order of magnitude), and gives ≈ 20% improvement on σ 8 . Further details are given in the main text, with corresponding results for the redshift-space bispectrum shown in Tab. 2.
It is instructive to compare the parameter constraints extracted using the one-loop bispectrum to those from the power spectrum multipoles alone, i.e. P (k) at the baseline k max = 0.16h Mpc −1 . These are shown in the left panel of Tab. 1. We find an improvement of 31% in σ 8 , whilst the error-bars on bias parameters tighten by an order of magnitude in some cases. Despite the noticeable increase in signal-to-noise of the data-set, we find a modest improvement in σ 8 : this is linked to the proliferation of bias, counterterm, and stochasticity parameters needed to describe the one-loop bispectrum in an unbiased manner.
It is also useful to compare our one-loop bispectrum results with those from the treelevel bispectrum. We cannot directly use results from [80] since the former work also varied other cosmological parameters such as H 0 and Ω m . To obtain a cleaner comparison, we repeat the tree-level analysis of [80] with the same analysis settings as here, using k B max = 0.08h Mpc −1 for the tree-level bispectrum. We find ∆σ 8 /σ 8 = 0.002 ± 0.0053, i.e. an 21% improvement over the power spectrum only result. Comparing this the present analysis, we see that the addition of the one-loop bispectrum likelihood yields an extra 10% improvement over the tree-level bispectrum likelihood. Finally, we compare our results with those of [81]. Unlike our work, [81] used the power spectrum of halos and galaxies in real space, leading to the notorious b 1 − σ 8 degeneracy being largely unbroken. This explains why our results on σ 8 are much better -most of the constraining power comes from redshift-space distortions omitted in [81]. Despite this difference, our analysis does confirm a general trend pointed out in [81] -the returns from the one-loop bispectrum are limited by the large number of nuisance parameters. As such, it will be important to obtain better priors on them in the future, for example using hydrodynamical simulations.

Primordial non-Gaussianity
It is interesting to study to what extent the one-loop bispectrum model can help improve constraints on primordial non-Gaussianity (PNG), following constraints from the tree-level bispectrum in [96,97]. We consider here the case of equilateral PNG, which induces the following three-point correlation of the linear density field (see [96] for further details), where ζ is the primordial curvature fluctuation with dimensionless amplitude ∆ ζ , and we have introduced the transfer function T (k) = (P 11 (k)/P ζ (k)) 1/2 . Non-Gaussianity in the initial conditions generates three main effects [96]: (1) an additional contribution B 111 to the tree-level bispectrum, (2) an extra one-loop power spectrum correction P 12 , and (3) further contributions in the galaxy bias expansion, which modifies the tree-level expressions by introducing the so-called "scale-dependent" bias. The latter stems from the following expression, where "..." denote non-linear PNG corrections which can be ignored for the purposes of this paper. The term ∇ 2 ζ generates tree-level "scale-dependent" bias corrections to the power spectrum. Note that these corrections are suppressed in the equilateral case w.r.t. the case of local PNG, where the scale-dependent bias is a leading effect on the galaxy power spectrum, and thus the power spectrum dominates the constraining power on f loc NL . As shown in [96], for the one-loop power spectrum and tree-level bispectrum, we must include all three of the P 12 , B 111 and b ζ -related terms in our model. In this paper, we consider the one-loop bispectrum, which technically requires additional non-linear f equil NL corrections to the galaxy bispectrum, such as B 113 [110,122]. However, given that f eq NL ∆ ζ is small, these will be suppressed, thus we leave their systematic calculation for future work, focussing only on the leading terms, similar to [96,98]. For f equil NL 500, the next-to-leading order contributions are subdominant to two-loop matter corrections for k 0.2h Mpc −1 .
Concerning scale-cuts, we find that use of the P and Q 0 statistics at high k max can lead to biases in the recovered values of f equil NL . This is consistent with the estimates of [96], which showed that the two-loop corrections can actually be larger than the non-Gaussian P 12 contribution at small scales. Thus, we choose k max = 0.2h Mpc −1 for the Q 0 statistics and k max = 0.14h Mpc −1 for P in the PNG analysis of this section. For B real we use the baseline data cut k B max = 0.21h Mpc −1 , motivated by the discussion above. To perform the analyses including PNG, we fit the parameter f equil NL in addition to σ 8 and nuisance parameters. Since the PT challenge simulations were run using purely Gaussian initial conditions, we expect to find f equil NL consistent with zero. Indeed, our nominal constraint on the amplitude of the equilateral shape is given by We stress that these results are obtained without any external priors on σ 8 or the non-linear bias coefficients. This constraint can be compared with that obtained from the tree-level real-space bispectrum likelihood, At face value, our results imply that the addition of the one-loop bispectrum can improve constraints on f equil NL by ∼ 30%.

Results: Redshift-Space
In this section we present the analysis of the data-vector [P (k), Q 0 (k), B 0 (k 1 , k 2 , k 3 )], where all statistics are in redshift space and include projection and coordinate-distortion effects. This set-up thus fully matches an analysis of a realistic galaxy survey such as BOSS [2]. As for the power spectrum, we expect that the addition of redshift-space distortions (particularly the fingers-of-God effect [118], hereafter FoG), will reduce the non-linear scale, thus it is likely that k B max , and the constraining power of the bispectrum monopoole, will decrease.

Clustering amplitude and bias parameters
Let us discuss the recovery of the mass clustering amplitude σ 8 and bias parameters. As a point of comparison, we fix the fiducial bias parameters to those extracted from the real space analysis before with k B max = 0.21h Mpc −1 . We fit the redshift-space bispectrum data for the following choices of scale-cut: P +Q 0 +B 0 (k B max = 0.15 hMpc 1 ) P +Q 0 +B 0 (k B max = 0.17 hMpc 1 ) P +Q 0 +B 0 (k B max = 0.20 hMpc 1 ) P +Q 0 +B 0 (k B max = 0.22 hMpc 1 ) P +Q 0 +B real (k B max = 0.21 hMpc 1 ) Figure 3. Posterior distributions of the clustering amplitude and certain nuisance parameters obtained from MCMC analyses of the power spectra and one-loop redshift-space bispectrum monopole B 0 . The power spectrum likelihood is the same for all cases. We show results for different values of the bispectrum data cut k B max , as indicated by the caption. This is analogous to Fig. 2 (whose optimal constraint is shown by the purple curve), but utilizes the redshift-space bispectrum. Corresponding marginalized parameter constraints with k B max = 0.15h Mpc −1 are given in Tab. 2.
the bispectrum data is not sufficient to constrain all the nuisance parameters entering the theory model. This gives rise to significant marginalization projection effects, which can be naïvely interpreted as a bias in our model. We study these effects in Section E and show that the measurements at k B max < 0.15h Mpc −1 are consistent with our baseline choice k B max = 0.15h Mpc −1 once projection effects are taken into account.
0.09 ± 0.10 Table 2. One-dimensional marginalized constraints on low-order bias parameters and the clustering amplitude σ 8 extracted from the PT Challenge data-set. We show the fit from the combined likelihood including power spectrum multipoles, Q 0 , the tree-level redshift-space bispectrum at k B max = 0.06h Mpc −1 (left table) and the one-loop redshift-space bispectrum at k B max = 0.15h Mpc −1 (right table). The parameters are normalized relative to their true values. Whilst we find significant enhancements in the bias parameter constraints compared to the power spectrum alone (cf. Tab. 1), the constraint on σ 8 does not improve appreciably.
At k B max = 0.17h Mpc −1 the clustering amplitude σ 8 , the rescaled linear bias b 1 σ 8 , b 2 and b Γ 3 become biased w.r.t. their optimal values coming from the real-space bispectrum analysis. These biases are accompanied with a significant increase in the χ 2 statistics. Thus, we conclude that the two-loop bispectrum corrections are not negligible at this scale. This is further supported by the bias growing with k B max : in particular, at k B max = 0.22h Mpc −1 the bias on σ 8 reaches 2%, which is significant in the context of the PT Challenge simulation volume. In conclusion, we find that the one-loop galaxy bispectrum model in redshift space works well up to k B max = 0.15h Mpc −1 for the precision that corresponds to the total volume of the PT Challenge simulation (which we recall is significantly larger than current and forthcoming datasets). The optimal values of cosmological and bias parameters for this choice are presented in Tab. 2. For comparison, we also show the results from the tree-level bispectrum analysis akin to [80]. 8 That k B max is lower in redshift-space than real-space is no surprise: this indicates that the characteristic scale of FoG effects (σ FoG ) is smaller than that of non-linearities (k −1 NL ). For the power spectrum in redshift-space, higher-order counterterms were important to model FoG, scaling as k 4 z . An analogous set of nuisance parameters may be included here, but we caution that their number is large due to the higher dimensionality of the bispectrum.
Considering the marginalized posteriors directly, we find that the one-loop bispectrum likelihood (at k B max = 0.15h Mpc −1 ) yields only a 12% improvement on σ 8 compared to the power spectrum alone, though the constraints on bias parameters (and thus astrophysics) improve markedly. Comparing this with the tree-level case, we see that the inclusion of oneloop corrections actually lead to a somewhat worse result than for tree-level bispectrum,  Posterior distributions of the clustering amplitude and low-order nuisance parameters from MCMC analyses of the power spectrum data (P + Q 0 , in green), and the combination of the power spectrum and redshift-space bispectrum monopole (P +Q 0 +B 0 ), using tree-level (blue, from [80]) and one-loop (red) theory. The covariance is rescaled to match the volume of the BOSS survey, and we assume k max = 0.2 h −1 Mpc for both the one-loop power spectrum and bispectrum.
which tightens the σ 8 constraint by 18% for the analysis settings adopted in this work. This is consistent with previous studies considering the real-space bispectrum [81], and arises primarily due to the large number of nuisance parameters appearing in the one-loop calculation, especially in redshift-space. A similar situation takes place in the context of the one-loop redshift-space power spectrum, whose information content is limited by marginalization over nuisance parameters [24]. We will discuss this issue in detail later.

Implications for the BOSS survey
In this section we estimate the potential performance of the one-loop bispectrum model applied to the BOSS survey data [2], which is the largest publicly-available spectroscopic galaxy clustering dataset. This survey has significantly smaller volume than our mock simulation data, so one can expect that the analysis can be pushed to smaller scales [16,80]. Indeed, the relevant parameter in this problem is the ratio of the theory systematic bias in a certain parameter to the statistical error on that parameter. For the BOSS volume the statistical errors are significantly larger than the PT Challenge simulation volume, due to a ratio of volumes of ≈ 100. As the theoretical errors do not depend on the volume, the ratio between the theoretical error and statistical errors thus becomes smaller, and hence any residual theoretical systematics becomes less sizable in relative terms.
To demonstrate this, we repeat the likelihood analysis above for the redshift-space data vector [P , Q 0 , B 0 ], but rescale the covariance to match the BOSS volume V BOSS = V PT Challenge /100 6 (h −1 Gpc) 3 . We select k max = 0.20h Mpc −1 for power spectrum multipoles P and k B max = 0.20h Mpc −1 for the bispectrum monopole B 0 ; significantly larger than that found in §6. The results of this analysis are presented in Fig. 4 Table 3. One-dimensional marginalized constraints on low-order bias parameters and the clustering amplitude σ 8 extracted from the PT Challenge dataset, with covariance adjusted to match the volume of the BOSS survey. We show results for the P + Q 0 only analysis (top left), the tree-level P + Q 0 + B 0 likelihood (top right), and the one-loop P + Q 0 + B 0 likelihood (bottom). The inclusion of the bispectrum sharpens constraints on σ 8 by ≈ 24%, with some ≈ 10% improvement arising from the addition of the one-loop contributions.
We see that in the context of BOSS, the addition of the one-loop bispectrum yields an ≈ 24% improvement over the power spectrum-only result and an ≈ 10% improvement over the tree-level bispectrum likelihood result. However, this leads to a noticeable shift in nuisance parameters, with b G 2 approximately 1.7σ from its fiducial value. This could simply be a prior-volume effect however (since the effect of the priors becomes more important at lower simulation volume), especially given that b G 2 departs from its fiducial value at 1.1σ already for the power-spectrum alone. The tree-level bispectrum analysis, however, results in an unbiased recovery of all nuisance and cosmological parameters.
It is also instructive to study whether the one-loop bispectrum can improve constraints on equilateral PNG. Incorporating this parameter in the analysis as before (varying both f equil NL and σ 8 ) we find f equil NL = 197 ± 350. For the clustering amplitude we find ∆σ 8 /σ 8 = −0.026±0.035, with a slight 0.6σ shift w.r.t. the ground truth. For comparison, we have also run an analysis using the tree-level bispectrum at k B max = 0.08h Mpc −1 instead of the oneloop bispectrum, and found f equil NL = 420±440, ∆σ 8 /σ 8 = −0.025±0.040. First, we see some bias in f equil NL , which can be attributed to prior volume effects and somewhat more optimistic data cuts for the power spectrum that we use in our analysis here. Indeed, in [96] it was shown that the tree-level model yields unbiased results on f equil NL for k B max = 0.08h Mpc −1 and k P max = 0.17h Mpc −1 . Second, we notice that the bound on f equil NL in the one-loop case is 30% better than that of the tree-level analysis. The improvement is quite modest as a consequence of the fact that the one-loop model introduces many nuisance parameters, which cannot be constrained by the data. In our analysis we use highly conservative but still physically-motivated priors; if more aggressive priors on nuisance parameters are used, the constraints are likely to improve further.
In conclusion, we note that the addition of the one-loop bispectrum may yield some ≈ 30% improvement on the amplitude of equilateral PNG. We stress, however, that this comes with two important caveats. First, the k B max used for this study (0.2h Mpc −1 ) results in noticeable biases on the nuisance parameters, suggesting that the errorbar on f equil NL may be underestimated due to over-fitting. Whether this induces a bias on f equil NL is unclear; such an error would likely show up only in the analysis of simulations containing f equil NL = 0. Secondly, we have neglected the PNG-induced one-loop corrections to the bispectrum (as in [98]), which can be marginally important for the scales of interest (particularly in the tails of the f equil NL posterior), as can be easily estimated with the scaling universe approximation outlined in [122].
Finally, we note that our analysis indicates a more modest improvement on f equil NL than that reported in [98]. [98] suggest that the one-loop bispectrum improves f equil NL constraints over the tree-level result by a factor of few. This follows from a comparison with the treelevel bispectrum analysis of [96]. This comparison is misleading, however, since the baseline analysis of [96] varies cosmology whilst [98] always keeps cosmological parameters fixed. We have checked that this accounts for most of the difference between [96] and [98]. A more detailed comparison with [98] is not currently possible because the former work has not yet presented sufficient details about their analysis and theory model. It will be interesting to compare our results with their analysis in the future.

Conclusions and Discussion
In this work, we have presented and validated a complete calculation of the galaxy bispectrum monopole in redshift space at one-loop order in effective field theory. Our model includes one-loop corrections due to mode coupling, as well as the full set of EFT counterterms that are needed to regulate the UV behavior of loop integrals and capture the physical effects of backreaction of short scales onto the large-scale modes. Furthermore, we incorporate a bias expansion up to fourth order (noting that many operators van-ish after renormalization of the power spectrum and bispectrum) as well as fourth order redshift-space distortions. In addition, our calculation includes IR resummation to capture the non-linear evolution of baryon acoustic oscillations (both for the power spectrum and bispectrum), as well as projection and binning effects. In short, we include all relevant ingredients needed to compare theory with observational galaxy clustering data.
We have studied the performance of the one-loop bispectrum model in terms of cosmological parameter constraints, focusing primarily on the mass fluctuation amplitude σ 8 . To validate our model we use the PT Challenge simulation suite [100], which are equivalent to a BOSS-like survey with a hundred times larger volume, thus allowing for high-precision tests. We analyze a data vector that consists of the standard redshift-space power spectrum multipoles, the real-space power spectrum proxy Q 0 [116], and the redshift-space bispectrum monopole. In this setup, we have found that the inclusion of the one-loop corrections allows us to extend the agreement between bispectrum theory and data up to k max = 0.15h Mpc −1 , or k max = 0.21h Mpc −1 in real-space. This can be contrasted with the tree-level model bispectrum model, which works only up to k max = 0.08h Mpc −1 [80]. We caution that these scale-cuts depend on both the survey volume and galaxy type: for BOSS, we can use k max = 0.20h Mpc −1 , and it is likely that the wavenumber reach is larger for emission line galaxies, which boast smaller FoG effects [12,64]. Further, one might hope to extend the k-reach by specializing to some real-space bispectrum analog (similar to Q 0 ) at high-k: this will be considered in future work.
Despite a significant extension of the k-space reach, we have not found the bispectrum to lead to noticeable improvements in the σ 8 constraints compared to those with obtained from tree-level theory when applied to the PT Challenge simulations. This is a consequence of the large number of the EFT nuisance parameters that appear in the one-loop calculation (particularly in redshift-space), and must be marginalized over in our analysis. For a BOSSvolume survey (and accompanying systematic error thresholds), we find greater utility, with the one-loop bispectrum improving constraints by ∼ 10% over the tree-level case, though it remains to be seen whether any accompanying shifts in nuisance parameters are real (and malignant) or just prior volume effects.
We have additionally studied whether the one-loop bispectrum can help constrain equilateral primordial non-Gaussianity (and thus single-field inflation), finding that, for the BOSS survey, the one-loop bispectrum may improve constraints on the non-Gaussianity parameter f equil NL by ≈ 30% compared to the tree-level theory. Achieving this, however, requires pushing the bispectrum analysis to k B max = 0.2h Mpc −1 , where the shifts in the bias parameters become evident. It remains to be seen if this problem can be alleviated with better priors on nuisance parameters or with one-loop PNG-induced corrections to the bispectrum, which were omitted in this study. If one is interested in astrophysics, the one-loop bispectrum is much more useful: we find a significant tightening in the posteriors of parameters such as linear and tidal bias compared to those with only tree-level theory.
An important conclusion from our study is that we need better knowledge of the EFT nuisance parameters if we wish to extract more cosmological information from the bispectrum. This can be done in several ways. First, one can include data from the higher-order angular moments of the redshift-space bispectrum [37]. Since these moments depend on the same set of parameters, their inclusion should tighten the EFT nuisance parameters posteriors, aiding determination of the cosmological parameters of interest. Second, one can constrain the EFT nuisance parameters with higher order statistics, such as the trispectrum, see e.g. [102,123] for work in this direction. Finally, one can obtain better priors on the extra nuisance parameters using high fidelity N-body or hydrodynamical simulations. A powerful route by which to acheive this involves EFT field level techniques, see e.g. [85,[124][125][126][127][128][129][130][131][132]. We plan to investigate these options in future work.
Though the one-loop bispectrum analysis of this work was limited only to two cosmological parameters, the mass fluctuation amplitude and the equilateral non-Gaussianity parameter, it may be similarly extended to other parameters such as local primordial non-Gaussianity, or the neutrino mass. The improvement on other parameters, especially those beyond the minimal ΛCDM model, could be significantly larger, particularly when some new feature is introduced that is not degenerate with the smooth loop corrections. If the parameter of interest enters the theoretical model linearly, the analysis can proceed as above; if this is not the case, one would require an optimization of our one-loop bispectrum pipeline, since the FFTLog-based approach does not currently allow for a fast re-calculation of the theoretical template as the power spectrum is varied. If only the α , α ⊥ parameters are varied however, the templates do not need to be recomputed, only rebinned (via 2.21). Analyses including such effects will be natural next steps in our research program.
Finally, we note that the bispectrum data offers novel probes of new physics. In particular, constructing the bispectrum from different tracers will allow one to probe the equivalence principle [53][54][55][56][57][58][59]. Such an analysis is complicated if one considers only the power spectrum since the effects sensitive to the equivalence principle appear there only at the one-loop order. In contrast, the cross-bispectrum of different kinds of tracers can be a sensible probe of the equivalence principle, whose violation would generate new bispectrum shapes that are not present in the ΛCDM model. This, in particular, will help one derive new constraints on the violation of Lorentz symmetry in the dark matter sector [133,134]. We leave this and other tests of new physics with the bispectrum for future work.
primarily performed on the Helios cluster at the Institute for Advanced Study, Princeton, and additional computations were carried out using the Princeton Research Computing resources at Princeton University, which is a consortium of groups led by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology's Research Computing Division.

A.2 Redshift-Space
The redshift-space kernels are obtained by expanding the RSD mapping of (2.7) and expanding all fields in terms of the linear density field. Following a lengthy computation, we find the following forms in terms of the real-space kernels: writing µ i···j ≡ µ q i +···+q j and dropping the argument of K 1 for clarity.

B Computation of the One-Loop Bispectrum with FFTLog
In this appendix, we discuss practical computation of the loop integrals given in (2.9). Before considering the redshift-space case, we will first examine how to compute the realspace integrals, which follow a similar logic, but are significantly simpler. Our approach follows [99], but is extended to the case of biased tracers and redshift-space.

B.1.1 Formalism
As noted in §2, the first step in the bispectrum computation is the expansion of the perturbation theory kernels (Appendix A) as polynomials in q 2 , |k 1 − q| 2 and |k 2 + q| 2 or their reciprocals (utilizing permutation symmetries). In practice, this results in a sum over many thousands of terms, once the relevant symmetries have been imposed, and is automated using mathematica. For B 222 , each term is proportional to q P L (q)P L (|k 1 − q|)P L (|k 2 + q|) q α 1 |k 1 − q| α 2 |k 2 + q| α 3 , (B.1) for integer α i , with a similar form found for the other loop integrals except with fewer factors of P L . Expanding the linear power spectrum as a sum of complex polynomials in k, i.e. P L (k) = m c m k ν+iηm for frequencies η m , coefficients c m , and (real) FFTLog bias ν (which sets the eventual integral convergence properties), we can rewrite (B.1) in the form for 2ν j = α j −ν −jη m i , where all the cosmology dependence (encoded in c m ) is now outside the integral. The remaining integral can be computed using path integral methods as where x 2 = k 3 /k 1 , y 2 = k 2 /k 1 and J (with complex arguments ν i ) can be expressed as a sum of hypergeometric functions and Gamma functions [99]. This reduces the computation of bispectrum templates to a set of matrix multiplications and function evaluations, as noted in §3. For B I 321 we find a similar form to (3.1), except with rank-two matrices, whilst B II 321 and B 411 involve only a one-dimensional sum (and one set of c m i coefficients).

B.1.2 Limiting Behavior
When computing spectra via FFTLog, it is important to verify whether the relevant loop integrals actually converge. This is achieved by taking the UV and IR limits of the integration kernels and assessing the dependence on the ν i parameters appearing in (B.2).
As an example, we consider the contribution of three δ operators to B 222 (involving three copies of b 1 F 2 (k − q, q)). This has the following limits in the equilateral configuration k 1 ∼ k 2 ∼ k 3 ∼ k: For P L (q) ∼ q ν , the integral is UV convergent for ν < 1 and IR convergent for ν > −1.
By choosing the bias in this range, FFTLog will give accurate values for the integrals. In contrast, if ν is chosen to be outside this range, we must add the relevant UV or IR limits by hand (taking care to include subleading divergences if necessary).
Considering all bias terms, the limits of B 222 take the following schematic form for equilateral triangles: where {f i } are some polynomials, and we consider only the leading-order contribution for each bias parameter. Inserting P L (q) ∼ q ν as before shows that a term containing K powers of b 2 is UV convergent for ν < 1 − 2K/3, implying ν < −1 for b 3 2 , significantly tighter than the ν < 1 limit for matter (i.e. b 3 1 ). However, the UV limit of b 3 2 is fully degenerate with the bispectrum shot-noise ( 0 in 2.13), and should be subtracted off in practice, as for the b 2 2 contribution to P 22 . If we adopt ν > −1, this term will not be captured by the FFTLog formalism, thus the subtraction becomes implicit. In this case, we require ν < −1/3 to avoid the b 2 2 divergence (and the second-order b 3 2 divergence, both of which are degenerate with the 2 stochasticity in 2.13). In the IR, (B.5) shows that the integral is convergent for ν > −1 for terms involving two or more powers of b 1 , and ν > −3 else. To satisfy all the conditions simultaneously, we may take −1 < ν < −1/3, dropping the shot-noise piece.
For B I 321 , the limiting UV and IR form is given by where X ∈ {b 2 , b 3 , γ × 2 } and ellipses are taken to mean bias operators excluding X (in the UV) or b 1 (in the IR). (B.6) implies that UV divergences can be avoided if we take ν < (1 − 2K)/2 when the term involves K powers of b 2 , b 3 , or γ × 2 ; these are all the composite operators appearing at third order. Furthermore, as in B 222 , the UV limits of the terms involving two powers of b 2 , b 3 and γ × 2 are proportional to shot-noise (this time of the η 0 variety in 2.13), 9 and should be subtracted off in practice (or dropped implicitly by fixing ν > −3/2). In the IR, divergences vanish for ν > −1 for terms involving b 3 1 or b 2 1 and ν > −3 else. Overall, we require a bias of −1 < ν < −1/2 to satisfy all conditions, assuming subtraction of the η 0 shot-noise contributions. For B II 321 , we require the UV and IR convergence properties of q Z 3 (k, q, −q)P L (q), which we labelP 13 (k) by analogy with the galaxy power spectrum (2.9). This natively involves all bias operators in (2.2) up to third order; however, this set is reduced to just {δ, G 2 (Φ v ), G 2 (ϕ 2 , ϕ 1 )} when the renormalization conditions are applied. These conditions demand lim for renormalized operator X and linear density field δ L (k), i.e. there can be no loop contributions which do not decay in the UV limit [106]. The contribution of all composite operators (e.g., δ 2 ) toP 13 is exactly that of a non-decaying loop diagram (since there is no suppression by the F 3 kernel), thus must vanish when the operators are properly renormalized. This leaves only G 3 , which evaluates to zero after averaging over the angular part of q. Following these redefinitions (which do not affect B 222 and B I 321 ), we find the UV and IR limits:P UV divergences occur unless ν < −1, and IR divergences occur unless ν > −1 (for b 1 ) or ν > −3 (else). As for P 13 [99], there is no range of biases which satisfy all the conditions; in this case, we can choose −1 < ν < 1 (satisfying the IR limits, and avoiding subleading UV divergences at ν > 1), and correct the UV part by adding the relevant limit by hand, which takes the following explicit form in real-space: proportional to the velocity divergence σ 2 v . Finally, we consider B 411 . This contains the fourth-order bias operators, and involves Wick contractions of linear density fields within the same operator, permitting simplification via the renormalization condition: which is proportional to the UV limit of B 411 . The first effect of this is to remove contributions from any fourth-order composite local evolution operator (such as δ 4 or δ G 3 ); these operators were already dropped from the bias expansion in (2.2). Secondly, this will remove a number of UV divergences in the below. Before bias renormalization, the UV and IR limits of the remaining terms take the form: where ellipses represent additional bias terms which impatience lead us to ignore. The first line is UV convergent for ν < −3 (first term, involving composite operators) or ν < −1 (second term, no composite operators). However, the first term possesses a UV limit that does not decay as an (negative) integer power of q 2 , violating the renormalization condition (B.10). The precise action of bias operator renormalization is to remove such terms (and only these, as far as this diagram is concerned). By evaluating the diagram with ν > −1, such contributions will be avoided, i.e. the operators will be correctly renormalized. In the IR, we find that divergences can be avoided by setting ν > −1 (for terms involving the first and second order operators proportional to b 1 , γ 2 or b 2 ), or ν > −3 (for the remaining terms). As forP 13 , there is no single bias that will simultaneously remove all the UV and IR divergences in B 411 , even after bias renormalization. Fixing −1 < ν < 1, we may compute the full expression by manually adding the appropriate UV limit to the FFTLog result. These limits can be computed straightforwardly from the kernels in mathematica and are omitted from this publication to avoid unnecessary tedium.

B.2 Redshift-Space
In redshift-space the perturbation theory kernels depend not only on the lengths q, |k 1 − q| and |k 2 + q| but also the LoS angles µ i ≡k i ·n andq ·n. 10 Although we are primarily interested only in the bispectrum monopole (i.e. that integrated over µ 1,2 , with a suitable Lebesque measure, as in 2.20), the full dependence on µ i is necessary for accurate calculation of coordinate distortions (2.21), thus we cannot simply average over µ i before computing the loop integrals; furthermore, this is difficult to perform analytically due to the presence of high powers ofq ·n. After expanding the kernels as polynomials, we will find loop integrals of the form: for n ∈ {0, 1, · · · , 6} (cf. 3.1), with prefactors depending on µ i , k i , biases and f (z). Below, we consider how to compute this utilizing the FFTLog procedure, generalizing the approach of [112] for the power spectrum. First, we expand theq ·n angles as Cartesian sums, i.e. 3 i=1q in i , and pull out the LoS vectors from the integral. The remaining function is a fully symmetric rank-n tensor, given by This has dependence only on k 1 and k 2 ; as such, its tensorial dependence can be written in terms of the components of k 1 , k 2 , and any isotropic tensors of relevance, i.e. the Kronecker delta. 11 Explicitly, this takes the form: where {O k } is the set of all independent symmetric rank-n combinations ofk i 1 ,k i 2 and δ ij K . As an example, the n = 2 operators are {δ ij We then define an "overlap matrix", giving the correlation between basis elements: (assuming Einstein summation conventions); this allows extraction of the A k coefficients via A k = [I] −1 kk O k i 1 ···in F i 1 ···in , where the second term is just the contraction of (B.14) with various powers ofk 1 andk 2 . Finally, we contract (B.14) with n copies ofn i to yield where the first set of parentheses contains a set of (k 1,2 ·q coefficients inside the q integral, and the second contains powers of µ 1,2 . To make this explicit, we give the n = 1 case: writing ν 12 ≡k 1 ·k 2 . In this manner, the FFTLog integral can be performed for arbitrarily large n. We adopt this method to compute the bispectrum templates in redshift space, applying it as a simplification step before the loop integrals are computed as in Appendix B.1. Notably, the above decomposition breaks down in the limit ofk 1 ·k 2 → −1, i.e. for k 1 k 2 , whence there is only one angle in the problem. The corresponds to flattened triangles (with k 1 = k 2 + k 3 or √ x + √ y = 1), which contain the divergence 1 − ν 2 12 → 0. Strictly speaking, this divergence is cancelled by the numerators, oncek 1 = −k 2 is identified; however, if one separately computes the loop integral coefficients proportional to powers of µ 1 and µ 2 , numerical issues will arise. In this limit, we adopt a different angular decomposition, noting that (B.13) can depend only onk 1 and the Kronecker delta. The basis tensors are much simpler in this case, for example, with {δ ij K ,k i 1k j 1 } for n = 2, and facilitate computation in an analogous manner to the above. For n = 1, (B.17) becomes which does not diverge. This divergence also illustrates the importance of expanding the bispectrum templates in the {µ, χ} basis (with µ ≡ µ 1 , χ ≡ 1 − µ 2 cos φ) rather than {µ 1 , µ 2 } (cf. §3): the former is undefined for flattened triangles, whence µ 1 = −µ 2 , whilst the latter simply has dependence only on µ in this limit (noting that µ 2 = µ ν 12 − 1 − ν 2 12 χ → −µ as ν 12 → −1).

B.3 Implementation
The above tricks allow us to efficiently compute the one-loop bispectrum in redshift space. A rough overview of the computation is the following: 1. Expand the relevant (symmetrized) perturbation theory kernels as polynomials in q, |k i ± q|, µ i andq ·n.
5. Compute the bispectrum templates for each of the 47 combinations of µ i χ j and the relevant combinations of bias parameters using the FFTLog algorithm. This is performed for a grid of values of k 1 , x, y, with flattened templates (obeying √ x+ √ y = 1) computed separately, using the alternate angular decomposition given in Appendix B.2, and involving only 7 non-trivial powers of µ.
6. Create a three-dimensional linear interpolator for each template using the precomputed bispectrum shapes (combining full and flattened configurations).
8. Compute the full bispectrum as a sum over templates, weighted by the bias configurations and any necessary discreteness weights.
Notably, only steps (5) and beyond depend on the power spectrum template, and thus the cosmological survey in question. We note one further subtlety: computing the bispectrum templates near (but not at) the flattened limit of √ x + √ y = 1 can lead to numerical issues due to large values of 1/(1 − ν 2 12 ), which appear in the angular decompositions of (q ·n) raised to the n-th power (cf. B.17). To counter this, when the templates are being computed and the condition √ x + √ y < 1.1 is met, we replace the FFTLog prefactor by its Taylor series in (1 + ν 12 ), artificially removing the divergent terms (which are present only due to numerical inaccuracies). Mathematica and Python code implementing all of the above steps is publicly available at GitHub.com/OliverPhilcox/OneLoopBispectrum. Following initial testing of the FFTLog routines against explicit numerical integration for a small number of bins, we use the following choices of FFTLog bias: ν = −0. 6   √ y, subject to the triangle conditions. This is sufficient to ensure that the spectra are subpercent accurate in the regime of interest; the results are largely unchanged if the number of FFTLog frequencies is reduced by a factor of two. The computation requires ∼ 10 4 CPU-hours to compute all templates (entirely performed within Mathematica), with the majority of time devoted to B 222 , and could certainly be optimized further. Calculations have been compared against explicit numerical integration of the (unsimplified) bispectrum kernels, and we find excellent sub-percent agreement in all cases.

C Counterterms from Redshift-Space Distortions
The RSD mapping to O(δ 4 1 ) can be obtained by expanding (2.7) to fourth order: To facilitate renormalization, we must smooth this expansion with a low-pass filter of some size R = Λ −1 . Products of fields at the same point (contact terms) are sensitive to short-scaled modes and hence must to be smoothed and renormalized. We denote these operations by square brackets, [...] R . Galilean symmetry implies the following schematic structure of the renormalized correlators (see [135][136][137] for the first order results), where u i , δ are the smoothed long-wavelength velocity and density fields (for clarity, we will drop the subscript in the below). To preserve Galilean symmetry, the operators O should not depend on the smoothed velocity field.
Note that the velocity field scales like k −1 δ k at the linear order, i.e.

D Stochastic terms
In this section we discuss stochastic contributions to the one-loop bispectrum.

D.1 Real space
The stochastic contributions to the galaxy density field in real space are given in terms of by (some operators are present in [137]): (D.1) Note that δ = 0 by definition. There are two non-trivial possibilities to contract operators in δ to obtain the tree-level bispectrum contributions (with free coefficients shown in color): where primes denote that we drop the Dirac delta function. These match the operators present in the tree-level bispectrum model [80]. At the one-loop order we find three distinct contractions: (k 2 2 + k 2 3 )P (k 1 ) + cyc. n , O k 2n−2 : d 1 k 2 1 k 1 k 2 k 3 =d 1 A shot k 2 1 + cyc. n 2 .
In addition, there are purely stochastic terms that generate the bispectrum of the order k 2n−2 . These arise from the following combinations: After performing angular integration (in the absence of coordinate-distortion effects), this term takes the same form as the real-space term ∼ k 2n−2 . Thus, it does not produce a new contribution to the bispectrum monopole, though is important if higher-order multipoles are also considered.

E Prior volume effects
In this section we study the prior volume effects present in our posteriors when the one-loop bispectrum likelihood is analyzed with small data cuts, such as k B max = 0.12h Mpc −1 . At face value, the posterior distributions from this analysis are several σ away from the true values. However, here we show that as much as half of this shift can be explained by prior volume (marginalization projection) effects. Indeed, such effects are expected to be present when the data volume is not sufficient to tightly constrain model parameters, which is the case for analyses with low k B max . We performed the following test: rerunning our full analysis on the mock data generated by our fitting pipeline for the best-fit cosmology at k B max = 0.15 h Mpc −1 . This mock data is simply a theory curve without any statistical scatter. In the absence of prior volume effects our pipeline must exactly recover the input parameters. However, when we fit this mock bispectrum data at k B max = 0.12 h Mpc −1 , we find that the mean values recovered from our pipeline are shifted relative to the input values at the (1 − 1.5)σ level, as shown in Fig. 5. This is evidence of prior volume effects. Furthermore, the shifts are in the directions of the apparent biases observed in the actual data ( §6). Thus, if we subtract these shifts from the actual posteriors at k B max = 0.12, the mean posterior values would match the true input parameter values at least within the 99% CL. Finally, we note that the parameters b 2 , b G 2 and b Γ 3 are highly correlated; this means that a shift in one would induce a shift in both.
As an additional check, we repeat our mock analysis for k B max = 0.15h Mpc −1 . Overall, we find much improved agreement between the mock and actual analyses. The posteriors for σ 8 and b Γ 3 are still shifted with respect to the ground truth by 1σ (which is smaller than 1.5σ shifts in the k B max = 0.12h Mpc −1 case), but all other parameters are recovered without noticeable bias. P +Q 0 +B 0 (k B max = 0.12 hMpc 1 ) Mock, P +Q 0 +B 0 (k B max = 0.12 hMpc 1 ) Figure 5.
Posterior distributions of the clustering amplitude and low-order nuisance parameters from MCMC analyses of the power spectrum and bispectrum likelihoods from the redshift-space analysis at k B max = 0.12h Mpc −1 for the PT Challenge simulation data (in green) and for the mock bispectrum data vector (in gray) computed with our pipeline for the bestfit cosmology at k B max = 0.15h Mpc −1 . Dashed lines show the input values for the mock data-vector, whose discrepancies with the grey posteriors indicate clear evidence for prior volume effects.