{'fraction_non_alphanumeric': 0.03649283591506035, 'fraction_numerical': 0.012983347445667513, 'mean_word_length': 4.376896982681664, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Time-sliced perturbation theory with primordial non-Gaussianity and effects of large bulk flows on inflationary oscillating features

We extend the formalism of time-sliced perturbation theory (TSPT) for cosmological large-scale structure to include non-Gaussian initial conditions. We show that in such a case the TSPT interaction vertices acquire new contributions whose time-dependence factorizes for the Einstein-de Sitter cosmology. The new formulation is free from spurious infrared (IR) enhancements and reveals a clear IR structure of non-Gaussian vertices. We use the new technique to study the evolution of oscillating features in primordial statistics and show that they are damped due to non-linear effects of large bulk flows. We derive the damping factors for the oscillating primordial power spectrum and bispectrum by means of a systematic IR resummation of relevant Feynman diagrams.


Introduction
Large scale structure (LSS) provides us with a powerful tool to study the dynamical and statistical properties of our universe. This structure has been formed during the epoch of matter domination and its dynamics on large scales is governed by dark matter and baryons acting as single pressureless perfect fluids. The three-dimensional distribution of LSS can potentially supersede the cosmic microwave background (CMB) measurements in the number of Fourier modes available for observations. On the other hand, the relevant cosmological information encoded in LSS is concealed by the non-linear effects related to gravitational clustering. Disentangling the cosmological information from these effects presents a highly non-trivial challenge that must be addressed in order to efficiently use the data of the ongoing and future surveys.
The standard tool to model dark matter clustering in the non-linear regime is N-body simulations, which are progressively becoming more robust and accessible. On the other hand, reaching the accuracy level required by future surveys still remains computationally expensive [1]. Besides, implementing non-minimal cosmological scenarios or initial conditions (e.g. primordial non-Gaussianity) still requires non-trivial extensions of the standard N-body codes.
An alternative way is the development of analytic methods based on perturbation theory [2]. Although limited to a range of large scales, this approach has a number of advantages. First, analytic methods allow for a shorter computational time. Second, they give valuable insights into the physics of non-linear clustering. Third, they are flexible and can be easily extended to models beyond the standard ΛCDM cosmology.
One of the goals of future LSS surveys is probing primordial non-Gaussianity (PNG), which can shed light on the dynamics of inflation or its alternatives. This task is complicated by a number of factors. First, generic non-Gaussian phenomena are encoded in higher-order statistics, which are hard to measure and interpret. Second, the primordial non-Gaussian signal must be disentangled from that generated by non-linear gravitational instability even from initially Gaussian fields. Perturbation theory provides a natural and computationally efficient way to do this [3]. Perturbative calculations including primordial non-Gaussianity were performed for the power spectrum [4] and bispectrum [5] of matter and biased tracers. The analysis of galaxy bias was later refined in [6], see [7] for a review. The Effective Field Theory (EFT) of LSS [8,9] has been applied to compute the power spectrum and bispectrum with non-Gaussian initial conditions [10]. Based on these developments, an analysis of the ability to extract the PNG signal from LSS was carried out in [11][12][13][14][15].
A particularly interesting question is existence of oscillatory features in the statistics of primordial density fluctuations. Such resonance features naturally appear in models inspired by the axion monodromy inflation [16] where the inflaton potential is modulated by small sinusoidal oscillations [17,18]. Alternatively, they can be generated by the interaction between the inflaton and massive particles, thereby encoding information about the particle spectrum at inflation [19][20][21]. The latter type of non-Gaussianity has got the name of cosmological collider signatures.
Oscillating features in the primordial power spectrum have been constrained from the CMB data by the Planck collaboration [22,23]. The prospects to improve these constraints using future LSS and 21 cm surveys were studied in [24][25][26][27][28], see also the recent white paper [29] and references therein. The Planck collaboration has also performed searches for oscillating primordial non-Gaussianity in CMB with the null results [30]. The imprint of resonant non-Gaussianity on the scale dependent halo bias was discussed in [31], whereas Ref. [28] analyzed the sensitivity of 21 cm intensity mapping to this type of PNG. Refs. [32,33] presented forecasts for the reach of LSS and 21 cm surveys in probing cosmological collider signatures. As pointed out in Ref. [34], one expects primordial oscillations in LSS statistics to get damped at late times due to the non-linear effects of large bulk flows, similarly to the damping of baryon acoustic oscillations (BAO). Indeed, the displacement of modes with Fourier momentum k caused by non-linear coupling to long modes with momentum q k washes out the features with a characteristic period of oscillations in Fourier space k osc q smaller that the long-mode momentum. However, this effect was not discussed beyond the power spectrum.
In this paper we study the non-linear damping of primordial oscillations in the inflationary 3-point function (bispectrum). We work out a systematic procedure of infrared (IR) resummation that captures the effect of large bulk flows. Our method is inspired by techniques developed in the context of the BAO and is based on Time-Sliced Perturbation Theory (TSPT) [35]. The latter has been proven to be a powerful tool for IR resummation of the power spectrum and higher-order statistics [36,37].
TSPT is a novel way of performing cosmological perturbation theory utilizing the machinery of quantum field theory and statistical physics. The key idea of this method is to study the time-dependent probability distribution function of density and velocity fields, instead of these fields themselves as done in the standard cosmological perturbations theory (SPT). The central object is the generating functional for cosmological correlators whose time evolution is governed by a Liouville equation of statistical mechanics. A perturbative expansion of the generating functional leads to Feynman rules similar to that of a 3-dimensional Euclidean quantum field theory, in which time plays a role of an external parameter. This expansion produces equal-time correlation functions of density and velocity at fixed time (redshift) slices, which explains the name to the method.
TSPT deals directly with the statistical quantities, such as equal time n-point correlation functions, and thus provides a natural framework for studying non-Gaussian initial conditions. We will show that the incorporation of such conditions is straightforward and results in new contributions to the TSPT vertices. The new contributions take particularly simple form for the Einstein-de Sitter (EdS) cosmology, where time dependence of the TSPT vertices factors out. We introduce a new convenient diagrammatic technique for the primordial non-Gaussian contributions. As suggested by observations, the latter should be treated as small perturbations on top of the standard correlation function corresponding to Gaussian initial conditions. The eventual expressions for TSPT correlation functions at fixed order of perturbation theory coincide with the SPT result, although the individual diagrams are quite different. In particular, we will show that the non-Gaussian TSPT vertices do not contain any singularities originating from dynamical non-linear mode-coupling. In this sense they are IR-safe. Of course, they may still contain physical IR singularities which may be contained in the primordial statistics themselves.
Finally, we show how the new reformulation allows for a systematic IR resummation of the large bulk flow effects on the evolution of oscillatory primordial statistics. We show that IR resummation in this case closely resembles that of the BAO for the Gaussian initial statistics and leads to a smearing of oscillating features. In particular, the relevant contributions to be resummed have a graphical representation of daisy diagrams. We derive explicit expressions for the non-linear damping factors of the power spectrum and bispectrum. In the latter case the damping depends on the full configuration of the wavevectors, including their relative angles.
The paper is organized as follows. In Sec. 2 we set up conventions and discuss the standard treatment of non-Gaussian initial conditions in cosmological perturbation theory. In Sec. 3 we recapitulate the Gaussian TSPT framework. In Sec. 4 we extend it by including the non-Gaussian initial conditions. Section 5 contains the proof of the IR safety of the non-Gaussian TSPT vertices. In Sec. 6 we present a framework for IR resummation of non-Gaussian vertices and apply it to the non-linear evolution of the resonant power spectrum and bispectrum from axion monodromy inflation. We conclude in Sec. 7. Some additional material is given in Appendices.
In the numerical calculations throughout the paper we use the Planck 2018 cosmology [38] with massless neutrinos. The linear power spectrum is computed with the Boltzmann code CLASS [39].

Preliminaries
In this section we set up the conventions and discuss the standard treatment of non-Gaussian initial conditions in LSS. Throughout the paper we will mostly focus on the nontrivial initial 3-point function (bispectrum). However, our results can be straightforwardly generalized to the case of higher non-Gaussian statistics. We will comment on this point at the relevant places in the text.

Primordial and linear statistics
We will distinguish between primordial and linear bispectra. We adopt the inflationary paradigm and call primordial the quantities referring to the curvature perturbation ζ generated during inflation. Its statistical properties are encoded in the power spectrum and the bispectrum, D is the 3-dimensional Dirac delta-function. The power spectrum of the primordial curvature perturbations can be written as, where k * is a pivot scale, A 2 ζ is the amplitude and (n s − 1) is the tilt of the spectrum. Their mean values from the latest Planck CMB temperature and polarization measurements [38] for k * = 0.05 Mpc −1 are 1 Let us outline some properties of B ζ (k 1 , k 2 , k 3 ) in typical inflationary models, see Ref. [40] for more details. Due to momentum conservation and rotational invariance the bispectrum 1 Note that A 2 ζ is related to the scalar amplitude As used in the Planck analysis by A 2 ζ = 2π 2 As.
is a function of three independent variables only. These can be chosen to be either the three momenta norms k 1 , k 2 , k 3 , or the norms of a pair of momenta, say k 1 , k 2 , and the cosine µ 12 of the angle between them. Approximate scale invariance implies that B ζ is a homogeneous function of degree −6, symmetric in its arguments. The shapes commonly used in data analysis are 2 local [41], equilateral [42], and orthogonal [43], where in the last two expressions we for simplicity neglected the departures from scale invariance (the tilt). The current observational constraints on f N L (at 68% CL) from the CMB measurements are [30], The key objects for the description of LSS are the matter density contrast δ = (ρ−ρ)/ρ and velocity divergence Θ ∝ ∇ · u fields. Hereρ is the average density of the Universe and u is the peculiar flow velocity. Initially small, these fields start growing after recombination and become non-linear. This non-linear evolution is captured by the cosmological perturbation theory, which uses as a seed the linear fields δ L , Θ L evolved up to the present epoch as if the perturbations were always in the linear regime. In the case of adiabatic initial conditions corresponding to a growing mode the two linear fields are identical. They are related to the curvature perturbation ζ by a transfer function T (t 0 ; k) which encodes the evolution of perturbations from inflation, through recombination, up to the present time t 0 . Thus we write, We will understand by linear statistics the properties of the field δ L . Defining the linear matter power spectrum, we see that From this we immediately deduce the relation between the linear bispectrum of LSS and the primordial one, (2.9) 2 We use the Planck conventions [30] for the bispectrum amplitude fNL. Note that the Planck collaboration uses the variable Φ ≡ 3ζ/5, instead of ζ, in the analysis of non-Gaussianity.
Note that the linear bispectrum is parametrically suppressed by the amplitude of the primordial fluctuations. Indeed, neglecting the shape of the bispectrum one has, schematically, The linear power spectrum (2.8) and bispectrum (2.9) serve as the input for the equations describing the non-linear evolution of LSS statistics at times after recombination.

Standard perturbation theory with non-Gaussian initial conditions
We are interested in the correlation functions of the overdensity and the velocity divergence fields, whose time-evolution is governed by the following pressureless perfect fluid equations: where φ is the gravitational potential obeying the Poisson equation, Here t is the conformal time, the spatial derivatives are taken with respect to the comoving coordinates, H = aH is the rescaled Hubble parameter, a is the scale factor and Ω m is the matter density fraction. The perfect fluid approximation breaks down at short scales due to free streaming of dark matter particles. In the EFT of LSS the influence of short scale nonlinearities on physics at large scales is encapsulated by an effective stress tensor added to the r.h.s. of (2.11) [8,9]. For simplicity, we neglect the EFT corrections in this paper, keeping in mind that eventually they will have to be properly included. It is well-known [2] that in the case of an EdS universe dominated by non-relativistic matter (Ω m = 1) the above equations can be cast in a form free from any explicit time dependence. This is achieved by introducing the time parameter η = ln D(t) , (2.13) where D is the linear growth factor 3 , and appropriately rescaling the velocity divergence, (2.14) For the realistic ΛCDM cosmology, the above substitution leaves a small residual time dependence which, however, has little effect on the dynamics [44]. Following conventional practice we will neglect this explicit time dependence in the equations, even though our analysis does not crucially depend on this restriction.
Within the above approximation Eqs. (2.11) can be rewritten in Fourier space as 4 where we used the notations, and introduced the non-linear kernels (2.16) These equations can be solved perturbatively by using the following power series Ansatz: where δ L is the linear density field and the non-linear kernels F n and G n are recursively derived from (2.15). The correlation functions of the fields δ k , Θ k are computed by averaging the expressions (2.17) using the statistics of the linear field δ L . The resulting expressions are conveniently represented as diagrammatic expansions in the number of loops. This method is known as Eulerian standard perturbation theory (SPT) [2]. Typically, the linear statistics are assumed to be Gaussian and adiabatic, i.e. they are fully characterized by the power spectrum (2.7). Recall that the initial conditions for structure formation are set shortly after recombination when baryons and dark matter started to behave as a single pressureless fluid, but the gravitational instability did not yet have enough time to form non-linear structures. For simplicity, we neglect non-Gaussianity generated by the second order effects at radiation domination and recombination [45][46][47]. Suppose now that in addition to the power spectrum the linear distribution of the density contrast is also characterized by a non-trivial 3-point function, This generates additional terms in the SPT perturbative expansion, which originate from averages of odd number of fields. For instance, the leading-order (tree-level) matter bis-pectrum and the 1-loop correction to the matter power spectrum read, Note that the primordial non-Gaussian contributions have one less power of the growth factor D(η) compared to the terms coming from the non-linear evolution. Let us discuss an important point. The kernels α and β entering the Euler equations have poles if one of their arguments vanishes, see Eqs. (2.16). These IR singularities are inherited by the non-linear kernels F n , G n ; for example, This leads, in turn, to IR poles in the integrands of the loop expressions. The origin of the IR singularities can be traced back to large displacements of short-wavelength density perturbations by long-wavelength bulk flows. This effect should, however, cancel in the correlation functions of fields taken at the same moment of time. Indeed, the cancellation of IR singularities in equal-time correlators is a well-known property in Gaussian SPT, see [48] and references therein. For non-Gaussian initial conditions, the cancellation of contributions due to the IR poles can be tracked explicitly at the lowest orders of the perturbation theory; however, we are not aware of a general proof in the literature. Clearly, it is desirable to have a framework where spurious IR singularities are absent altogether. Such framework is provided by the time-sliced perturbation theory [35] described in the subsequent sections.
It is important to stress that the IR poles discussed above should be distinguished from those that may be present in the primordial bispectrum itself and thus are physical. Throughout the paper 'IR poles' or 'IR singularities' will only refer to the singularities originated from the non-linear mode coupling.

Review of Gaussian TSPT
The main idea of the TSPT approach is to substitute the time evolution of the overdensity and velocity divergence fields, δ and Θ, by that of the their time dependent probability distribution functional (PDF). This idea is natural when one is only interested in equal time correlation functions. For adiabatic initial conditions only one of the two fields is statistically independent. It is convenient to choose the velocity divergence field Θ as an independent variable and its PDF will be denoted by P[Θ; η]. At any moment in time, the field δ can be expressed in terms of Θ as Note that, in contrast to the SPT formulae (2.17), we have here on the r.h.s. the full non-linear field Θ(q). Recursion relations for the kernels K n are given in Appendix A. Equation (3.1) can be used to eliminate the density field from Eq. (2.15b) and obtain the following equation for the velocity divergence, with I 1 ≡ 1 corresponding to the adiabatic mode in the perfect fluid approximation. The other kernels I n can be derived from the fluid equations (2.15), see Appendix A. In the EdS approximation both kernels K n and I n are time independent. One introduces the generating functional, Equal-time correlation functions for Θ and δ are obtained by taking functional derivatives of Z with respect to the external sources J or J δ , respectively. For example, the matter power spectrum is given by The conservation of probability implies the Liouville equation for the probability density functional, It is useful to expand ln P[Θ; η] as a power series in Θ, where N is a normalization factor. Then, using Eq. (3.2) we obtain the following chain of equations for the vertices Γ tot n , where the sum in the second term on the l.h.s. is taken over all permutations σ of n indices. It is convenient to decompose the solution of this equation into two pieces, where Γ n is the solution of Eq. (3.7) with vanishing r.h.s. and with the initial conditions reflecting the initial statistical distribution, whereas C n is the solution of the inhomogeneous equation with vanishing initial conditions. The vertices Γ n thus have a physical meaning of 1-particle irreducible contributions to the tree-level correlators with amputated external propagators; the vanishing of the average velocity dispersion, Θ(k) = 0, implies Γ 1 = 0. On the other hand, C n are counterterms that cancel certain ultraviolet divergences appearing in the loop corrections [35]. The vertices Γ n and C n satisfy a hierarchy of equations which replaces the dynamical equations of SPT. From now on we specify to the EdS approximation where the kernels I n are constant in time. In order to find Γ n , we use an Ansatz which separates time and momentum dependence, This leads to the following recursion relations for Γ (l) n with l < n, (3.10) Note that the r.h.s. contains only vertices Γ (l) n with n less than n, implying that all Γ (l) n , l < n, are uniquely determined from vertices of lower orders. The remaining function Γ (n) n must be fixed from the initial conditions. Not to overload the paper with unnecessary formulae, we do not reproduce here the equations for the counterterms C n which can be found in [35]. The upshot is that C n are completely fixed in terms of the kernels I n and are time independent.
To sum up, the PDF (3.6) is fully specified by providing the 'diagonal' vertices Γ (n) n , n ≥ 2. These define the early-time asymptotics of the full vertex functions, where we have formally sent the initial time η i to −∞.
For the Gaussian initial conditions the solution to (3.10) is greatly simplified. In that case the PDF P[Θ; η] should reduce to a Gaussian weight at early times, upon rescaling the field with the linear growth factor, Θ →Θ = Θ/D(η). In other words, one requires, where P L (k) is the linear power spectrum. This implies the following initial conditions for the vertices, Then all Γ (l) n with l > 2 vanish, leaving the solution, where we have introduced the notationsΓ n ≡ Γ (2) n and Remarkably, in the case of Gaussian initial conditions and the EdS background, all Γ n vertices have universal dependence on time through the factor g −2 (η). As will be discussed shortly, g(η) plays the role of expansion parameter in TSPT. Due to momentum conservation, the vertices are proportional to a δ-function of the sum of their arguments. In what follows we use primes to denote the quantities stripped off such δ-functions, Note that the recursion relations (3.14b) and the initial conditions (3.13a) imply that all Γ n vertices are functionals of the linear power spectrum P L (k). The TSPT perturbative series is obtained upon expanding the generating functional (3.3) over the Gaussian part of P. Following the standard rules of quantum field theory, this expansion can be represented as a sum of Feynman diagrams. These are built from vertices corresponding to Γ n , n ≥ 3, and lines corresponding to the 'propagator' g 2 P L (k), see Fig. 1. One should also include vertices corresponding to counterterms C n , n ≥ 1, in order to subtract certain ultraviolet divergences of loop diagrams. To compute an n-point correlation function of the velocity divergence Θ one needs to draw all diagrams with n external legs. It is straightforward to see that diagrams with higher number of loops are proportional to higher powers of g(η). Hence, g(η) should be interpreted as the coupling constant of the theory. For the correlators of the density field δ one uses the expression (3.1) which is reminiscent of expressions for composite operators in quantum field theory. It gives rise to additional vertices proportional to the kernels K n ; these are denoted by an external arrow.

TSPT with primordial non-Gaussianity and Feynman rules
In the presence of primordial non-Gaussianity the initial conditions (3.13) get modified.
One expects now a non-vanishing initial value of the 3-point vertex, Γ 3 = 0. The earlytime asymptotics of the PDF become, (4.1) Substituting this expression into the generating functional (3.3) and taking variational derivatives with respect to the external sources one derives the correlation functions of the Θ-field at early times. The latter are to be matched to the linear statistics. Assuming that the non-Gaussian contribution is small, as will be confirmed shortly, we can restrict to linear order in the cubic vertex Γ 3 and obtain the matching condition, We observe that Γ 3 is proportional to the linear bispectrum. Let us estimate the size of the cubic term in (4.1). Taking for the characteristic amplitude of the modesΘ k ∼ P L (k) we obtain Γ As discussed before, the latter quantity is of order f N L A ζ which is much smaller than unity, given the existing bounds on non-Gaussianity. This justifies our expansion to linear order in Γ 3 . The TSPT vertices receive a contribution seeded by the primordial non-Gaussianity, Similarly to the Gaussian verticesΓ n , the non-Gaussian termsΓ N G n can be systematically constructed by induction. An explicit expression forΓ N G 4 is given in Appendix B. A comment is in order. The non-Gaussian vertices are proportional to an extra inverse power of the coupling g compared to the Gaussian part. This may give a wrong impression that they are enhanced. Actually, the contrary is true: they are suppressed due to the small amplitude of the primordial non-Gaussianity. Indeed, let us estimate the ratio between the non-Gaussian and Gaussian vertices. As usual, we neglect the momentum dependence and focus only on the overall amplitude. We have, where we have used thatΓ N G n andΓ n are proportional toΓ N G 3 andΓ 2 given by Eqs. (4.4a), (3.13a) respectively. In the last expression we also used the estimate (2.10). The ratio g(η) √ P L /A ζ is nothing, but the transfer function T (η), which is much bigger than unity at all relevant times and momenta. So, unless f N L is very large, the r.h.s. in (4.5) is much smaller than one. This suggests that the primordial non-Gaussianity should be treated perturbatively, on top of the standard vertices generated by non-linear clustering. It is convenient to consider f N L as a bookkeeping parameter which in all expressions is associated to the suppression (4.5). In other words, in the case of primordial non-Gaussianity we should also make an expansion in f N L on top of the TSPT expansion in powers of g(η). For practical purposes, it is enough to keep terms up to linear order in f N L .
We now discuss the diagrammatic technique. In principle, one could just modify the explicit structure of vertices Γ n and keep the same Feynman rules as in the Gaussian case. The propagator is still determined by the Gaussian part and is given by the linear power spectrum. The other building blocks of TSPT, such as the counterterms C n and kernels K n , are defined by the dynamical equations of motion (I n kernels) and thus are not modified when including non-Gaussianity. However, since the non-Gaussian vertices are suppressed as (4.5) and have a different scaling with the TSPT coupling constant, it is useful to single them out by introducing new graphical notations. We will denote them with an open star, e.g. for the non-Gaussian 3-point vertex we will use: Each non-Gaussian vertex contains one power of f N L , thus in the diagrammatic expansion we should keep only the graphs containing at most one non-Gaussian element. For instance, the Feynman diagrams contributing to the 1-loop non-Gaussian matter power spectrum are given by One can show that, similarly to the Gaussian case, the TSPT result for non-Gaussian contributions at linear order in f N L and a given order in g coincides with the result of SPT at the corresponding loop order. We verify this explicitly for the one-loop corrections to the power spectrum (4.7) in Appendix B; verification for the one-loop bispectrum was performed in [49]. At first glance, the new diagrammatic technique may seem somewhat excessive, as the calculation of the non-Gaussian contribution to the power spectrum in TSPT involves three graphs instead of only one in SPT. However, the new representation has an important advantage: as will be shown in the next section, all TSPT graphs are explicitly IR safe. This property allows, in particular, for a resummation of physical IR contributions, which are enhanced if the initial statistics have oscillating features.
It is straightforward to generalize the non-Gaussian TSPT formalism to include the redshift space distortions along the lines of Ref. [37]. To this end, one describes the mapping from the real to redshift space by means of a fictitious free motion of matter particles along the z-direction. This evolution proceeds in fictitious 'redshift time' F that varies from 0 to f , where f is the logarithmic growth rate defined in (2.14) (not to be confused with f N L ). The fictitious flow transforms the PDF of the cosmological fields. Applying the TSPT technique to this flow produces a set of equations for the redshift space vertices, which we will distinguish by the superscript '(s)'. These equations are linear and therefore do not mix the Gaussian and non-Gaussian parts of the PDF. Thus, we can directly use the formulas from [37] and write, ∂ ∂FΓ where is the redshift space kernel and the notationǩ i means that the momentum k i is absent from the arguments of the corresponding function. These equations should be integrated from F = 0 with the 'initial conditions' set by the TSPT vertices in real space, The final redshift space vertices are read off at F = f and have the form of polynomial expressions in the logarithmic growth rate f ; explicitly, the lowest vertices are, The expressions for the other building blocks of TSPT in redshift space, such as the Gaussian verticesΓ n , are not modified by the primordial non-Gaussianity and can be found in [37].
Before closing this section, let us briefly discuss the case of an arbitrary primordial non-Gaussianity encoded by 1-particle irreducible primordial l-point functions, (4.12) These give rise to linear l-point functions, The TSPT vertices are then found from the recursion relations (3.10) where now the initial conditions for the diagonal elements Γ  . (4.14) The full TSPT vertices can be written in the form, where the l-th term in the sum is generated from (4.14) through the chain equations (3.10).
Note that a primordial l-point function contributes only into the TSPT vertices with n ≥ l.
Similarly to the bispectrum case, one can introduce new diagrammatic notations for the non-Gaussian vertices with different l, whose magnitude will be controlled by where f (l) N L stands for the overall normalization of the primordial l-point function.

IR safety
The TSPT framework is naturally adapted to produce the equal-time correlation functions of the cosmological fields. The latter are protected from IR singularities by the equivalence principle [50][51][52][53][54]. Thus, one expects the TSPT formalism to be manifestly IR safe. Indeed, it was shown in Ref. [35] that in the Gaussian TSPT the kernels K n , counterterms C n and verticesΓ n are free from IR singularities. We now extend the proof to the non-Gaussian case. As K n , C n ,Γ n are not affected by the presence of non-Gaussianity, we need to analyze only the new non-Gaussian verticesΓ N G n . We are interested only in the IR singularities that may appear due to the dynamics of gravitational clustering and do not consider the poles that can be present in the initial bispectrum itself 5 . Thus, for the sake of the argument, we will assume that the linear bispectrum is free from singularities. Then the seed non-Gaussian vertexΓ N G 3 is IR safe, see Eq. (4.4a). In order to show IR safety for an arbitrary vertex, we proceed by induction. Consider a non-Gaussian (n − 1)-point vertex whose momentum arguments can be either soft (denoted by q) or hard (denoted by k) satisfying q k .
Suppose thatΓ N G n−1 is IR safe when any subset of its arguments uniformly go to zero, Consider now the recursion relation (4.4b) determining the vertexΓ N G n . According to the formulas from the Appendix A, the kernels I n with n ≥ 3 are proportional to K n and hence are IR safe. Further, the IR singularity of I 2 comes entirely from the kernel β defined in (2.16). Therefore, to prove the IR safety ofΓ N G n it is sufficient to focus on the part of Eq. The sum in the square brackets is Due to momentum conservation the sum of all the arguments entering the vertexΓ N G n is zero, which implies We conclude that the poles in the second line of (5.3) cancel, and the vertexΓ N G n is IR safe.

IR resummation of primordial oscillatory features
A well-known effect in Gaussian cosmological perturbation theory is the suppression of BAO due to non-linear mode coupling. Physically, it originates from large bulk flows that tend to smear short-scale features in the correlation functions. TSPT provides an accurate description of this effect via a systematic resummation of IR enhanced contributions both in real [36] and redshift [37] spaces. In this section we show that similar resummation technique can be developed within TSPT for the case of oscillating features in primordial statistics.
As a concrete example we consider resonant power spectrum and bispectrum produced in the axion monodromy inflation [18,31], where k t ≡ k 1 + k 2 + k 3 is the sum of wavenumbers, ϕ * is a constant phase whose value is not relevant for us, and The parameter γ is related to the inflationary slow-roll parameter 1 = −Ḣ inf /H 2 inf , the Planck mass M P and the axion decay constant f a as γ = √ 2 1 M P /f a 1; it is typically large, γ 1. On the other hand, the second parameter of the model δn s is assumed to be small. Planck constraints on oscillating features in the power spectrum [22] set an approximate bound [31] δn s 2 · 10 −3 γ 0.63 in the range 5 γ 10 3 . An interesting feature of the resonant non-Gaussianity is that its amplitude can be naturally large even when the corrections to the power spectrum are small: the value of f res N L as high as 2 · 10 4 is allowed for γ ∼ 10 3 . Note that though the second term in (6.1b) is suppressed by 1/γ, it is important to recover the correct squeezed limit of the bispectrum [18], so we keep it in our analysis.
Linear matter power spectrum and bispectrum are obtained by multiplying (6.1) with the transfer function, see Eqs. (2.8), (2.9). This superimposes on top of the primordial oscillations the standard BAO. Due the smallness of the BAO component in the transfer function, its consequences can be studied separately from the primordial features. The BAO contribution to various cosmological correlation functions, including damping by large bulk flows in Gaussian perturbation theory, was analyzed using TSPT approach in [36,37]. In principle, in the non-Gaussian case BAO will be imprinted in all linear statistics. For example, the linear bispectrum will contain BAO wiggles even if the primordial bispectrum is smooth. Similarly to the Gaussian case, the amplitude of these wiggles will be damped by large bulk flows. For completeness, we study this effect in Appendix D. In the main text we focus on oscillations that are already present in the primordial statistics. We work in real space, leaving the generalization to redshift space for future.

Power spectrum
To warm up, we first neglect the bispectrum and study the effect of IR modes on the primordial oscillations in the power spectrum. The analysis closely parallels the BAO resummation in Ref. [36]. We decompose the linear power spectrum into oscillatory ('wiggly') and smooth ('non-wiggly') parts, where, according to (6.1a), The Gaussian verticesΓ n are functionals of the linear power spectrum. By expanding them to the first order in P w L we separate them into non-wiggly and wiggly parts, Γ n (k 1 , ..., k n ) =Γ nw n (k 1 , ..., k n ) +Γ w n (k 1 , ..., k n ) . (6.5) As in Sec. 5, let us split the momenta into hard {k i } and soft {q j } and analyze the structure of the wiggly vertices in the limit ε ∼ q/k 1. Consider first the cubic vertex. The leading contribution in the soft limit is Consistently with the results of Sec. 5, the pole of the first factor cancels with the numerator of the second factor at q → 0. However, the cancellation does not occur if q is bigger than the period of oscillations, q k/γ, in which regime the contribution is enhanced by 1/ε. To keep the track of this enhancement, we introduce the linear operator D q acting on the wiggly power spectrum, Each insertion of this operator will be treated as a quantity of order 1/ε. Note that this operator is exactly the same as in the case of BAO resummation [36]. Next, one can show that the leading contribution in an n-point vertex with n − 2 soft momenta is of order (1/ε) n−2 and has the form, The proof of this formula is given in Appendix C.
We are now ready to identify the most IR enhanced diagrams among the loop corrections to the wiggly power spectrum. These diagrams should contain vertices with the maximal number of legs carrying soft momenta. They correspond to daisy graphs, obtained by dressing the wiggly vertices with soft loops ('petals'). Thus, in the leading IR order we have, where we denoted by a wavy line and shaded circles the wiggly linear power spectrum and the wiggly vertices respectively. The smooth lines represent the non-wiggly propagator g 2 P nw L (k). The term with loops in this expression is of the order g 2 (g 2 /ε 2 ) . We see that the loop suppression represented by g 2 is partially compensated by the IR enhancement, so these contributions must be resummed. The leading part of the diagram with loops reads, where in passing to the second line we used the formula (6.8). We have also restricted the loop integrals to the IR domain q ≤ k S . The separation scale k S defining this domain must belong to the range k/γ < k S < k. Otherwise its choice is arbitrary and represents an intrinsic freedom in the resummation scheme. We take k S to be a constant multiple of the momentum k, k S = κ · k , (6.11) and vary the proportionality coefficient κ between 0.3 and 0.7. By choosing κ rather close to 1 we expect to take into account all relevant IR modes. The sensitivity of the final result to the precise value of κ provides an estimate of the theoretical uncertainty which is due to the fact that one considers only the IR part of the loops and drops the integrals over the hard momenta q > k S . This sensitivity will decrease when one includes in the calculation higher orders in the hard loops (see below). Adding up the contributions (6.10) with all we obtain the leading-order IR-resummed wiggly power spectrum in a closed form, where we have introduced a new operator, The evaluation of its action at leading order in γ and ε gives (see Appendix C for details), where j 0 , j 2 are the spherical Bessel functions. Substituting this into (6.12) we see that the wiggly power spectrum is damped. The situation is very similar to the case of BAO [36,55], with the only difference that the dependence of the damping factor Σ P on k is now more complicated than quadratic. We show this dependence in the left panel of Fig. 2; the right panel illustrates the suppression of the wiggly power spectrum. The total IR-resummed power spectrum at the leading order is obtained by adding the non-wiggly contribution that remains unmodified, The resummation procedure can be extended to include the corrections from the hard loops. We do not repeat the derivation here referring the interested reader to [36]. At the next-to-leading order in the hard loops the result reads, where the last term is the usual one-loop correction to the power spectrum evaluated using the linear spectrum with damped oscillations. Note that an extra contribution g 2 e −g 2 Σ P P w L in the first brackets is precisely what is needed to avoid double counting the soft part of the one-loop daisy diagram.

Bispectrum
We now include into consideration the primordial bispectrum. We will assume that it is purely oscillatory, as in our working example of the axion monodromy inflation (6.1b). In a more general case the derivation below would apply to the wiggly part of the bispectrum. In this subsection we neglect the oscillations in the primordial power spectrum and focus on the features coming from primordial non-Gaussianity.

Generic momenta
We first consider the case when all three momenta in the bispectrum are of the same order; the squeezed case with one soft momentum will be studied afterwards. We start by writing the IR enhanced terms in the quartic non-Gaussian vertex with three hard and one soft momentum. Using the relations (4.4) and the form of the I 2 kernel (A.1a) we obtain, where we temporarily adopted the representation of the linear bispectrum as a function of two wavevectors (more precisely, of their lengths and the relative angle). Though not symmetric with respect to permutations in the three-point correlator, this representation turns out to be convenient for the derivation of the resummation formula. Equation (6.17) suggests to extend the action of operator D q introduced in Sec. 6.1 to the linear bispectrum, Clearly, this operator is of order 1/ε in the IR power counting. By induction one can prove an analog of Eq. (6.8) for the leading IR-enhanced part of the non-Gaussian vertices (see Appendix C), With this result at hand, it is straightforward to identify and resum the leading IR corrections to the non-Gaussian part of the bispectrum. Repeating the reasoning of Sec. 6.1 we obtain, Evaluating the diagrams and summing them up we arrive at, where the operator S is still given by Eq. (6.13), with D q acting on the bispectrum as in (6.18). Note that until now we have not used any specific expression of the linear bispectrum, so the result of IR resummation in the operator form (6.21) is valid for any oscillating non-Gaussianity. Evaluation of the action of S leads to the damping of the oscillations in the bispectrum, where we have switched back to the symmetric representation of B as the function of three wavenumbers. The concrete expression for the damping factor depends on the model. In Appendix C we perform the calculation for the case of the resonant bispectrum (6.1b) with the result,

(6.23)
Herek i is a unit vector directed along the momentum k i . One observes that, apart from the dependence on the absolute values of the momenta, the damping has a non-trivial dependence on the shape of the momentum triangle. This is illustrated in Fig. 3 where we show the dependence of Σ B on the momentum for several shapes. The IR separation scale is varied in the range 0.15k t ≤ k S ≤ 0.35k t . Note that with this choice Σ B reduces to the power spectrum damping Σ P (k 1 ) in the squeezed limit k 3 k 1 ∼ k 2 (more on the squeezed limit below). We see that the dependence on k S is rather mild. As noted before, it provides an estimate of the uncertainty due to the neglect of the higher loop corrections. To illustrate the effect of damping, we split the bispectrum into the envelope 6 B env and In Fig. 4 we show the effect of IR resummation on b osc for the value of the frequency parameter γ = 30. We observe that the oscillations are significantly damped at k 0.07 h/Mpc and get essentially washed out at k 0.15 h/Mpc.

The squeezed limit
So far we have been assuming that all three arguments of the bispectrum are of the same order, k 1 ∼ k 2 ∼ k 3 . This assumption does not cover the squeezed limit when one of the momenta is soft. The latter limit is relevant for the scale-dependent halo bias. Besides, in the cosmological collider models of [19][20][21] this is precisely the limit that encodes the information about the masses and spins of heavy particles. Thus, it deserves a separate attention. We are going to see that it leads to new IR-enhanced contributions. Consider again the non-Gaussian four-point vertex where now only two momenta are hard. A brief calculation gives its leading IR-enhanced part, 25) where the operator D q acts only on the hard momentum in the linear bispectrum, In Appendix C we generalize this formula to an arbitrary n-point vertex, We observe that the leading IR term is of order 1/ε n−3 , as in (6.19), but its structure is different. Therefore, the leading contributions to the squeezed bispectrum are still given by the diagrams (6.20), which, however, need to be re-evaluated. Denoting the soft external momentum by q 0 we obtain for the -loop term in the series, (6.28) Summation over yields the leading order squeezed bispectrum in the operator form, (6.29) The first term here is the same as Eq. (6.22). To see this we evaluate the action of the operator S, where Σ P is the power-spectrum damping factor (6.14). As discussed before, the latter coincides with the squeezed limit of the bispectrum damping Σ B . The second contribution in (6.29) is new. It describes the distortion of the oscillatory pattern produced by the long mode q 0 . Of course, this effect is non-vanishing only if q 0 k/γ. In real space this corresponds to the requirement that the length of the long mode should be still shorter than the correlation length l corr ∼ γ/k set up by the momentumspace oscillations. Spelling out the operators in (6.29) we obtain an explicit expression, where Note that we have shifted the hard momentum k → k − q 0 2 to somewhat simplify the subsequent formulas.
Evaluation of the resonant bispectrum in the squeezed limit is subtle. Taken at face value, the expression (6.1b) will be dominated at q k/γ by the second term in the square brackets. The latter can be written as 7 which is precisely of the form dictated by the Maldacena consistency condition [56]. It is known that such terms drop from physical observables. In particular, they cancel out if one transforms the bispectrum from the comoving coordinates used to derive (6.1b) to conformal Fermi coordinates (CFC) [57]. In the latter frame the leading contribution into squeezed primordial bispectrum of a single-field slow-roll inflation scales as O(q 2 /k 2 ). Thus, to avoid unphysical contributions, the term (6.33) must be subtracted. On the other hand, the construction of the CFC frame is possible only if the wavelength of the long mode is larger than all other relevant scales. It breaks down when the length of the long mode becomes comparable to the correlation length 8 l corr ∼ γ/k. In other words, while we can easily obtain an unambiguous physical squeezed limit of the bispectrum at q k/γ using CFC, it is no longer possible at k/γ q k. A proper definition of the physical bispectrum in the intermediate range presents an open issue and may require taking into account the projection effects pertaining to the actual physical quantities measured by an observer on Earth. Such study is beyond the scope of this paper. Instead, we adopt a simplified approach and model the expected form of the physical bispectrum by subtracting from (6.1b) the term (6.33) weighted with a function w γq 2k that interpolates between w = 1 + O(x 2 ) at x → 0 and w = 0 at x → ∞; in the practical calculations we will use w = (1 + x 2 ) −1 . This model cannot be judged fully satisfactory. In particular, it likely introduces an order-one error compared to the true bispectrum at γq 2k ∼ 1. Nevertheless, it is sufficient for our purposes of illustrating the dynamical effects of large bulk flows.
With the above caveat in mind, we will use the following form of the squeezed linear bispectrum, where the smooth enveloping function B env is defined in (6.24). This form directly determines the contribution B (1) of the IR resummed bispectrum, which differs from the linear bispectrum only by the damping of oscillations. To find the second contribution B (2) one substitutes (6.34) into (6.32b). Integration over the directions of q yields, 8 Indeed, it can be shown using the results of [58] that the terms introduced into the squeezed bispectrum by the transformation from the comoving frame to CFC are organized as an expansion in powers of (qlcorr) 2 . with Note the peculiar dependence of the expression (6.35) on the angle between the soft and hard momenta: for q 0 k/γ it rapidly oscillates as a function of the angle, whereas at q 0 k/γ it reduces to the sum of monopole and quadrupole. We plot the contributions B (1) , B (2) (for collinear k and q 0 ) as functions of the hard momentum k in Fig. 5. We observe the damping of oscillations in B (1) above k ∼ 0.15h/Mpc, though the effect is somewhat weaker than in the equilateral and flattened cases. The behavior of B (2) is quite different: its amplitude grows from small k up to k ∼ 0.2h/Mpc and only then starts decaying exponentially. In Fig. 6 we plot the components of the squeezed bispectrum as functions of the soft momentum at fixed k.
Finally, let us notice that the squeezed limit of the resonant bispectrum, through the mode coupling, produces an additional oscillating contribution in the power spectrum on top of the primordial one (6.4). Its leading part can be found using again the method of IR resummation. It is given by the sum of daisy diagrams, 38) where the coefficients C 1 , C 2 are given in (6.36). Given the measured value of A ζ , this contribution is always small. Further, the relation (6.2) implies that it is subdominant compared to the primordial wiggly power spectrum as long as γ < 100. Thus, in our working example of the primordial non-Gaussianity produced in the axion monodromy inflation this contribution appears uninteresting. However, it may be relevant in other models where the oscillating component in the primordial power spectrum is absent or strongly suppressed. This situation occurs, for example, in the context of cosmological colliders [19][20][21].

Conclusion
In this paper we worked out the TSPT framework for non-Gaussian primordial statistics. They are included by setting appropriate initial conditions for the TSPT statistical weight, which eventually modifies only the vertices Γ n . Other building blocks of TSPT, such as the composite operator kernels K n and counterterms C n remain intact. In this way TSPT can easily accommodate any sort of non-Gaussian initial conditions without a significant modification of the formalism. Treating initial non-Gaussian contributions as small deformations of the Gaussian statistical weight generates a new perturbative expansion on top of the usual TSPT perturbation series. We have shown that the new non-Gaussian TSPT vertices are manifestly IR safe, implying IR safety of the non-Gaussian TSPT Feynman diagrams.
TSPT with non-Gaussian initial conditions opens a new way to study the effects of primordial non-Gaussianity on LSS clustering. As an application, we developed a new procedure to capture the non-linear evolution of primordial oscillating features that may be imprinted in the power spectrum and higher statistics during inflation. We showed that such oscillations are suppressed by nonlinear effects of bulk motions, just like BAO. We demonstrated how to calculate the damping factor by resummation of IR enhanced diagrams and computed it explicitly for the resonant power spectrum and bispectrum from axion monodromy inflation. We point out that the damping of primordial oscillating features due to large bulk flows has been overlooked in many studies of future survey sensitivities. Taking it into account can significantly affect the forecasts. It will also be interesting to study to what extent the damping can be removed using the density field reconstruction technique [59,60].
We have found that in the case of squeezed bispectrum the IR resummation, apart from damping, produces an additional contribution depending on the angle between the hard and soft momenta. It will be interesting to study observable signatures of this term, in particular, to what extent it can be discriminated from the primordial angular dependence of the squeezed bispectrum in the cosmological collider models.
The results of this work reinforce the status of IR resummation in TSPT as a universal systematic procedure, applicable beyond the case of Gaussian initial conditions. They can be extended in several directions. While in the current paper we focused on the power spectrum and bispectrum, the IR resummation can be straightforwardly applied to higher primordial statistics, in particular, to the trispectrum which may possess oscillating features due to an exchange of massive particles [19][20][21]. Next, it would be interesting to incorporate the effects of redshift space distortions in the IR-resummed non-Gaussian statistics; this can be done along the lines of [37]. One would also like to extend the non-Gaussian IR resummation beyond the leading order to include the contributions of loop diagrams with hard momenta and subleading IR corrections, as done for the Gaussian case in Ref. [36].
Last but not least, we adopted in this paper the pressureless perfect fluid approximation. To properly renormalize the contributions of short scales, the non-Gaussian TSPT should be embedded in the framework of EFT of LSS [8][9][10]. We leave this task for future.
Note added. When this paper was completed, we became aware of the work [61] that also discusses the damping of the primordial oscillation features in the power spectrum. When overlap, our results agree. We thank the authors of [61] for sharing with us the preliminary version of the paper.
(4.4b) and read, Now we are ready to write the expressions corresponding to the three diagrams in (4.7). For the first diagram, which we call 'daisy', we obtain, .

(B.5)
Collecting the three contributions together we arrive at The kernels I 2 , K 2 are to be taken from Eqs. (A.1a), (A.2a). Using the relation between the TSPT and SPT kernels established in [35], ..,q i , ...,q j ..., q n−2 ) + . . . , where we have denoted by q tot the sum of all soft momenta. Dots in the last line stand for the terms involving the verticesΓ w n with n ≤ n − 2 and kernels I m with m ≥ 3. All the latter kernels are IR safe, hence the omitted terms are at most of order 1/ε n−4 . At the leading IR order they can be neglected. The kernel I 2 has a singularity only if one of its momenta is hard and the other is soft. Otherwise it is of order one. This implies that the terms in the first and last lines of (C.1) are at most of order 1/ε n−3 and can be neglected as well. One is left with the contributions of the second and third lines. The next step is to use the IR asymptotics and the expression (6.8) forΓ w n−1 . This yields, Γ w n (k, −k − q tot , q 1 , ..., q n−2 ) The terms in the square brackets combine to produce the operator D q i . Using the commutativity of operators D q with different q's one arrives at Eq. (6.8).
We now turn to the proof of Eq. (6.19). We again assume that it holds for n ≤ n − 1 and consider the relation (4.4b) defining the n-th vertex. Keeping only the terms containing the IR-enhanced kernels I 2 we obtain, where dots correspond to the terms that are at most of order 1/ε n−4 . Using next Eq. (C.2) and substituting the leading IR part ofΓ N G n−1 we get, Recalling that k 12 = k 1 + k 2 we recognize in the square brackets the operator D q i which completes the product of operators to n−3 j=1 D q j . Summation over i then cancels the factor (n − 3) in the denominator and we obtain (6.19).
Finally, we give the proof of Eq. (6.27). The strategy is the same as in the two previous cases: we assume that it is valid for n ≤ n − 1 and study the vertex Γ N G n . The leading contribution containing the IR-enhanced I 2 kernels is ..q n−2 ) + . . . .

(C.6)
Substituting the IR asymptotics of I 2 andΓ N G n−1 we obtain, The sum over i cancels the factor (n − 3) from the denominator which yields Eq. (6.27).

C.2 Damping factors
We first compute the action of the operator S given by Eq. (6.13) on the wiggly power spectrum (6.4). Using q · ∇ k e iγ ln k/k * = iγ (q · k) k 2 e iγ ln k/k * (C. 8) we have in the leading order in γ, Evaluating the integral over directions of the momentum q we arrive at Eq. (6.14).

D BAO in the non-Gaussian contribution to the bispectrum
Even if the primordial bispectrum is smooth, the linear one will acquire an oscillatory component due to the BAO wiggles in the transfer function, see Eq. (2.9). In this Appendix we study these wiggles using the frameworks of TSPT and discuss their damping by large bulk flows. Following [36], where such study was performed for Gaussian initial conditions, we separate the linear power spectrum into the non-wiggly and wiggly parts, where the latter encapsulates the BAO. It is characterized by oscillations in k with the period 2π/r BAO , where r BAO 110h −1 Mpc is the BAO scale. As the BAO component is small, we can expand all quantities to linear order in P w L . Thus, the linear bispectrum becomes, where B nw L is given by Eq. (2.9) with P L replaced by P nw L and Correspondingly, the TSPT vertices, both Gaussian and non-Gaussian, split into smooth and wiggly parts. In particular, for the cubic non-Gaussian vertex we have, .

(D.5)
Note that the sign in the last formula is different from (4.4a) because of contributions coming from the expansion of the power spectra in the denominator. Other wiggly non-Gaussian vertices are generated fromΓ N G, w 3 by the recursion relation (4.4b). We can now readily use the results of Sec. 6 to work out the effect of IR resummation. Thus, the IR enhanced contributions to wiggly bispectrum with three hard momenta are obtained by dressing the wiggly elements in the two tree-level diagrams + (D. 6) with daisy loops. Here the empty and shaded star stand for the non-wiggly and wiggly non-Gaussian vertices, respectively; recall also that the wiggly line denotes P w L , whereas the plain line stands for P nw L . We already know the result of this dressing from Eqs. (6.12), (6.21) that provide us with its operator form, It is straightforward to find the action of the operator S on P w L and B w L . The result amounts to an exponential damping of the wiggly power spectra that enter in (D.7) both explicitly and implicitly through the expression (D.3), P w L (k) → e −g 2 σ 2 k 2 P w L (k) . (D.8) Here is the same factor that appears in the context of BAO damping in the Gaussian perturbation theory, see e.g. [34,36,55]. The choice of the separation scale k S delimiting the region of soft modes parameterizes the freedom in the resummation scheme. It is known that good results are obtained for k S in the range (0.05 ÷ 0.2)h/Mpc (see a detailed discussion in Ref. [36]). Substitution of (D.3), (D.8) into (D.7) yields, By adding the non-wiggly part we can cast the result in a compact form, B N G, LO δδδ (η; k 1 , k 2 , k 3 ) = g 3 (η) P res L (η; k 1 )P res L (η; k 2 )P res L (η; k 3 ) P ζ (k 1 )P ζ (k 2 )P ζ (k 3 ) B ζ (k 1 , k 2 , k 3 ) , (D.11) where the resummed linear power spectrum is P res L (η; k) = P nw L (k) + e −g 2 (η)σ 2 k 2 P w L (k) .

(D.12)
This expression has a clear intuitive interpretation. It tells us that to obtain the leadingorder IR resummed bispectrum all we have to do is to take the usual expression, Eq. (2.9), and replace all linear matter power spectra in it with their resummed counterparts. One can show that the final formula (D.11) is also valid in the squeezed limit. An extra contribution coming from the daisy dressing of the non-Gaussian wiggly vertexΓ N G, w n and corresponding to the second term in Eq. (6.29) in fact cancels against the contribution of the fish diagram composed of the wiggly Gaussian and non-wiggly non-Gaussian vertices, We leave the proof of this cancellation to the reader. The effect of IR resummation on the BAO feature in the non-Gaussian bispectrum is illustrated in Fig. 7. We see a clear damping of the oscillations due to large bulk flows. Of course, even apart from the damping, the contribution of BAO in the linear bispectrum is at most several per cent. Given the current constraints on non-Gaussianity, measuring this feature appears unfeasible even with futuristic surveys.