# Perturbation Theory vs Simulation: Quasi-linear Scale, Binning Effect, and
Visualization of Bispectrum

Joseph Tomlinson1,2 and Donghui Jeong1,2,3
1 Department of Astronomy and Astrophysics, The Pennsylvania State University,
University Park, PA 16802, USA
3 Institute for Gravitation and the Cosmos, The Pennsylvania State University,
University Park, PA 16802, USA
3 School of Physics, Korea Institute for Advanced Study, 85 Hoegiro,
Dongdaemun-gu, Seoul, 02455, Republic of Korea E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Recently, Perturbation Theory (PT), specifically the Effective Field Theory of
Large Scale Structure (EFTofLSS) and its equivalents, have proven powerful in
analyzing observational data. To further this pursuit, we present a
quantitative analysis for the accuracy of PT modeling by comparing its
analytical prediction to the result from a suite of Quijote simulations.
Specifically, we determine $k_{\rm NL}$, the wavenunmber below which the
analytical prediction matches well with the N-body result, for both leading
order (LO) and next-to-leading order (NLO) power spectrum and bispectrum at
redshifts $z=0$, $0.5$, $1$, $2$, $3$. We also quantify the binning effect in
Fourier space and show that an appropriate correction must be applied to the
analytic predictions in order to compare them with the discrete Fourier
transform results obtained from N-body-simulation or real data. Finally, we
have devised a novel spherical-Bispectrum visualization scheme fully capturing
the scale and configuration dependences. The new scheme facilitates
bispectrum-amplitude comparison, for example, between theory and N-body
results.

###### keywords:

cosmology: large-scale structure of Universe - method: data analysis

††pubyear: 2022††pagerange: Perturbation Theory vs Simulation: Quasi-linear
Scale, Binning Effect, and Visualization of Bispectrum–References

##  1 Introduction

Upcoming large-scale structure (LSS) surveys (Doré et al., 2015; Hill et al.,
2008; Laureijs et al., 2011; Levi et al., 2013; Collaboration, 2012; Maartens
et al., 2015; Spergel et al., 2015; Takada et al., 2014) will allow us to
probe the clustering of galaxies on ever larger scales while measuring them on
existing scales with high precision. This increase in precision on small-scale
galaxy clustering measurements calls for more accurate models that have well
characterized systematics. There has been significant development of the
theoretical models in recent years, from fast methods for computing next-to-
leading and higher-order corrections (Schmittfull et al., 2016; McEwen et al.,
2016; Slepian, 2018; Schmittfull & Vlah, 2016; Simonović et al., 2018), to
developments of the galaxy bias models to bridge theory and observations
(Scoccimarro, 2004; Taruya et al., 2010; Zheng et al., 2010; Hirata, 2009;
Blazek et al., 2015; Eggemeier et al., 2019). See Desjacques et al. (2018a)
for a review, and Desjacques et al. (2018b); Tomlinson et al. (2020); Ivanov
et al. (2020) for models combining both.

While most of this work has been done on the power spectrum, this increase in
precision will also allow for higher-order statistics to be used, typically
the three-point function or the bispectrum. Some significant progress has been
made in modelling the bispectrum, the Fourier transform of the three-point
function, in a similar way to the power spectrum (Eggemeier et al., 2019;
Smith et al., 2008; Rampf & Wong, 2012; Assassi et al., 2014; Bernardeau et
al., 2012; Lazanu et al., 2016; Angulo et al., 2015; Eggemeier et al., 2021;
McCullagh et al., 2015). Using the bispectrum is important for properly
breaking degeneracies that are present in the two-point function or power
spectrum. The use of the bispectrum has been shown to increase constraints on
various cosmological parameters of importance, dark energy (Sefusatti et al.,
2006; Song et al., 2015; Byun et al., 2017), primordial non-Gaussianity
(Sefusatti et al., 2012; Tellarini et al., 2016; Karagiannis et al., 2020;
Dizgah et al., 2020), and neutrino mass (Chudaykin & Ivanov, 2019; Hahn et
al., 2020; Hahn & Villaescusa-Navarro, 2020; Kamalinejad & Slepian, 2020). The
bispectrum has been previously applied to data (Gil-Marín et al., 2016;
Slepian et al., 2017a, b; Cabass et al., 2022), but only very recently has
this been done on anything but large scales where it most well modeled. Though
recently Eggemeier et al. (2021) has shown results past the linear regime.

To exploit the full potential of the dataset, these models often go beyond
leading order (LO) to the next-to-leading order (NLO) in perturbation theory,
allowing for accurate modelling of smaller-scale clustering signatures; for a
review see Bernardeau et al. (2002). The exact scale at which NLO perturbation
theory is required is extremely important not only for accurately estimating
the information gained by including NLO corrections, but also for properly
keeping theoretical systematics under control. This has been studied to some
extent for the power spectrum (Jeong & Komatsu, 2006b; Nishimichi et al.,
2009; Taruya et al., 2009; Osato et al., 2019), but extensive measurements for
a variety of redshifts are required, since the range of validity is sensitive
to the redshift. This is especially true for the bispectrum. With accurate
knowledge of the difference between the applicable scales of the LO and NLO
models, along with knowledge of covariance matrices (e.g. Villaescusa-Navarro
et al. (2020)) to measure mode coupling at small scales which erases some of
the information gained from the extra modes at those scales, one could then
accurately predict the constraining power gained from using a NLO model over a
LO model.

This is the main purpose of this paper. In this paper, we shall measure
$k_{\rm NL}$, the maximum wavenumber below which perturbation theory provides
an accurate modeling of the nonlinear density field, by comparing the matter
power spectrum and bispectrum obtained from a suite of Quijote N-body
simulation Villaescusa-Navarro et al. (2020). We have defined $k_{\rm NL}$ by
requiring the cumulative deviation of perturbation theory prediction to stay
below 1 % and 2 % for the power spectrum, and 2 % and 5 % for the bispectrum.
We repeat the analysis for the LO prediction and NLO prediction to obtain
$k_{\rm NL}$ for each case. We also study the triangular-configuration
dependencies of $k_{\rm NL}$ for the nonlinear bispectrum.

This paper is organized as follows. We start in Sec. 2 with providing details
on the data we used from the Quijote N-body simulations and the basics of the
theory we used to match them, including binning effects. We then present our
new method for visualizing the bispectrum, which we denote the spherical
bispectrum, in Sec. 3. Finally our results for the range of validity for both
the power spectrum and the bispectrum, including configurations, at all
redshifts is included in Sec. 4. We conclude in Sec. 5.

Throughout the paper, we use the following conventions and shorthand
notations. Our definition of the bispectrum is

$\left<\delta(\bm{k}_{1})\delta(\bm{k}_{2})\delta(\bm{k}_{3})\right>=(2\pi)^{3}B(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})\delta^{D}(\bm{k}_{123})\,.$
(1)

Here, $\delta^{D}$ is the Dirac-delta operator, and we use the shorthanded
notation of $\bm{k}_{1\cdots n}\equiv\bm{k}_{1}+\cdots\bm{k}_{n}$. We denote
the amplitude of a vector $\bm{k}_{i}$ as $k_{i}$.

##  2 Data & Theory

###  2.1 Quijote Simulations

For our N-body data set we use a subset of the Quijote simulations
(Villaescusa-Navarro et al., 2020), a set of 43,100 N-body simulations with a
wide variety of cosmological parameters and results available at five
different redshifts ($z=0$, $0.5$, $1$, $2$, and $3$). For this study, we use
the high-resolution, fiducial-cosmology sample, a set of 100 simulations run
with 8 times as many dark matter particles. In detail, this sample of 100
simulations each has $1024^{3}$ dark matter particles distributed over a
volume of 1 Gpc3/$h^{3}$, giving a fundamental wavenumber $k_{f}\approx
0.0063$ $h$/Mpc, with a cosmology in agreement with Planck Collaboration et
al. (2020a). For more details about the simulations themselves, see
Villaescusa-Navarro et al. (2020).

We choose the high resolution simulations in particular because the higher
resolution gives better accuracy at small scales, where we need it to
determine the range of validity, at the cost of large scale error bars due to
only having 100 simulations. We find this tradeoff to be worthwhile since it
is already well known that LO and NLO perturbation theory are consistent with
the nonlinear density field on large scales, so large error bars there do not
impact the results of our analysis. In this work, we have devised a criteria
for defining $k_{\rm NL}$ without being affected by the deviation due to
large-scale cosmic variance. From these simulations we use the precomputed
power spectra (Villaescusa-Navarro et al., 2020) and then use the full dark-
matter particle distribution to measure the bispectrum, all of which exists
for five redshifts. The reason we remeasure the bispectrum is because the
precomputed bispectra are not minimally binned, with $\Delta k=2k_{f}$, so we
remeasure them with minimal binning, $\Delta k=k_{f}$ to ensure the most
accurate representation of the underlying data we can get and to minimize
binning effects.

An important part of measuring the bispectrum is distributing the N-body
particles onto a grid. For our grid decomposition we use nbodykit(Hand et al.,
2018) with a custom Julia wrapper. We also follow the recommendations in
Sefusatti et al. (2016) and use interlaced grids along with the third-order
triangular-shaped clouds (TSC) method to dramatically reduce aliasing effects.
For measuring the bispectrum from the decomposed grid we use the fast
Scoccimarro estimator (Scoccimarro, 2015; Sefusatti et al., 2016; Tomlinson et
al., 2019), specifically we use the code presented in Tomlinson et al. (2019)
with some efficiency improvements. We measure the bispectrum up to $k_{\rm
NL}=80k_{f}\approx 0.503$ $h$/Mpc to ensure we have sufficiently small enough
scales to properly capture the deviation from perturbation theory at high
redshift.

###  2.2 Eulerian Perturbation Theory

As for the theory model, we use the standard Eulerian perturbation theory(SPT)
model, see Bernardeau et al. (2002) for a review. In brief, SPT factors the
full nonlinear matter density into different parts

$\delta(\bm{k})=\delta_{L}(\bm{k})+\delta^{(2)}(\bm{k})+\delta^{(3)}(\bm{k})+\cdots\,,$
(2)

where $\delta^{(n)}(\bm{k})$ stands for the $n$-th order density contrast
proportional to the $n$ linear quantities such as density field, velocity
field, and tidal field. We can then use Wick’s theorem to find the leading
order (LO) and the next-to-leading order (NLO) contribution to the power
spectrum straightforwardly

$\displaystyle\left<\delta(\bm{k})\delta(\bm{k}^{\prime})\right>$
$\displaystyle=\left<\delta_{L}(\bm{k})\delta_{L}(\bm{k}^{\prime})\right>$
$\displaystyle+2\left<\delta_{L}(\bm{k})\delta^{(3)}(\bm{k}^{\prime})\right>+\left<\delta^{(2)}(\bm{k})\delta^{(2)}(\bm{k}^{\prime})\right>$
(3)

We write each term in the above equation as

$P_{\rm LO+NLO}(k)=P_{L}(k)+2P_{13}(k)+P_{22}(k)\,,$ (4)

where the leading order power spectrum is simply the linear power spectrum
$P_{\rm LO}(k)=P_{L}(k)$, and the NLO terms consist of

$\displaystyle P_{22}(k)$ $\displaystyle=2\int\frac{{\rm
d}^{3}\bm{q}}{(2\pi)^{3}}F_{2}(\bm{q},\bm{k}-\bm{q})^{2}P_{L}(q)P_{L}(|\bm{k}-\bm{q}|)$
$\displaystyle P_{13}(k)$ $\displaystyle=3P_{L}(k)\int\frac{{\rm
d}^{3}\bm{q}}{(2\pi)^{3}}F_{3}(\bm{q},-\bm{q},\bm{k})P_{L}(q)\,.$ (5)

Here, $F_{n}(\bm{q}_{1},\cdots,\bm{q}_{n})$ is the $n$-th order integration
kernel, see Bernardeau et al. (2002). For the form suitable for a faster
implementation of these integrals using the FFTlog-based method, see
Schmittfull et al. (2016); McEwen et al. (2016).

Since the initial density field is close to Gaussian (Bennett et al., 2013;
Planck Collaboration et al., 2020b), we assume that the linear-order density
field follows Gaussian statistics, and there is no linear expression for the
bispectrum. The leading-order bispectrum expression is

$\displaystyle B_{{\rm LO}}(k_{1},k_{2},k_{3})=$ $\displaystyle
B_{211}(k_{1},k_{2},k_{3})$ $\displaystyle=$
$\displaystyle\,2F_{2}(\bm{k}_{1},\bm{k}_{2})P_{L}(k_{1})P_{L}(k_{2})+(2~{}{\rm
cyclic})\,.$ (6)

Like the NLO power spectrum calculation, we can use Wick’s theorem to find the
NLO contribution to the bispectrum

$\displaystyle B_{\rm LO+NLO}(k_{1},k_{2},k_{3})=B_{\rm
LO}(k_{1},k_{2},k_{3})+B_{222}(k_{1},k_{2},k_{3})$
$\displaystyle+B_{411}(k_{1},k_{2},k_{3})+B_{123}^{I}(k_{1},k_{2},k_{3})+B_{123}^{II}(k_{1},k_{2},k_{3})\,,$
(7)

where

$\displaystyle B_{222}(k_{1},k_{2},k_{3})=\,8\int\frac{{\rm
d}^{3}\bm{q}}{(2\pi)^{3}}F_{2}(\bm{q},\bm{k}_{1}-\bm{q})F_{2}(-\bm{q},\bm{k}_{2}+\bm{q})$
$\displaystyle\quad\times
F_{2}(\bm{q}-\bm{k}_{1},-\bm{k}_{2}-\bm{q})P_{L}(q)P_{L}(|\bm{k}_{1}-\bm{q}|)P_{L}(|\bm{k}_{2}+\bm{q}|)\,,$
(8) $\displaystyle
B_{411}(k_{1},k_{2},k_{3})=\,12P_{L}(k_{1})P_{L}(k_{2})\int\frac{{\rm
d}^{3}\bm{q}}{(2\pi)^{3}}$ $\displaystyle\quad\times
F_{4}(\bm{q},-\bm{q},-\bm{k}_{1},-\bm{k}_{2})P_{L}(q)+({\rm 2\,cyclic})\,,$
(9) $\displaystyle
B_{123}^{I}(k_{1},k_{2},k_{3})=\,6P_{L}(k_{1})\int\frac{{\rm
d}^{3}\bm{q}}{(2\pi)^{3}}F_{2}(\bm{q},\bm{k}_{2}-\bm{q})$
$\displaystyle\quad\times
F_{3}(-\bm{q},\bm{q}-\bm{k}_{2},-\bm{k}_{1})P_{L}(q)P_{L}(|\bm{k}_{2}-\bm{q}|)+({\rm
5\,cyclic})\,,$ (10) $\displaystyle
B_{123}^{II}(k_{1},k_{2},k_{3})=\,6P_{L}(k_{1})P_{L}(k_{2})F_{2}(\bm{k}_{1},\bm{k}_{2})$
$\displaystyle\quad\times\int\frac{{\rm
d}^{3}\bm{q}}{(2\pi)^{3}}F_{3}(\bm{q},-\bm{q},\bm{k}_{2})P_{L}(q)+(\rm{5\,cyclic})\,.$
(11)

A fast way to compute these integrals exists (Simonović et al., 2018), but not
one using the standard formalism and has not been extended to redshift-space
or biased tracers.

These perturbation theory expressions are what we use throughout this work.
For the linear power spectrum we take the one provided with the Quijote
simulations generated using CAMB (Lewis et al., 2000). To calculate the NLO
power spectrum model we use an adaptation of the methods described in
Schmittfull et al. (2016). The LO bispectrum model is generated using the
provided linear power spectrum and Eq. (2.2). We generate our NLO bispectrum
model through direct integration.

###  2.3 Binned Theory Models

Figure 1: (Left): The binning correction on the LO + NLO power spectrum. The
binning correction is plotted for three plotting redshifts ($z=0$, $1$, $3$)
to emphasize the contributions originating from the NLO part of the model,
since only the NLO part has any redshift evolution. We see that binning
correction matters only for very large scales, as it stays below 1% in all
bins except the largest scale bin. (Right): The binning correction on the LO +
NLO bispectrum, as a function of $k_{\rm sph}$ that we define later in Sec. 3.
The visible binning correction is independent of redshift, signifying that the
NLO contribution is much smaller than the LO contribution, so we only plot one
redshift at $z=0$. We see significantly larger errors than the case of the
power spectrum with the largest bin reaching an error around 16 %.

Fourier-space clustering statistics measured from a density field defined at
grid points are inherently binned with the fundamental frequency
$k_{F}=2\pi/L$. On the other hand, analytic theory models, like those in Sec.
2.2, can be calculated at each wavenumber without binning. This causes a
discrepancy between the two even if the underlying statistic is identical, so
a correction is essential when comparing analytic theory models and binned
measurements. We refer to this as a binning correction. Note that theory
models that are inherently gridded, such as the result from GridSPT (Taruya et
al., 2018), do not suffer from this problem but instead they are statistical
and so require many realizations to reduce statistical errors.

To compute the binning correction we need to bin the analytic models the exact
same way that we bin the statistical measures from simulations or real data.
This is typically done by using an estimator, but, instead of taking data as
input, the theory-model calculation is used as input. We then quantify the
binning correction as the difference between the binned and unbinned theory
relative to the unbinned theory. In this section, we shall measure the binning
correction for the power spectrum and bispectrum.

For the power spectrum, the binning correction is typically very small except
on scales near the fundamental wavenumber. To bin the theoretical power
spectrum, we first generate a three-dimensional grid in Fourier space and
assign the LO and NLO power spectrum calculated by using Eq. (4) to each grid
point. We then measure the binned theory model for both LO and NLO power
spectrum the same binning scheme as the direct measurement:

$P_{\rm
bin}(k_{i})=\frac{1}{N_{k}}\sum\limits_{|\bm{k}-k_{i}|\leq\frac{1}{2}k_{F}}P(\bm{k})\,.$
(12)

Note that this is just the traditional direct power spectrum estimator but
with $|\delta(\bm{k})|^{2}$ replaced with $P(\bm{k})$. We show the binning-
correction error for the LO+NLO power spectrum in the left pane of Fig. 1. We
have plotted the errors for our three plotting redshifts ($z=0$, $1$, $3$) to
emphasize the contributions coming from the NLO part, those that evolve with
redshift, while the contributions that do not evolve with redshift are from
the LO part, the linear power spectrum. We see that the binning effect
contributes overwhelmingly more to the LO model, likely because it is only
significant for the very largest scales. Our results are broadly consistent
with the folklore that the binning correction only matters on very large
scales, with only the very largest bin having an error larger than 1%.
However, to use the large-scale bins for precision work, for example for
measuring the local-type primordial non-Gaussianity (Slosar et al., 2008;
Hamaus et al., 2011; Giannantonio et al., 2014; Agarwal et al., 2014; Barreira
et al., 2020; Moradinezhad Dizgah et al., 2021; Rezaie et al., 2021), an
accurate modeling of the binning correction is necessary.

For the bispectrum the binning effect is significantly more important than the
power spectrum. Bernardeau et al. (2012) and more recently Eggemeier et al.
(2021) have reported significant binning effects between a few and up to 10
percent and its triangular-configuration dependencies. The binning effect on
the bispectrum is computationally expensive to study, because we must sum over
all theory contribution on the grid:

$\displaystyle B_{\rm bin}(k_{1},k_{2},k_{3})=$
$\displaystyle\frac{1}{N_{123}}\sum\limits_{|\bm{k}_{i}-k_{1}|\leq\frac{1}{2}k_{F}}\sum\limits_{|\bm{k}_{j}-k_{2}|\leq\frac{1}{2}k_{F}}$
$\displaystyle\,\qquad\times\sum\limits_{|\bm{k}_{k}-k_{3}|\leq\frac{1}{2}k_{F}}B(\bm{k}_{i},\bm{k}_{j},\bm{k}_{k})\,.$
(13)

For example, for the bispectra considered in this work, it takes around a week
for each of LO and NLO bispectrum models to get the binned form.

The right panel of Fig. 1 shows the binning correction for the NLO bispectrum,
as a function of the $k_{\rm sph}$ that we define in Sec. 3. Unlike the power
spectrum we do not plot this for multiple redshifts, since the contribution
from the NLO part of the bispectrum, and hence the redshift evolution, is
indistinguishably small on this plot. This means that nearly all of the
binning effect, at least in $k_{\rm sph}$-space, is due to the LO part. When
plotted in $k_{\rm sph}$-space, we see the bispectrum display similar behavior
to the power spectrum, with significant binning errors on very large scale
modes, reaching a maximum error of 16% in the largest bin.

Although we use the computationally slow exact binning in this work, as we
only needed to compute it a few times, the long computation time of binning
the bispectrum is a clear obstacle for precision bispectrum analysis given how
significant the effect can be. A few alternative techniques have been
developed. Sefusatti et al. (2010) have developed a method by averaging over
the wavenumbers in a bin to generate an effective wavenumber

$\displaystyle k_{{\rm eff},i}(k_{1},k_{2},k_{3})=$
$\displaystyle\frac{1}{N_{123}}\sum\limits_{|\bm{k}_{i}-k_{1}|\leq\frac{1}{2}k_{F}}\sum\limits_{|\bm{k}_{j}-k_{2}|\leq\frac{1}{2}k_{F}}\sum\limits_{|\bm{k}_{k}-k_{3}|\leq\frac{1}{2}k_{F}}\bm{k}_{i}\,,$
(14)

and evaluating the theoretical-model bispectrum at that effective wavenumber.
Following Yankelevich & Porciani (2018), Oddo et al. (2020) improves the
method by sorting the wavenumbers, achieving a less than 5% binning error.
Recently, Eggemeier et al. (2021) has introduced a scheme which uses Delaunay
interpolation and constructs tetrahedra to reduce the computation time,
achieving around 1% accuracy when using full Delaunay binning.

Finally, a comment is in order for the interpolation scheme that we use to
compute the binned power spectrum and bispectrum. When computing the power
spectrum and the NLO bispectrum at the grid points, we use cubic-spline
interpolation to alleviate needing to compute the model at each grid point.
For the power spectrum, which is a very smooth function, this is unlikely to
cause any significant effects. For the bispectrum, however, our testing shows
this can have a significant effect. For example, we find $\sim$ 10% error on
binning correction at the largest scales when computing Eq. (13) with the
interpolated LO bispectrum instead of the exact LO bispectrum. Note, however,
that we only use interpolation for the NLO-bispectrum binning correction,
which is not very significant on any scales. Therefore, the interpolation
should not affect our results.

##  3 The Spherical Bispectrum

Figure 2: (Left): The LO and NLO bispectrum models plotted in the usual
flattened index space. The flattened index is generated from creating an
ordered set of tuples ($k_{1},k_{2},k_{3}$) with $k_{1}\geq k_{2}\geq k_{3}$
and incrementing in row-major order, that is incrementing $k_{3}$ first.
Binned with bin size $k_{f}$ (Right): The same two models plotted in spherical
bispectrum space. It is both much easier to see the differences between the
models and also more clear at exactly what scale they begin to diverge. Both
plots are at $z=0$.

The Bispectrum is a function of three semi-independent variables, $k_{1}$,
$k_{2}$, and $k_{3}$, connected by the triangle condition: $k_{a}+k_{b}\leq
k_{c}$ for all combination of $a\neq b\neq c$. While this does not hamper the
bispectrum’s ability as a cosmological probe, it is desirable to have a
visualization tool showing the full scale- and configuration- dependence and
facilitating the comparison among different models.

In literature, the most common method of visualizing the bispectrum is by
flattening the three input wavenumbers into a contiguous one-dimensional
flattened index each of which maps onto the triplet ($k_{1}$, $k_{2}$,
$k_{3}$). The flattened index plot is not smooth when $k_{1}$ and $k_{2}$
values vary so the plot features are dominated by indexing effects rather than
the nonlinearities. Also, it is non-trivial to assign a numeric scale to a
specific bispectrum element, since there are three separate scales that
contribute towards it. As a result, it is difficult to assign a single $k_{\rm
NL}$ by examining the flattened index plot.

Other methods, such as heatmap triangle plots (Jeong & Komatsu, 2009;
McCullagh et al., 2015), also exist although they can only be done by holding
one of the wavenumbers constant. Thus, while highlighting the configuration-
dependence of the bispectrum, the latter visualization methods are not
convenient for direct comparison among different bispectra, for example, from
theory calculation and simulation measurement.

To facilitate the comparison between different models with each other and
against data visually, here we introduce the spherical bispectrum, which is
essentially a coordinate transformation from the wavenumber triplet ($k_{1}$,
$k_{2}$, $k_{3}$) to the spherical coordinate system ($k_{\rm sph}$,
$\theta_{\rm sph}$, $\phi_{\rm sph}$). That is, by treating the wavenumbers of
the bispectrum as equivalent to a Cartesian coordinate system (i.e.
$k_{1}=k_{x}$, $k_{2}=k_{y}$, $k_{3}=k_{z}$), we can then convert this three-
dimensional coordinate system into spherical coordinates like so

$\displaystyle k_{\rm sph}=$
$\displaystyle\sqrt{k_{1}^{2}+k_{2}^{2}+k_{3}^{2}}/\sqrt{3}$
$\displaystyle\phi_{\rm sph}=$ $\displaystyle\arctan\left(k_{2}/k_{1}\right)$
$\displaystyle\theta_{\rm sph}=$
$\displaystyle\arctan\left(\sqrt{k_{1}^{2}+k_{2}^{2}}/k_{3}\right)\,.$ (15)

We then bin the bispectrum results $B(k_{\rm sph},\phi_{\rm sph},\theta_{\rm
sph})$ based on $k_{\rm sph}$ and compute a weighted average off all the modes
that fall into the same bin, weighted by the number of triangles that
contribute to each mode. Formally, we may write the spherical bispectrum as

$B_{\rm sph}(k_{\rm sph})=\frac{\sum\limits_{|q_{\rm sph}-k_{\rm
sph}|<\Delta/2}B(q_{1},q_{2},q_{3})N_{\rm
tri}(q_{1},q_{2},q_{3})}{\sum\limits_{|q_{\rm sph}-k_{\rm
sph}|<\Delta/2}N_{\rm tri}(q_{1},q_{2},q_{3})}\,,$ (16)

with the bin size $\Delta$. The spherical bispectrum allows us to get a good
measure of the bispectrum at a given scale, which is also easily plot-able in
a standard form identical to the power spectrum. The $k_{\rm sph}$ can be
thought of as a geometrically motivated way of defining a single effective
wavenumber for each wavenumber triplet $(k_{1},k_{2},k_{3})$.

In Fig. 2, we plot two theory curves, the leading order (LO) and next-to-
leading order (NLO) bispectrum in two different visualization schemes: the
traditional flattened index style (Left) and the new spherical bispectrum
style (Right). As shown in Fig. 2, we see significantly enhanced clarity with
regards to where the two models deviate when looking at them in $k_{\rm
sph}$-space compared to flattened-index space.

Note that if we only use the radial binning, the spherical bispectrum does
lose some information compared to using the full bispectrum.

We can differentiate, however, the bispectrum’s configuration dependence by
taking advantage of the angular information in the spherical coordinates. In
Sec. 3.1, we use six triangular configurations (five used in the previous
bispectrum analysis along with a new one), which are differentiated based on
angular cuts in the spherical coordinates, giving us easy visual access to
different parts of the angular domain.

###  3.1 Spherical Triangle Configurations

Figure 3: On the left is the formal definintions of each config in
$\theta_{\rm sph}-\phi_{\rm sph}$ space. On the right is our expanded
definitions for “-like" configs we use in this work. Here we have relaxed
equality conditions by 5% and consider any mode with $k_{1}/k_{3}>3$ to be
squeezed.

Probing specific range of angles in spherical coordinates allows us to look at
the bispectrum’s specific triangular configurations and to regain the
information that is lost in the radial binning procedure. These configuration
dependencies have been studied as ways to probe different physical mechanism
of generating the bispectrum because the bispectrum at each configuration
responds differently for different physical origins such as early Universe
physics, non-linear growth of structure, and nonlinear bias (Desjacques et
al., 2018a).

There are six primary configurations we study here, and they are formally
defined as follows:

  * •

Equilateral: $k_{1}=k_{2}=k_{3}$

  * •

Isosceles: $k_{1}=k_{2}>k_{3}$

  * •

Folded: $k_{1}=2k_{2}=2k_{3}$

  * •

Elongated: $k_{1}=k_{2}+k_{3}$

  * •

Squeezed: $k_{1}=k_{2}\gg k_{3}$

  * •

Obtuse-Isosceles (O-Isosceles): $k_{1}>k_{2}=k_{3}$

However, some of these configurations contain too few modes per bin and, even
with reasonable bin sizes, using their strict definition makes them very
stochastic. The problem is particularly severe for the Squeezed where the
equality holds up for only limited number of configurations.

We alleviate this problem, and also include more information, by expanding
their definitions while still only including modes that capture the same type
of mode ratios. When used in this way we refer to them as “-like"
configurations, for the remainder of this section, to differentiate from the
formal definitions given above. In Fig. 3, we show the different
configurations in $\theta_{\rm sph}-\phi_{\rm sph}$ space, both formally
(Left) and with our expanded definitions (Right). For the rest of the work we
drop the “-like" moniker and simply refer to the base name of the
configuration.

Specifically, we do not expand the definition of Equilateral and Folded modes,
as they are relatively smooth without the alteration. For every other mode, we
relax the condition of equality in their formal definition to instead be a
ratio lower limit. For example Isosceles modes formally impose $k_{1}=k_{2}$
but for our isosceles-like modes we impose $k_{2}/k_{1}>0.95$. We impose this
same relaxation on all “-like" modes. For the case of squeezed-like modes we
introduce one further relaxation, the minimum ratio of $k_{1}/k_{3}>3$. These
relaxations of the definitions not only smooth out the resulting functions but
also allow for more angular information to be captured while still having a
reasonable amount of overlap between different configurations.

###  3.2 Spherical Bispectrum Bin Completeness: $k_{\rm sph,max}$

An important consideration with the spherical bispectrum is that of bin
completeness. Since the spherical bispectrum is a combination of various
bispectrum elements with different ($k_{1}$, $k_{2}$, $k_{3}$) values, some
elements may be put into a bin where other elements would in theory go but
naively imposing a finite $k_{\rm max}$ could prevent those elements from
being added into the bin.

We can calculate the largest, in terms of scale, bin where this occurs by
finding the smallest bispectrum element which is not captured with a cut at
$k_{\rm max}$ and then compute which bin it would fall into. If we
characterize $k_{\rm max}$ with its integer form $n_{k,\rm max}=k_{\rm
max}/k_{F}$ ($k_{F}$ is the fundamental wavenumber) then depending on whether
it is even or odd there are different forms for $n_{\rm sph,max}$, the integer
form of the first incomplete bin, given by Tab. 1. There are a variety of
different limits depending on the configuration.

To summarize, with the maximum wavenumber of $n_{k,{\rm max}}$, the spherical
bispectrum is complete up to, but not including, $n_{\rm sph,max}$ indicated
in Tab. 1.

##  4 Result: The Range of Validity

|  $n_{\rm sph,max}$ for even values of $n_{k,\rm max}$ |  $n_{\rm sph,max}$ for odd values of $n_{k,\rm max}$
---|---|---
General | $\sqrt{n_{k,\rm max}(n_{k,\rm max}+2)+4/3}/\sqrt{2}$ | $(n_{k,\rm max}+1)/\sqrt{2}$
Folded | $(n_{k,\rm max}+2)/\sqrt{2}$ | $(n_{k,\rm max}+1)/\sqrt{2}$
Elongated | $\sqrt{n_{k,\rm max}(n_{k,\rm max}+2)+4/3}/\sqrt{2}$ | $(n_{k,\rm max}+1)/\sqrt{2}$
O-Isosceles | $\sqrt{n_{k,\rm max}(n_{k,\rm max}+8/3)+2}/\sqrt{2}$ | $(n_{k,\rm max}+1)/\sqrt{2}$
Equilateral | $(n_{k,\rm max}+1)$ | $(n_{k,\rm max}+1)$
Isosceles | $\sqrt{2n_{k,\rm max}(n_{k,\rm max}+2)+3}/\sqrt{3}$ | $\sqrt{2n_{k,\rm max}(n_{k,\rm max}+2)+3}/\sqrt{3}$
Squeezed | $\sqrt{2n_{k,\rm max}(n_{k,\rm max}+2)+3}/\sqrt{3}$ | $\sqrt{2n_{k,\rm max}(n_{k,\rm max}+2)+3}/\sqrt{3}$
Table 1: The value of $n_{\rm sph,max}=k_{\rm sph,max}/k_{f}$ for each
configuration when $n_{k,\rm max}$ is even/odd. This is calculated by finding
the smallest uncaptured mode for each configuration and then calculating its
$k_{\rm sph}$.

In this section, we present out results for the range of validity for both
power and bispectrum at both LO and NLO. First, let us define the range of
validity.

While the general consensus in literature is that the nonlinear scale, or
range of validity for the linear model, is where the NLO correction is a
sizable fraction of the LO contribution, the exact definition of the nonlinear
scale can vary throughout the literature. In this work, we use a more data-
driven approach and compute a $\chi^{2}$-like statistic to determine $k_{\rm
NL}$ where the simulation and theory models diverge beyond our accuracy
threshold. Specifically, for a clustering statistic $F$, e.g. the power
spectrum or bispectrum, we compute $n_{NL}=k_{NL}/k_{f}$ through

$\frac{1}{n_{NL}}\sum_{i=1}^{n_{NL}}\frac{(F_{{\rm theory},i}-F_{{\rm
sim},i})^{2}}{\sigma_{F_{i},\rm sim}^{2}+\sigma_{F_{i},\rm
theory}^{2}}<(1.5)^{2}\,,$ (17)

where $\sigma_{F_{i},\rm sim}$ is the statistical error on the mean of the
clustering statistic measured from the simulations at point $i$. The
$\sigma_{F_{i},\rm theory}$ in denominator is an induced statistical error on
the theory to take our accuracy threshold into account, defined as

$\sigma_{F_{i},\rm theory}=\frac{{\rm acc}}{2}F_{\rm theory}(k_{i})\,,$ (18)

where acc is our desired accuracy threshold, e.g. $1\%$ for the power
spectrum. This theoretical error term, while only contributing acc/2
correction to cosmic variance, dominates over the shot-noise error near
$k_{\rm NL}$. The contribution from the first five modes is suppressed by the
usual cosmic variance and stays less than 1% of the summation, which is
dominated by contributions near $k_{\rm NL}$. In effect, what this computes is
the average deviation from our accuracy threshold in units of the error, i.e.
in units of sigma, and finds where that crosses 1.5 sigma.

A final note about our choice of 1.5 sigma is in order. Here, the 1.5 sigma
constraint is somewhat arbitrary and was determined by empirical inspection,
to match the $k_{\rm NL}$ we determine in Eq. (17) and the values from the
comparison plots such as Fig. 5 and Fig. 7: using 1 sigma is too strict and 2
sigma is too lenient. The choice of threshold sigma, however, does not alter
the results in any significant way and merely shifts the specific values of
$k_{NL}$ up or down while maintaining the same general relationships.
Typically the difference from 1 to 1.5 and 1.5 to 2 is about a $10\%$
difference at each step. For values of $k_{\rm NL}$ that do not fall exactly
on a simulation bin center we use cubic spline interpolation to increase our
precision.

###  4.1 Power Spectrum

Figure 4: (Left): The power spectrum for all three of our models at all three
plotting redshifts (0, 1, 3). We plot all three redshifts on the same plot
since the scaling makes them mostly distinct. The largest valued set of power
spectra is the $z=0$ set, while the smallest is the $z=3$ set. It can clearly
be seen by eye how much better the high redshift curves match at smaller
scales compared to the $z=0$ set.
(Right): The various power spectrum models divided by the no-wiggle power
spectrum from Eisenstein & Hu (1998). The redshifts follow the same pattern as
on the left. Note that for both plots the errors are actually five-sigma
errors to enhance visibility.

To visualize the differences between the two (LO and LO+NLO) models and the
N-body result we present a plot of all three curves at all three plotting
redshifts ($z=0$, $1$, $3$) in Fig. 4. To facilitate visualization of the
nonlinearities in BAO, we also divide the curves by the no-wiggle power
spectrum from Eisenstein & Hu (1998) on the right panel of Fig. 4.

Here, we see excellent agreement between everything on large scales, with the
small-scale agreement being highly dependent on redshift, as expected (Jeong &
Komatsu, 2006a). The non-wiggle plots emphasize the baryon acoustic
oscillation (BAO) feature and show that LO+NLO perturbation theory prediction
can accurately capture the nonlinearities in BAO at $z=3$, while failing at
lower redshifts ($z=0$, $1$).

We also show the error plots for our three plotting redshifts ($z=0$, $1$,
$3$) in the left panel of Fig. 5 along with overplotted lines to highlight 1%
and 2% residuals. We see the expected redshift evolution that both LO and NLO
perturbation theory predictions become more accurate at higher redshifts.
Corresponding $k_{\rm NL}$ values as we compute following Eq. (17) can be seen
for all redshifts and models in the right panel of Fig. 5. Note that there is
a sharp discontinuity in the “Linear 1%” model. This sudden increase of
$k_{\rm NL}$ is due to the first BAO wiggle being inconsistent to 1% at low
redshift but consistent at high redshift. Besides that, the scaling shows the
expected behavior with slow increases on either side of the linear
discontinuity and steep increases at each redshift step for the LO+NLO model.

Figure 5: (Left): The relative error for each model compared to the simulation
at three different redshifts, from top to bottom $z=0$, $z=1$, $z=3$. Left:
The error for the three linear models. Here we see significant redshift
evolution, with the first point inconsistent with 1% going from $\sim$ 0.05
$h$/Mpc at $z=0$ to $\sim$ 0.15 $h$/Mpc at $z=3$. Right: The error for the
three NLO models sees a similar redshift evolution, going from diverging at
$\sim$ 0.09 $h$/Mpc at $z=0$ to $\sim$ 0.3 $h$/Mpc at $z=3$. It also remains
consistently better than just the linear model. Overall we see significant
accuracy improvement from both increasing the linear model to NLO and from
increasing redshift.
(Right): The value of $k_{\rm NL}$ for each redshift using both linear and
L+NLO models. These values were computed using Eq. (17). We see a clear
increase with redshift for both linear and NLO models, expected for
perturbative methods. We fit a simple model to our results, given by Eq. (19)
that we plot alongside the model given by Nishimichi et al. (2009). Note that
although the form of the model is cosmology independent it was calibrated
using a significantly different cosmology so differences of these levels are
not unexpected. The same caution should be taken when applying our fit model.

We fit a simple model to the results of our calculation for the LO+NLO results
and get an expression

$\frac{k_{\rm NL}(z)}{h/{\rm
Mpc}}=0.0502+0.0251\left(\frac{D(z)}{D(0)}\right)^{-1.939}$ (19)

where $D(z)$ is the linear growth factor. We do not fit a similar curve for
the linear model due to the discontinuity caused by the first BAO wiggle.
Along with our basic fit model we also look at the model described in
Nishimichi et al. (2009). This takes the form

$\frac{k_{\rm NL}(z)^{2}}{6\pi^{2}}\int_{0}^{k_{NL}(z)}P_{L}(q,z)\,{\rm
d}q<0.18$ (20)

Where the constant on the right hand side varies depending on the model, here
we use their $C_{1\%}^{\rm SPT}$. We find these $k_{NL}$ values for our L+NLO
model and plot it alongside our data and model in the right panel of Fig. 5.
Our fitted model is quite consistent with our results, particularly at higher
redshifts. The Nishimichi et al. (2009) model has similar behaviour to our
model, albeit with different scaling at high redshift, not unexpected for a
model that was calibrated with a significantly different cosmology, WMAP3
(Spergel et al., 2007) ($\Omega_{m}=0.234$ $\sigma_{8}=0.76$) vs Planck18
(Planck Collaboration et al., 2020a) ($\Omega_{m}=0.315$ $\sigma_{8}=0.811$).
The difference suggests that care should be taken before applying Eq. (19)
that, even though it is written in a form to be generalizable to other
cosmologies, is likely to have errors when applied to moderately different
cosmologies.

###  4.2 Bispectrum

Similarly to the power-spectrum case, here we present spherical bispectrum
plots for all redshifts ($z=0$, $1$, $3$) and models (LO and LO+NLO) in Fig.
6. We see very similar behavior to the power spectrum, one of our motivations
for introducing the spherical-bispectrum visualization technique, with the
LO+NLO model having better overall agreement and the agreement becomes better
with redshift. Since the full bispectrum does not have a single scale
associated with each element, we define the range of validity $k_{\rm sph,NL}$
based off the spherical bispectrum. The use of $k_{\rm sph,NL}$ allows us to
clearly analyze a single effective scale for deviation of the two models from
N-body data.

We essentially use Eq. (17) but with $B(k_{\rm sph})$ as $F$, and with a few
different cutoff errors, 2-5%, instead of 1% as for the power spectrum since
both observational and statistical errors are larger for the bispectrum. This
is because the bispectrum is a higher-order statistic, involving a larger
number of individual measurements. We present relative error plots of both LO
and LO+NLO perturbative bispectrum models at our three plotting redshifts
($z=0$, $1$, $3$) in the top panels of Fig. 7 with overplotted lines to
highlight 2% and 5% errors, the bounds of our analysis. We again see a similar
pattern with LO+NLO having better accuracy than LO, and the accuracy also
increasing with redshift for both models. Our final results for the total
bispectrum $k_{\rm sph,nl}$ are presented in the Fig. 7. Due to the difference
between $k_{\rm sph}$ and $k$ and the nature of $B_{\rm sph}$ itself, it is
hard to make a direct comparison between the bottom panel of Fig. 7 and the
right panel of Fig. 5. Nevertheless, we can read off the trend that $k_{\rm
NL}$ with 1% accuracy in power spectrum corresponds to $k_{\rm sph,NL}$ with
3% accuracy in bispectrum.

We can also perform this comparison between the theory calculation and N-body
results for specific configurations to disentangle the averaging effect in
$B_{\rm sph}$. This study can also clarify that the accuracy of PT’s
bispectrum modeling depends on the configuration. To compute $k_{\rm
sph,NL}^{\rm(config)}$ at each cutoff percentage, we again use Eq. (17) but
with the bispectrum filtered to a specific configuration, instead of the
bispectrum of all modes. We repeat the calculation for each model and with all
six configurations.

Figure 6: Left: The bispectrum for all three of our models at all three
plotting redshifts (0, 1, 3). We plot all three redshifts on the same plot
since the scaling makes them mostly distinct. The largest valued set of
bispectra is the z=0 set, while the smallest is the z=3 set. It can clearly be
seen by eye how much better the high redshift curves match at smaller scales
compared to the z=0 set. Right: The various bispectrum models divided by the
no-wiggle bispectrum, the LO bispectrum generated using the no-wiggle power
spectrum from Eisenstein & Hu (1998). The redshifts follow the same pattern as
on the left.

Figure 7: (Top): The relative error for each model compared to the simulation
at three different redshifts, from top to bottom $z=0$, $z=1$, $z=3$. Left:
The error for the three LO models. Here we see significant redshift evolution,
with the first point inconsistent with 5% going from $\sim$ 0.2 $h$/Mpc at
$z=0$ to $\sim$ 0.375 $h$/Mpc at $z=3$. Right: The error for the three LO+NLO
models sees a similar redshift evolution, going from diverging at $\sim$ 0.2
$h$/Mpc at $z=0$ to $\sim$ 0.5 $h$/Mpc at $z=3$. It also remains consistently
better than just the linear model. Overall we see significant accuracy
improvement from both increasing the LO model to LO+NLO and from increasing
redshift.
(Bottom): The value of $k_{\rm sph,nl}$ found using Eq. (17) with the
spherical bispectrum at four different threshold values. Left: The LO
bispectrum results at four different accuracy thresholds. In general we see a
smooth increase in $k_{\rm sph,nl}$ with both redshift and accuracy threshold.
Right: The results for the LO+NLO model at four accuracy thresholds. Generally
the results are similar to the LO results, except for a slight
flattening/decrease between redshifts 2 and 3. The 5% results also hit the
maximum value for which $k_{\rm sph}$ is complete, implying it is a lower
bound. Figure 8: The value of $k_{\rm sph,nl}$ found using Eq. (17) with the
spherical bispectrum, filtered to a specific configuration, at four different
threshold values. Left: The results using the LO bispectrum model. Inset in
each plot is a representative triangle showing the definition of the
configuration. The equilateral results show that the threshold accuracy,
within our bounds, does not change $k_{\rm sph,nl}$, signifying that the loss
of accuracy is abrupt, but has a reasonable redshift scaling. Squeezed
configurations seem to behave as expected, although the dependence on accuracy
threshold is quite weak. Folded configurations seems to have little accuracy
dependence, and experiences a significant jump from redshift 0 to 0.5. Right:
The results for the LO+NLO model. At NLO Equilateral configurations retain
their lack of dependence on accuracy threshold. It also hits the maximum of
the data at both redshift 2 and 3, representing lower bounds. NLO Squeezed
configurations have a much stronger accuracy dependence, with an added
flattening between redshift 2 and 3. At NLO Folded configurations have more of
a dependence on accuracy, hitting the maximum value at redshift 2, for some
thresholds, and redshift 3, for all thresholds. Figure 9: The value of $k_{\rm
sph,nl}$ found using Eq. (17) with the spherical bispectrum, filtered to a
specific configuration, at four different threshold values. Left: The results
using the LO bispectrum model. Inset in each plot is a representative triangle
showing the definition of the configuration. The elongated results show
expected redshift and accuracy scaling, except between redshift 0 and 1 where
the higher thresholds scale rapidly while the lower thresholds more slowly.
The Isosceles configuration results show reasonable scaling with redshift,
with a weak dependence on accuracy threshold. Similarly the Obtuse-Isosceles
results are very similar to the Isosceles results. Right: The results for the
LO+NLO model. At NLO Elongated configurations display the expected scaling,
with some of the redshift 2 and 3 results hitting the maximum value, implying
a lower bound. The NLO Isosceles configurations have a strange feature where
there is a significant drop from redshift 2 to 3, at least for 5% accuracy.
There are many potential causes for this, an unlikely statistical fluctuation,
an unexpected binning error, or many other possibilities. At NLO Obtuse-
Isosceles behaves similarly to Isosceles, with a drop from redshift 2 to 3,
although it does hit the maximum value at redshift 2.

We show the results for the each configuration in Fig. 8 and Fig. 9. Here we
see that most configurations behave in a similar way, LO+NLO is more accurate
than just LO, often by a significant amount, and accuracy improves
dramatically at higher redshifts with our accuracy threshold having the
expected minor scaling. There are a number of mildly anomalous results in that
there are a number of points where LO is equally as accurate as LO+NLO, in the
All, Isosceles, and Obtuse-Isosceles plots we see a slight decrease in
accuracy going from $z=2$ to $z=3$. We attribute these features to either be
statistical fluctuations or a consequence of imperfect binning corrections,
some configurations are significantly more dependent on the binning effect
than others (Bernardeau et al., 2012; Eggemeier et al., 2021). Due to the
inaccuracy of our large-scale bispectrum measurements we leave the
investigation of this problem for future work. A few other features of
interest in the plot are the Equilateral curves having little variation with
threshold. This represents the Equilateral models being very accurate with a
sudden divergence from the model causing all the thresholds to be passed at
the same point. There are also a few discontinuous points with the same cause
as for the LO power spectrum (the blue line in the right panel of Fig. 5), a
feature that isn’t captured at lower redshifts accurately that falls under the
threshold as the redshift increases.

##  5 Conclusion

The galaxy power spectrum and bispectrum are some of the main summary
statistics for studying large-scale structure (LSS) cosmology. Unlike the CMB
anisotropies, where $\Delta T/T\simeq 10^{-5}$ is small, linear perturbation
theory fails to describe late-time cosmic density field whose r.m.s.
fluctuation is of order unity on relevant scales. For LSS cosmology,
therefore, it is essential to model the nonlinear evolution of the power
spectrum and bispectrum, and standard perturbation theory (SPT) provides an
analytical expression for the nonlinear contributions, we call NLO (next-to-
leading order) in this paper. To properly exploit the power of an analytic
theory a firm grasp on the systematics and region of validity are required. In
this paper, we analyze the range of validity of SPT by comparing the
analytical calculation to a suite of N-body simulations.

We first compare the SPT power spectrum at both leading order (LO) and NLO to
the results of one hundred high-resolution run of Quijote simulation and
derive a range of validity, or $k_{\rm NL}$ the maximum wavenumber for which
the SPT theory is accurate to 1% in comparison to the N-body result. For the
LO+NLO SPT matter power spectrum the we find a fitting formula for $k_{\rm
NL}[D(z)]$ in Eq. (19).

To carry out a similar analysis for the matter bispectrum, we have introduced
a new visualization technique, the spherical bispectrum. The spherical
bispectrum allows for plotting and analysis of the bispectra as a one-
dimensional function of of $k_{\rm sph}$ that show the length-scale
dependence. At the same time, we can show the triangular-configuration
dependence by fixing the range of $\theta_{\rm sph}$ and $\phi_{\rm sph}$ to
define the configuration. This scheme, therefore, retains all information in
the bispectrum. By using this technique we find a good proxy for the maximum
wavenumber for the bispectrum in Fig. 7.

These results allow for more accurate assessment for the analytic theory’s
constraining power by properly limiting the smallest scale above which each
theory can accurately model the nonlinearities. That is, combined with an
accurate covariance matrix (Villaescusa-Navarro et al., 2020), since on small
scales there is often significant mixing between different wave modes which
reduces the information gain by including more modes, one can achieve very
accurate predicted constraints. We leave that as a future work.

One issue that the bispectrum analysis must resolve, but we have only
discussed briefly in this paper, is a more detailed analysis for proper
binning of the continuous theory models, to confront the theory calculation
with the data. There are many recent techniques attempting to alleviate the
computational cost while maintaining good accuracy (Bernardeau et al., 2012;
Sefusatti et al., 2010; Oddo et al., 2020; Eggemeier et al., 2021), but to our
knowledge no comprehensive study has been done to compare the various methods
while dealing with other potential issues such as interpolation scheme and bin
size. Unlike for the power-spectrum analysis where binning the theory
prediction causes only a minor problem for the first few modes, this is a key
issue for the bispectrum that still requires more systematic study for
accurate data analysis.

Of course, the range of validity analysis must also be done for the next-order
clustering statistic, the trispectrum or four point function. The challenge
with this type of extension is the significant increase in estimation cost,
estimating a trispectrum with similar bin size and max wavenumber as a
bispectrum takes $\sim$10 times as long (Tomlinson et al., 2019). Similarly,
because the cost of binning a continuous theory is proportional to the cost of
estimation the binning problem is even larger for analytic trispectrum
theories. All of these issues compound with the high computational cost of
evaluating the NLO trispectrum, likely limiting study to only the cheap to
compute LO part of the trispectrum. On the other hand, the trispectrum is
expected to be an important probe of primordial non-Gaussanities (Gualdi et
al., 2021a; Gualdi et al., 2021b), as well as parity-violating new physics
(Cahn et al., 2021). The LSS analysis with trispectrum has just begun recently
(Philcox et al., 2021; Gualdi & Verde, 2022), and the high cost of including
the trispectrum in future analysis demands a detailed study of the information
contained in the trispectrum, similar to what Hahn et al. (2020) did for the
bispectrum, is needed to assess how practical the trispectrum is to use.

In this work, we have focused solely on the analytical SPT models that do not
include any free parameters. Most extended models are built off a foundation
of SPT with a couple of free parameters. Therefore, more advanced models such
as Effective Field Theory (EFT) (Baumann et al., 2012; Carrasco et al., 2012;
Carrasco et al., 2014a, b; Hertzberg, 2014; Senatore & Zaldarriaga, 2015),
Regularized Perturbation Theory (RegPT) (Taruya et al., 2012), Renormalized
Lagrangian Perturbation Theory (RLPT) (Matsubara, 2008; Rampf & Wong, 2012),
Galilean-invariant Renormalized Perturbation Theory (gRPT) (Eggemeier et al.,
2020) or the general bias expansion (Desjacques et al., 2018a, b; Tomlinson et
al., 2020) can use this type of study as the basis for either their own model
limitations or to seed a more detailed study of each model. Generically, the
bounds $k_{\rm NL}$ in this work should serve as lower bounds for these more
advanced methods which often have different corrections to improve accuracy.
For example, EFT often uses a counterterm like $\Delta P\propto k^{2}P_{L}$
which based off the results in this work, Fig. 5, should be a good choice to
improve accuracy. Note that the same argument applies to the stochastic bias
parameter in the general bias method. For a comparison of the range of
validity of a few different methods at $z=1$, see Alkhanishvili et al. (2022).

## Acknowledgements

The authors would like to thank Francisco Villaescusa-Navarro for assistance
with data and computing access. We would also like to thank the whole Quijote
team for making the simulation suite available and well documented. The author
acknowledge that the work reported on in this paper was substantially
performed using the Princeton Research Computing resources at Princeton
University which is consortium of groups led by the Princeton Institute for
Computational Science and Engineering (PICSciE) and Office of Information
Technology’s Research Computing. This work was supported at Pennsylvania State
University by NASA ATP program (80NSSC18K1103).

## Data Availability

The underlying data for this work is available at https://quijote-
simulations.readthedocs.io/en/latest/data.html. The 100 high resolution
bispectra we computed for this work are available either by request from the
authors or on the Princeton Tiger cluster.

## References

  * Agarwal et al. (2014) Agarwal N., Ho S., Shandera S., 2014, Journal of Cosmology and Astroparticle Physics, 2014, 038
  * Alkhanishvili et al. (2022) Alkhanishvili D., Porciani C., Sefusatti E., Biagetti M., Lazanu A., Oddo A., Yankelevich V., 2022, Monthly Notices of the Royal Astronomical Society
  * Angulo et al. (2015) Angulo R. E., Foreman S., Schmittfull M., Senatore L., 2015, Journal of Cosmology and Astroparticle Physics, 2015, 039–039
  * Assassi et al. (2014) Assassi V., Baumann D., Green D., Zaldarriaga M., 2014, Journal of Cosmology and Astroparticle Physics, 2014, 056–056
  * Barreira et al. (2020) Barreira A., Cabass G., Schmidt F., Pillepich A., Nelson D., 2020, J. Cosmology Astropart. Phys., 2020, 013
  * Baumann et al. (2012) Baumann D., Nicolis A., Senatore L., Zaldarriaga M., 2012, J. Cosmology Astropart. Phys., 2012, 051
  * Bennett et al. (2013) Bennett C. L., et al., 2013, ApJS, 208, 20
  * Bernardeau et al. (2002) Bernardeau F., Colombi S., Gaztañaga E., Scoccimarro R., 2002, Physics Reports, 367, 1–248
  * Bernardeau et al. (2012) Bernardeau F., Crocce M., Scoccimarro R., 2012, Physical Review D, 85
  * Blazek et al. (2015) Blazek J., Vlah Z., Seljak U., 2015, Journal of Cosmology and Astroparticle Physics, 2015, 015–015
  * Byun et al. (2017) Byun J., Eggemeier A., Regan D., Seery D., Smith R. E., 2017, Monthly Notices of the Royal Astronomical Society, 471, 1581–1618
  * Cabass et al. (2022) Cabass G., Ivanov M. M., Philcox O. H. E., Simonović M., Zaldarriaga M., 2022, arXiv e-prints, p. arXiv:2201.07238
  * Cahn et al. (2021) Cahn R. N., Slepian Z., Hou J., 2021, arXiv e-prints, p. arXiv:2110.12004
  * Carrasco et al. (2012) Carrasco J. J. M., Hertzberg M. P., Senatore L., 2012, Journal of High Energy Physics, 2012
  * Carrasco et al. (2014a) Carrasco J. J. M., Foreman S., Green D., Senatore L., 2014a, Journal of Cosmology and Astroparticle Physics, 2014, 056
  * Carrasco et al. (2014b) Carrasco J. J. M., Foreman S., Green D., Senatore L., 2014b, Journal of Cosmology and Astroparticle Physics, 2014, 057
  * Chudaykin & Ivanov (2019) Chudaykin A., Ivanov M. M., 2019, Journal of Cosmology and Astroparticle Physics, 2019, 034–034
  * Collaboration (2012) Collaboration L. D. E. S., 2012, Large Synoptic Survey Telescope: Dark Energy Science Collaboration (arXiv:1211.0310)
  * Desjacques et al. (2018a) Desjacques V., Jeong D., Schmidt F., 2018a, Physics Reports, 733, 1–193
  * Desjacques et al. (2018b) Desjacques V., Jeong D., Schmidt F., 2018b, Journal of Cosmology and Astroparticle Physics, 2018, 035–035
  * Dizgah et al. (2020) Dizgah A. M., Biagetti M., Sefusatti E., Desjacques V., Noreña J., 2020, Primordial Non-Gaussianity from Biased Tracers: Likelihood Analysis of Real-Space Power Spectrum and Bispectrum (arXiv:2010.14523)
  * Doré et al. (2015) Doré O., et al., 2015, Cosmology with the SPHEREX All-Sky Spectral Survey (arXiv:1412.4872)
  * Eggemeier et al. (2019) Eggemeier A., Scoccimarro R., Smith R. E., 2019, Physical Review D, 99
  * Eggemeier et al. (2020) Eggemeier A., Scoccimarro R., Crocce M., Pezzotta A., Sánchez A. G., 2020, Phys. Rev. D, 102, 103530
  * Eggemeier et al. (2021) Eggemeier A., Scoccimarro R., Smith R. E., Crocce M., Pezzotta A., Sánchez A. G., 2021, Testing one-loop galaxy bias: joint analysis of power spectrum and bispectrum (arXiv:2102.06902)
  * Eisenstein & Hu (1998) Eisenstein D. J., Hu W., 1998, The Astrophysical Journal, 496, 605–614
  * Giannantonio et al. (2014) Giannantonio T., Ross A. J., Percival W. J., Crittenden R., Bacher D., Kilbinger M., Nichol R., Weller J., 2014, Phys. Rev. D, 89, 023511
  * Gil-Marín et al. (2016) Gil-Marín H., Percival W. J., Verde L., Brownstein J. R., Chuang C.-H., Kitaura F.-S., Rodríguez-Torres S. A., Olmstead M. D., 2016, Monthly Notices of the Royal Astronomical Society, 465, 1757–1788
  * Gualdi & Verde (2022) Gualdi D., Verde L., 2022, arXiv e-prints, p. arXiv:2201.06932
  * Gualdi et al. (2021a) Gualdi D., Novell S., Gil-Marín H., Verde L., 2021a, Journal of Cosmology and Astroparticle Physics, 2021, 015
  * Gualdi et al. (2021b) Gualdi D., Gil-Marín H., Verde L., 2021b, Journal of Cosmology and Astroparticle Physics, 2021, 008
  * Hahn & Villaescusa-Navarro (2020) Hahn C., Villaescusa-Navarro F., 2020, Constraining $M_{\nu}$ with the Bispectrum II: The Total Information Content of the Galaxy Bispectrum (arXiv:2012.02200)
  * Hahn et al. (2020) Hahn C., Villaescusa-Navarro F., Castorina E., Scoccimarro R., 2020, Journal of Cosmology and Astroparticle Physics, 2020, 040–040
  * Hamaus et al. (2011) Hamaus N., Seljak U., Desjacques V., 2011, Phys. Rev. D, 84, 083509
  * Hand et al. (2018) Hand N., Feng Y., Beutler F., Li Y., Modi C., Seljak U., Slepian Z., 2018, The Astronomical Journal, 156, 160
  * Hertzberg (2014) Hertzberg M. P., 2014, Phys. Rev. D, 89, 043521
  * Hill et al. (2008) Hill G. J., et al., 2008, The Hobby-Eberly Telescope Dark Energy Experiment (HETDEX): Description and Early Pilot Survey Results (arXiv:0806.0183)
  * Hirata (2009) Hirata C. M., 2009, Monthly Notices of the Royal Astronomical Society, 399, 1074–1087
  * Ivanov et al. (2020) Ivanov M. M., Simonović M., Zaldarriaga M., 2020, Journal of Cosmology and Astroparticle Physics, 2020, 042–042
  * Jeong & Komatsu (2006a) Jeong D., Komatsu E., 2006a, ApJ, 651, 619
  * Jeong & Komatsu (2006b) Jeong D., Komatsu E., 2006b, The Astrophysical Journal, 651, 619–626
  * Jeong & Komatsu (2009) Jeong D., Komatsu E., 2009, ApJ, 703, 1230
  * Kamalinejad & Slepian (2020) Kamalinejad F., Slepian Z., 2020, A Non-Degenerate Neutrino Mass Signature in the Galaxy Bispectrum (arXiv:2011.00899)
  * Karagiannis et al. (2020) Karagiannis D., Fonseca J., Maartens R., Camera S., 2020, Probing primordial non-Gaussianity with the bispectrum of future 21cm intensity maps (arXiv:2010.07034)
  * Laureijs et al. (2011) Laureijs R., et al., 2011, arXiv e-prints, p. arXiv:1110.3193
  * Lazanu et al. (2016) Lazanu A., Giannantonio T., Schmittfull M., Shellard E., 2016, Physical Review D, 93
  * Levi et al. (2013) Levi M., et al., 2013, The DESI Experiment, a whitepaper for Snowmass 2013 (arXiv:1308.0847)
  * Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, The Astrophysical Journal, 538, 473–476
  * Maartens et al. (2015) Maartens R., Abdalla F. B., Jarvis M., Santos M. G., 2015, Cosmology with the SKA – overview (arXiv:1501.04076)
  * Matsubara (2008) Matsubara T., 2008, Phys. Rev. D, 77, 063530
  * McCullagh et al. (2015) McCullagh N., Jeong D., Szalay A. S., 2015, Monthly Notices of the Royal Astronomical Society, 455, 2945–2958
  * McEwen et al. (2016) McEwen J. E., Fang X., Hirata C. M., Blazek J. A., 2016, Journal of Cosmology and Astroparticle Physics, 2016, 015–015
  * Moradinezhad Dizgah et al. (2021) Moradinezhad Dizgah A., Biagetti M., Sefusatti E., Desjacques V., Noreña J., 2021, J. Cosmology Astropart. Phys., 2021, 015
  * Nishimichi et al. (2009) Nishimichi T., et al., 2009, Publications of the Astronomical Society of Japan, 61, 321–332
  * Oddo et al. (2020) Oddo A., Sefusatti E., Porciani C., Monaco P., Sánchez A. G., 2020, Journal of Cosmology and Astroparticle Physics, 2020, 056–056
  * Osato et al. (2019) Osato K., Nishimichi T., Bernardeau F., Taruya A., 2019, Physical Review D, 99
  * Philcox et al. (2021) Philcox O. H. E., Hou J., Slepian Z., 2021, arXiv e-prints, p. arXiv:2108.01670
  * Planck Collaboration et al. (2020a) Planck Collaboration et al., 2020a, A&A, 641, A6
  * Planck Collaboration et al. (2020b) Planck Collaboration et al., 2020b, A&A, 641, A9
  * Rampf & Wong (2012) Rampf C., Wong Y. Y., 2012, Journal of Cosmology and Astroparticle Physics, 2012, 018–018
  * Rezaie et al. (2021) Rezaie M., et al., 2021, Monthly Notices of the Royal Astronomical Society, 506, 3439
  * Schmittfull & Vlah (2016) Schmittfull M., Vlah Z., 2016, Physical Review D, 94
  * Schmittfull et al. (2016) Schmittfull M., Vlah Z., McDonald P., 2016, Physical Review D, 93
  * Scoccimarro (2004) Scoccimarro R., 2004, Physical Review D, 70
  * Scoccimarro (2015) Scoccimarro R., 2015, Physical Review D, 92
  * Sefusatti et al. (2006) Sefusatti E., Crocce M., Pueblas S., Scoccimarro R., 2006, Physical Review D, 74
  * Sefusatti et al. (2010) Sefusatti E., Crocce M., Desjacques V., 2010, Monthly Notices of the Royal Astronomical Society, p. no–no
  * Sefusatti et al. (2012) Sefusatti E., Crocce M., Desjacques V., 2012, Monthly Notices of the Royal Astronomical Society, 425, 2903–2930
  * Sefusatti et al. (2016) Sefusatti E., Crocce M., Scoccimarro R., Couchman H. M. P., 2016, Monthly Notices of the Royal Astronomical Society, 460, 3624–3636
  * Senatore & Zaldarriaga (2015) Senatore L., Zaldarriaga M., 2015, Journal of Cosmology and Astroparticle Physics, 2015, 013
  * Simonović et al. (2018) Simonović M., Baldauf T., Zaldarriaga M., Carrasco J. J., Kollmeier J. A., 2018, J. Cosmology Astropart. Phys., 2018, 030
  * Simonović et al. (2018) Simonović M., Baldauf T., Zaldarriaga M., Carrasco J. J., Kollmeier J. A., 2018, Journal of Cosmology and Astroparticle Physics, 2018, 030–030
  * Slepian (2018) Slepian Z., 2018, On Decoupling the Integrals of Cosmological Perturbation Theory (arXiv:1812.02728)
  * Slepian et al. (2017a) Slepian Z., et al., 2017a, Mon. Not. Roy. Astron. Soc., 468, 1070
  * Slepian et al. (2017b) Slepian Z., et al., 2017b, Monthly Notices of the Royal Astronomical Society, 469, 1738–1751
  * Slosar et al. (2008) Slosar A., Hirata C., Seljak U., Ho S., Padmanabhan N., 2008, Journal of Cosmology and Astroparticle Physics, 2008, 031
  * Smith et al. (2008) Smith R. E., Sheth R. K., Scoccimarro R., 2008, Physical Review D, 78
  * Song et al. (2015) Song Y.-S., Taruya A., Oka A., 2015, Journal of Cosmology and Astroparticle Physics, 2015, 007–007
  * Spergel et al. (2007) Spergel D. N., et al., 2007, ApJS, 170, 377
  * Spergel et al. (2015) Spergel D., et al., 2015, Wide-Field InfrarRed Survey Telescope-Astrophysics Focused Telescope Assets WFIRST-AFTA 2015 Report (arXiv:1503.03757)
  * Takada et al. (2014) Takada M., et al., 2014, Publications of the Astronomical Society of Japan, 66, R1
  * Taruya et al. (2009) Taruya A., Nishimichi T., Saito S., Hiramatsu T., 2009, Physical Review D, 80
  * Taruya et al. (2010) Taruya A., Nishimichi T., Saito S., 2010, Physical Review D, 82
  * Taruya et al. (2012) Taruya A., Bernardeau F., Nishimichi T., Codis S., 2012, Phys. Rev. D, 86, 103528
  * Taruya et al. (2018) Taruya A., Nishimichi T., Jeong D., 2018, Physical Review D, 98
  * Tellarini et al. (2016) Tellarini M., Ross A. J., Tasinato G., Wands D., 2016, Journal of Cosmology and Astroparticle Physics, 2016, 014–014
  * Tomlinson et al. (2019) Tomlinson J., Jeong D., Kim J., 2019, The Astronomical Journal, 158, 116
  * Tomlinson et al. (2020) Tomlinson J., Gebhardt H. S. G., Jeong D., 2020, Physical Review D, 101
  * Villaescusa-Navarro et al. (2020) Villaescusa-Navarro F., et al., 2020, The Astrophysical Journal Supplement Series, 250, 2
  * Yankelevich & Porciani (2018) Yankelevich V., Porciani C., 2018, Monthly Notices of the Royal Astronomical Society, 483, 2078–2099
  * Zheng et al. (2010) Zheng Z., Cen R., Trac H., Miralda-Escudé J., 2010, The Astrophysical Journal, 726, 38