# The Bispectrum in Lagrangian Perturbation Theory

Shi-Fan Chen  Zvonimir Vlah  Martin White

###### Abstract

We study the bispectrum in Lagrangian perturbation theory. Extending past
results for the power spectrum, we describe a method to efficiently compute
the bispectrum in LPT, focusing on the Zeldovich approximation, in which
contributions due to linear displacements are captured to all orders in a
manifestly infrared (IR) safe way. We then isolate the effects of these linear
displacements on oscillatory components of the power spectrum like baryon
acoustic oscillations or inflationary primordial features and show that the
Eulerian perturbation theory (EPT) prescription wherein their effects are
resummed by a Gaussian damping of the oscillatory component in the linear
power spectrum arise as a saddle-point approximation of our calculation. These
two methods of IR resummation are in excellent agreement at 1-loop in the
bispectrum, though at tree level resummed EPT does less well to capture the
nonlinear damping of the oscillations, and the LPT calculation does not
require an artificial split of the power spectrum into smooth and oscillatory
components, making it particularly useful for modeling exotic features in
higher-order N-point functions. We finish by extending our analysis of IR
resummation in LPT to N-point functions of arbitrary order.

##  1 Introduction

The observed large-scale structure is thought to arise through the action of
gravitational instability in a cold dark matter dominated, expanding Universe
from (very nearly) Gaussian initial conditions [1]. The pattern of non-
Gaussianity thus imparted is of a very particular form, and the tower of
cumulants of the density and velocity fields on large scales can be
systematically computed within perturbation theory [2]. While for a Gaussian
field the power spectrum is a sufficient statistic, and thus contains all of
the primordial information, once the field becomes non-Gaussian there is
information to be gained about gravitational interactions, cosmological
parameters and the bias relation from measurements beyond the power spectrum.
For a given set of galaxy bias parameters cosmological perturbation theory
makes consistent predictions for all of these n-point functions, which can be
combined with the power spectrum in data analyses to either improve power-
spectrum-only constraints or access new physical signals not present in the
2-point function. Recent years have in particular seen rapid advances in the
use of the 3-point bispectrum to measure primordial non-Gaussianity [3, 4, 5]
and tighten constraints on cosmological parameters [6]. N-point functions are
particularly appealing because, on scales where the non-Gaussianity is weak
and well-controlled, the available cosmological signal is predominantly at low
N such that the tower of cumulants can be truncated at some finite order. In
fact, recent work suggests that measurement of the low-order cumulants
contains effectively all of the perturbative information that is present in
the galaxy density field if attention is restricted to the larger scales where
the modeling is most robust, unless the contributions of certain modes is
accidentally parametrically large compared to expectations at a given order in
perturbation theory [5].

The presence of a feature in the linear theory power spectrum due to baryon
acoustic oscillations [7, 8] is one such complication in the perturbative
calculation of N-point functions, since it leads to a “small” parameter
related to the large-scale displacements of galaxies that is numerically quite
large. In order to handle this it is necessary to sum a subset of the full
perturbative terms to high order, a procedure that goes by the name of “IR
resummation.” Since it fundamentally traces the advection of matter due to
long-wavelength perturbations, IR resummation is very naturally handled within
the Lagrangian formulation of perturbation theory (LPT; [9, 10, 11, 12, 13,
14, 15, 16]), and indeed most implementations of IR resummation use LPT at
some point in their formulation (see e.g. [17, 18] and the discussion in
[19]).

Beyond the direct Lagrangian resummation of long-wavelength displacements it
is also possible to specifically isolate the BAO feature using a “wiggle no-
wiggle split”, estimate the effect of long-displacements on it, and resum only
those particular effects. Since this procedure resums the linear displacements
only in the wiggle component we will call it RWiggle for short throughout this
paper. The result is effectively a Gaussian damping of BAO wiggles in the
linear power spectrum [20] and can be shown using a simple saddle-point
approximation of the full Lagrangian calculation [21] or by diagramatically
selecting enhanced diagrams in time-sliced perturbation theory (TSPT) [22].
The agreement between these two procedures has been very well studied for the
2-point function, i.e. the power spectrum, and the level of agreement of
different schemes is well understood [23]. A similar study has not yet been
performed for cumulants beyond second order, with essentially all calculations
to date using the Gaussian damping ansatz due in part to its simplicity
compared to the full Lagrangian calculation for statistics beyond the power
spectrum.

Beyond the theoretical advantage of being able to compare different IR
resummation schemes, the full Lagrangian calculation also has the advantage
that it does not require that oscillatory features in the initial conditions
be artificially separated. This separation is not inherent to the underlying
physical theory and, while well-studied in the case of the BAO, will need to
be done on a case-by-case basis if, for example, searching for oscillations in
the power spectrum due to exotic inflationary models. The IR resummation of
the wiggles can still be performed in this case [24], but depends on more
conditions such as distinguishing between smooth and oscillatory components
when there is no single BAO-like scale in the problem, making the out-of-the
box nature of the Lagrangian calculation an appealing alternative. This will
be particularly useful for upcoming spectroscopic surveys, which will enable
us to place extremely tight bounds on these inflationary features [25, 26, 27,
28, 29], and for extending the treatment of IR resummation to oscillations
with complicated shapes in the bispectrum to constrain e.g. interactions of
the inflaton with massive particles with spin (also known as the “cosmological
collider”) [30, 31, 5] .

The purpose of this paper is to investigate how the bispectrum and higher
order functions can be efficiently computed within LPT while fully resumming
the effect of long-wavelength displacements, and what we learn about IR
resummation from considering the LPT bispectrum. We are not the first to
compute n-point functions in LPT for this purpose: ref. [32] for example
considered the configuration-space 3-point function in the Zeldovich
approximation by brute force, multi-dimensional integration while ref. [33]
formulated the LPT bispectrum in the spirit of integrated perturbation theory
(iPT) [12]. However, we update upon earlier work and extend the calculation of
the 3-point function to Fourier space using fast FFTLog-based numerical
techniques and study the phenomenology of oscillations in the bispectrum in
detail, fully taking into account mode-coupling terms and elucidating the
connection between the full Lagrangian calculation and the RWiggle
prescription. The structure of the paper is as follows: After describing the
general structure and expression for the bispectrum in LPT in §2, we show how
it can be efficiently evaluated in an IR-safe way cancelling unwanted
divergences in §3. In §4 we describe the behavior of BAO wiggles in some
detail, comparing RWiggle with the direct Lagrangian calculation via numerical
techniques derived in the previous section, as well as making contact with
other IR resummation schemes in the literature. This discussion is extended to
n-point functions of arbitrary order in §5. We conclude in §6 and discuss the
straightforward extension to nonlinear bias and redshift-space distortions as
well as an example calculation of the damping of inflationary features with a
logarithmic dependence on wavenumber. Various technical details are derived in
the Appendices. Numerical calculations throughout this paper will be done
assuming the Planck 2018 best-fit flat $\Lambda$CDM cosmology [27]
($\Omega_{m}=0.3111$, $h=0.6766$, $\omega_{b}=0.02242$, $n_{s}=0.9665$,
$\sigma_{8}=0.8102$) at redshift $z=0$ unless otherwise indicated.

##  2 The Bispectrum in Lagrangian Perturbation Theory

Lagrangian perturbation theory (LPT; [9, 10, 11, 12, 13, 14, 15]) is by now
well established and efficient codes for numerical computation of power
spectra exist [23]. We will not repeat all of the details here, instead
referring the reader to the above-cited articles; briefly LPT models the
growth of structure by following the displacements of fluid elements from
their initial, Lagrangian positions ($\mathbf{q}$) to their final, Eulerian
positions ($\mathbf{x}$) via a displacement:
$\mathbf{x}=\mathbf{q}+\mathbf{\Psi}(\mathbf{q})$. The displacement,
$\mathbf{\Psi}$, is then systematically expanded order by order in the initial
density perturbation. Biased tracers are modeled by introducing a functional,
$F$, of the initial density, velocity and tidal shear and their low-order
derivatives. Within this formalism the galaxy density $\delta_{g}(\textbf{x})$
has a Fourier transform given by [12, 13]

$\delta_{g}(\textbf{k})=\int d^{3}\textbf{q}\
e^{-i\textbf{k}\cdot\textbf{q}}\Big{(}F(\textbf{q})e^{-i\textbf{k}\cdot\Psi(\textbf{q})}-1\Big{)}\,.$
(2.1)

where $F(\mathbf{q})$ indicates the bias functional evaluated at the
Lagrangian coordinate $\mathbf{q}$. For convenience we will ignore the $-1$
piece below, since it is nonzero only for vanishing wavenumbers and is anyway
trivial to evaluate. As cosmological perturbation theory is an effective
theory modeling the large-scale fluid limit of structure formation, both the
bias functional $F(\textbf{q})$ and the perturbative expansion of the
displacement $\Psi(\textbf{q})$ have contributions whose sizes are free
parameters, allowing us to marginalize over unknown small scale effects like
galaxy formation (see e.g. [34, 19, 35]); since we focus on the role of long-
wavelength linear displacements in this work we will not further discuss this
aspect of LPT, and PT in genereal, except where relevant, leaving interested
readers to consult the literature on this well-studied topic (see e.g. refs.
[36, 37] for thorough treatments of the matter bispectrum at 1- and 2-loops).

From the above the bispectrum can then be written as

$\left\langle\delta(\textbf{k}_{1})\delta(\textbf{k}_{2})\delta(\textbf{k}_{3})\right\rangle=\int_{\textbf{q}_{1},\textbf{q}_{2},\textbf{q}_{3}}e^{-i\textbf{k}_{1}\cdot\textbf{q}_{1}-i\textbf{k}_{2}\cdot\textbf{q}_{2}-i\textbf{k}_{3}\cdot\textbf{q}_{3}}\left\langle
F_{1}F_{2}F_{3}\
e^{-i\textbf{k}_{1}\cdot\Psi_{1}-i\textbf{k}_{2}\cdot\Psi_{2}-i\textbf{k}_{3}\cdot\Psi_{3}}\right\rangle\,,$

where we have used the shorthand $f_{n}=f(\textbf{q}_{n}).$111For clarity of
presentation, we will adopt the shorthands  $\int_{\textbf{q}_{n}}\equiv\int
d^{3}\textbf{q}_{n}\,,\quad\int_{\textbf{p}_{n}}\equiv\int\frac{d^{3}\textbf{p}_{n}}{(2\pi)^{3}}\,,$
(2.2) for configuration and Fourier space integrals, respectively. We also
also label Lagrangian quantities by their coordinates, e.g.
$F_{n}=F(\textbf{q}_{n})$.  This expression can be simplified by noting that
the expectation value is translation invariant, i.e. does not depend on all
three Lagrangian coordinates $\textbf{q}_{1,2,3}$ independently. Adopting the
coordinate system
$(\textbf{q},\textbf{r},\textbf{Q})=(\textbf{q}_{1}-\textbf{q}_{3},\textbf{q}_{2}-\textbf{q}_{3},\textbf{q}_{3})$,
we can therefore eliminate the Q dependence in the bracketed mean and write

$\displaystyle\left\langle\delta(\textbf{k}_{1})\delta(\textbf{k}_{2})\delta(\textbf{k}_{3})\right\rangle$
$\displaystyle=\int_{\textbf{q},\textbf{r},\textbf{Q}}e^{-i\textbf{k}_{1}\cdot\textbf{q}-i\textbf{k}_{2}\cdot\textbf{r}-i\textbf{k}_{123}\cdot\textbf{Q}}\left\langle
F(\textbf{q})F(\textbf{r})F(\textbf{0})\
e^{-i\textbf{k}_{1}\cdot\Psi(\textbf{q})-i\textbf{k}_{2}\cdot\Psi(\textbf{r})-i\textbf{k}_{3}\cdot\Psi(\textbf{0})}\right\rangle$
$\displaystyle=(2\pi)^{3}\delta_{D}(\textbf{k}_{123})\
\int_{\textbf{q},\textbf{r}}e^{-i\textbf{k}_{1}\cdot\textbf{q}-i\textbf{k}_{2}\cdot\textbf{r}}\left\langle
F_{1}F_{2}F_{3}\
e^{-i\textbf{k}_{1}\cdot\Psi_{1}-i\textbf{k}_{2}\cdot\Psi_{2}+i(\textbf{k}_{1}+\textbf{k}_{2})\cdot\Psi_{3}}\right\rangle\,,$

where $\textbf{k}_{123}=\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3}$, and in
the last line, we have used that this correlator is zero unless
$\textbf{k}_{123}=0$. The bispectrum is defined to be this correlator without
the $\delta$-function, i.e. [33]

$B(\textbf{k}_{1},\textbf{k}_{2})=\int_{\textbf{q},\textbf{r}}e^{-i\textbf{k}_{1}\cdot\textbf{q}-i\textbf{k}_{2}\cdot\textbf{r}}\left\langle
F_{1}F_{2}F_{3}\
e^{-i\textbf{k}_{1}\cdot\Delta_{13}-i\textbf{k}_{2}\cdot\Delta_{23}}\right\rangle_{(\textbf{q},\textbf{r})=(\textbf{q}_{1}-\textbf{q}_{3}\textbf{q}_{2}-\textbf{q}_{3})}\,,$
(2.3)

where the pairwise displacement is $\Delta_{nm}=\Psi_{n}-\Psi_{m}$.
Importantly, the right-hand side of Equation 2.3 is invariant under uniform
shifts of the displacement, i.e. satisfies Galilean symmetry.

###  2.1 Matter Bispectrum in the Zeldovich Approximation

Let us start in the case of matter in the Zeldovich approximation [9], i.e.
lowest order LPT. Then $F_{n}=1$ and $\Psi=\Psi^{(1)}$ with
$\Psi^{(1)}(\textbf{k})=i(\textbf{k}/k^{2})\delta_{\rm lin}$. In the Zeldovich
approximation displacements are Gaussian, since they are first order in
$\delta_{\rm lin}$, so the cumulant theorem allows us to write

$\left\langle
e^{-i\textbf{k}_{1}\cdot\Delta_{13}-i\textbf{k}_{2}\cdot\Delta_{23}}\right\rangle=\exp\\{-\frac{1}{2}\Big{(}k_{1,i}k_{1,j}\left\langle\Delta_{13}^{i}\Delta_{13}^{j}\right\rangle+k_{2,i}k_{2,j}\left\langle\Delta_{23}^{i}\Delta_{23}^{j}\right\rangle+2k_{1,i}k_{2,j}\left\langle\Delta_{13}^{i}\Delta_{23}^{j}\right\rangle\Big{)}\Big{\\}}\,.$

The moments of the pairwise displacements can be written as functions of the
pairwise separations q, r by translation invariance [14]

$\left\langle\Delta_{13}^{i}\Delta_{13}^{j}\right\rangle=A_{ij}(\textbf{q})\,,\quad\left\langle\Delta_{23}^{i}\Delta_{23}^{j}\right\rangle=A_{ij}(\textbf{r})\,,$

and can be decomposed into tensor comonents as
$A_{ij}(\textbf{q})=X(q)\delta_{ij}+Y(q)\hat{q}_{i}\hat{q}_{j}$ [12]. The
“mixed”, third term is slightly more complicated; it can be written as

$\displaystyle\left\langle\Delta_{13}^{i}\Delta_{23}^{j}\right\rangle$
$\displaystyle=\left\langle(\Psi_{1}-\Psi_{3})^{i}(\Psi_{2}-\Psi_{3})^{j}\right\rangle$
$\displaystyle=\left\langle\Psi_{1}^{i}\Psi_{2}^{j}\right\rangle-\left\langle\Psi_{1}^{i}\Psi_{3}^{j}\right\rangle-\left\langle\Psi_{3}^{i}\Psi_{2}^{j}\right\rangle+\left\langle\Psi_{3}^{i}\Psi_{3}^{j}\right\rangle$
$\displaystyle=\left(\left\langle\Psi_{1}^{i}\Psi_{2}^{j}\right\rangle-\left\langle\Psi^{i}\Psi^{j}\right\rangle\right)+\left(\left\langle\Psi^{i}\Psi^{j}\right\rangle-\left\langle\Psi_{1}^{i}\Psi_{3}^{j}\right\rangle\right)+\left(\left\langle\Psi^{i}\Psi^{j}\right\rangle-\left\langle\Psi_{3}^{i}\Psi_{2}^{j}\right\rangle\right)$
$\displaystyle=\frac{1}{2}\left[A_{ij}(\textbf{r})+A_{ij}(\textbf{q})-A_{ij}(\textbf{q}_{12})\right];\quad\textbf{q}_{12}=\textbf{q}_{1}-\textbf{q}_{2}=\textbf{q}-\textbf{r}\,.$
(2.4)

The full expectation value can thus be written as

$\displaystyle\left\langle
e^{-i\textbf{k}_{1}\cdot\Delta_{13}-i\textbf{k}_{2}\cdot\Delta_{23}}\right\rangle$
$\displaystyle=\exp\\{-\frac{1}{2}\Big{(}k_{1,i}k_{1,j}A_{ij}(\textbf{q})+k_{2,i}k_{2,j}A_{ij}(\textbf{r})+k_{1,i}k_{2,j}\big{(}A_{ij}(\textbf{q})+A_{ij}(\textbf{r})-A_{ij}(\textbf{q}-\textbf{r})\big{)}\Big{)}\Big{\\}}$
$\displaystyle=\exp\\{\frac{1}{2}\Big{(}k_{1,i}k_{3,j}A_{ij}(\textbf{q})+k_{2,i}k_{3,j}A_{ij}(\textbf{r})+k_{1,i}k_{2,j}A_{ij}(\textbf{q}-\textbf{r})\Big{)}\Big{\\}}\,,$
(2.5)

where in the last line, we have used momentum conservation
$\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3}=0$ to restore the symmetry
between the sides of the triangle. Note that this form is symmetric in both
the wavenumbers $\textbf{k}_{1,2,3}$ and the three vertices of the triangle
$\textbf{q}_{1,2,3}$.

###  2.2 Tree Level

It will be quite instructive to check equation 2.5 at the lowest order, i.e.
tree level, in such a way as to make contact with the more familiar Eulerian
perturbation theory (EPT) calculation. In the Zeldovich approximation,
expanding order-by-order in the initial conditions, the density can be written
as

$\delta(\textbf{k})=\sum_{n=0}^{\infty}\int_{\textbf{p}_{1},...,\textbf{p}_{n}}\
Z_{n}(\textbf{p}_{1},...,\textbf{p}_{n})\
\delta_{0}(\textbf{p}_{1})...\delta_{0}(\textbf{p}_{n})\
(2\pi)^{3}\delta_{D}(\textbf{k}-\textbf{p}_{\rm tot})\,,$

where [38]

$Z_{n}(\textbf{p}_{1},...,\textbf{p}_{n})=\frac{1}{n!}\left(\frac{\textbf{p}_{\textrm{tot}}\cdot\textbf{p}_{1}}{p_{1}^{2}}\right)\cdots\left(\frac{\textbf{p}_{\textrm{tot}}\cdot\textbf{p}_{n}}{p_{n}^{2}}\right)\,.$
(2.6)

In particular, we have that

$Z_{2}(\textbf{p}_{1},\textbf{p}_{2})=\frac{1}{2}\Big{(}1+\frac{\textbf{p}_{1}\cdot\textbf{p}_{2}}{p_{1}^{2}}\Big{)}\Big{(}1+\frac{\textbf{p}_{1}\cdot\textbf{p}_{2}}{p_{2}^{2}}\Big{)}\,,$
(2.7)

such that the tree-level bispectrum of the Zeldovich matter field is given by

$B^{\rm tree}_{\rm
Zel}(\textbf{k}_{1},\textbf{k}_{2},\textbf{k}_{3})=2\,Z_{2}(\textbf{k}_{1},\textbf{k}_{2})\,P_{\rm
L}(k_{1})P_{\rm L}(k_{2})+2\,\textrm{cycle}\,.$ (2.8)

Let us compare this to Equation 2.5. To get a nonzero Fourier transform, we
have to start at the second order term when Taylor expanding the exponential

$\displaystyle B^{\rm tree}_{\rm Zel}=$
$\displaystyle\int_{\textbf{q},\textbf{r}}e^{-i\textbf{k}_{1}\cdot\textbf{q}-i\textbf{k}_{2}\cdot\textbf{r}}\
\Big{[}\Big{(}\frac{1}{2}k_{1,i}k_{3,j}A_{ij}(\textbf{q})\Big{)}\Big{(}\frac{1}{2}k_{2,i}k_{3,j}A_{ij}(\textbf{r})\Big{)}$
(2.9)
$\displaystyle+\Big{(}\frac{1}{2}k_{1,i}k_{3,j}A_{ij}(\textbf{q})\Big{)}\Big{(}\frac{1}{2}k_{1,i}k_{2,j}A_{ij}(\textbf{q}-\textbf{r})\Big{)}+\Big{(}\frac{1}{2}k_{2,i}k_{3,j}A_{ij}(\textbf{r})\Big{)}\Big{(}\frac{1}{2}k_{1,i}k_{2,j}A_{ij}(\textbf{q}-\textbf{r})\Big{)}\Big{]}\,.$

The quadratic terms shown here are symmetric under permutations of the
triangle coordinates in coordinate and Fourier space, so we need only look at
the first term, which is the simplest due to being factorizable in the
coordinates we have chosen with $\textbf{q}_{3}$ ‘at the origin’. In fact,
this first piece gives

$\int_{\textbf{q}}e^{-i\textbf{k}_{1}\cdot\textbf{q}}\
\Big{(}\frac{1}{2}k_{1,i}k_{3,j}A_{ij}(\textbf{q})\Big{)}\int_{\textbf{r}}e^{-i\textbf{k}_{2}\cdot\textbf{r}}\
\Big{(}\frac{1}{2}k_{2,i}k_{3,j}A_{ij}(\textbf{r})\Big{)}=2\,Z_{2}(\textbf{k}_{1},\textbf{k}_{2})\,P_{\rm
L}(k_{1})P_{\rm L}(k_{2})\,,$ (2.10)

where in obtaining this result we have used that
$\textbf{k}_{3}=-(\textbf{k}_{1}+\textbf{k}_{2})$. Transforming the latter two
pieces of $B^{\rm tree}_{\rm Zel}$ give the cyclically permuted results so the
full result is in agreement with equation 2.8. Note that evaluating the
Fourier transform above technically also incurs terms that are linear in the
power spectrum, e.g.

$\int_{\textbf{q},\textbf{r}}e^{-i\textbf{k}_{1}\cdot\textbf{q}-i\textbf{k}_{2}\cdot\textbf{r}}\
\left(\frac{1}{2}k_{1,i,}k_{3,j}A_{ij}(\textbf{q})\right)=\left(\frac{\textbf{k}_{1}\cdot\textbf{k}_{3}}{k_{1}^{2}}\right)P_{\rm
L}(k_{1})(2\pi)^{3}\delta_{D}(\textbf{k}_{2})\,,$

as well as a zeroth order term equal to two $\delta$ functions; this is
because we have ignored the $-1$ in Equation 2.3 which cancels out these
disconnected pieces.

##  3 Direct IR-Safe Resummation of Zeldovich Displacements

While LPT and EPT agree order by order, as we showed in the previous section
at tree level, the nonlinear structure of the exponentiated displacements LPT
has physical consequences beyond this level, particularly for capturing the
nonlinear damping of the BAO, which we will return to in §4. It is thus of
great interest to develop efficient numerical techniques to evaluate the
double integral in Equation 2.3.

###  3.1 Preliminary Steps: Generalizing the Zeldovich Power Spectrum

Motivated by previously developed techniques to efficiently compute the power
spectrum within LPT, we introduce

$\mathcal{E}_{12}(\textbf{p})\equiv\mathcal{E}(\textbf{k}_{1},\textbf{k}_{2},\textbf{p})=\int_{\textbf{q}}e^{-i\textbf{p}\cdot\textbf{q}}\,e^{+\frac{1}{2}\textbf{k}_{1,i}\textbf{k}_{2,j}A_{ij}(\textbf{q})}\,,~{}~{}{\rm
i.e.}~{}~{}e^{+\frac{1}{2}\textbf{k}_{1,i}\textbf{k}_{2,j}A_{ij}(\textbf{q})}=\int_{\textbf{p}}\,e^{i\textbf{p}\cdot\textbf{q}}\,\mathcal{E}_{12}(\textbf{p})\,,$
(3.1)

which is equal to the Zeldovich power spectrum when
$\textbf{p}=\textbf{k}_{1}=-\textbf{k}_{2}$ (i.e. for
$\textbf{k}_{1}+\textbf{k}_{2}=0$). Note the $+$ sign in the exponent
containing $A_{ij}$. Techniques for evaluating this Fourier transform in that
special case were developed in ref. [15], and we extend them to allow for
arbitrary $\textbf{k}_{1,2}$ and p in Appedix A; roughly, while in the case of
the power spectrum all three vectors $\textbf{k}_{1,2}$ and p were co-linear,
here they can form arbitrary spatial configurations, such that the result is a
function of both the individual lengths $k_{1},k_{2},p$ as well as the angles
between the vectors specified by their dot products
$\hat{k}_{1}\cdot\hat{k}_{2},\hat{k}_{1}\cdot\hat{p},\hat{k}_{2}\cdot\hat{p}$.
We find that accounting for these new degrees of freedom, for a given
bispectrum configuration specified by $\textbf{k}_{1,2}$, the integral can be
performed for all p in terms of the product of two matrices, one consisting of
radial Hankel transforms, which can be efficiently computed using the FFTLog
algorithm [39], and another depending solely on angular polynomial
coefficients.

An important point is that since the integrand in Equation 3.1 does not
asymptote to zero as $\textbf{q}\to\infty$, $\mathcal{E}$ will have both a
finite part and one proportional to a $\delta$-function, i.e.

$\mathcal{E}_{12}(\textbf{p})=\mathcal{E}_{12}^{\rm
fin}(\textbf{p})+(2\pi)^{3}\delta_{D}(\textbf{p})\
e^{\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})\Sigma^{2}}\,,$ (3.2)

where $3\Sigma^{2}$ is the mean pairwise displacement of two points separated
by $\textbf{q}\rightarrow\infty$. While we focus on computing the numerically
tricky finite piece in Appendix A, keeping track of both pieces is important
for the stability of the integral, as we will discuss below. A particularly
useful way to re-express Equation 3.2 is

$\displaystyle\mathcal{E}_{12}(\textbf{p})=\mathcal{E}^{\rm
fin}_{12}(\textbf{p})+\left(1-\int_{\textbf{p}^{\prime}}\mathcal{E}^{\rm
fin}_{12}(\textbf{p}^{\prime})\right)(2\pi)^{3}\delta_{D}(\textbf{p})\,.$
(3.3)

This can be seen by noting that, since $A_{ij}(\textbf{q}\rightarrow 0)=0$,
the $\textbf{q}=0$ Fourier transform of $\mathcal{E}$ must be equal to unity
and is again a consequence of Galilean invariance.

###  3.2 An IR-safe Integrand for the Bispectrum

Since the Zeldovich bispectrum is the Fourier transform of three exponentiated
$A_{ij}$’s, it can be written using Equation 2.5 in terms of their Fourier
transforms $\mathcal{E}$ as

$\displaystyle
B(\textbf{k}_{1},\textbf{k}_{2})=\int_{\textbf{q}_{1},\textbf{q}_{2}}\,e^{-i\textbf{k}_{1}\cdot\textbf{q}_{1}-i\textbf{k}_{2}\cdot\textbf{q}_{2}}\left\langle
e^{-i\textbf{k}_{1}\cdot\Delta_{1}-i\textbf{k}_{2}\cdot\Delta_{2}}\right\rangle=\int_{\textbf{p}}\,\mathcal{E}_{13}(\textbf{k}_{1}-\textbf{p})\mathcal{E}_{23}(\textbf{k}_{2}+\textbf{p})\mathcal{E}_{12}(\textbf{p})\,.$
(3.4)

The configuration space integral above represents a double Fourier transform,
which is difficult to evaluate numerically and thus has represented an
obstacle in any practical implementation of bispectrum LPT results.
Conversely, recasting it in the momentum space, the integral becomes a
straightforward convolution amenable to direct integration, given the
techniques to compute $\mathcal{E}$ in the previous subsection. Indeed, by
splitting up the finite and $\delta$-function pieces, we obtain

$\displaystyle\int_{\textbf{p}}\,\mathcal{E}_{13}(\textbf{k}_{1}-\textbf{p})\mathcal{E}_{23}(\textbf{k}_{2}+\textbf{p})\mathcal{E}_{12}(\textbf{p})=\,$
$\displaystyle
e^{\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})\Sigma^{2}}\mathcal{E}^{\rm
fin}_{13}(\textbf{k}_{1})\mathcal{E}^{\rm fin}_{23}(\textbf{k}_{2})+2\,{\rm
cycle}$ $\displaystyle+\int_{\textbf{p}}\,\mathcal{E}^{\rm
fin}_{13}(\textbf{k}_{1}-\textbf{p})\mathcal{E}^{\rm
fin}_{23}(\textbf{k}_{2}+\textbf{p})\mathcal{E}^{\rm fin}_{12}(\textbf{p})\,,$
(3.5)

where we have dropped any disconnected pieces requiring one of the external
momenta $\textbf{k}_{i}$ to be zero. However, as we shall show now, it is
possible to re-absorb these $\delta$-funtion contributions into the integral
to obtain a numerically better-behaved integral.

The above integral, if implemented numerically as written, is rather sensitive
to the behavior of the integrand terms when the argument of $\mathcal{E}$’s
approaches zero, where they have a pole. As a particularly simple example, let
us first look at the linearized result, where

$\mathcal{E}^{\rm
lin}_{12}(\textbf{p})=-\frac{(\textbf{k}_{1}\cdot\textbf{p})(\textbf{k}_{2}\cdot\textbf{p})}{p^{4}}P_{\rm
L}(p)+\left(1+\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})\Sigma^{2}\right)(2\pi)^{3}\delta_{D}(\textbf{p})\,,$
(3.6)

that has a pole at $p=0$. In the long-wavelength (infrared) region where
$p<k_{\rm IR}\ll k_{1,2}$ we get the contribution

$\left[1+k_{1,i}k_{2,j}\left(\frac{1}{2}\Sigma^{2}\delta_{ij}-\int_{|\textbf{p}|<k_{\rm
IR}}\frac{p_{i}p_{j}}{p^{4}}P_{\rm L}(p)\right)\right]\ \mathcal{E}_{\rm
lin}(\textbf{k}_{1},\textbf{k}_{3},\textbf{k}_{1})\mathcal{E}_{\rm
lin}(\textbf{k}_{2},\textbf{k}_{3},\textbf{k}_{2})\,,$ (3.7)

where “1” piece is simply the tree level result derived in §2.2. More
interestingly, in the parentheses we see that the infrared contribution to
Zeldovich square displacement $\Sigma^{2}$ from the $\delta$-function piece is
exactly cancelled by the equivalent infrared contribution by convolving the
finite pieces, leading to no net effect. However, as written in Equation 3.5,
this cancellation depends on the agreement between a convolution for the
finite piece and direct integration of the linear power spectrum for the
$\delta$-function piece, which is not guaranteed. While the p integral in
Equation 3.7 is not formally divergent in $\Lambda$CDM, its contribution on
typical scales for LSS analyses typically has $k^{2}\Sigma^{2}\gtrsim 1$,
requiring us to rely on the cancellation of two numbers in Equation 3.5.

Figure 1: The bispectrum integral separated into three domains according to
Equation 3.8. Each domain is defined to by proximity to one of the infrared
poles of the integral ${\bf 0},\textbf{k}_{1},-\textbf{k}_{2}$, such that they
are separated by the perpendicular bisectors of a triangle with sidelengths
$k_{1,2,3}$ in the plane formed by the poles, and join at the circumcenter. By
performing coordinate transformations separately in each domain such that the
pole is sent to the origin, we can rewrite the bispectrum as a sum of three
integrals, each with a domain defined by the $k_{1,2,3}$ triangle shifted to a
different position such that a different vertex is at the origin. We see also
that in addition to $\pm\textbf{k}_{1,2}$ some of these triangles also develop
poles at $\pm\textbf{k}_{3}$, restoring the permutation symmetry of the
bispectrum.

In order to build the infrared cancellation more explicitly into the
bispectrum, it will be convenient to slightly rearrange the arguments in
Equation 3.4. This is because the $\mathcal{E}$’s in the integrand have
arguments $\textbf{p},\textbf{k}_{1}-\textbf{p},\textbf{k}_{2}+\textbf{p}$,
such that the integrand has “infrared” regions about three points ${\bf
0},\textbf{k}_{1}$ and $-\textbf{k}_{2}$. However, much like in the case of
the power spectrum [40] and 1-loop bispectrum [36] it is possible to isolate
each singularity in the integrand and map them onto the origin to yield a more
symmetric expression. Specifically, we make use of the fact that the locus of
points closest to each pole partitions the space into three non-overlapping
regions (Fig. 1) to split the integral as

$\displaystyle\int_{\textbf{p}}$ $\displaystyle=\int_{p\ <\
|\textbf{k}_{1}-\textbf{p}|,\
|\textbf{k}_{2}+\textbf{p}|}+\int_{|\textbf{k}_{1}-\textbf{p}|\ <\ p,\
|\textbf{k}_{2}+\textbf{p}|}+\int_{|\textbf{k}_{2}+\textbf{p}|\ <\
|\textbf{k}_{1}-\textbf{p}|,\ p}$ $\displaystyle=\int_{p\ <\
|\textbf{k}_{1}-\textbf{p}|,\ |\textbf{k}_{2}+\textbf{p}|}+\int_{p^{\prime}\
<\ |\textbf{k}_{1}-\textbf{p}^{\prime}|,\
|\textbf{k}_{3}+\textbf{p}^{\prime}|}+\int_{p^{\prime\prime}\ <\
|\textbf{p}^{\prime\prime}-\textbf{k}_{2}|,\
|\textbf{k}_{3}+\textbf{p}^{\prime\prime}|}\,.$ (3.8)

In the second line we have for the latter two integrals performed the
coordinate transformations $\textbf{p}^{\prime}=\textbf{k}_{1}-\textbf{p}$ and
$\textbf{p}^{\prime\prime}=\textbf{k}_{2}+\textbf{p}$. This corresponds to a
Voronoi tessellation of the poles in the integrand and can be extended also to
higher N-point functions (§5); the above decomposition was also used in EPT by
ref. [36] for the $B_{222}$ contribution to the 1-loop bispectrum.222See also
ref. [41] who use a slightly different domain decomposition to remove the same
IR poles. We see that each of the singularities in the bispectrum integrand
exists in one of these partitions, which we have mapped onto
$\textbf{p},\textbf{p}^{\prime},\textbf{p}^{\prime\prime}=\bf{0}$, such that
in each domain the closest singularity is the one at the origin. These domains
have boundaries set by the perpendicular bisectors formed by the triangle
formed by the poles with side lenghts $k_{1,2,3}$ and meet at the
circumcenter. Rewriting the above in terms of $\Theta$ functions, the
bispectrum integral is now

$\displaystyle B(\textbf{k}_{1},\textbf{k}_{2})=\int_{\textbf{p}}$
$\displaystyle\mathcal{E}_{12}(\textbf{p})\
\mathcal{E}_{13}(\textbf{k}_{1}-\textbf{p})\
\mathcal{E}_{23}(\textbf{k}_{2}+\textbf{p})\
\Theta(|\textbf{k}_{1}-\textbf{p}|-p)\ \Theta(|\textbf{k}_{2}+\textbf{p}|-p)$
$\displaystyle+\ \ $
$\displaystyle\mathcal{E}_{12}(\textbf{k}_{1}-\textbf{p})\
\mathcal{E}_{13}(\textbf{p})\ \mathcal{E}_{23}(-\textbf{k}_{3}-\textbf{p})\
\Theta(|\textbf{k}_{1}-\textbf{p}|-p)\ \Theta(|\textbf{k}_{3}+\textbf{p}|-p)$
$\displaystyle+\ \ $
$\displaystyle\mathcal{E}_{12}(-\textbf{k}_{2}+\textbf{p})\
\mathcal{E}_{13}(-\textbf{k}_{3}-\textbf{p})\ \mathcal{E}_{23}(\textbf{p})\
\Theta(|\textbf{k}_{2}-\textbf{p}|-p)\
\Theta(|\textbf{k}_{3}+\textbf{p}|-p)\,.$ (3.9)

This integral can be made manifestly symmetric in $\textbf{k}_{1,2,3}$ by
using that $\mathcal{E}_{nm}$ is an even function of p, but we will not
perform this additional step here for the sake of brevity.

We can now write this integrand in an IR safe way by explicitly re-introducing
the $\delta$-function piece. Plugging Equation 3.3 back into the bispectrum
integral (Eq. 3.4) we have

$\displaystyle B(\textbf{k}_{1},\textbf{k}_{2})$
$\displaystyle=\int_{\textbf{p}}\mathcal{E}_{12}^{\rm
fin}(\textbf{p})\mathcal{E}_{13}(\textbf{k}_{1}-\textbf{p})\mathcal{E}_{23}(\textbf{k}_{2}+\textbf{p})+\left(1-\int_{\textbf{p}}\mathcal{E}_{12}^{\rm
fin}(\textbf{p})\right)\mathcal{E}_{13}(\textbf{k}_{1})\mathcal{E}_{23}(\textbf{k}_{2})$
$\displaystyle=\mathcal{E}_{13}(\textbf{k}_{1})\mathcal{E}_{23}(\textbf{k}_{2})+\int_{\textbf{p}}\mathcal{E}_{12}^{\rm
fin}(\textbf{p})\left(\mathcal{E}_{13}(\textbf{k}_{1}-\textbf{p})\mathcal{E}_{23}(\textbf{k}_{2}+\textbf{p})-\mathcal{E}_{13}(\textbf{k}_{1})\mathcal{E}_{23}(\textbf{k}_{2})\right)\,.$
(3.10)

Both pieces of the integral are now IR safe: the product in front is since it
is derived purely from pairwise displacements, while the integrand in the
second piece vanishes as $p\rightarrow 0$. From this we can see that the full
bispectrum integrand is given by

$\displaystyle B(\textbf{k}_{1},\textbf{k}_{2})$ $\displaystyle=\
\mathcal{E}_{13}^{\rm fin}(\textbf{k}_{1})\mathcal{E}_{23}^{\rm
fin}(\textbf{k}_{2})+\mathcal{E}_{13}^{\rm
fin}(\textbf{k}_{3})\mathcal{E}_{12}^{\rm
fin}(\textbf{k}_{2})+\mathcal{E}_{12}^{\rm
fin}(\textbf{k}_{1})\mathcal{E}_{23}^{\rm fin}(\textbf{k}_{3})$
$\displaystyle+\int_{\textbf{p}}\Big{\\{}\mathcal{E}_{12}^{\rm
fin}(\textbf{p})\mathcal{E}_{13}^{\rm
fin}(\textbf{k}_{1}-\textbf{p})\mathcal{E}_{23}^{\rm
fin}(\textbf{k}_{2}+\textbf{p})\
\Theta(|\textbf{k}_{1}-\textbf{p}|-p)\Theta(|\textbf{k}_{2}+\textbf{p}|-p)$
$\displaystyle+\ \ \mathcal{E}^{\rm
fin}_{12}(\textbf{k}_{1}-\textbf{p})\mathcal{E}^{\rm
fin}_{13}(\textbf{p})\mathcal{E}_{23}^{\rm fin}(-\textbf{k}_{3}-\textbf{p})\
\Theta(|\textbf{k}_{1}-\textbf{p}|-p)\Theta(|\textbf{k}_{3}+\textbf{p}|-p)$
$\displaystyle+\ \ \mathcal{E}_{12}^{\rm
fin}(-\textbf{k}_{2}+\textbf{p})\mathcal{E}_{13}^{\rm
fin}(-\textbf{k}_{3}-\textbf{p})\mathcal{E}_{23}^{\rm fin}(\textbf{p})\
\Theta(|\textbf{k}_{2}-\textbf{p}|-p)\Theta(|\textbf{k}_{3}+\textbf{p}|-p)$
$\displaystyle-\ \ \mathcal{E}_{12}^{\rm fin}(\textbf{p})\mathcal{E}_{13}^{\rm
fin}(\textbf{k}_{1})\mathcal{E}_{23}^{\rm
fin}(\textbf{k}_{2})-\mathcal{E}_{23}^{\rm
fin}(\textbf{p})\mathcal{E}_{13}^{\rm
fin}(\textbf{k}_{3})\mathcal{E}_{12}^{\rm
fin}(\textbf{k}_{2})-\mathcal{E}_{13}^{\rm
fin}(\textbf{p})\mathcal{E}_{12}^{\rm
fin}(\textbf{k}_{1})\mathcal{E}_{23}^{\rm fin}(\textbf{k}_{3})\Big{\\}}\,,$
(3.11)

such that the poles in each region picked out by the $\Theta$ functions are
explicitly cancelled by the poles at the origin in the last line.

The above also enables us to write down an infrared-safe expression for the
1-loop Zeldovich bispectrum. Splitting $\mathcal{E}^{\rm
fin}=\mathcal{E}^{(1)}+\mathcal{E}^{(2)}+...$ into its $\mathcal{O}(P_{\rm
L}^{n})$ contributions we have that the tree level bispectrum is obtained by
the replacement to the first line of Equation 3.11 with $B_{\rm
tree}\sim\mathcal{E}^{(1)}\mathcal{E}^{(1)}$ while the 1-loop contributions
are given by the contributions $B_{\rm
1-loop}\sim\mathcal{E}^{(1)}\mathcal{E}^{(2)}$ and
$\int_{\textbf{p}}\mathcal{E}^{(1)}\mathcal{E}^{(1)}\mathcal{E}^{(1)}$.
Appendix A details how the order $n$ contributions to $\mathcal{E}$ can be
isolated in numerical calculations and efficiently calculated. Since the
infrared safety of Equation 3.11 holds order by order, the equivalent
expressions here at 1-loop are manifestly infrared safe as well. This is in
contrast to the IR-safe integrands for the 1-loop EPT bispectrum developed in
ref. [36, 41] where the integrals of all the loop diagrams have to be combined
after different remappings to cancel IR divergences. In Appendix B we
diagrammatically review the 1-loop contributions to the bispectrum, showing
why unlike in our method where each piece is individually IR safe each of the
diagrams shown in Figure 2 contain IR contributions which have to be combined
together to cancel in an IR-safe way.

Figure 2: Contributions to the 1-loop bispectrum. Each of the first three
diagrams is required to cancel the IR contribution to the fourth.

###  3.3 Numerical Implementation

We are now in a position to numerically investigate the Zeldovich bispectrum.
In order to do so we have developed
triceratops333https://github.com/sfschen/triceratops a Python code that
computes the Zeldovich bispectrum in various approximations (tree-level,
1-loop, full resummation) using the methods described in this section.
triceratops computes $\mathcal{E}$ efficiently using the FFTLog algorithm [39]
and performs the convolution integral using a three-dimensional, spherical
grid in p that is logarithmically spaced in $p$ and use Gaussian weights for
$\theta$ and $\phi$. As an additional check, we also developed in parallel a
Mathematica code to perform the same calculations using both the in-built
integration routines as well as the vegas algorithm [42]—while not as fast as
triceratops, the two codes are in excellent agreement, showing that the
underlying integrals are well-behaved.

Figures 3 and 4 show the bispectrum computed using the methods described above
in the squeezed and equilateral configurations. Reassuringly, both the full
and 1-loop Zeldovich calculations approach the leading-order tree-level
calculation at low wavenumbers, while differing at higher ones, with the
1-loop and full Zeldovich calculations agreeing with each other to slightly
higher wavenumbers than the tree-level ones.

Figure 4 also shows the full Zeldovich calculation using both the naive (Eqn.
3.5) and IR safe (Eqn. 3.11) integrands. Both are in excellent agreement for
the settings shown in the plot. However, in order to compare their behavior
when the sampling of small wavenumbers is varied, we compute the 1-loop
contribution to the bispectrum in Figure 5 using both methods—here it is quite
clear that the IR-safe integrand is far less sensitive to the scales sampled,
recovering the correct result even when the smallest wavenumber sampled
$k_{\rm min}$ is quite close to the wavenumber arguments in the bispectrum
itself. This is of course as expected since the IR-safe integrand has
explicitly vanishing contributions from long-wavelength modes. Finally, we
have also tested the sensitivity of our bispectrum calculation to the maximum
wavenumber $k_{\rm max}$ sampled in the integral, as well as the number of
Hankel transforms required when calculation $\mathcal{E}$—we find that the
algorithm is relatively lenient in these requirements, since the dependence on
$\mathcal{E}$ is not very spread out in Fourier space but rather localized to
physically-relevant wavenumbers. In particular, each $\mathcal{E}$ is
sufficiently captured by including up to five spherical-bessel transforms,
much like is typically needed for the power spectrum on perturbative scales,
as long as the bispectrum is evaluated also on these scales.

Figure 3: The tree level, 1-loop Zeldovich and full Zeldovich matter bispectra
for triangles with longest sidelengths of $k=0.10$ (left) and $k=0.20\,h\,{\rm
Mpc}^{-1}$ (right) at $z=0$. Each row shows the logarithm of the tree level
bispectrum in $\,h^{-3}\,{\rm Mpc}^{3}$ units (top) and the ratios of the loop
(middle) and full Zeldovich (bottom) bispectra to the tree-level prediction.
Figure 4: The $z=0$ matter Zeldovich bispectrum in the squeezed (left) and
equilateral configurations (right), computed using the methods described in
§3, as a function of side length. The blue and orange lines show the
calculation performed using Equation 3.5 and the IR-safe Equation 3.11,
respectively, which are in excellent agreement when the numerics are properly
converged. The dashed line shows the 1-loop bispectrum, which agrees with the
full Zeldovich calculation at low $k$. The squeezed configuration is taken to
be one where the short side is $5\%$ that of the longer sides. Figure 5:
Numerical stability of the bispectrum calculation to the minimum wavenumber
sampled in the integral $k_{\rm min}$, in the squeezed (left) and equilateral
(right) configurations. Dashed and solid lines show the 1-loop contribution to
the Zeldovich bispectrum using the naive and IR-safe integrands, respectively,
with the latter converging to the correct result faster at higher $k_{\rm
min}$ as expected.

##  4 Infrared Resummation of BAO and the Saddle-Point Approximation

###  4.1 Review: Resummation for the Power Spectrum

In the preceding sections we developed the formalism to compute the bispectrum
of galaxies in Lagrangian perturbation theory with the linear displacements
fully resummed. In practice, it is often desirable to look specifically at the
effect of these displacements on the BAO feature, where their effects are most
enhanced. In order to do so it is useful to isolate the BAO feature in the
linear spectrum by performing a “wiggle no-wiggle split” $P_{\rm L}=P_{\rm
L}^{nw}+f_{b}P_{\rm L}^{w}$ [20] where $f_{b}P_{\rm L}^{w}$ is a purely
oscillatory component carrying the BAO while $P_{\rm L}^{nw}$ is the smooth
component of the power spectrum. We have included the baryon fraction $f_{b}$
as a counting parameter since we will be interested primarily in working at
first order in the BAO, but we will keep it only implicitly below, with the
understanding that each superscript $w$ carries a factor of $f_{b}$. The
definition of $P_{\rm L}^{nw}$ is inexact and many algorithms to compute it
exist in the literature with differing results (see e.g. [43] for a review of
commonly used methods). Throughout the numerical portions of this work we will
adopt the B-spline method in ref. [21], which has the nice property that the
total variance in linear displacements is preserved, though we note that the
exact choice is somewhat a matter of taste and an inherent systematic of
calculations involving this split.

The most prominent effect of the exponentiated linear displacements in the
power spectrum is to produce a roughly Gaussian damping of the BAO feature,
which can therefore be resummed by including this Gaussian factor whenever the
wiggle piece of the power spectrum appears (RWiggle). This effect has long
been known empirically in the large-scale structure literature [20, 44, 12]
and was derived analytically in perturbation theory as an infrared resummation
of long-wavelength modes enhanced by the BAO by refs. [17, 18, 21, 22]. An
equivalent derivation, presented in ref. [21], shows that this form can be
simply derived by performing a saddle-point approximation about the BAO in
Lagrangian perturbation theory, and we briefly review it below.

Using the wiggle no-wiggle split we can separate
$A_{ij}=A^{nw}_{ij}+A^{w}_{ij}$ into two pieces due to the wiggle and no-
wiggle power spectra

$A_{ij}^{w,nw}(\textbf{q})=2\int_{\textbf{k}}\left(1-e^{i\textbf{k}\cdot\textbf{q}}\right)\left(\frac{k_{i}k_{j}}{k^{4}}\right)P_{\rm
L}^{w,nw}(k)\,.$ (4.1)

The former will be strongly peaked at the BAO scale $q=r_{s}$ while the latter
will be essentially a smooth version of $A_{ij}$ without the BAO feature.
Performing this split we have that the Zeldovich matter power spectrum for
matter, at first order in the wiggles, is

$\displaystyle P^{\rm Zel}_{w}(k)$
$\displaystyle=-\frac{1}{2}k_{i}k_{j}\int_{\textbf{q}}e^{i\textbf{k}\cdot\textbf{q}-\frac{1}{2}k_{i}k_{j}A^{nw}_{ij}(\textbf{q})}A^{w}_{ij}(\textbf{q})$
$\displaystyle\approx-\frac{1}{2}k_{i}k_{j}e^{-\frac{1}{2}k^{2}X_{s}}\int_{\textbf{q}}e^{i\textbf{k}\cdot\textbf{q}-\frac{1}{2}k_{i}k_{j}Y_{s}\hat{q}_{i}\hat{q}_{j}}A^{w}_{ij}(\textbf{q})$
$\displaystyle\approx e^{-\frac{1}{2}k^{2}\Sigma^{2}_{s}}P^{w}_{\rm
L}(k),\quad kr_{s}\gg 1\,.$ (4.2)

In going from the first to second lines, we have used the saddle-point444Here
by saddle point we refer to the approximation where when integrating the
product of two functions $f,g$, if $f$ is narrowly peaked at $x_{0}$ we can
have  $\int dx\ f(x)g(x)=g(x_{0})\int dx\ f(x)+\mathcal{O}(\sigma)$ (4.3)
where $\sigma$ is the width of $f$ (Appendix C).  approximation to substitute
the values of $X^{nw}(q),\ Y^{nw}(q)$ at the BAO scale $r_{s}$, which we have
denoted using subscripts $X_{s}$, $Y_{s}$, rather than integrate over their
full $q$ dependence. Corrections to this approximation can be computed by
expanding the smooth component about $r_{s}$, i.e. by Taylor expanding
$e^{-\frac{1}{2}k_{i}k_{j}A^{nw}_{ij}}$ in $\Delta q=q-r_{s}$, with the
correction scaling as powers of the BAO width $r_{D}$ multiplying the mean
square density contrast on BAO scales. We provide further details in Appendix
C.

The remaining subtlety is the angular integral in the second line of Equation
4.2. In Appendix D, we show that the angular integral in the penultimate line
can be performed exactly in the limit that $kr_{s}\gg 1$ to yield the final
expression with

$\Sigma^{2}_{s}=X_{s}+Y_{s}=\int\frac{dk}{3\pi^{2}}\ P_{\rm L}(k)\
\left(1-j_{0}(kr_{s})+2j_{2}(kr_{s})\right)\,.$ (4.4)

At intermediate $k$ this approximation breaks down, leading to corrections to
the simple exponential damping form. The $1-j_{0}+2j_{2}$ in the parentheses
cuts off contributions from modes with wavelengths longer than the BAO scale,
such that coherent displacements that move the BAO, i.e. act locally as a
Galilean transformation, cannot damp it. As shown in Figure 6, RWiggle
describes the wiggles in the power spectrum almost perfectly even when
everything else is kept to linear order. It is important to note that the
above is a description of the parametrically large effects of linear
displacements on the BAO feature—beyond them nonlinear effects and effective-
theory corrections, which do not have to be limited to the Gaussian form
above, can also play a role. These effects have been extensively studied in
the perturbation-theory literature so we will not describe them in any detail
in what follows, except to emphasize that even a complete model of the linear
displacements does not fully specify the nonlinear shape of the wiggles.

We can also connect the saddle-point calculation above to the resummation
developed in §3. Specifically, we can rewrite the integral in the first line
of Equation 4.2 as

$\displaystyle P^{NL}_{w}(k)$
$\displaystyle=\int_{\textbf{p}}\left\\{-\frac{1}{2}k_{i}k_{j}\tilde{A}^{w}_{ij}(\textbf{k}-\textbf{p})\times\text{FT}\left[\exp\left(-\frac{1}{2}k_{i}k_{j}A^{nw}_{ij}(\textbf{q})\right)\right](\textbf{p})\right\\}$
$\displaystyle=\int_{\textbf{p}}\left\\{-\frac{1}{2}k_{i}k_{j}\tilde{A}^{w}_{ij}(\textbf{k}-\textbf{p})\times\mathcal{E}^{nw}(\textbf{k},-\textbf{k},\textbf{p})\right\\}$
$\displaystyle=e^{-\frac{1}{2}k^{2}\Sigma^{2}}P^{w}_{\rm
L}(k)+\int_{\textbf{p}}\left\\{-\frac{1}{2}k_{i}k_{j}\tilde{A}^{w}_{ij}(\textbf{k}-\textbf{p})\times\mathcal{E}^{nw}_{\rm
fin}(\textbf{k},-\textbf{k},\textbf{p})\right\\}\,,$ (4.5)

where $\mathcal{E}^{nw}$ is computed using the no-wiggle spectrum only and
where in the last line we have used the finite-infinite split in Equation 3.2
as required when numerically evaluating $\mathcal{E}$. The first term looks
superficially similar to the damping form derived in Equation 4.2, with
$\Sigma^{2}$ replacing $\Sigma^{2}_{s}$. However, it would be incorrect and
unphysical to assume that $P_{w}$ is damped by the exponential in the
$\delta$-function piece, which is what we would get if we dropped the finite
piece.

It is instructive to recover the result for the dampened wiggle obtained in
Equation 4.2 starting from the momentum space representation given in Equation
4.5. We first note that there are two potentially IR divergent contributions
arising from the poles $\textbf{p}\to 0$ and $\textbf{p}\to\textbf{k}$. The
latter does not, in fact, contribute given that the wiggle power spectrum has
no support in that regime. Thus, only the $\textbf{p}\to 0$ pole remains and
we can simply regularise the integral

$\displaystyle P^{NL}_{w}(k)$
$\displaystyle=e^{-\frac{1}{2}k^{2}\Sigma^{2}}P^{w}_{\rm
L}(k)+\frac{1}{2}k_{i}k_{j}\tilde{A}^{w}_{ij}(\textbf{k})\int_{\textbf{p}}\mathcal{E}^{nw}_{\rm
fin}(\textbf{k},-\textbf{k},\textbf{p})$ $\displaystyle\hskip
73.97733pt-\frac{1}{2}k_{i}k_{j}\int_{\textbf{p}}\left(\tilde{A}^{w}_{ij}(\textbf{k}-\textbf{p})-\tilde{A}^{w}_{ij}(\textbf{k})\right)\mathcal{E}^{nw}_{\rm
fin}(\textbf{k},-\textbf{k},\textbf{p})$ $\displaystyle=P^{w}_{\rm
L}(k)-\frac{1}{2}k_{i}k_{j}\int_{\textbf{p}}\left(\tilde{A}^{w}_{ij}(\textbf{k}-\textbf{p})-\tilde{A}^{w}_{ij}(\textbf{k})\right)\mathcal{E}^{nw}_{\rm
fin}(\textbf{k},-\textbf{k},\textbf{p})\,,$ (4.6)

where we used that $\int_{\textbf{p}}\mathcal{E}^{nw}_{\rm
fin}(\textbf{k},-\textbf{k},\textbf{p})=1-\exp\left(-k^{2}\Sigma^{2}\right)$.
This allows us to write an IR safe expression where we have explicitly
cancelled the large bulk displacements, which do not therefore affect the
shape of the wiggle power spectrum, as is expected from the equivalence
principle and Galilean invariance. To derive the physical effect of long
displacements on wiggles, we use the fact that the BAO bump is localized in
the configuration space and write

$\displaystyle-\frac{1}{2}k_{i}k_{j}$
$\displaystyle\int_{\textbf{p}}\left(\tilde{A}^{w}_{ij}(\textbf{k}-\textbf{p})-\tilde{A}^{w}_{ij}(\textbf{k})\right)\
\mathcal{E}^{nw}_{\rm fin}(\textbf{k},-\textbf{k},\textbf{p})$
$\displaystyle=-\frac{1}{2}k_{i}k_{j}\int_{\textbf{q},\textbf{p}}e^{i\textbf{k}\cdot\textbf{q}}\left(e^{i\textbf{p}\cdot\textbf{q}}-1\right)A^{w}_{ij}(\textbf{q})\mathcal{E}^{nw}_{\rm
fin}(\textbf{k},-\textbf{k},\textbf{p})\approx\left(e^{-\frac{1}{2}k^{2}\Sigma^{2}_{s}}-1\right)P^{w}_{\rm
L}(k)\,.$ (4.7)

In order to resum only IR modes, we could have also explicitly restricted the
p integral, i.e. modes for which we have $p\ll\Lambda\lesssim k$, which would
have further restricted the contributions to $\Sigma_{s}$. As can be seen in
the configuration and Fourier-space derivations above, this is not required at
this stage to derive RWiggle for the power spectrum. However, these IR modes
come into play if we want to show that the resummation is under control and
that the subleading contributions can indeed be neglected (see Appendixes C
and D).

Figure 6: The wiggle component of the Zeldovich power spectrum, computed in
full (blue), at linear order (orange) and in the RWiggle scheme. The Zeldovich
prediction is suppressed relative to the linear one in a way that is very well
described by a Gaussian damping.

###  4.2 Saddle-Point IR Resummation of the Bispectrum BAO

Figure 7: Triangle configuration relevant for the IR resummation of the BAO
feature in Lagrangian space, where two points $\textbf{q}_{1,3}$ are separated
by $q=r_{s}$ equal to the BAO radius.

Let us now extend the above analysis and extract the IR-resummed wiggle
contributions to the bispectrum (RWiggle). The final result shows that the
largest effects, due to long displacements, can be captured by linearized
wiggle contributions to the linear power spectrum. We can derive this result,
as in the case of the 2-point function, using a saddle-point approximation. We
do so in two related ways, first (more restrictively) by focusing on the
contributions of strictly long-wavelength modes in configuration space and
subsequently relaxing this condition and showing how the full linear Zeldovich
displacements can be resummed and the residuals computed in an IR-safe way.
Our result, up to the linear terms in $P_{w}$, corresponds to the case where
the effects of the long displacements are captured by substituting the damped
linear power spectrum (Eqn. 4.2) for the linear power spectrum. The latter
result was derived diagrammatically by ref. [22] using time-sliced
perturbation theory (TSPT).

Let us begin, however, by considering the leading contribution expanding in
the wiggle power spectrum $P_{w}$ to $\mathcal{E}$. Since each $A_{ij}$ in
Equation 2.5, or each $\mathcal{E}$ in Equation 3.4, can have a power of
$A_{ij}^{w}$ expanded out, we have three distinct contributions in total
proportional to $P_{w}(k_{1,2,3})$. In particular, we can use the same logic
as in the previous section to write each as

$\mathcal{E}_{w}(\textbf{k}_{1},\textbf{k}_{2},\textbf{p})=\int_{\textbf{q}}e^{i\textbf{p}\cdot\textbf{q}+\frac{1}{2}k_{1,i}k_{2,j}A^{nw}_{ij}(\textbf{q})}\left(\frac{1}{2}k_{1,i}k_{2,j}A^{w}_{ij}(\textbf{q})\right)\approx
e^{\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})\Sigma^{2}_{s}}A(\textbf{k}_{1},\textbf{k}_{2},\textbf{p})P_{\rm
L}^{w}(p)\,,$ (4.8)

where we have defined the angular factor
$A(\textbf{k}_{n},\textbf{k}_{m},\textbf{p})=-(\textbf{k}_{n}\cdot\textbf{p})(\textbf{k}_{m}\cdot\textbf{p})/p^{4}$.
Note that the damping of the BAO in $\mathcal{E}$ is itself IR safe, i.e. does
not depend on displacements from modes longer than the BAO radius, since it is
defined from pairwise displacements.

We can now proceed to consider BAO oscillations in the bispectrum itself. We
will do this in two ways—first, by considering explicitly the effect of long-
wavelength modes in (Lagrangian) configuration space and deriving the
bispectrum version of RWiggle through a saddle-point calculation therein, and
then by using an IR-safe integrand in Fourier space to extract the effect of
all linear displacements and derive a convenient expression for evaluating the
residual term.

####  4.2.1 Infrared Modes and Configuration-Space Saddle Point

We first consider the effect of purely infrared modes—that is, ones shorter
than the wavenumbers entering a given bispectrum configuration, starting from
Equation 3.11. Specifically, in each of the three pieces of the integral we
can extract contributions due to p with $p$ smaller than $k_{1,2}/2$,
$k_{1,3}/2$ or $k_{2,3}/2$. These regions are naturally set by the angular
bisectors defining the IR regions. Focusing on one of the three domains
without lack of generality, we have in the first one

$\displaystyle(\text{IR\ Region of }k_{1,2})$
$\displaystyle=\int_{|\textbf{p}|<k_{1,2}/2}\mathcal{E}_{12}^{\rm
fin}(\textbf{p})\left(\mathcal{E}_{13}^{\rm
fin}(\textbf{k}_{1}-\textbf{p})\mathcal{E}_{23}^{\rm
fin}(\textbf{k}_{2}+\textbf{p})-\mathcal{E}_{13}^{\rm
fin}(\textbf{k}_{1})\mathcal{E}_{23}^{\rm fin}(\textbf{k}_{2})\right)$
$\displaystyle=\int_{\textbf{q}_{1,2}}\
e^{-i\textbf{k}_{1}\cdot\textbf{q}_{1}-i\textbf{k}_{2}\cdot\textbf{q}_{2}}\
\left(\tilde{\mathcal{E}}^{\rm
fin}_{12}(\textbf{q}_{1}-\textbf{q}_{2})-\tilde{\mathcal{E}}^{\rm
fin}_{12}(\bf{0})\right)_{\rm IR}\tilde{\mathcal{E}}^{\rm
fin}_{13}(\textbf{q}_{1})\tilde{\mathcal{E}}^{\rm fin}_{23}(\textbf{q}_{2})$
$\displaystyle=\int_{\textbf{q}_{1,2}}\
e^{-i\textbf{k}_{1}\cdot\textbf{q}_{1}-i\textbf{k}_{2}\cdot\textbf{q}_{2}}\
\left(e^{\frac{1}{2}k_{1,i}k_{2,j}A_{ij}(\textbf{q}_{1}-\textbf{q}_{2})}-1\right)_{\rm
IR}\tilde{\mathcal{E}}^{\rm fin}_{13}(\textbf{q}_{1})\tilde{\mathcal{E}}^{\rm
fin}_{23}(\textbf{q}_{2})\,,$ (4.9)

where the subscript denotes that all the non-IR modes are removed, and in
going from the second to third lines, we have used the integral relation in
Equation 3.3. We can now add in the “tree-level” piece, which cancels the
“minus one”:

$\displaystyle\mathcal{E}^{\rm fin}_{13}(\textbf{k}_{1})\mathcal{E}^{\rm
fin}_{23}(\textbf{k}_{2})$ $\displaystyle+(\text{IR\ Region of }k_{1,2})$
$\displaystyle=\int_{\textbf{q}_{1,2}}\
e^{-i\textbf{k}_{1}\cdot\textbf{q}_{1}-i\textbf{k}_{2}\cdot\textbf{q}_{2}}\
\left(e^{\frac{1}{2}k_{1,i}k_{2,j}A_{ij}(\textbf{q}_{1}-\textbf{q}_{2})}\right)_{\rm
IR}\tilde{\mathcal{E}}^{\rm fin}_{13}(\textbf{q}_{1})\tilde{\mathcal{E}}^{\rm
fin}_{23}(\textbf{q}_{2})\,.$ (4.10)

By construction, the BAO wiggle piece needs to come from the last two factors
in the integrand and not the infrared first piece; splitting the two
$\mathcal{E}$ into a wiggle and no-wiggle piece yields the tree level power
spectrum contributions

$\displaystyle\mathcal{E}_{13}^{w,\rm
fin}(\textbf{k}_{1})\mathcal{E}_{23}^{nw,\rm fin}(\textbf{k}_{2})+$
$\displaystyle\mathcal{E}_{13}^{nw,\rm
fin}(\textbf{k}_{1})\mathcal{E}_{23}^{w,\rm fin}(\textbf{k}_{2})$
$\displaystyle=Z_{2}(\textbf{k}_{1},\textbf{k}_{2})\left(P^{w}_{\rm
L}(k_{1})P^{nw}_{\rm L}(k_{2})+P^{nw}_{\rm L}(k_{1})P^{w}_{\rm
L}(k_{2})\right)+\mathcal{O}(P_{\rm L}^{3})\,,$ (4.11)

yielding six total such contributions when all three IR domains of Equation
3.11 are taken into account.

We now consider the BAO saddle-point of Equation 4.10. Isolating out the
contribution containing $\mathcal{E}^{w}_{13}(\textbf{q}_{1})$, we see that
its support is sharply peaked around $x_{1}=r_{d}$, as shown in Figure 7.
Unlike in the power spectrum case, however, we have to contend with the
exponentiated factor running over $\textbf{q}_{1}-\textbf{q}_{2}$, the latter
of which is not fixed; however, if we fix $\textbf{q}_{2}=0$ in the infrared
piece, we get

$\displaystyle\mathcal{E}_{23}^{nw,\rm
fin}(\textbf{k}_{2})\int_{\textbf{q}_{1}}\
e^{-i\textbf{k}_{1}\cdot\textbf{q}_{1}}\left(e^{\frac{1}{2}k_{1,i}k_{2,j}A_{ij}(\textbf{q}_{1})}\right)_{\rm
IR}\tilde{\mathcal{E}}^{w,\rm fin}_{13}(\textbf{q}_{1})$ $\displaystyle\approx
e^{\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})\Sigma^{2}_{s}}\mathcal{E}_{13}^{w,\rm
fin}(\textbf{k}_{1})\mathcal{E}_{23}^{nw,\rm fin}(\textbf{k}_{2})$ (4.12)
$\displaystyle\approx
e^{-\frac{1}{2}k_{1}^{2}\Sigma^{2}_{s}}A(\textbf{k}_{1},\textbf{k}_{2},\textbf{p})P_{\rm
L}^{w}(p)\mathcal{E}_{23}^{nw,\rm fin}(\textbf{k}_{2})\,,$

i.e. that the tree-level bispectrum is damped by IR displacements in a
Gaussian way depending only on the wavenumber of the BAO wiggle itself.

The corrections due to this approximation come from the statistics of the
other legs of the 3-point function triangle and can be computed as terms in a
Taylor series

$\int_{\textbf{q}_{2}}e^{-i\textbf{k}_{2}\cdot\textbf{q}_{2}}(\textbf{q}_{2,i_{1}}...\textbf{q}_{2,i_{n}})\
\tilde{\mathcal{E}}^{nw,\rm fin}_{23}(\textbf{q}_{2})\int_{\textbf{q}_{1}}\
e^{-i\textbf{k}_{1}\cdot\textbf{q}_{1}}\
\left(\frac{\partial^{n}}{\partial\textbf{q}_{1,i_{1}}...\partial\textbf{q}_{1,i_{n}}}\right)\left(e^{\frac{1}{2}k_{1,i}k_{2,j}A_{ij}(\textbf{q}_{1})}\right)_{\rm
IR}\tilde{\mathcal{E}}^{w,\rm fin}_{13}(\textbf{q}_{1})\,,$

where the first factor is simply the $\textbf{k}_{2}$ derivative of
$\mathcal{E}^{nw,\rm fin}_{23}$. The second factor involves derivatives
evaluated about $x_{1}=r_{d}$; these take the form

$\left(\frac{\partial^{n}}{\partial\textbf{q}_{1,i_{1}}...\partial\textbf{q}_{1,i_{n}}}\right)\left(e^{\frac{1}{2}k_{1,i}k_{2,j}A_{ij}(\textbf{q}_{1})}\right)_{x_{1}=r_{d}}\sim\frac{1}{2}k_{1,i}k_{2,j}\int_{\textbf{p},\
p<1/r_{d}}\frac{p_{i}p_{j}p_{i_{1}}...p_{i_{n}}}{p^{4}}P_{\rm L}(p)\,.$ (4.13)

In a power-law universe where $P_{\rm L}=Ap^{n_{s}}$ it is easy to see that
this integral scales as $r_{d}^{-n}$ with $n$, such that the nth order
contribution is suppressed by $(k_{2}r_{d})^{-n}$. We thus see that the
damping can be captured by a simple damping of the input linear power
spectrum, with perturbatively small corrections due to modes with wavelengths
on the BAO scale.

####  4.2.2 Fourier Space Resummation of Linear Displacements

Here, we show that the equivalent result as derived above can be obtained from
the Fourier space representation of the bispectrum, allowing us to relax the
requirement of restricting to the modes smaller than $k_{1,2}/2$. We thus
start with the expression for the bispectrum (still working in Zeldovich
approximation) given in Equation 3.4. Since we are free to shift the
integration variable, we can recover any relevant momenta combination by
considering only the case when the wiggle contribution of $\mathcal{E}^{\rm
w}$ appears as

$\displaystyle\int_{\textbf{p}}\,\mathcal{E}^{\rm
nw}_{13}(\textbf{k}_{1}-\textbf{p})\mathcal{E}^{\rm
nw}_{23}(\textbf{k}_{2}+\textbf{p})\mathcal{E}^{\rm w}_{12}(\textbf{p})\,.$
(4.14)

Given that the remaining two contributions can be recovered by transforming
$\textbf{p}\to-\textbf{k}_{2}+\textbf{p}$, and
$\textbf{p}\to\textbf{k}_{1}-\textbf{p}$, we can write the linearized wiggle
contribution to the bispectrum as

$\displaystyle B^{w}$
$\displaystyle\approx\int_{\textbf{p}}\,\mathcal{E}^{nw}_{13}(\textbf{k}_{1}-\textbf{p})\mathcal{E}^{nw}_{23}(\textbf{k}_{2}+\textbf{p})\mathcal{E}^{w}_{12}(\textbf{p})+(2~{}{\rm
cycle})$
$\displaystyle=e^{\frac{1}{2}(\textbf{k}_{3}\cdot\textbf{k}_{1})\Sigma^{2}}\Big{(}\mathcal{E}^{nw,{\rm
fin}}_{32}(\textbf{k}_{3})\mathcal{E}_{12}(\textbf{k}_{1})+\mathcal{E}^{w}_{32}(\textbf{k}_{3})\mathcal{E}^{nw,{\rm
fin}}_{12}(\textbf{k}_{1})\Big{)}$ $\displaystyle\hskip
128.0374pt+\int_{\textbf{p}}\,\mathcal{E}^{nw,{\rm
fin}}_{12}(\textbf{k}_{1}-\textbf{p})\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{k}_{2}+\textbf{p})\mathcal{E}^{w}(\textbf{k}_{1},\textbf{k}_{2},\textbf{p})+2~{}{\rm
cycle}$ $\displaystyle=\mathcal{E}^{nw,{\rm
fin}}_{32}(\textbf{k}_{3})\mathcal{E}^{w}_{12}(\textbf{k}_{1})+\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{k}_{2})\mathcal{E}^{\rm w}_{13}(\textbf{k}_{1})$
$\displaystyle\hskip 128.0374pt+\int_{\textbf{p},{\rm
IR}}\,\mathcal{E}^{nw,{\rm
fin}}_{13}(\textbf{k}_{1}-\textbf{p})\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{k}_{2}+\textbf{p})\mathcal{E}^{w}_{12}(\textbf{p})+2~{}{\rm
cycle}\,,$ (4.15)

where we have again used that $\int_{\textbf{p}}\mathcal{E}^{\rm
fin}_{12}(\textbf{p})=1-\exp(\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})\,\Sigma^{2})$.
In the above we also introduced the IR-regularized mode coupling integral

$\displaystyle\int_{\textbf{p},{\rm IR}}\,\mathcal{E}^{nw,{\rm
fin}}_{13}(\textbf{k}_{1}-\textbf{p})\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{k}_{2}+\textbf{p})\mathcal{E}^{w}_{12}(\textbf{p})$ (4.16)
$\displaystyle~{}~{}\equiv\int_{\textbf{p}}\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{p})\Big{(}\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3}+\textbf{p})\mathcal{E}^{w}_{21}(\textbf{k}_{2}-\textbf{p})\Theta\left(|\textbf{k}_{3}+\textbf{p}|-p\right)-\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3})\mathcal{E}^{w}_{21}(\textbf{k}_{2})\Big{)}+(\textbf{k}_{1}\leftrightarrow\textbf{k}_{2})\,.$

In the above regularization we used the fact that we can split the integral
domain into two parts, each characterized by the poles at
$\textbf{p}\to-\textbf{k}_{2}$ and $\textbf{p}\to\textbf{k}_{1}$, such that we
can write

$\int_{\textbf{p}}\,\mathcal{E}^{nw,{\rm
fin}}_{13}(\textbf{k}_{1}-\textbf{p})\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{k}_{2}+\textbf{p})\mathcal{E}^{w}_{12}(\textbf{p})=\int_{|\textbf{k}_{1}-\textbf{p}|<|\textbf{k}_{2}+\textbf{p}|}+\int_{|\textbf{k}_{1}-\textbf{p}|>|\textbf{k}_{2}+\textbf{p}|}\,.$
(4.17)

Furthermore, each contribution can be rewritten as

$\displaystyle\int_{|\textbf{k}_{1}-\textbf{p}|<|\textbf{k}_{2}+\textbf{p}|}\mathcal{E}^{nw,{\rm
fin}}_{13}(\textbf{k}_{1}-\textbf{p})$ $\displaystyle\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{k}_{2}+\textbf{p})\mathcal{E}^{w}_{12}(\textbf{p})$
$\displaystyle=\int_{p<|\textbf{k}_{3}+\textbf{p}|}\mathcal{E}^{nw,{\rm
fin}}_{13}(\textbf{p})\mathcal{E}^{nw,{\rm
fin}}_{32}(\textbf{k}_{3}+\textbf{p})\mathcal{E}^{w}_{12}(\textbf{k}_{1}-\textbf{p})\,.$
(4.18)

At this point there are still remaining IR contributions arising from taking
$\textbf{p}\to\textbf{k}_{3}$. We can further extract these IR modes by
isolating the $p<k_{3}$ region from the integral above so that

$\displaystyle\int_{\textbf{p},{\rm IR}}\,\mathcal{E}^{nw,{\rm fin}}_{13}$
$\displaystyle(\textbf{k}_{1}-\textbf{p})\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{k}_{2}+\textbf{p})\mathcal{E}^{w}_{12}(\textbf{p})$
$\displaystyle=\int_{\textbf{p}}\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{p})\Big{(}\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3}+\textbf{p})\Theta\left(|\textbf{k}_{3}+\textbf{p}|-p\right)-\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3})\Big{)}\mathcal{E}^{\rm
w}_{21}(\textbf{k}_{2}-\textbf{p})$
$\displaystyle~{}~{}~{}~{}+\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3})\int_{\textbf{p}}\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{p})\Big{(}\mathcal{E}^{w}_{21}(\textbf{k}_{2}-\textbf{p})-\mathcal{E}^{w}_{21}(\textbf{k}_{2})\Big{)}+(\textbf{k}_{1}\leftrightarrow\textbf{k}_{2})$
$\displaystyle\approx\int_{\textbf{p}}\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{p})\Big{(}\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3}+\textbf{p})\Theta\left(|\textbf{k}_{3}+\textbf{p}|-p\right)-\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3})\Big{)}\mathcal{E}^{w}_{21}(\textbf{k}_{2}-\textbf{p})$
$\displaystyle~{}~{}~{}~{}+\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3})\left(e^{\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})\Sigma^{2}_{s}}-e^{-\frac{1}{2}k^{2}_{2}\Sigma^{2}_{s}}\right)\frac{(\textbf{k}_{1}\cdot\textbf{k}_{2})}{k_{2}^{2}}P^{w}_{\rm
L}(k_{2})+(\textbf{k}_{1}\leftrightarrow\textbf{k}_{2})\,.$ (4.19)

In the second line above, we used the approximation that the wiggles are
localized in the configuration space, corresponding to the configuration-space
saddle point, i.e.

$\displaystyle\int_{\textbf{p}}\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{p})\Big{(}$
$\displaystyle\mathcal{E}^{w}_{21}(\textbf{k}_{2}-\textbf{p})-\mathcal{E}^{w}_{21}(\textbf{k}_{2})\Big{)}$
$\displaystyle=e^{\frac{1}{2}\textbf{k}_{1,i}\textbf{k}_{2,j}A^{nw}_{ij}(\textbf{q}_{*})}\tfrac{1}{2}\textbf{k}_{1,i}\textbf{k}_{2,j}\int_{\textbf{q}}\textbf{q}\
e^{-i\textbf{k}_{2}\cdot\textbf{q}}A^{w}_{ij}(\textbf{q})\int_{\textbf{p}}\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{p})\left(e^{i\textbf{p}\cdot\textbf{q}}-1\right)$
$\displaystyle=e^{\frac{1}{2}\textbf{k}_{1,i}\textbf{k}_{2,j}A^{nw}_{ij}(\textbf{q}_{*})}\left(e^{\frac{1}{2}\textbf{k}_{2,i}\textbf{k}_{3,j}A^{nw}_{ij}(\textbf{q}_{*})}-1\right)\tfrac{1}{2}\textbf{k}_{1,i}\textbf{k}_{2,j}\int
d^{3}\textbf{q}\ e^{-i\textbf{k}_{2}\cdot\textbf{q}}\ A^{w}_{ij}(\textbf{q})$
$\displaystyle\approx-\left(e^{-\frac{1}{2}\textbf{k}_{2,i}\textbf{k}_{2,j}A^{nw}_{ij}(\textbf{q}_{*})}-e^{\frac{1}{2}\textbf{k}_{1,i}\textbf{k}_{2,j}A^{\rm
nw}_{ij}(\textbf{q}_{*})}\right)\frac{(\textbf{k}_{1}\cdot\textbf{k}_{2})}{k_{2}^{2}}P_{\rm
L}^{w}(k_{2})$
$\displaystyle\approx\left(e^{\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})\Sigma^{2}_{s}}-e^{-\frac{1}{2}k^{2}_{2}\Sigma^{2}_{s}}\right)\frac{(\textbf{k}_{1}\cdot\textbf{k}_{2})}{k_{2}^{2}}P_{L}^{w}(k_{2})\,.$
(4.20)

We see that we obtain two exponential contributions, one containing the
angular dependence of the two modes $\textbf{k}_{1}\cdot\textbf{k}_{2}$ and
the second one depending only on $k_{2}^{2}$. One might wonder if the first
piece can give rise to bispectrum configurations where the wiggles are
enhanced. However, this does not happen as these contributions are exactly
equal in magnitude but opposite in sign from the ones in the tree level part
of Equation 4.15. Thus, combining all these results and cancelling all
isolated IR contributions, we get

$\displaystyle B^{w}$ $\displaystyle=-\Big{(}\mathcal{E}^{nw,{\rm
fin}}_{32}(\textbf{k}_{3})(\textbf{k}_{1}\cdot\textbf{k}_{2})+\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{k}_{2})(\textbf{k}_{1}\cdot\textbf{k}_{3})\Big{)}e^{-\frac{1}{2}k_{1}^{2}\Sigma^{2}_{s}}P_{\rm
L}^{w}(k_{1})/k_{1}^{2}$
$\displaystyle~{}~{}~{}+\int_{\textbf{p}}\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{p})\Big{(}\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3}+\textbf{p})\Theta\left(|\textbf{k}_{3}+\textbf{p}|-p\right)-\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3})\Big{)}\mathcal{E}^{w}_{21}(\textbf{k}_{2}-\textbf{p})+(\textbf{k}_{1}\leftrightarrow\textbf{k}_{2})$
$\displaystyle~{}~{}~{}+\left({\rm 2cycle}\right)\,.$ (4.21)

We have thus obtained the wiggle bispectrum in an explicitly IR-safe form,
with the resummed long displacements acting on the linearized wiggle
contribution. The mode coupling integral above can also be rewritten in terms
of the $\textbf{k}_{1}$ and $\textbf{k}_{2}$ modes as

$\displaystyle\int_{\textbf{p}}\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{p})\Big{(}$ $\displaystyle\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3}+\textbf{p})\Theta\left(|\textbf{k}_{3}+\textbf{p}|-p\right)-\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3})\Big{)}\mathcal{E}^{w}_{21}(\textbf{k}_{2}-\textbf{p})$
(4.22) $\displaystyle=\int_{\textbf{p}}\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{k}_{2}+\textbf{p})\Big{(}\mathcal{E}^{nw,{\rm
fin}}_{13}(\textbf{k}_{1}-\textbf{p})\Theta\left(|\textbf{k}_{1}-\textbf{p}|-|\textbf{k}_{2}+\textbf{p}|\right)-\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3})\Big{)}\mathcal{E}^{\rm w}_{12}(\textbf{p})\,.$

Let us consider the possible corrections to this IR resummed result. As stated
above, we have linearized the result in $P_{\rm L}^{w}$, and one can, of
course, extend the analysis to include quadratic and higher-order
contributions. Moreover, the localized-wiggle (saddle-point) approximation,
used in Equation 4.19, gives rise to correction terms discussed in Appendix C.
On top of these, however, we can also consider the residual contributions of
the IR modes in the mode coupling integral in Equation 4.21. Isolating the
modes $p<\Lambda<k_{2,3}/2$ we can expand the integrand in the IR-safe mode
coupling integral to obtain the correction due to these modes

$\displaystyle\int_{\textbf{p}}\mathcal{E}^{nw,{\rm fin}}_{23}(\textbf{p})$
$\displaystyle\Big{(}\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3}+\textbf{p})-\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3})\Big{)}\mathcal{E}^{w}_{21}(\textbf{k}_{2}-\textbf{p})$
$\displaystyle=\sum_{n=1}^{\infty}\frac{1}{n!}\nabla_{i_{1}..i_{n}}^{n}\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3})\
\int_{\textbf{p}}\textbf{p}_{i_{1}}...\textbf{p}_{i_{n}}\mathcal{E}^{nw,{\rm
fin}}_{23}(\textbf{p})\ \mathcal{E}^{w}_{21}(\textbf{k}_{2}-\textbf{p})$
$\displaystyle=\sum_{n=1}^{\infty}\frac{1}{n!}\nabla_{i_{1}..i_{n}}^{n}\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3})\
\int_{\textbf{q}}e^{-i\textbf{k}_{2}\cdot\textbf{q}}\left(\frac{\partial^{n}\tilde{\mathcal{E}}_{23}^{nw,\rm
fin}(\textbf{q})}{\partial\textbf{q}_{i_{1}}...\partial\textbf{q}_{i_{n}}}\right)\tilde{\mathcal{E}}^{w}_{12}(\textbf{q})$
$\displaystyle\approx\sum_{n=1}^{\infty}\frac{1}{n!}\nabla_{i_{1}..i_{n}}^{n}\mathcal{E}^{nw,{\rm
fin}}_{31}(\textbf{k}_{3})\
\int_{\textbf{q}}e^{-i\textbf{k}_{2}\cdot\textbf{q}}\left(\frac{\partial^{n}\tilde{\mathcal{E}}_{23}^{nw,\rm
fin}(\textbf{q})}{\partial\textbf{q}_{i_{1}}...\partial\textbf{q}_{i_{n}}}\right)_{q=r_{s}}\tilde{\mathcal{E}}^{w}_{12}(\textbf{q})\,,$
(4.23)

where the $\nabla$ derivatives are with respect to $\textbf{k}_{3}$. Already
in the second line we see that the smoothness of the no-wiggle component leads
to a suppresssion, since the $\nabla$ derivatives multiplied by the integral
over $\textbf{p}_{i}$ lead to a supppression of $(\Lambda/k_{3})^{n}$ at each
order in the Taylor series. However, the final correction is even further
suppressed than this estimate, since by re-writing the correction as the
Fourier transform in the final line we can see that the correction can again
be evaluated using a a saddle-point approximation owing to the localization of
the wiggles in configuration space. In Fourier-space this corresponds to the
wiggle at $\textbf{k}_{2}-\textbf{p}$ losing coherence with $\textbf{k}_{2}$
for $p$ wider than the wavelength of the oscillation. This leads to a
correction of the form in Equation 4.13 due to modes at the BAO scale, which
is indeed not surprising since the two derivations use the same IR domain in
different spaces. In the next subsection, where we explicitly show and compare
various IR resummation schemes, we also explicitly evaluate the residual IR
and UV corrections due to the mode coupling integral in Equation 4.21, showing
that the latter are numerically highly suppressed as argued here.

###  4.3 Numerical Comparison of Resummation Schemes

Figure 8: As in Figure 3 but with the BAO wiggle component isolated and
normalized by the no-wiggle tree-level prediction. Significant deviations in
the wiggle phase and amplitude are evident in the $k=0.20\,h\,{\rm Mpc}^{-1}$
triangles. Figure 9: As in Figure 4, BAO in the squeezed (left) and
equilateral (right) bispectrum configurations. The full Zeldovich calculation
(blue) approaches the 1-loop (orange) and tree-level (dotted) predictions
towards low $k$ but exhibit dramatically different oscillatory behavior
towards smaller scales. The Zeldovich prediction is well approximated by an
exponential damping of BAO wiggles in the linear power spectrum, particularly
in the squeezed configuration.

We can now numerically investigate the BAO signal in the bispectrum as
predicted by various perturbation theory schemes at tree-level and 1-loop,
pre- and post-IR resummation, using the full Zeldovich calculation as a
benchmark where all the effects of long-wavelength linear displacements are
manifestly taken into account. In order to do so it is convenient to define
the wiggle and no-wiggle components of bispectra as

$B_{w}=B[P_{\rm L}=P_{nw}+P_{w}]-B[P_{nw}],\ B_{nw}=B[P_{nw}]\,,$ (4.24)

where the bispectra on the right-hand side are evaluated assuming the linear
power spectra in the square brackets.

Figures 8 and 9 show the BAO components of the tree level, 1-loop and full
Zeldovich bispectra, so defined, for bispectra of various side lengths and
geometries. The 1-loop and full Zeldovich power spectra begin to behave quite
differently at large wavenumbers, with the latter remaining mostly in phase
with the tree-level predictions while the former veers off significantly and
even changes sign. This is as expected because the damping effect by long
displacements is controlled by a rather larger parameter $k^{2}\Sigma^{2}_{s}$
and its effects become rather non-perturbative at the higher wavenumbers show
in these plots. The enhancement of this large parameter at 1-loop in
perturbation theory is described explicitly at the end of Appendix B.

In comparison to the Zeldovich 1-loop prediction, the dashed lines in Figure 9
show the prediction of the resummed tree-level bispectrum using the RWiggle
scheme. Evidently, this scheme captures the damping of the BAO component quite
well, with particularly excellent agreement in the case of the squeezed
configuration, though it slightly underpredicts the damping in the equilateral
triangle. The agreement in the former case can be understood by noting that in
the squeezed configuration the contribution due to the integral piece is quite
small, so that the wiggles are essentially due to the $\mathcal{E}\mathcal{E}$
terms for which RWiggle works almost perfectly.

Figure 10: The 1-loop bispectrum resummed using the RWiggle scheme (orange),
compared to the full Zeldovich bispectrum (black dashed) and the resummed
tree-level bispectrum (blue). All three are in reasonably good agreement in
squeeed triangles (left) other than a slight phase shift not captured by the
tree-level spectrum, but the tree-level RWiggle visibly underestimates the
damping of the equilateral configuration even above $0.05\,h\,{\rm Mpc}^{-1}$.

We can also compare the prediction of the Zeldovich calculation to the RWiggle
at 1-loop order. In this case we adopt the scheme from ref. [22] (see also §5)

$B^{\tt RWiggle}_{\rm 1-loop}=B_{\rm
tree}\big{[}P_{nw}+(1+\frac{1}{2}k^{2}\Sigma_{s}^{2})\
e^{-\frac{1}{2}k^{2}\Sigma_{s}^{2}}P_{w}\big{]}+B_{\rm
loop}[P_{nw}+e^{-\frac{1}{2}k^{2}\Sigma_{s}^{2}}P_{w}]\,,$ (4.25)

where the brackets indicate the bispectrum evaluated substituting the argument
for the linear power spectrum. This comparison is shown in Figure 10: we see
that the 1-loop corrections to RWiggle almost entirely corrects for the
discrepancy with the full Zeldovich calculation, including the small phase-
shift visible in the squeezed configuration and the underpredicted damping in
the equilateral one. This was not entirely unexpected; as shown in §4.2
corrections to RWiggle are due to parametrically small density modes at the
BAO scale, whose contributions are expected to be perturbative.

Figure 11: Contributions to the mode coupling integral when resumming the
effects of linear displacements on the wiggle component in Equation 4.21.
Contributions are split into IR and UV regions longer and shorter than the
sides of the equilateral triangle, with the effect on the resummed wiggles
almost entirely due to the IR modes.

To further check this intuition, Figure 11 shows the correction to RWiggle
derived in Equation 4.21 split into contributions from large and small
wavenumbers. Here we have evaluated the $\mathcal{E}$ in the mode-coupling
integrals at tree-level for simplicity, since this captures the leading 1-loop
corrections to the tree-level wiggles. The full mode-coupling integral is
shown as a black-dashed line and is evidently composed almost entirely of the
IR contribution (blue). The UV portion (orange), while nonzero, is close to an
order of magnitude smaller than the IR one and out-of-phase with the tree-
level oscillations. This numerically validates the argument in §4.2.1 and
4.2.2 that only specific IR modes close to the BAO scale contribute correcting
the saddle-point approximation. Both contributions are in addition very small
($\approx$ 20x) compared to the tree-level prediction, showing that RWiggle
captures the dominant effect of the linear displacements.

###  4.4 Further IR-Resummation Schemes

In addition to the direct Zeldovich calculation and RWiggle a few other
schemes for IR resummation exist in the literature, particularly as applied to
the power spectrum. The goal of this subsection is to briefly sketch their
relation to our calculations and how they may be extended to the case of the
bispectrum.

Ref. [17] constructed an IR resummation scheme of the two-point function with
the goal of IR resummation operating on the usual EPT loop terms. This ensures
that the IR resummed theory has approximately the same UV properties
(counterterm structure and values) as the canonical, un-resummed, EPT. We
follow the presentation from Appendix B of ref. [19], which re-derived this
result starting from the canonical Zeldovich power spectrum. By defining
$A_{ij}(\textbf{q})=\left\langle\Delta_{i}\Delta_{j}\right\rangle=\Sigma^{2}\delta^{\rm
K}_{ij}-2\eta_{ij}(\textbf{q})$, we can write

$\displaystyle\mathcal{E}^{\rm fin}_{12}(\textbf{p})$
$\displaystyle=\int_{\textbf{q}}\,e^{i\textbf{p}\cdot\textbf{q}}\left(e^{\frac{1}{2}k_{1i}k_{2j}A_{ij}(\textbf{q})}-e^{\frac{1}{2}k_{1i}k_{2j}A_{ij}(\infty)}\right)$
$\displaystyle=\int_{\textbf{q}}\,e^{i\textbf{p}\cdot\textbf{q}}\,e^{\frac{1}{2}k_{1i}k_{2j}A_{ij}(\textbf{q})}\left(1-e^{k_{1i}k_{2j}\eta_{ij}}\right)$
$\displaystyle=-\sum_{n=1}^{\infty}\frac{1}{n!}\int_{\textbf{q}}\,e^{i\textbf{p}\cdot\textbf{q}}\,K_{12}(\textbf{q})\left(k_{1i}k_{2j}\eta_{ij}\right)^{n}\,,$
(4.26)

where we introduce the Gaussian kernel
$K_{12}(\textbf{q})=\exp\left\\{\frac{1}{2}k_{1i}k_{2j}A_{ij}(\textbf{q})\right\\}$.
On the other hand, $\mathcal{E}^{\rm fin}$ can also be expanded in the
conventional EPT fashion in powers of linear power spectra. This leads to the
form

$\displaystyle\mathcal{E}^{\rm fin}_{12}(\textbf{p})$
$\displaystyle=\sum_{n=0}^{\infty}\mathcal{E}_{12}^{(n)-{\rm
loop}}(\textbf{p})\,,$ (4.27)

where for $n=0$ we have the linear theory results, while the higher $n$ give
us higher loop results

$\mathcal{E}_{12}^{(n)-{\rm
loop}}(\textbf{p})=\frac{1}{(n+1)!}\sum_{m=1}^{n+1}(-1)^{m}\binom{n+1}{m}\left\\{\left(\frac{1}{2}\left(\textbf{k}_{1}\cdot\textbf{k}_{2}\right)\Sigma^{2}\right)^{n-m+1}\big{(}k_{1i}k_{2j}\eta_{ij}\big{)}^{m}\right\\}_{\rm
FT}\,.$ (4.28)

Up to one-loop order, this gives us

$\displaystyle\mathcal{E}_{12}^{\rm lin}$
$\displaystyle=-\Big{\\{}k_{1i}k_{2j}\eta_{ij}\Big{\\}}_{\rm FT}\,,$
$\displaystyle\mathcal{E}_{12}^{1-{\rm loop}}$
$\displaystyle=\frac{1}{2}\Big{\\{}\left(k_{1i}k_{2j}\eta_{ij}\right)^{2}\Big{\\}}_{\rm
FT}-\frac{1}{2}\big{(}\textbf{k}_{1}\cdot\textbf{k}_{2}\Sigma^{2}\big{)}\Big{\\{}k_{1i}k_{2j}\eta_{ij}\Big{\\}}_{\rm
FT}\,.$ (4.29)

where in the one-loop term we recognise the familiar $P_{22}+P_{13}$
structure. Also, note that these are different from the usual power spectrum
loops (in the Zeldovich approximation) since they depend on three wave vectors
$\textbf{k}_{1}$, $\textbf{k}_{2}$ and p. However, aligning these as
$\textbf{k}_{1}=-\textbf{k}_{2}=\textbf{p}$ reproduces the power spectrum
results, as was noted below Equation 3.1.

Utilising this EPT style expansion, it can be shown that Equation 4.26 can be
rewritten as

$\displaystyle-\sum_{n=1}^{\infty}\frac{1}{n!}\left(k_{1i}k_{2j}\eta_{ij}\right)^{n}=\sum_{n=0}^{\infty}\sum_{m=0}^{n}\left[K^{-1}_{12}\right]_{n-m}\left\\{\mathcal{E}_{12}^{(m)-{\rm
loop}}\right\\}_{\rm FT}$ (4.30)

where we have introduced the inverse of the $K$ kernel and its expansion

$\displaystyle K^{-1}_{12}(\textbf{q})$
$\displaystyle=\exp\left\\{-\frac{1}{2}k_{1i}k_{2j}A_{ij}(\textbf{q})\right\\}\,,$
$\displaystyle\left[K^{-1}_{12}\right]_{n}(\textbf{q})$
$\displaystyle=\frac{1}{n!}\left(k_{1i}k_{2j}\eta_{ij}-\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})\Sigma^{2}\right)^{n}\,.$
(4.31)

Combining these results, we can write

$\displaystyle\mathcal{E}^{\rm fin}_{12}(\textbf{p})$
$\displaystyle=\sum_{n=0}^{\infty}\sum_{m=0}^{n}\int_{\textbf{p}^{\prime}}\mathcal{M}^{(n-m)}_{12}\left({\textbf{p}-\textbf{p}^{\prime}}\right)\mathcal{E}_{12}^{(m)-{\rm
loop}}(\textbf{p}^{\prime})\,,$ (4.32)

where the momentum space kernel $\mathcal{M}$ contains only effects of the
long displacements

$\mathcal{M}^{(n)}_{12}(\textbf{k})=\int_{\textbf{q}}\,e^{i\textbf{k}\cdot\textbf{q}}\,K_{12}(\textbf{q})\left[K^{-1}_{12}\right]_{n}(\textbf{q})\,.$
(4.33)

Note that, strictly speaking, even in the ordinary power spectrum case, the
integral above is not a simple convolution of $\mathcal{M}^{(n)}$ kernels with
the EPT loops, given that there are always external modes present in the loops
of $\mathcal{E}_{12}$. Nonetheless, it is easy to see that this result reduces
to EPT in the case that $\mathcal{M}^{(0)}$ is a delta function and the rest
of $\mathcal{M}^{(m)}$ vanish. Moreover, ref. [17] shows that a good
approximation of this result can also be achieved when the external
$\textbf{k}_{1,2,3}$ and internal p modes in $\mathcal{E}_{12}$ are
identified, and thus the latter term indeed corresponds to the EPT result.
This analysis can be extended to the bispectrum, for example, using Equation
3.4, or by generalizing $K_{12}(\textbf{q})$ to the bispectrum exponent in
Equation 2.5.

Another related method to resum Zeldovich displacements was described in ref.
[45] using the so-called shifted operators

$\tilde{O}(k)=\int d^{3}\textbf{q}\
e^{-i\textbf{k}\cdot(\textbf{q}+\Psi(\textbf{q}))}\ O(\textbf{q}).$ (4.34)

In this scheme the matter clustering ($F(\textbf{q})=1$) term in LPT (Eqn.
2.3) is re-written in terms of bias operators, such that for example the
Zeldovich matter density is given by [45]

$\displaystyle\delta_{\rm Zel}(\textbf{k})$ $\displaystyle=\int
d^{3}\textbf{q}\
e^{-i\textbf{k}\cdot(\textbf{q}+\Psi)}\Big{\\{}\delta(\textbf{q})+\frac{1}{2}\big{(}s^{2}(\textbf{q})-\frac{2}{3}\delta^{2}(\textbf{q})\Big{)}+\ldots\big{\\}}\equiv\tilde{\delta}(k)+\frac{1}{2}\tilde{\mathcal{G}}_{2}(\textbf{k})+\ldots$
(4.35)

up to quadratic order. The shifted operator basis is an intermediate between
Lagrangian and Eulerian perturbation theory, almost keeping the form of the
latter but with bulk displacements kept to arbitrary order in the exponent
rather than expanded perturbatively. A useful fact is that the nonlinear
matter field itself can be written in terms of the shifted operators above
and, conveniently, the bispectrum for the shifted operators up to the order
shown above can be easily computed using generalizations of $\mathcal{E}$, as
we describe in more detail in Appendix E.

##  5 IR Resummation for N-point Functions

We can also extend our analysis to the IR resummation of n-point functions
beyond the bispectrum. To do so, we need to consider the generalised version
of the cumulant in Equation 2.5

$\left\langle
e^{-i\sum_{n=1}^{N}k_{n,i}\Psi^{(1)}_{i}(\textbf{q}_{i})}\right\rangle_{\textbf{k}_{\rm
tot}=0}=\exp\left\\{\frac{1}{2}\sum_{n=1}^{N}\sum_{m>n}k_{n,i}k_{m,j}A_{ij}(\textbf{q}_{nm})\right\\}\,.$
(5.1)

Here, the sum runs over all distinct pairs $(n,m)$, and we have used the
identity Equation 2.4. In the above we have defined
$\textbf{q}_{nm}=\textbf{q}_{n}-\textbf{q}_{m}$. As before it is useful to
shift to a coordinate system
$\textbf{q}_{i<N}\rightarrow\textbf{q}_{i<N}-\textbf{Q}$, with
$\textbf{Q}=\textbf{q}_{N}$, such that the expectation values can be written
independently of Q by translation invariance. The Fourier-space N-point
function is then given by

$C_{N}(\textbf{k}_{1},...,\textbf{k}_{N-1})=\left(\prod_{a=1}^{N-1}\int_{\textbf{q}_{a}}e^{-i\textbf{k}_{a}\cdot\textbf{q}_{a}}\right)\
\exp\left\\{\frac{1}{2}\sum_{n\neq
m}k_{n,i}k_{m,j}A_{ij}(\textbf{q}_{nm})\right\\}\
f_{N}(\textbf{q}_{1},...,\textbf{q}_{N-1})\,,$ (5.2)

where we have used “$n\neq m$” as the shorthand for the sum over distinct
pairs. Here, $f_{N}$ captures terms due to galaxy bias and higher-order
displacements and is equal to unity for matter in the Zeldovich approximation.
We can equivalently express the above in terms of $\mathcal{E}$ as

$C(\textbf{k}_{1},...,\textbf{k}_{N-1})=\prod_{n\neq m}\int_{\textbf{p}_{nm}}\
\mathcal{E}(\textbf{k}_{n},\textbf{k}_{m},\textbf{p}_{nm})\prod_{a=1}^{N-1}(2\pi)^{3}\delta_{D}\big{(}\textbf{k}_{a}-\sum_{n>a}p_{an}+\sum_{n<a}p_{na}\big{)}\,,$
(5.3)

where the momentum conservation comes from integrating over the vertices
$\textbf{q}_{a}$, and we remind the reader that there is no sum over the Nth
vertex. Note that the bispectrum and power spectrum are special cases in that
there is only one fewer vertex (2 and 1) than pairwise separations
$\textbf{q}_{nm}$ (3 and 2), leading to a reduction in the number of integrals
in Fourier space. With the trispectrum, for example, one has to do an equal
number of integrals in each space. For higher N-point functions there is a
rapid scaling in the number of internal momenta $\textbf{p}_{nm}$ with $N$,
though the Fourier space integral still has the advantage of being smooth
compared to the Fourier transforms required in configuration space.

Separating the contributions to the above integral via a wiggle-no wiggle
split it is easy to work out the effect of long-wavelenght displacements from
on BAO wiggles in the N-point functions. In particular, expanding the wiggles
to linear order we will have contributions that are, schematically,

$\displaystyle
C_{w}\sim\left(\prod_{m=2}^{N-1}\int_{\textbf{q}_{m}}e^{-i\textbf{k}_{m}\cdot\textbf{q}_{m}}\right)$
$\displaystyle\exp\left\\{\frac{1}{2}\sum_{n,m>1}^{n\neq
m}k_{n,i}k_{m,j}A_{ij}(\textbf{q}_{nm})\right\\}$
$\displaystyle\int_{\textbf{q}_{1}}e^{-i\textbf{k}_{1}\cdot\textbf{q}_{1}}\exp\left\\{\frac{1}{2}\sum_{m\neq
1}k_{1,i}k_{m,j}A_{ij}(\textbf{q}_{1m})\right\\}\xi_{w}(\textbf{q}_{1})\,.$
(5.4)

Here we have chosen a BAO features in the $\textbf{q}_{1}=\textbf{q}_{1N}$ leg
of the N-point function $f_{N}$ in Lagrangian space but, since the final
expression is invariant under permutations of the vertices we can do so
without loss of generality. The top line of Equation 5.4 is smooth in
$\textbf{q}_{1}$ so we can ignore it. If we make the same approximation
$\textbf{q}_{1m}\rightarrow\textbf{q}$ as in the bispectrum in the previous
section then the exponent in the second line becomes at the saddle point

$\frac{1}{2}k_{1,i}\left(\sum_{m\neq
n}k_{m,j}\right)A_{ij}(\textbf{q}_{s})=-\frac{1}{2}k_{1}^{2}\Sigma^{2}_{s}\,,$
(5.5)

by momentum conservation. This implies that the effects of long-wavelength
displacements on wiggles in N-point functions can be well-described, at
leading order in the wiggles, by the substitution

$P_{w}(k)\rightarrow e^{-\frac{1}{2}k^{2}\Sigma^{2}_{s}}P_{w}(k)\,.$ (5.6)

This prescription is a consequence of translation invariance and momentum
conservation; it is interesting to observe that this condition is also
satisfied for N-point functions with collapsed legs, such that the loops can
be computed with this substitution as well to approximate the effects of IR
displacements.

##  6 Conclusions

Figure 12: Squeezed (left) and equilateral (right) Zeldovich bispectra when
the linear power spectrum has logarithmic wiggles on top of our fiducial
$\Lambda$CDM one. The dotted and dashed lines show the tree level and RWiggle
predictions—the latter captures some of the damping due to the linear
displacements, as in the case of the BAO, but with still clearly visible
differences in the amplitude.

The presence of sharp, localized features in the linear correlation
function—or, oscillatory “wiggles” in the power spectrum—lead large-scale
structure observables to be sensitive to the bulk motion of matter on the
large-scales at which these features reside. In $\Lambda$CDM universes close
to our own, the contributions of these infrared (IR) modes to galaxy N-point
statistics can be numerically quite large, such that they need to be resummed
in order for perturbative calculations to be well-behaved. Indeed, this is the
case both for features that are the bread and butter of modern galaxy surveys,
i.e. the baryon acoustic oscillations standard ruler, and proposed targets of
future experiments, e.g. features in the primordial power spectrum that probe
the particular shape of the inflaton potential. A natural way to formulate
this procedure, known as IR resummation, is within the Lagrangian picture of
fluid mechanics, wherein fluids are described via the displacements of
individual fluid elements which give rise to these bulk flows on large scales.
Within Lagrangian perturbation theory (LPT), the effect of the bulk
displacements can be understood as the smearing of the correlation function
due to pairwise displacements of fluid elements separated by the scale of the
feature (e.g. the BAO radius $r_{d}$).

In this paper, we build on previous works on the galaxy 2-point function to
study the galaxy 3-point function, or bispectrum, in Lagrangian perturbation
theory. In the first part, we extend past techniques used to efficiently
compute the resummed LPT power spectrum to the resummed bispectrum, focusing
specifically on the case of the bispectrum within the Zeldovich approximation
which contains all the long-wavelength displacements we wish to include. Our
method involves computing the Fourier transform of the variance of the
pairwise displacement, projected along arbitrary wavenumbers, and convolving
the resulting kernels in Fourier space. These kernels contain an IR pole when
the wavenumber $\textbf{p}=\textbf{0}$ whose effects cancel in the full
calculation, and in order to properly treat their convolution numerically we
split the integral into Voronoi regions closest to each pole and re-organize
the expression for the integrand in each such that IR divergences explicitly
vanish. This is a generalization of the so-called IR-safe integrands in
Eulerian perturbation theory, with the effects of bulk displacements now
carried to arbitrary order. We make the code to perform these calculations,
triceratops, publicly available.

In the second part of the paper, we focus on computing the effect of the
resummed displacements on features in the linear 2-point function, taking the
BAO peak as a particular example. In the 2-point function, the presence of
such a feature allows for a simple saddle-point approximation of these
displacements at the peak, leading to a simple Gaussian damping of the power
spectrum wiggles in Fourier space. The situation in the 3-point function is
not so simple: here, requiring that two points in a triangle be separated by
the BAO radius nonetheless leaves substantial freedom for the remaining point
and two legs of the triangle to move around. However, we can nonetheless
perform a saddle point approximation of the Zeldovich displacements in the
bispectrum—we show that this leads to the familiar form where the effect of
the displacements is a Gaussian smoothing of the linear power spectrum in each
leg $k$, which we call RWiggle, with corrections from the extra degrees of
freedom due to perturbatively small density modes on the BAO scale, suppressed
by inverse powers of $kr_{d}$. This conclusion holds for arbitrary $N$-point
functions due to momentum conservation. We numerically compare the wiggle
component computed in this approximation to that in the full Lagrangian
calculation, finding close agreement, especially at 1-loop order after the IR
resummation of the wiggle component. We also make contact with other IR
resummation schemes in the literature, deriving an extension of the scheme
developed for the power spectrum in ref. [17] from the Lagrangian calculation
and showing how the bispctrum in the shifted-operator formalism developed for
field-level analyses in ref. [45] can be computed uinsg the techniques in this
paper.

We have thus established that resummed 1-loop perturbation theory very-well
captures the effect of long-wavelength displacements on the BAO signal in
cosmological N-point functions. However, it is worth noting that the direct
Lagrangian treatment has distinct advantages over the RWiggle prescription, in
that the wiggle (feature) component does not have to be distinctly treated.
One advantage of this is that in the Lagrangian method one does not need to
perform the wiggle/no-wiggle split, which is somewhat arbitrary, since there
is no clear criterion that the smooth and wiggle components must satisfy, and
indeed likely need to be adapted individually for each inflationary model
treated. Moreover, features in the linear power spectrum beyond the BAO, e.g.
inflationary signatures in the primordial spectrum, may have nonlinear
dispersion relations not localized to any particular scale in configuration
space—such a case is handled out-of-the-box in LPT but needs to be handled on
a case-by-case basis in RWiggle, since each frequency in Fourier space
corresponds to a distinct IR scale. This makes the Lagrangian treatment
indispensable for future spectroscopic searches for exotic physics.

As a simple example let us consider oscillations in the bispectrum due to
logarithmic oscillations in the primordial power spectrum

$P_{\rm L}(k)=P_{\rm L}^{\Lambda\rm
CDM}(k)\left(1+A\sin(\omega\log\frac{k}{k_{\ast}})\
e^{-(kr_{d})^{2}/2}\right)\,,$ (6.1)

which can appear in certain models of inflation [46, 28]. Previous works [24,
47] have shown that, in the limit of large log frequency, these logarithmic
oscillations are approximately damped in a wavenumber-dependent fashion555In
LPT this corresponds to the limit in which the oscillations can be locally
well-approximated by a Taylor expansion  $\phi(k)=\omega f(k)=\omega
f(k_{0})+\omega f^{\prime}(k_{0})(k-k_{0})+...$ (6.2) Around $k_{0}$ the
oscillations correspond to a physical scale of $q(k)=\omega
f^{\prime}(k_{0})=\omega/k_{0}$ for logarithmic oscillations, though we note
that this analysis in principle applies to any nonlinear dispersions.  by
$\Sigma^{2}(k)=X(q(k))+Y(q(k))$, where $q(k)=\omega/k$ while also exactly
computed in the full LPT calculation of the power spectrum. Using the results
in this paper, we can straightforwardly extend these calculations to the
bispectrum, as we have done in Figure 12. Here we have followed ref. [47] and
chosen $A=0.05$, $\omega=10$, $k_{\ast}=0.05\,h\,{\rm Mpc}^{-1}$ and
$r_{d}=2.5\ h^{-1}\text{Mpc}$. As in the case of the BAO the RWiggle
prescription at tree level qualitatively describes the damping of the wiggles,
though unlike there deviations can be seen even in the squeezed configuration.
For current and future surveys, the galaxy bispectrum carries substantial
promise as a probe of primordial features on a comparable footing with the
power spectrum [48]; however, since searches for primordial features need to
recover or constrain not just the frequency but the amplitude of power
spectrum oscillations it will be important to consider these differences in IR
resummation schemes when using the bispectrum as a probe of inflationary
features. For that, a full treatment of the bispectrum in Lagrangian
perturbation theory beyond what was developed here, or Eulerian calculations
at 1-loop order, will be essential.

Indeed, while we have focused mostly on the matter power spectrum in the
Zeldovich approximation in this work, our calculations have immediate
implications for galaxy n-point functions in real and redshift space.
Extending to redshift-space is a simple matter of coordinate transformations.
The linear displacements are converted to their redshfit-space counterparts
via a matrix transformation boosting the line-of-sight ($\hat{n}$) components,
$\Psi^{(1)}_{s,i}=R_{ij}\Psi^{(1)}_{j}$, where
$R_{ij}=\delta_{ij}+f(z)\hat{n}_{i}\hat{n}_{j}$ where $f(z)$ is the linear
growth rate. This is equivalent to a remapping of the wavenumbers, since
displacements always appear in the combination
$k_{i}\Psi_{s,i}=(R.\textbf{k})_{i}\Psi_{i}$. In RWiggle this is equivalent to
setting $k^{2}\Sigma_{s}^{2}\rightarrow k^{2}(1+f(2+f)\mu^{2})\Sigma_{s}^{2}$,
where $\mu=\hat{n}\cdot\hat{k}$. The extension to biased tracers can be
similarly accomplished by letting the bias functionals in Equation 2.3 depend
on local observables like the tidal tensor. Since this extension is not
expected to change any conclusions about the role of infrared displacements,
which instead live in the exponent of Equation 2.3, we will leave detailed
calculations to future work. The interested reader is, however, directed to
Appendix E for expressions in real-space up to quadratic order in the bias
expansion.

## Acknowledgments

SC would like to thank the CECC in Taiwan, which generously hosted him in the
early stages of this project. SC acknowledges the support of the National
Science Foundation at the Institute for Advanced Study. Z.V. acknowledges the
support of the Kavli Foundation. MW is supported by NASA and the DOE.

This research has made use of NASA’s Astrophysics Data System and the arXiv
preprint server. This research is supported by the Director, Office of
Science, Office of High Energy Physics of the U.S. Department of Energy under
Contract No. DE-AC02-05CH11231, and by the National Energy Research Scientific
Computing Center, a DOE Office of Science User Facility under the same
contract.

##  Appendix A Numerical Evaluation of $\mathcal{E}$

In §3 we introduced a generalized version of the Zeldovich power spectrum, the
function $\mathcal{E}(\textbf{k}_{1},\textbf{k}_{2},\textbf{p})$, through
which the effects of large linear displacements can be captured. The goal in
this appendix is to develop numerical techniques to efficiently evaluate this
function. In §A.1 we outline one way to perform this calculation by performing
an angular decomposition of the integral into parts that can be rapidly
computed by Hankel transform. §A.2 shows how the first few terms of this
decomposition reduce to linear theory at leading order. Finally, §A.3 presents
an alternative method to perform this calculation whose form will also be
useful in computing the angular structure of the BAO resummation in §D.

###  A.1 Angular Decomposition and Fast Evaluation Using Hankel Transforms

Let us now examine $\mathcal{E}$ in more detail. By symmetry, we have the
angular dependence

$\mathcal{E}(\textbf{k}_{1},\textbf{k}_{2},\textbf{p})=\mathcal{E}\left(k_{1},k_{2},p,\hat{k}_{1}\cdot\hat{k}_{2},\hat{k}_{1}\cdot\hat{p},\hat{k}_{2}\cdot\hat{p}\right)=\int
d^{3}q\,e^{-i\textbf{p}\cdot\textbf{q}}e^{+\frac{1}{2}\textbf{k}_{1,i}\textbf{k}_{2,j}A_{ij}(\textbf{q})}\,,$
(A.1)

where $A_{ij}=X\delta^{K}_{ij}+Y\hat{q}_{i}\hat{q}_{j}$ [14]. This gives us

$\frac{1}{2}\textbf{k}_{1,i}\textbf{k}_{2,j}A_{ij}=\frac{1}{2}k_{1}k_{2}\left[\hat{k}_{1}\cdot\hat{k}_{2}X(q)+(\hat{k}_{1}\cdot\hat{q})(\hat{k}_{2}\cdot\hat{q})Y(q)\right]\,.$
(A.2)

and thus

$\displaystyle\mathcal{E}$
$\displaystyle\left(k_{1},k_{2},p,\hat{k}_{1}\cdot\hat{k}_{2},\hat{k}_{1}\cdot\hat{p},\hat{k}_{2}\cdot\hat{p}\right)$
$\displaystyle\hskip 71.13188pt=\int
q^{2}dq\,e^{\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})\left[X(q)+Y(q)\right]}\mathcal{I}\left(qp,\tfrac{1}{2}k_{1}k_{2}Y(q),\hat{k}_{1}\cdot\hat{k}_{2},\hat{k}_{1}\cdot\hat{p},\hat{k}_{2}\cdot\hat{p}\right)\,.$
(A.3)

We are thus faced with solving integrals of type

$\mathcal{I}\left(\alpha,\beta,\hat{k}_{1}\cdot\hat{k}_{2},\hat{k}_{1}\cdot\hat{p},\hat{k}_{2}\cdot\hat{p}\right)=\int
d\Omega_{\hat{q}}\,e^{-i\alpha\hat{p}\cdot\hat{q}+\beta\hat{k}_{1,i}\hat{k}_{2,j}\left(\hat{q}_{i}\hat{q}_{j}-\delta^{K}_{ij}\right)}=\sum_{n=0}^{\infty}\frac{\beta^{n}}{n!}I_{n}\,.$
(A.4)

The key part to compute is the integral

$I_{n}=\int
d\Omega_{\hat{q}}\,\left[\hat{k}_{1,i}\hat{k}_{2,j}\left(\hat{q}_{i}\hat{q}_{j}-\delta^{\rm
K}_{ij}\right)\right]^{n}e^{i\alpha\hat{p}\cdot\hat{q}}\,.$ (A.5)

We can be a bit more general by replacing $\hat{p}$ with p. We thus have

$I_{n}=\int
d\Omega_{\hat{q}}\,\left[\hat{k}_{1,i}\hat{k}_{2,j}\left(\hat{q}_{i}\hat{q}_{j}-\delta^{\rm
K}_{ij}\right)\right]^{n}e^{i\alpha\textbf{p}\cdot\hat{q}}\,.$ (A.6)

Notice that these integrals form the recursion

$I_{n+1}=-\hat{k}_{1,i}\hat{k}_{2,j}\left(\frac{1}{\alpha^{2}}\frac{\partial^{2}}{\partial
p_{i}\partial p_{j}}+\delta^{\rm K}_{ij}\right)I_{n}\,.$ (A.7)

For $n=0$ we have

$\displaystyle I_{0}=4\pi j_{0}(\alpha p)\,,$ (A.8)

and thus

$I_{n}=4\pi\bigg{[}-\hat{k}_{1,i}\hat{k}_{2,j}\left(\frac{1}{\alpha^{2}}\frac{\partial^{2}}{\partial
p_{i}\partial p_{j}}+\delta^{\rm K}_{ij}\right)\bigg{]}^{n}j_{0}(\alpha p)\,.$
(A.9)

We note that in case when $\hat{k}_{1}=-\hat{k}_{2}=\hat{p}$ we have

$I_{n}=4\pi\Big{(}1+\partial_{\alpha p}^{2}\Big{)}^{n}j_{0}(\alpha
p)\,=4\pi~{}n!\left(\frac{2}{\alpha p}\right)^{n}j_{n}(\alpha p)\,,$ (A.10)

recovering the results for the power spectrum in ref. [15].

In order to find a closed expression of the solution we try the ansatz

$I_{n}=(-1)^{n}4\pi\left(F_{n}(a,b,c,\alpha p)\frac{j_{n}(\alpha p)}{(\alpha
p)^{n}}+(\alpha p)^{2}G_{n}(a,b,c,\alpha p)\frac{j_{n+1}(\alpha p)}{(\alpha
p)^{n+1}}\right)\,,$ (A.11)

with $a=\hat{k}_{1}\cdot\hat{k}_{2}$, $b=\hat{k}_{2}\cdot\hat{p}$ and
$c=\hat{k}_{1}\cdot\hat{p}$. This gives us the following recursive solutions

$\displaystyle F_{n+1}$
$\displaystyle=\hat{k}_{1,i}\hat{k}_{2,j}\Big{(}\left(2G_{n,j}-F_{n,j}\right)p_{i}+(i\leftrightarrow
j)+(2n+3)\tfrac{F_{n,ij}}{\alpha^{2}}+(\alpha
p)^{2}\tfrac{G_{n,ij}}{\alpha^{2}}\Big{)}$ (A.12) $\displaystyle\hskip
180.67517pt+2a(n+1)F_{n}+\left(2a-(\alpha p)^{2}(bc-a)\right)G_{n}\,,$
$\displaystyle G_{n+1}$
$\displaystyle=\hat{k}_{1,i}\hat{k}_{2,j}\Big{(}-\tfrac{F_{n,ij}}{\alpha^{2}}-G_{n,j}p_{i}+(i\leftrightarrow
j)\Big{)}+(bc-a)F_{n}+\big{(}(2n+1)bc-a\big{)}G_{n}\,,$

where ${}_{,i}\equiv\partial/\partial p_{i}$. Given that $F_{0}=1$ and
$G_{0}=0$ we can run the recursion and get for the first few orders

$\displaystyle F_{1}$ $\displaystyle=2a\,,$ (A.13) $\displaystyle G_{1}$
$\displaystyle=bc-a\,,$

and

$\displaystyle F_{2}$ $\displaystyle=1+7a^{2}-(\alpha p)^{2}(bc-a)^{2}\,,$
(A.14) $\displaystyle G_{2}$
$\displaystyle=-a^{2}-b^{2}-c^{2}-4abc+7b^{2}c^{2}\,.$

This motivates the further ansatz

$\displaystyle F_{n}(a,b,c,z)$ $\displaystyle=\sum_{m=0}^{\lfloor
n/2\rfloor}(-1)^{m}f_{n,m}(a,b,c)z^{2m}\,,$ (A.15) $\displaystyle
G_{n}(a,b,c,z)$
$\displaystyle=\sum_{m=0}^{\lfloor(n-1)/2\rfloor}(-1)^{m}g_{n,m}(a,b,c)z^{2m}\,.$

Introducing operators

$\displaystyle\hat{O}_{1}(m)$
$\displaystyle=\left(a-bc\right)\partial_{b}+\left(1-c^{2}\right)\partial_{c}+2mc\,,$
(A.16) $\displaystyle\hat{O}_{2}(m)$
$\displaystyle=\left(1-b^{2}\right)\partial_{b}+\left(a-bc\right)\partial_{c}+2mb\,,$
$\displaystyle\hat{O}_{12}(m)$
$\displaystyle=\left(a-bc\right)\left(1-b^{2}\right)\partial^{2}_{b}+\left(a-bc\right)\left(1-c^{2}\right)\partial^{2}_{c}+\left(\left(1-c^{2}\right)\left(1-b^{2}\right)+\left(a-bc\right)^{2}\right)\partial_{b}\partial_{c}$
$\displaystyle~{}~{}+\left(2m(c+ab-2b^{2}c)+3b^{2}c-2ab-c\right)\partial_{b}+\left(2m(b+ac-2bc^{2})+3bc^{2}-2ac-b\right)\partial_{c}$
$\displaystyle~{}~{}+2m\big{(}a+2(m-1)cb\big{)}\,,$

and since $\hat{O}_{12}(0)\tilde{f}_{n,0}=0$ we end up with recursion
relations

$\displaystyle f_{n+1,m}$
$\displaystyle=\left(2c\hat{O}_{2}(m)+2b\hat{O}_{1}(m)+\hat{O}_{12}(m)+2a\right)g_{n,m}$
(A.17)
$\displaystyle~{}~{}~{}+\left(-c\hat{O}_{2}(m)-b\hat{O}_{1}(m)+2a(n+1)\right)f_{n,m}-(2n+3)\hat{O}_{12}(m+1)f_{n,m+1}+(bc-a)g_{n,m-1}\,,$
$\displaystyle g_{n+1,m}$
$\displaystyle=\hat{O}_{12}(m+1)f_{n,m+1}+(bc-a)f_{n,m}+\Big{(}-c\hat{O}_{2}(m)-b\hat{O}_{1}(m)+(2n+1)bc-a\Big{)}g_{n,m}\,.$

These reproduce the same results as obtained before. One can imagine pushing
the recursion one layer further, assuming a polynomial form for
$f_{n,m}(a,b,c)$ and $g_{n,m}(a,b,c)$. However, this does not seem practical,
given that the resulting expression for coefficients requires three additional
indices for powers of $a$, $b$ and $c$.

We can use the form above to factor out the $p$ dependence of these integrals.

$\displaystyle\mathcal{E}$ $\displaystyle=4\pi\sum_{n}\frac{(-1)^{n}}{n!}\int
dq\ q^{2}\
e^{\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})(X+Y)}\beta^{n}(pq)^{-n}\Big{(}F_{n}-(pq)^{2}\big{(}\frac{n\Theta_{n>0}}{\beta}\big{)}G_{n-1}\Big{)}j_{n}(pq)$
$\displaystyle=\sum_{n,m}4\pi\frac{(-1)^{n+m}}{n!}\int dq\ q^{2}\
e^{\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})(X+Y)}\beta^{n}(pq)^{2m-n}\Big{(}f_{n,m}+\big{(}\frac{n\Theta_{n>0,m>0}}{\beta}\big{)}g_{n-1,m-1}\Big{)}j_{n}(pq)$
$\displaystyle\equiv\sum_{n,m}\frac{(-1)^{n+m}}{n!}p^{2m-n}\epsilon_{n,m}(p)$
(A.18)

where as a reminder $F_{n}$ and $G_{n}$ are functions of $a,b,c$ and $pq$. In
the above we have defined

$\epsilon_{n,m}(p)=4\pi\int dq\ q^{2}\
e^{\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})(X+Y)}\beta^{n}q^{2m-n}\Big{(}f_{n,m}+\big{(}\frac{n\Theta_{n>0,m>0}}{\beta}\big{)}g_{n-1,m-1}\Big{)}j_{n}(pq)$

which can be computed once for all $p$ via Hankel transform. Furthermore,
since the only dependence on $\hat{p}$ comes in via the $f_{n,m},g_{n,m}$ we
can precompute the scale dependence of these integrals as

$\epsilon_{n,m}(p)=f_{n,m}\Gamma_{n,m}+ng_{n-1,m-1}\Gamma^{(-)}_{n,m}$

where we have defined

$\displaystyle\Gamma_{n,m}$ $\displaystyle=4\pi\int dq\ q^{2}\
e^{\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})(X+Y)}\beta^{n}q^{2m-n}j_{n}(pq)$
(A.19) $\displaystyle\Gamma^{(-)}_{n,m}$
$\displaystyle=4\pi\Theta_{n>0,m>0}\int dq\ q^{2}\
e^{\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})(X+Y)}\beta^{n-1}q^{2m-n}j_{n}(pq).$
(A.20)

For any given triangle $\textbf{k}_{1,2,3}$ these components can be pre-
computed for all $p$ such that for a given p evaluating
$\mathcal{E}(\textbf{k}_{1},\textbf{k}_{2},\textbf{p})$ reduces to simple
matrix multiplication.

###  A.2 Recovering Tree Level Results

The various kernels and infinite sums in the above being somewhat opaque, in
this section we show how to recover the expected results at tree level by
evaluating the first few summands at lowest order. We can begin by writing

$\displaystyle\mathcal{E}=4\pi\sum_{n}\frac{(-1)^{n}}{n!}\int dq\
q^{2}\Big{(}1$
$\displaystyle+\frac{\textbf{k}_{1}\cdot\textbf{k}_{2}}{2}(X+Y)\Big{)}\Big{(}\frac{k_{1}k_{2}Y}{2}\Big{)}^{n}$
$\displaystyle(pq)^{-n}\Big{(}F_{n}-\frac{2n\Theta_{n>0}}{k_{1}k_{2}Y}(pq)^{2}G_{n-1}\Big{)}j_{n}(pq)+\mathcal{O}(P_{\rm
L}^{2}).$ (A.21)

Since each $X,Y$ carries one power of $P_{\rm L}$, we can stop the sum at
$n=2$ in order to recover linear theory:

$\displaystyle\mathcal{E}_{\rm lin}=4\pi\int$ $\displaystyle dq\
q^{2}\Bigg{\\{}\Big{(}\frac{\textbf{k}_{1}\cdot\textbf{k}_{2}}{2}\Big{)}(X+Y)\
F_{0}\ j_{0}(pq)$ $\displaystyle-(pq)^{-1}\Big{(}\frac{k_{1}k_{2}Y}{2}\Big{)}\
F_{1}\
j_{1}(pq)+(pq)^{-2}\Big{(}\frac{k_{1}k_{2}Y}{2}\Big{)}^{2}\Big{(}-\frac{4}{k_{1}k_{2}Y}G_{1}(pq)^{2}\Big{)}j_{2}(pq)\Bigg{\\}}.$

Using the definitions of $F_{n}$, $G_{n}$ and $X$ and $Y$ as well as the
identity $j_{1}(x)/x=(j_{0}+j_{2})/3$ it is straightforward to recover

$\mathcal{E}_{\rm lin}^{\rm
fin}(\textbf{k}_{1},\textbf{k}_{2},\textbf{p})=-\frac{(\textbf{k}_{1}\cdot\textbf{p})(\textbf{k}_{2}\cdot\textbf{p})}{p^{4}}P_{\rm
L}(p).$ (A.22)

It is interesting to note that unlike in the case of the Zeldovich power
spectrum, where $\textbf{k}_{1}=-\textbf{k}_{2}$, recovering linear theory
requires going to $n=2$. In general, in order to isolate contributions to
$\mathcal{E}^{\rm fin}$ at a given order we note that we can write

$\Gamma^{(-)}_{n,m}=\sum_{r}\frac{4\pi}{r!}\int dq\ q^{2}\
\left(\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})(X+Y)\right)^{r}\beta^{n(-1)}q^{2m-n}j_{n}(pq)$

such that the order-$N$ contributions to $\Gamma^{(-1)}_{n,m}$ come from
$N=r+n(-1)$, and the corresponding contributions to $\mathcal{E}$ come from
setting all other terms in the sum to zero.

We can now proceed to obtain the tree-level bispectrum in our numerical
implementation. Pluggin in the linear contribution to $\mathcal{E}$ as derived
above into Equation 3.11 we get to leading order

$\displaystyle B^{\rm tree}(\textbf{k}_{1},\textbf{k}_{2})$
$\displaystyle=\mathcal{E}_{\rm lin}^{\rm
fin}(\textbf{k}_{1},\textbf{k}_{1}+\textbf{k}_{2},-\textbf{k}_{1}-\textbf{p})\mathcal{E}_{\rm
lin}^{\rm
fin}(\textbf{k}_{2},\textbf{k}_{1}+\textbf{k}_{2},-\textbf{k}_{2}-\textbf{p})+...$
$\displaystyle=\Big{(}1+\frac{\textbf{k}_{1}\cdot\textbf{k}_{2}}{k_{1}^{2}}\Big{)}\Big{(}1+\frac{\textbf{k}_{1}\cdot\textbf{k}_{2}}{k_{2}^{2}}\Big{)}P_{\rm
L}(\textbf{k}_{1})P_{\rm L}(\textbf{k}_{2})+\textrm{cycl.}$ (A.23)

as expected.

###  A.3 Alternative Angular Decomposition

Figure 13: Convergence of the Hankel-transform series in Equation A.26
compared to the fiducial method described in §A.1 in the case of the Zeldovich
power spectrum. The series both converge to the correct values as the number
$N_{j}$ of Hankel transforms increases, though Equation A.26 does so
substantially faster at almost all wavenumbers.

It is also possible to perform the angular integral in $\mathcal{E}$ by
directly Taylor-expanding the angular dependence in $Y(q)$. In this case we
have

$\displaystyle\int d^{3}\textbf{q}\
e^{-i\textbf{p}\cdot\textbf{q}+\frac{1}{2}k_{1,i}k_{2,j}A_{ij}(\textbf{q})}$
$\displaystyle=\int d^{3}\textbf{q}\
e^{-i\textbf{p}\cdot\textbf{q}+\frac{1}{2}((\textbf{k}_{1}\cdot\textbf{k}_{2})X(q)+(\textbf{k}_{1}\cdot\hat{q})(\textbf{k}_{2}\cdot\hat{q})Y(q))}$
$\displaystyle=\sum_{n=0}^{\infty}\frac{1}{n!}\int d^{3}\textbf{q}\
e^{-i\textbf{p}\cdot\textbf{q}+\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})X(q)}k_{1,i_{1}}...k_{1,i_{n}}k_{2,j_{1}}...k_{2,j_{n}}\left(\frac{Y(q)}{2}\right)^{n}\hat{q}_{i_{1}}...\hat{q}_{j_{n}}$
$\displaystyle=\sum_{n=0}^{\infty}\frac{1}{n!}\int d^{3}\textbf{q}\
e^{-i\textbf{p}\cdot\textbf{q}+\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})X(q)}k_{1,i_{1}}...k_{1,i_{n}}k_{2,j_{1}}...k_{2,j_{n}}\left(\frac{Y(q)}{2}\right)^{n}\sum_{m=0}^{n}C_{n,m}[\mathcal{L}_{2m}(\hat{q})]_{i_{1}...j_{n}}$
$\displaystyle=\sum_{n=0}^{\infty}\sum_{m=0}^{n}4\pi\left(k_{1,i_{1}}...k_{1,i_{n}}k_{2,j_{1}}...k_{2,j_{n}}[\mathcal{L}_{2m}(\hat{p})]_{i_{1}...j_{n}}\right)\frac{(-1)^{m}}{n!}\int
dq\ q^{2}\
e^{\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})X(q)}C_{n,m}\left(\frac{Y(q)}{2}\right)^{n}j_{2m}(pq)$
$\displaystyle=\sum_{m=0}^{\infty}\sum_{n=m}^{\infty}\frac{C_{n,m}}{n!}\left(k_{1,i_{1}}...k_{1,i_{n}}k_{2,j_{1}}...k_{2,j_{n}}[\mathcal{L}_{2m}(\hat{p})]_{i_{1}...j_{n}}\right)4\pi(-1)^{m}\int
dq\ q^{2}\
e^{\frac{1}{2}(\textbf{k}_{1}\cdot\textbf{k}_{2})X(q)}\left(\frac{Y(q)}{2}\right)^{n}j_{2m}(pq)$

where we have defined the coefficients

$\mu^{2n}=\sum_{m=0}^{n}C_{n,m}\mathcal{L}_{2m}(\mu).$ (A.24)

As in the previous method the final expression is a series of Hankel
transforms whose integrands are independent of $\hat{p}$ and proportional to
powers of $Y(q)$, multiplied by the angular coefficient

$\mathcal{D}_{2m}(\hat{k}_{1},\hat{k}_{2},\hat{p})=\hat{k}_{1,i_{1}}...\hat{k}_{1,i_{n}}\hat{k}_{2,j_{1}}...\hat{k}_{2,j_{n}}[\mathcal{L}_{2m}(\hat{p})]_{i_{1}...j_{n}}.$
(A.25)

Here $\mathcal{L}_{2m}$ are the tensor equivalents of the Legendre polynomials
(see e.g. ref. [49]). The $\mathcal{D}_{2m}$ are polynomials in the dot
products of the three unit vector arguments and can be worked out recursively
using the recurrence relations of the Legendre polynomials.

Let us consider the special case of the power spectrum
$\textbf{k}=\textbf{k}_{1}=-\textbf{k}_{2}=\textbf{p}$. In this case the
angular contraction is

$\mathcal{D}(\hat{k},-\hat{k},\hat{k})=(-1)^{n}\hat{k}_{i_{1}}...\hat{k}_{j_{n}}[\mathcal{L}_{2m}(\hat{k})]_{i_{1}...j_{n}}=(-1)^{n}$

such that we can write

$P_{\rm Zel}(k)=\sum_{m=0}^{\infty}4\pi(-1)^{m}\int dq\ q^{2}\
e^{-\frac{1}{2}k^{2}X(q)}j_{2m}(pq)F_{m}(q)$

where

$\displaystyle F_{m}(q)$
$\displaystyle=\sum_{n=m}^{\infty}\frac{C_{n,m}}{n!}\left(-\frac{k^{2}Y(q)}{2}\right)^{n}$
$\displaystyle=\frac{(4m+1)\Gamma\left(m+\frac{3}{2}\right)\,}{(2m+1)\Gamma(2m+\frac{3}{2})}\
_{1}F_{1}\left(m+\frac{1}{2};2m+\frac{3}{2};-\frac{k^{2}Y(q)}{2}\right)\left(-\frac{k^{2}Y(q)}{2}\right)^{m}.$
(A.26)

Note that $F_{m}$ is proportional to $Y^{m}$ to leading order, so this formula
is a re-summation of the method presented above where terms with different
powers were separated. Figure 13 shows the convergence of this sum as
$m_{max}$ is increased—it is in excellent agreement with the fiducial method
and indeed exhibits significantly faster converge at large wavenumbers.

##  Appendix B IR Cancellation and Enhancement by BAO at 1-loop Order

Let us first consider the infrared properties of the kernels $Z_{n}$, i.e.
their asymptotic limits when one or more momenta are small compared to the
total momentum. In particular, if the momenta $\textbf{p}_{m}$ satisfy
$p_{m}\ll k$ we have that

$Z_{n}(\textbf{p}_{1},...,\textbf{p}_{m},\textbf{k}_{1},...,\textbf{k}_{n-m})\rightarrow\frac{(n-m)!}{n!}\left(\frac{\textbf{p}_{1}\cdot\textbf{k}}{p_{1}^{2}}\right)...\left(\frac{\textbf{p}_{m}\cdot\textbf{k}}{p_{m}^{2}}\right)Z_{n-m}(\textbf{k}_{1},...,\textbf{k}_{n-m}),$
(B.1)

i.e. the kernels factor into “shifts” and lower-order kernels. Here
$\textbf{k}=\textbf{k}_{1}+...+\textbf{k}_{n-m}\approx\textbf{p}_{1}+...+\textbf{p}_{m}+\textbf{k}_{1}+...+\textbf{k}_{n-m}$.

Let us now consider each contributing diagram to the 1-loop bispectrum. Since
the IR limits of the kernels $Z_{n}$ reduce to shift terms involving long-
wavelength modes $\textbf{p}_{m}$ we integrate over and lower order kernels,
let us look specifically for diagrams whose IR limits approach piece of the
tree level bispectrum proportional to
$Z_{2}(\textbf{k}_{1},\textbf{k}_{2})P_{\rm L}(k_{1})P_{\rm L}(k_{2})$. All
other IR contributions will be equivalent up to permutations of the external
wavenumbers.

###  B.1 The 114 contribution

In this case we want to consider the diagram

$\displaystyle B_{114}$
$\displaystyle=\left\langle\delta^{(1)}(\textbf{k}_{1})\delta^{(1)}(\textbf{k}_{2})\delta^{(4)}(\textbf{k}_{3})\right\rangle^{\prime}$
$\displaystyle=12\int_{\textbf{p}}Z_{4}(\textbf{p},-\textbf{p},\textbf{k}_{1},\textbf{k}_{2})\
P_{\rm L}(p)P_{\rm L}(k_{1})P_{\rm L}(k_{2})$
$\displaystyle\supset\int_{\textbf{p}}^{k_{\rm
IR}}\left(\frac{\textbf{p}\cdot\textbf{k}_{12}}{p^{2}}\right)\left(\frac{-\textbf{p}\cdot\textbf{k}_{12}}{p^{2}}\right)Z_{2}(\textbf{k}_{1},\textbf{k}_{2})\
P_{\rm L}(p)P_{\rm L}(k_{1})P_{\rm L}(k_{2})$
$\displaystyle=\left(-k_{3,i}k_{3,j}\int_{\textbf{p}}^{k_{\rm
IR}}\frac{p_{i}p_{j}}{p^{4}}P_{\rm
L}(p)\right)Z_{2}(\textbf{k}_{1},\textbf{k}_{2})\ P_{\rm L}(k_{1})P_{\rm
L}(k_{2})$ $\displaystyle=-\frac{1}{2}k_{3}^{2}\Sigma^{2}_{\rm IR}\
2Z_{2}(\textbf{k}_{1},\textbf{k}_{2})\ P_{\rm L}(k_{1})P_{\rm L}(k_{2})$

where we have defined

$\Sigma^{2}_{\rm IR}=\frac{1}{3}\int_{\textbf{p}}^{k_{\rm IR}}\frac{P_{\rm
L}(p)}{p^{2}}.$ (B.2)

and used that $\textbf{k}_{12}=-\textbf{k}_{3}$.

###  B.2 The 123 contribution

In this case we have the diagram

$\displaystyle B_{123}$
$\displaystyle=\left\langle\delta^{(1)}(\textbf{k}_{1})\delta^{(2)}(\textbf{k}_{2})\delta^{(3)}(\textbf{k}_{3})\right\rangle^{\prime}$
$\displaystyle\supset
6\int_{\textbf{p}}Z_{2}(\textbf{p},\textbf{k}_{2}-\textbf{p})Z_{3}(\textbf{p},\textbf{k}_{2}-\textbf{p},\textbf{k}_{1})\
P_{\rm L}(p)P_{\rm L}(|\textbf{k}_{2}-\textbf{p}|)P_{\rm L}(k_{1})$
$\displaystyle\supset 12\int_{\textbf{p}}^{k_{\rm
IR}}Z_{2}(\textbf{p},\textbf{k}_{2}-\textbf{p})Z_{3}(\textbf{p},\textbf{k}_{2}-\textbf{p},\textbf{k}_{1})\
P_{\rm L}(p)P_{\rm L}(|\textbf{k}_{2}-\textbf{p}|)P_{\rm L}(k_{1})$
$\displaystyle\rightarrow 2\int_{\textbf{p}}^{k_{\rm
IR}}\left(\frac{\textbf{p}\cdot\textbf{k}_{2}}{p^{2}}\right)\left(\frac{\textbf{p}\cdot\textbf{k}_{12}}{p^{2}}\right)Z_{2}(\textbf{k}_{1},\textbf{k}_{2})\
P_{\rm L}(p)P_{\rm L}(|\textbf{k}_{2}-\textbf{p}|)P_{\rm L}(k_{1})$
$\displaystyle=-2k_{2,i}k_{3,j}Z_{2}(\textbf{k}_{1},\textbf{k}_{2})\int_{\textbf{p}}^{k_{\rm
IR}}\frac{p_{i}p_{j}}{p^{4}}P_{\rm L}(p)P_{\rm
L}(|\textbf{k}_{2}-\textbf{p}|)P_{\rm L}(k_{1})$
$\displaystyle\rightarrow-(\textbf{k}_{2}\cdot\textbf{k}_{3})\Sigma^{2}_{\rm
IR}\ 2Z_{2}(\textbf{k}_{1},\textbf{k}_{2})P_{\rm L}(k_{1})P_{\rm L}(k_{2}).$

In the third line, we have used that there are two IR poles in the integrand
at $\textbf{p}=\bf{0},\textbf{k}_{2}$, and that the integrand is invariant
when shifting to the latter. There is an additional contribution to in
$B_{213}$ with $\textbf{k}_{1}\leftrightarrow\textbf{k}_{2}$, such that the
total IR contribution is

$-(\textbf{k}_{2}\cdot\textbf{k}_{3}+\textbf{k}_{1}\cdot\textbf{k}_{3})\Sigma^{2}_{\rm
IR}\ 2Z_{2}(\textbf{k}_{1},\textbf{k}_{2})P_{\rm L}(k_{1})P_{\rm
L}(k_{2})=+k_{3}^{2}\Sigma^{2}_{\rm IR}\
2Z_{2}(\textbf{k}_{1},\textbf{k}_{2})P_{\rm L}(k_{1})P_{\rm L}(k_{2}).$ (B.3)

###  B.3 The 132 contribution

Here we have the diagram

$\displaystyle B_{132}$
$\displaystyle=\left\langle\delta^{(1)}(\textbf{k}_{1})\delta^{(3)}(\textbf{k}_{2})\delta^{(2)}(\textbf{k}_{3})\right\rangle^{\prime}$
$\displaystyle\supset
6\int_{\textbf{p}}Z_{2}(\textbf{k}_{1},\textbf{k}_{2})Z_{3}(\textbf{p},-\textbf{p},\textbf{k}_{2})\
P_{\rm L}(p)P_{\rm L}(k_{1})P_{\rm L}(k_{2})$
$\displaystyle\rightarrow\int_{\textbf{p}}^{k_{\rm
IR}}\left(\frac{\textbf{p}\cdot\textbf{k}_{2}}{p^{2}}\right)\left(\frac{-\textbf{p}\cdot\textbf{k}_{2}}{p^{2}}\right)\
Z_{2}(\textbf{k}_{1},\textbf{k}_{2})\ P_{\rm L}(p)P_{\rm L}(k_{1})P_{\rm
L}(k_{2})$ $\displaystyle=\left(-\frac{1}{2}k_{2}^{2}\Sigma^{2}_{\rm
IR}\right)\ 2Z_{2}(\textbf{k}_{1},\textbf{k}_{2})P_{\rm L}(k_{1})P_{\rm
L}(k_{2}).$

Here again we can permute $\textbf{k}_{1,2}$ to get the total contribution

$\left(-\frac{1}{2}(k_{1}^{2}+k_{2}^{2})\Sigma^{2}_{\rm IR}\right)\
2Z_{2}(\textbf{k}_{1},\textbf{k}_{2})P_{\rm L}(k_{1})P_{\rm L}(k_{2}).$ (B.4)

###  B.4 The 222 contribution

Finally we have the diagram

$\displaystyle B_{222}$
$\displaystyle=\left\langle\delta^{(2)}(\textbf{k}_{1})\delta^{(2)}(\textbf{k}_{2})\delta^{(2)}(\textbf{k}_{2})\right\rangle^{\prime}$
$\displaystyle\supset
8\int_{\textbf{p}}Z_{2}(\textbf{p},\textbf{k}_{1}-\textbf{p})Z_{2}(-\textbf{p},\textbf{k}_{2}+\textbf{p})Z_{2}(\textbf{k}_{1}-\textbf{p},\textbf{k}_{2}+\textbf{p})\
P_{\rm L}(p)P_{\rm L}(|\textbf{k}_{1}-\textbf{p}|)P_{\rm
L}(|\textbf{k}_{2}+\textbf{p}|)$ $\displaystyle\rightarrow
2\int_{\textbf{p}}^{k_{\rm
IR}}\left(\frac{\textbf{p}\cdot\textbf{k}_{1}}{p^{2}}\right)\left(\frac{-\textbf{p}\cdot\textbf{k}_{2}}{p^{2}}\right)Z_{2}(\textbf{k}_{1},\textbf{k}_{2})\
P_{\rm L}(p)P_{\rm L}(|\textbf{k}_{1}-\textbf{p}|)P_{\rm
L}(|\textbf{k}_{2}+\textbf{p}|)$
$\displaystyle\rightarrow-(\textbf{k}_{1}\cdot\textbf{k}_{2})\Sigma^{2}_{\rm
IR}\ 2Z_{2}(\textbf{k}_{1},\textbf{k}_{2})\ P_{\rm L}(k_{1})P_{\rm L}(k_{2}).$

Note that this diagram has three IR regions, identical to the full resummation
discussed in §3.2–however, shifting to the other IR regions produces terms
proportional to $P(k_{3})$ instead, which are not degenerate with the terms
discussed here but rather permutations thereof, so we can ignore them.

Combining the previous results we have that the full IR contribution the
1-loop bispectrum is given by

$\frac{1}{2}\left(-k_{3}^{2}+2k_{3}^{2}-(k_{1}^{2}+k_{2}^{2})-2\textbf{k}_{1}\cdot\textbf{k}_{2}\right)\Sigma^{2}_{\rm
IR}\ 2Z_{2}(\textbf{k}_{1},\textbf{k}_{2})P_{\rm L}(k_{1})P_{\rm L}(k_{2})=0,$
(B.5)

that is long-wavelength displacements do not contribute to the 1-loop
bispectrum.

###  B.5 Enhancement by the BAO

In the above cancellations we have implicitly used that we can approximate

$P_{\rm L}(|\textbf{k}+\textbf{p}|)=P_{\rm
L}(k)\left(1+\mathcal{O}\left(\frac{p}{k}\right)\right).$ (B.6)

However, when the power spectrum contains an oscillatory component
$P_{w}\sim\sin(r_{s}k)$ we have that the corrections are instead of order
powers of $r_{s}$ times the signal, or enhanced by $r_{s}k$ compared to
expectations. In other words, a mode p which is in the infrared compared to k
can still have a large effect due to the size of $r_{s}$.

In the 1-loop bispectrum we can consider this enhancement by looking at the
contributions due to $B_{123}$ and $B_{222}$ for which we have the IR
contribution

$\displaystyle-2\left(k_{2,i}k_{3,j}+k_{2,i}k_{1,j}\right)$ $\displaystyle
Z_{2}(\textbf{k}_{1},\textbf{k}_{2})\int_{\textbf{p}}^{k_{\rm
IR}}\frac{p_{i}p_{j}}{p^{4}}P_{\rm L}(p)P_{\rm
L}(|\textbf{k}_{2}-\textbf{p}|)P_{\rm L}(k_{1})$
$\displaystyle=2k_{2,i}k_{2,j}Z_{2}(\textbf{k}_{1},\textbf{k}_{2})\int_{\textbf{p}}^{k_{\rm
IR}}\frac{p_{i}p_{j}}{p^{4}}P_{\rm L}(p)P_{\rm
L}(|\textbf{k}_{2}-\textbf{p}|)P_{\rm L}(k_{1})$

where we have dropped the p in $P_{\rm L}(|\textbf{k}_{1}-\textbf{p}|)$ to
focus on the effect on the BAO in $P_{\rm L}(k_{2})$. Isolating the BAO
component in $P_{\rm L}(\textbf{k}_{2}-\textbf{p})$ then gives us, subtracting
the $\textbf{p}\rightarrow 0$ limit by previous arguments

$\displaystyle k_{2,i}k_{2,j}$ $\displaystyle\
2Z_{2}(\textbf{k}_{1},\textbf{k}_{2})P_{\rm L}(k_{1})\int_{\textbf{p}}^{k_{\rm
IR}}\frac{p_{i}p_{j}}{p^{4}}P_{\rm
L}(p)(P_{w}(|\textbf{k}_{2}-\textbf{p}|)-P_{w}(k_{2}))$
$\displaystyle=\frac{1}{2}k_{2,i}k_{2,j}\
2Z_{2}(\textbf{k}_{1},\textbf{k}_{2})P_{\rm L}(k_{1})\int_{\textbf{p}}^{k_{\rm
IR}}\frac{p_{i}p_{j}}{p^{4}}P_{\rm
L}(p)(P_{w}(|\textbf{k}_{2}-\textbf{p}|)+P_{w}(|\textbf{k}_{2}+\textbf{p}|)-2P_{w}(k_{2}))$
$\displaystyle=\frac{1}{2}k_{2,i}k_{2,j}\
2Z_{2}(\textbf{k}_{1},\textbf{k}_{2})P_{\rm L}(k_{1})\int d^{3}\textbf{x}\
e^{-i\textbf{k}_{2}\cdot\textbf{x}}\xi_{w}(\textbf{x})\
\int_{\textbf{p}}^{k_{\rm IR}}\frac{p_{i}p_{j}}{p^{4}}P_{\rm L}(p)\
(2\cos(\textbf{p}\cdot\textbf{x})-2).$

Since the BAO component of the correlation function is peaked at $x=r_{s}$ we
see that any contribution with $p\ll 1/r_{s}$ is suppressed, but scales
between $1/r_{s}$ and $k_{2}$ can contribute to the BAO despite being in the
IR regime. Indeed, setting $\textbf{x}=r_{s}\hat{x}$ and performing the p
integral we can see that it yields precisely
$A_{ij}(\textbf{x})=X(r_{s})\delta_{ij}+Y(r_{s})\hat{x}_{i}\hat{x}_{j}$, such
that the total contribution to the bispectrum is

$-\frac{1}{2}k_{2}^{2}\Sigma^{2}_{s}P_{w}(k_{2})$ (B.7)

i.e. the leading order contribution to the damping of the BAO derived using
the wiggle no-wiggle split in §4.

##  Appendix C Corrections to Saddle-Point in Power Spectrum

Let us consider the integral

$P_{w,\rm Zel}(k)=\int d^{3}\textbf{q}\
e^{i\textbf{k}\cdot\textbf{q}-\frac{1}{2}k_{i}k_{j}A_{ij}^{nw}(\textbf{q})}\
\xi_{w}(q).$ (C.1)

To use the saddle-point approximation we make the Taylor expansion about
$\textbf{q}=r_{d}\hat{q}+(q-r_{d})\hat{q}$:

$A_{ij}^{nw}(\textbf{q})=A_{ij}^{nw}(r_{d}\hat{q})+(q-r_{d})\
\hat{q}_{a}\partial_{a}A_{ij}^{nw}(r_{d}\hat{q})+\ldots\ .$ (C.2)

The first term gives us the leading order saddle-point result, i.e.
exponential damping. The second term leads to a correction

$-\frac{1}{2}k_{i}k_{j}\int d^{3}\textbf{q}\
e^{i\textbf{k}\cdot\textbf{q}-\frac{1}{2}k_{i}k_{j}A_{ij}^{nw}(\textbf{q})}\
(q-r_{d})\ \hat{q}_{a}\partial_{a}A_{ij}^{nw}(r_{d}\hat{q})\ \xi_{w}(q).$
(C.3)

Defining

$\xi^{\ell}_{n}(q)=\int\frac{dk\ k^{2+n}}{2\pi^{2}}P(k)j_{\ell}(kq)$ (C.4)

we have that the derivative above is given by

$\displaystyle\frac{\partial A^{nw}_{ij}(\textbf{q})}{\partial
q_{k}}=\frac{2}{5}\xi^{1,nw}_{-1}(q)\hat{q}_{(i}\delta_{jk)}-\frac{4}{5}\xi^{3,nw}_{-1}(q)\left(\frac{5}{2}\hat{q}_{i}\hat{q}_{j}\hat{q}_{k}-\frac{1}{2}\hat{q}_{(i}\delta_{jk)}\right)$
(C.5)

and, working to leading order we get that the correction is

$-\frac{1}{2}k^{2}\Big{[}\frac{2}{3}\xi^{1}_{-1}(r_{d})P^{1}_{w,0}(k)-\frac{2}{3}\left(\frac{4}{5}\xi^{1}_{-1}(r_{d})-\frac{6}{5}\xi^{3}_{-1}(r_{d})\right)P^{1}_{w,2}(k)\Big{]}.$

In the above we have defined

$P^{n}_{w,\ell}(k)=4\pi\int dq\ q^{2}\ (q-r_{d})^{n}\ \xi_{w}(q)\
j_{\ell}(kq).$ (C.6)

Roughly speaking these should be order the width of the (linear) BAO
$r_{D}^{n}$. The cumulative correction to the damping is thus of order
$r_{D}r_{d}\sigma^{2}_{r_{d}}$, where we have used that
$\xi^{1}_{-1}(q)=q\sigma^{2}_{q}/3$ is related to the mean square density
$\sigma^{2}_{q}$ in spheres of radius $q$—in practice we find that the effect
is sufficiently small to be subdominant even to differences due to the choice
of the wiggle no-wiggle split for the BAO, though its relevance is
proportional to the width of the feature in the correlation function which can
vary for e.g. primordial features.

##  Appendix D Angular Integral in Wiggle No-Wiggle Split

In this section we perform the angular integral in Equation 4.2 to derive the
damping form for the BAO component of the power spectrum with
$\Sigma^{2}_{s}=X_{s}+Y_{s}$. In order to do so we will use the series and
angular decomposition in §A.3. Let us take the special case where the
Zeldovich integrand multiplies $\xi_{w}(q)$ such that $X,Y$ are isolated at
their values at the peak $X_{s},Y_{s}$. In the power spectrum for example this
corresponds to the BAO component of the $b_{1}^{2}$ term

$\displaystyle P_{w}^{NL}(k)$ $\displaystyle=\int d^{3}\textbf{q}\
e^{-i\textbf{k}\cdot\textbf{q}-\frac{1}{2}k_{i}k_{j}A_{ij}(\textbf{q})}\
\xi_{w}(q)$
$\displaystyle\approx\sum_{n=0}^{\infty}\sum_{m=0}^{n}4\pi\frac{(-1)^{m}}{n!}\int
dq\ q^{2}\
e^{-\frac{1}{2}k^{2}X_{s}}C_{n,m}\left(-\frac{k^{2}Y_{s}}{2}\right)^{n}\
\xi_{w}(q)\ j_{2m}(kq)$
$\displaystyle=e^{-\frac{1}{2}k^{2}X_{s}}\sum_{n=0}^{\infty}\frac{1}{n!}\
\left(-\frac{k^{2}Y_{s}}{2}\right)^{n}\ \sum_{m=0}^{n}C_{n,m}\ P_{w,2m}(k).$
(D.1)

Here we have defined

$P_{w,\ell}=4\pi(-i)^{\ell}\int dq\ q^{2}\ \xi_{w}(q)\ j_{\ell}(kq).$ (D.2)

Since $\xi_{w}$ is well localized at $r_{s}$, for any sufficiently large
$k>r_{s}^{-1}$ wherever the integral has support we can make the asymptotic
approximation that

$j_{2m}(kq)\approx\frac{1}{kq}\cos(kq-\pi/2-m\pi)=\frac{(-1)^{m}}{kq}\cos(kq-\pi/2)\approx(-1)^{m}j_{0}(kq)$
(D.3)

such that $P_{w,2m}\approx(-1)^{m}P_{w}$. Plugging this approximation back
into Equation D.1 and using $\sum_{m=0}^{n}C_{n,m}=1$ yields

$\displaystyle P_{w}^{NL}(k)=e^{-\frac{1}{2}k^{2}X_{s}}$
$\displaystyle\sum_{n=0}^{\infty}\frac{1}{n!}\left(-\frac{k^{2}Y_{s}}{2}\right)^{n}\left(\sum_{m=0}^{n}C_{n,m}\right)P_{w}(k)$
$\displaystyle=e^{-\frac{1}{2}k^{2}X_{s}}\sum_{n=0}^{\infty}\frac{1}{n!}\left(-\frac{k^{2}Y_{s}}{2}\right)^{n}P_{w}(k)=e^{-\frac{1}{2}k^{2}\Sigma^{2}_{s}}P_{w}(k)$
(D.4)

with $\Sigma^{2}_{s}=X_{s}+Y_{s}$ as desired.

We can also consider BAO contribtions with higher angular dependence in q
space. For example in the galaxy power spectrum we have linear terms
proportional to

$-\frac{1}{2}k_{i}k_{j}A^{w}_{ij}(\textbf{q})\sim
k^{2}\mathcal{L}_{0,2}(\mu)\int\frac{dp}{2\pi^{2}}\ j_{0,2}(pq)\
P_{w}(p),\quad ik_{i}U_{i}\sim k\mu\int\frac{dp\ p}{2\pi^{2}}\ j_{1}(pq)\
P_{w}(p)$

which Fourier transform into $P_{w}$ at leading order. For simplicity let us
consider the latter $ik_{i}U_{i}$ term: upon resumming the linear
displacements this yields

$\displaystyle P^{w}_{b_{1}}$ $\displaystyle=ik_{a}\int d^{3}\textbf{q}\
e^{-i\textbf{k}\cdot\textbf{q}-\frac{1}{2}k_{i}k_{j}A_{ij}(\textbf{q})}U_{w}(q)\hat{q}_{a}$
$\displaystyle=ik_{a}\sum_{n=0}^{\infty}\frac{\hat{k}_{i_{1}}...\hat{k}_{i_{2n}}}{n!}\int
d^{3}\textbf{q}\ e^{-i\textbf{k}\cdot\textbf{q}-\frac{1}{2}k^{2}X(q)}\
U_{w}(q)\left(\frac{-k^{2}Y(q)}{2}\right)^{n}\hat{q}_{a}\hat{q}_{i_{1}}...\hat{q}_{i_{2n}}$
$\displaystyle\approx
ik_{a}e^{-\frac{1}{2}k^{2}X_{s}}\sum_{n=0}^{\infty}\frac{\hat{k}_{i_{1}}...\hat{k}_{i_{2n}}}{n!}\left(\frac{-k^{2}Y_{s}}{2}\right)^{n}\int
d^{3}\textbf{q}\ e^{-i\textbf{k}\cdot\textbf{q}}\ U_{w}(q)\
\hat{q}_{a}\hat{q}_{i_{1}}...\hat{q}_{i_{2n}}$

From here we can again decompose

$\mu^{2n+1}=\sum_{m=0}^{n}C_{n,m}^{(1)}\mathcal{L}_{2m+1}(\mu)$ (D.5)

to write

$\int d^{3}\textbf{q}\ e^{-i\textbf{k}\cdot\textbf{q}}\ U_{w}(q)\