A Cosmological Bootstrap for Resonant Non-Gaussianity

Recent progress has revealed a number of constraints that cosmological correlators and the closely related field-theoretic wavefunction must obey as a consequence of unitarity, locality, causality and the choice of initial state. When combined with symmetries, namely homogeneity, isotropy and scale invariance, these constraints enable one to compute large classes of simple observables, an approach known as (boostless) cosmological bootstrap. Here we show that it is possible to relax the restriction of scale invariance, if one retains a discrete scaling subgroup. We find an infinite class of solutions to the weaker bootstrap constraints and show that they reproduce and extend resonant non-Gaussianity, which arises in well-motivated models such as axion monodromy inflation. We find no evidence of the new non-Gaussian shapes in the Planck data. Intriguingly, our results can be re-interpreted as a deformation of the scale-invariant case to include a complex order of the total energy pole, or more evocatively interactions with a complex number of derivatives. We also discuss for the first time IR-divergent resonant contributions and highlight an inconsequential inconsistency in the previous literature.


Introduction
In the study of a physical system, much of the initial effort is devoted to charting the landscape of theories that are both consistent and compatible with the data.In this endeavor, we can either sample the bulk of theories by building concrete models, or probe the boundaries of theory space using fundamental principles of physics.These two approaches, which can be roughly labelled as model building and bootstrapping, complement each other in very powerful ways and particle physics is a spectacular success story of this synergy in action.Much recent work aspires to achieve similar progress in the study of the primordial universe and the origin of initial conditions for our universe [1].Both model building and bootstrapping heavily leverage symmetries to reduce the complexity of the analysis and to sharpen physical predictions.For models of the early universe, a major role has been played by scale invariance, which in the cosmological context1 refers to a fixed scaling behaviour of n-point primordial correlation functions of curvature perturbations, B n (λk) ∼ λ 3(1−n) B n (k).Scale invariance has been observed to be an approximate symmetry of primordial perturbations, with percent deviations detected in the so-called spectral tilt of the scalar power spectrum.While it is natural to focus on models that explain such a surprising and unexpected symmetry, such as slow-roll inflation, it is also interesting to explore possible deviations from it that are still allowed by the data.At the phenomenological level, a minimal possibility is to reduce continue scale invariance to a discrete subgroup, by allowing for small periodic oscillations of correlators with scale, which average out over a large range of scales (for early proposals see e.g.[2]).Remarkably, this possibility appears to be also natural from the model building point of view, where the same setup that ensures continue scale invariance in perturbation theory displays small periodic deviations from it due to non-perturbative effects [3][4][5][6][7][8].
The goal of this work is to study cosmological correlators that enjoy only discrete scale invariance, as in models of resonant non-Gaussianity [2,9].We perform this analysis by adapting the tools of the boostless cosmological bootstrap, recently developed in a series of papers [10][11][12][13][14], building upon previous work [15][16][17], to the case where full discrete scale invariance is not a symmetry.Our analysis and results can be motivated by several points of views.
Resonant non-Gaussianity At the phenomenological level, we derive an infinite class of "resonant" scalar bispectra, which we then compare to CMB data by Planck to obtain observational constraints.The leading resonant non-Gaussianity had already been studied in a concrete model in [9] and within the framework of the effective field theory of inflation [18] in [19] (see also [20][21][22][23][24] for other studies).Here we find infinitely many new shapes that correspond to cubic scalar operators with higher derivative interactions.In a technically natural effective theory, higher-dimension operators are generically expected to produce smaller signals.While the calculation can be performed exactly, the regime of fast oscillations is particularly interesting because the result simplifies considerably and the signal is expected to be relatively large.This regime is controlled by the dimensionless parameter α, namely the frequency ω of background oscillations that are responsible for the resonance in units of the Hubble parameter, α = ω/H.For the bispectra that had already been discussed in the literature, our derivation reproduces the same leading order term in α ≫ 1. However we find a different result at the next-to-leading order in 1/α.Within our derivation the discrepancy is unavoidable because the resonant bispectra published in the literature do not satisfy the so-called manifestly local test, a necessary condition for the non-Gaussianity generated by local interactions [12].Upon further inspection we identify the source of the discrepancy: in previous calculations the stationary phase approximation was employed, but subleading corrections in α had not been correctly propagated to the final result, hence yielding an incorrect prediction of the next to leading and subsequent terms.Amending the mistake in the literature is in practice inconsequential for phenomenology, since it amounts to a small shift in the phase of oscillations, which is a free parameter of the model.However, the situation is a good reminder of the usefulness of model-independent constraints, both for bootstrapping and for checking the consistency of predictions.

The bootstrap without scale invariance
There is a second independent motivation for our analysis.So far all cosmological bootstrap approaches have made crucial use of scale invariance.Indeed, the de Sitter invariant approach of [16,[25][26][27] relied on the full group of dS isometries and the same applies to results obtained from "momentum-Mellin" space [28,29] or Anti-de Sitter [30,31].In a second approach, motivated by the no-go theorems of [17], invariance under dS boosts is relaxed [10], and additional constraints from locality [12], unitarity [11,14,32,33], and causality/analyticity [16,[34][35][36][37] are brought into the game, while still assuming scale invariance (see e.g.[38][39][40][41][42][43][44]).Several hybrid approaches have leveraged the dS invariance of a seed signal and combined it with the breaking of boosts in a systematic way (see e.g.[45][46][47][48]).It is therefore natural to ask how far one can go when scale invariance is relaxed.Here we show that the bootstrap approach can be successful in predicting a class of new non-Gaussian signals, even when the constraint of scale invariance has been weakened to a discrete subset.The price to pay is that we will have to rely more heavily on the structure of the initial ansatz.At a hand-wavy level, the ansatz can be justified by invoking the "simplest possible analytic structure", however we don't attempt to formalise this in a precise mathematical sense.An interesting byproduct of our analysis is that we highlight interesting modifications of well-known bootstrap constraints.For example, the Cosmological Optical Theorem (COT) does apply to the case of resonant non-Gaussianity (and indeed to any FLRW spacetime with time-dependent couplings [14]).However this constraint relates leading terms to exponentially suppressed terms and is therefore not very useful unless one attempts to bootstrap the full result.Instead, we show that the COT can be modified by taking advantage of the analytic dependence of the result on a non-kinematical parameter, namely the frequency of oscillations α.This leads to a modified cosmological optical theorem that constrains directly the leading non-Gaussianity, which in practice allows us to completely neglect the inconsequential exponentially suppressed terms.
A complex deformation Finally, we come to the third and last motivation.At face value, cosmological correlators are functions of the kinematics and of parameters in the theory.The kinematical dependence has been explored extensively, while the parameter dependence appears to be quite trivial in perturbation theory: each contribution is simply proportional to the product of the coupling constants involved in the corresponding diagram.However, there is one parameter that appears non-trivially.To see this, consider a contact interaction of a massless scalars in dS.At tree-level one expects a contribution to the wavefunction coefficients of the form where {k} collectively denotes the n three-momenta k a and the so-called total energy is the sum of their norms, k t ≡ a |k a |.Assuming scale invariance, it was shown that p is a positive integer2 given by [10] where the sum runs over each interaction of a given diagram and ∆ A is the mass dimension of the interaction at vertex A. This result applies also to the order of partial energy singularities and assumes the absence of IR divergences.For example, for n = 3 we need a cubic interaction and so p simply counts the total number of space and time derivatives.Surprisingly, the interesting parameter for us here is p.Indeed, the results in this paper show that resonant non-Gaussianity can be roughly thought of as a continuation to complex p of the standard scale-invariant result (albeit with different constraints on the polynomial).Indeed, taking p → p + iα, we obtain non-Gaussianity of the form where the scale that makes the logarithm well-defined corresponds to a model-dependent shift in the phase of oscillations, which is marginalized over when comparing to observations.From the perspective of the bulk time integral representation, adding an imaginary part to p corresponds to a complex twist of the integral where ω = αH and exp(iα log η) ∼ cos(ωt) represents a background oscillation that is periodic in cosmological time t.This interpretation is quite close in spirit to the proposal of [15] of relating what ultimately are flat space objects to cosmological wavefunctions by integrating over appropriate energy shifts of the interaction vertices.Indeed, resonant non-Gaussianity is produced on sub-Hubble scales, around kη ∼ α ≫ 1, where the flat-space approximation becomes increasingly good.As a consequence, the full correlator is specified by a corresponding flat-space amplitude A n , to leading order in α, It's important to contrast this with the usual statement that the UV-limit of a flat-space amplitude appears as the residue of the leading k T pole of a wavefunction coefficient [52,53].In the usual case the leading k T → 0 pole is leading only for unphysical complex momenta, while for physical momenta many other contributions arise, namely all subleading k T poles, that are of the same order as the leading k T -pole and are not directly related to flat space.In contrast, for resonant non-Gaussianity we have a physical parameter, namely α, that we can dial to make the flat-space amplitude contribution leading for all3 physical momenta.

Summary of the results
For the convenience of the busy modern reader, we summarize here the take-home messages of this work: • We have developed a (boostless) cosmological bootstrap that does not rely on scale invariance but only on a weaker discrete scaling subgroup, determined by the condition (2.8).We find an infinite class of non-Gaussianities (cubic wavefunction coefficients and bispectra) that satisfy this relaxed scaling condition as well as all other known physical constraints imposed by unitarity in the form of the cosmological optical theorem, locality in the form of the manifestly local test, and the choice of initial state in the form of the sufficiently constraining ansatz in (2.14).The result, summarized to all orders in derivatives in (3.5), turns out to be a generalization of the well-known resonant non-Gaussianity [9], which arises in models of axion monodromy inflation.Our results are possibly not exhaustive of all possibilities with the chosen symmetry, but are expected to capture the leading signal.We comment in Section 4 on additional contributions from the resonant excitation of the mode functions, which we neglect in the main derivation.
• As discussed above, our results are close to the scale invariant case and can be thought of as a deformation of the order p of the total energy pole to include an imaginary part, which in turn leads to oscillations in the total energy k T , and breaks scale invariance down to a discrete subgroup, as in (1.3).
• For the first time we discuss IR divergent resonant contributions (see Section 3.6).These appear in the wavefunction coefficients as well as in the correlators of a spectator scalar.However they disappear in the bispectrum of curvature perturbations in single field inflation, as expected on general grounds.
• We searched for several resonant non-Gaussian bispectra in the temperature anisotropies of the cosmic microwave background, as measured by the Planck satellite using the pipeline CMB-BEST [55].We find no evidence of a signal, as discussed in Section 5.
• Even though the cosmological optical theorem (COT) applies to time dependent couplings, it relates leading terms, which we bootstrap, to exponentially suppressed terms, which we wish to neglect and so is not practically useful.However, we derive a modified cosmological optical theorem that can be applied to the case of an imaginary and time-dependent coupling constant e ±iωt .This is useful because the modified COT relates the leading term in the wavefunction coefficient to itself and hence determines the reality of its coefficients.The new relation leverages the simple analytic dependence of the result on a parameter of the model, namely the frequency of oscillations ω or its dimensionless cousin α = ω/H.
We end this introduction with a guide to the organization of the rest of this paper.In Section 2 we discuss a set of bootstrap rules and put forward a bootstrap ansatz, which will be subsequently used to find an infinite class of resonant non-Gaussianities.In particular, we discuss the constraint from discrete scale invariance, the form of the ansatz appropriate for resonant production of non-gaussianity from a Bunch-Davies initial state, the modified cosmological optical theorem that applies to the leading resonant contributions and finally the manifestly local test.In Section 3 we use these constraints to derive an infinite class of cubic wavefunction coefficients and present some examples to the lowest orders in derivatives.We also study IR divergences in Subsection 3.6 and parity-odd interactions in Subsection 3.5.Then, in Section 4, we discuss subleading contributions coming from the resonant excitation of negative frequency modes in the free theory and their impact on non-Gaussianity.In Section 5, we present the observational constraints on the leading resonant bispectra that we have derived from the Planck data on CMB temperature and polarization anisotropies.We conclude in Section 6 with a discussion and an outlook.

Bootstrap ansatz and rules
In this section, we present a set of bootstrap rules and a bootstrap ansatz.Then we show how they lead to a class of solutions for the cubic wavefunction coefficient of a massless scalar field in de Sitter spacetime.Our derivation parallels that in [10,12], but we relax the requirement of exact continuous scale invariance in favour of discrete scale invariance.
Anticipating invariance under translations, we parameterize the field-theoretic wavefunction as where to simplify the notation we have introduced and η is conformal time for the flat (Poincaré) patch of de Sitter spacetime, Throughout we assume that ϕ is a spectator massless scalar field in de Sitter and we allow for boostbreaking interactions.

Homogeneity and isotropy
An n-point wavefunction coefficient ψ n ({k}; η) is in general a function of the n wavenumbers k a with a = 1, . . ., n and of time η.We will focus on the late time behaviour of ψ 3 , as η → 0, namely at the future conformal boundary of dS.
In writing (2.1) we already assumed invariance under spatial translations, which leads to the momentum conserving delta function.We furthermore impose invariance under rotations, which implies that ψ n is a function of rotation-invariant contractions of the momenta with the Kronecker delta δ ij or the Levi-Civita symbol e ijk .Invariance under rotations and translations reduces the number of variables on which ψ n depends from 3n to 3(n − 2).Here we will focus on n = 3.In this case the Levi-Civita symbol cannot appear because all the momenta lie in the same plane.Hence ψ 3 is a function of just three variables, which can be chosen to be the norms k a ≡ |k a | of the momenta Let us discuss the constraints on ψ 3 (k 1 , k 2 , k 3 ) from a particular form of discrete scale invariance.

Discrete scale invariance
Usually one is interested in imposing continuous scale invariance, namely the relation This transformation is the manifestation of the dilation isometry of dS for massless fields at the future conformal boundary, η → 0. To see this, notice that a dilation rescales x → λx and η → λη, which leaves the line element in (2.3) invariant.For a massless scalar we expect that the leading time dependence for η → 0 is a constant, ϕ ∼ η 0 .For the late-time wavefunction this in turn leads to (2.5). 4ere we are interested in relaxing this condition.If we made no assumption whatsoever about scale invariance, the family of possible solutions would be too large to make any useful predictions.Instead, we will assume that the continuous scale invariance in (2.5) is relaxed to a specific discrete form of scale invariance defined as follows: where λ is some real and positive number which is part of the data that specifies the theory.Notice that this then implies and consequently we can take λ ∈ (1, ∞) without loss of generality.With this information, we finally formulate our discrete scale invariance constraint as This condition is inspired by a class of models in which periodic oscillations appear in the inflaton potential, where V 0 is some slow-roll potential, and Λ and f are energy scales.This structure is a rather generic feature in models of axion inflation, where the periodic features are induced by non-perturbative effects in a gauge sector to which the axion couples via a shift-symmetric dimension five coupling (sometimes called "Chern-Simons" or "theta-term").A partial collection of references is [3,5,6,[56][57][58].One can generalize this setup by considering the effective field theory of inflation (EFToI) [18,59].In the most general case, the Wilsonian coupling constants appearing in the EFToI are arbitrary functions of time.Then one assumes that when they are approximately constant one recovers exact scale-invariant predictions.To generalize to oscillations in the inflaton potential, one demands that the Wilsonian couplings are not constants but are invariant under specific discrete cosmic time translations, t → t + const [19,60].The discrete shifts in time lead to a discrete rescaling in k as follows.First notice that in de Sitter t ∝ ln(−η) and so a shift in t corresponds to a rescaling of η.Second, recall that when the Wilsonian coefficients are constant the theory must be invariant under dilations and so only the combination kη can appear.Then k ∼ η ∼ e −Ht and the discrete scale invariance in (2.8) follows.Here we bypass this time evolution picture and simply use the constraint in (2.8) as a definition of the class of models we are considering.
Before proceeding we should mention that here we are assuming that there are no IR divergences and so we can safely evaluate ψ 3 (η) as η → 0. In the presence of IR divergences, we have to introduce a late-time regulator and discrete scale invariance needs to be formulated in a slightly modified way.In this section we will continue assuming no IR divergences, and we postpone a discussion of those to Section 3.6.

The bootstrap ansatz
Here we want to argue that a minimal ansatz that satisfies the constraint from discrete scale invariance in (2.8) is given by where is the so-called total energy, k * is some constant pivot scale, which is degenerate with the normalization of ψ3 , and ψ3 is a rational function of the energies k a that scales as (2.12) Since the bootstrap procedure cannot fix it, we will omit the overall phase k * throughout, and will only restore it when presenting the final results in Section 3. If we assume manifestly local interactions and a Bunch-Davies initial state, the only singularity that ψ3 can have for complex energies k a is at k T = 0 and hence it can be written as where p is the maximum order of the total-energy singularity, which equals the maximum number of space plus time derivatives one is considering, and the numerator is a homogeneous polynomial of degree p + 3. Notice that combining (2.10) and (2.13) we find simply This is almost identical to the usual ansatz for the scale invariant bispectrum [10], except for the crucial difference that the "order of the k T pole" is allowed to have an imaginary part.At a more abstract level, this arises by introducing a complex twist into the usual bulk time integrals, which is closely related to the ongoing work of [61].We take (2.14) to be part of the data that defines the theory we are studying.
Again, we postpone a discussion of possible IR divergences to Section 3.6.
In the following we will justify this ansatz in two ways.First by arguing that this ansatz possesses the "simplest" possible analytic structure.Second, we will see that this ansatz is precisely what is expected from the picture of bulk time evolution.
The simplest analytic structure An imprecise but general lesson of the bootstrap approach to amplitudes and cosmological correlators [1] is that the growing complexity of bulk perturbation theory manifest itself in the increasingly complex analytic structure of the corresponding "boundary" quantities as function of the kinematical variables k a .For example, contact diagrams of massless fields have only total energy poles, exchange diagrams have also partial energy poles and loop diagrams display a collection of branch points (see e.g.[15] and [35] for more systematic discussions in Minkowski).Hence we would like to start looking for solutions of the discrete scale invariance constraint that have the simplest possible analytic structure.While it would be nice to do this in a mathematically rigorous way by formally quantifying the "complexity" of a given analytic structure, here we proceed heuristically and consider first entire functions, then functions with a finite or infinite number of poles and finally functions with a finite or infinite number of branch points.A slightly more detailed discussion is given in Appendix A. First notice that the only way to have a finite number of poles is to have poles at k = 0 and k = ∞.Poles at any other location must have infinitely many images because of discrete scale invariance.In this simplest case ψ n is forced to also satisfy continuous scale invariance.This case has already been discussed in the literature, e.g. in [10,12], and so we set it aside.The next simplest possibility is to have an infinite number of poles, which necessarily accumulate at k = 0 and k = ∞.We weren't able to find any concrete function that satisfies this property and so we move on to the next possibility-two branch points at k = 0 and k = ∞ connected by a single branch cut (which in some sense can be thought of as a limiting case of an infinite number of poles).The rich class of discrete scale invariant functions with such an analytic structure is given by e −iα log (k T /k * ) times a homogeneous function of degree three, which finally brings us to (2.10), the ansatz that we will consider in the rest of this paper.
The bulk time evolution Here we will the use the representation of the cosmological wavefunction in terms of bulk time integrals to confirm our choice of the ansatz in (2.14).The bulk integrals that are commonly encountered in resonant models at tree level have the form 5 [62] where we are omitting any possible overall factors that might depend on H, ω, and the couplings.Several comments are in order.In this expression, the first factor comes from the background oscillations contained in the coupling constants, the second factor from the bulk-to-boundary propagators, and the function f (k a , η) contains the non-oscillatory part of the bulk-to-boundary propagators of a massless scalar field ϕ, its derivatives with respect to time, other polynomial factors of η coming from the vertices, and factors of momenta k a coming from spatial derivatives.In (2.15) we assumed Bunch-Davies mode functions (i.e.only positive frequencies e +ikη ), but resonant models excite a small component of the negative frequency solution.For the moment we neglect these subleading contributions but we will discuss them in Section 4.
Assuming the Bunch-Davies massless mode functions the bulk-to-boundary propagator is and the n-th time derivative takes the from (2.17) Hence the contributions to f (k a , η) from the mode functions are polynomial in k a η.If we further restrict to manifestly local interactions, namely products of fields and their derivatives at the same spacetime point, also time and space derivatives contribute to f (k a , η) additional monomials in because every factor of k is accompanied by a factor of η both in the massless dS mode functions and in the time and space derivatives, where it arises from the metric that contracts the indices.In summary, we can write where a q are possibly complex numerical coefficients and Poly q is a homogeneous polynomial of degree q.
In passing, we mention that we do allow for a speed of sound c s that is different from the speed of light, but we will assume throughout that c s is constant in time.To the order to which we are working, c s only appears as an overall factor normalizing the size of non-Gaussianity.Since the overall normalization is anyways model dependent, we simply drop all factors of c s throughout the calculation.
The integral in (2.15) gives a pair of gamma functions for each monomial in f (k a , η).To see this, let us use Ht = − log(−Hη), the complex exponential representation of the cosine, and let's introduce the dimensionless frequency of oscillations, Then, up to overall factors, we can write a master integral in the following form (here we set H = 1, which eventually leads to an inconsequential overall factor) for some integer m and where, as usual, η → −∞ is reached with a small positive imaginary part so that the integral converges.Under a Wick rotation, this integral can be related to two Euler gamma functions6 with a complex argument where the factors of e ±πα/2 result from the Wick rotation.This, and the fact that in the large-α limit both gamma functions behave as Γ(m ± iα) ∼ e −πα/2 , implies a relative exponential suppression between the two terms.We will be interested precisely in this regime α = ω/H ≫ 1 in which many background oscillations take place per Hubble time, so we will ignore the exponentially suppressed term of the solution.
In the following we call this the resonance approximation, where p is related to the highest power of η appearing in the sum (2.19).The coefficients of the homogeneous polynomial Poly p+3 (which should not be confused with the polynomials in (2.19)) are complex, and its degree p + 3 is fixed by the discrete scale invariance condition (2.8).
For all tree-level diagrams, p is given by [10] where ∆ A is the mass dimension of interaction A and the sum runs over all interactions involved in a given diagram.Notice that this formula also applies to the order of the singularity of partial energy poles, but we will not need this here.For the contact bispectrum that we are studying here, we can easily understand this formula as follows.There is an overall η −4 factor coming from the integral measure, an η nt+ni factor coming from the metric contractions of the n t time and n i spatial derivatives, and a maximum factor of η 3 coming from the three propagators (2.17).Hence, So p counts the total number of derivatives of the vertex.This agrees with (2.25) when we notice that there is a single cubic vertex of mass dimension ∆ = 3 + n t + n i .
Before concluding a comment is in order.In previous references [9], the integral in (2.21) was in practice approximated using the stationary phase approximation 7 (even though the exact gamma function expression was also derived).Here we do not invoke such approximation and our resonance approximation includes all power law suppressed terms in 1/α, which would be missed by the leading order stationary phase result.As we will see shortly, these subleading terms beyond the stationary phase approximation are actually fixed by locality in the form of the manifestly local test.
In summary, all wavefunction coefficients coming from tree-level bulk interactions in dS as in (2.15) take the form of our ansatz (2.14) up to exponentially suppressed terms in e −πα (resonance approximation).This ansatz provides a rich class of minimal solutions to the discrete scale invariance constraint in (2.8). 7The method of stationary phase says that as where the phase φ(t) has a stationary point φ ′ (c) = 0 at a < c < b.Note that if g(c) ∼ O(x n ) with n ̸ = 0, the error becomes O(x n−3/2 ).

Bose symmetry
We are computing wavefunction coefficients of a bosonic field ϕ, so the result must be symmetric under any permutation of the field labels.Since the oscillatory part of (2.10) is symmetric, this implies that ψ3 must also be so, which in turn implies that Poly p+3 (k a ) can be written in a unique way in terms of the three elementary symmetric polynomials and the whole wavefunction coefficient is (2.30)

Cosmological optical theorem
The Cosmological Optical Theorem (COT) [11] still applies in the presence of time dependent couplings such as the ones appearing in resonant models, as long as they are real (and hence the theory is unitary) and the state of the universe approaches the Bunch-Davies state in the infinite past.This implies that the exact wavefunction coefficients that we are interested in must satisfy the COT for contact diagrams, namely Cosmo optical theorem: Indeed the exact propagators appearing in the presence of an oscillatory background do obey hermitian analyticity [14] and this ensures that the proof of the COT proceeds in the same way as described in [13].Moreover, notice that this relation is valid at all times and hence also applies in the case of IR divergences when we need to introduce an IR regulator, see Section 3.6.
However there is a practical problem: in our approach we are not computing the exact wavefunction coefficient ψ exact 3 .Rather we aim to capture its resonance approximation in (2.23), by neglecting terms that are exponentially suppressed in α.The issue is that the analytic continuation to negative energies, k → −k, present in the COT, turns an exponentially suppressed term such as e −πα e iα log k T into an unsuppressed term.The resulting term ends up canceling out with the leading term to satisfy (2.31).In the resonance approximation we neglect the exponentially suppressed contribution and so this cancellation cannot be checked or used 8 .Another way to see the problem is to notice that the resonance approximation trades the real coupling constant cos(ωt) for the complex coupling constant e iωt , hence invalidating the assumptions of unitarity.Throughout this paper we denote the resonance approximation simply by (2.32) Having identified this problem, there are two possible solutions.One is to modify the ansatz including the exponentially suppressed terms that are being neglected, in such a way that the solution satisfies (2.31).
The other is to find a modified version of the COT that is satisfied by the resonance approximation.The latter turns out to be the most straightforward approach and the one we choose to follow.To find a modified COT notice that our effective coupling e −iαHt /2 remains unchanged if we take its complex conjugate and send α → −α.In other words, it behaves just as a real coupling would under complex conjugation.This in turn ensures that the derivation of the COT of [13] goes through 9 and applies to the resonance approximation.Hence, for contact diagrams we impose the following Modified cosmo optical theorem: where α ∈ R + and we have made the α-dependence explicit.It would be interesting to check this modified relation to all orders in perturbation theory, but we will leave that for the future.
To apply the modified COT to our solution (2.30) we need to take into account the overall factor of e πα/2 that results from the time integrals, cf.Equation (2.22), and that we are ignoring.We have absorbed this exponential factor into an overall normalization that we are not fixing, but it is essential for the COT to work and so we temporarily restore it here10 : (2.34) For real momenta the conjugate wavefunction is then and plugging this into the modified COT (2.33) we get At a given order p in derivatives, we can write the polynomial explicitly as follows, Poly p+3 (k T , e 2 , e 3 , α) = where ⌊. ..⌋ is the floor function and the coefficients C r,s can depend on any parameter of the theory but here only the α-dependence is relevant for our purposes.Since (2.36) must be satisfied for any real kinematic configuration, we can match both sides of the equation individually for every monomial in {k T , e 2 , e 3 } that appears in them.What we get is that coefficients C r,s must satisfy this is, any appearance of α in them must come with a factor of i. Remember that we are ignoring an overall complex factor in ψ 3 , whose phase we have absorbed in the definition of k * , so what unitarity really dictates is that there is a relative phase of π/2 between successive powers of α in the wavefunction coefficient.In the following, we will not write the α-dependence of the coefficients C r,s explicitly to avoid clutter.

Manifestly local test
To capture local models we are interested in wavefunction coefficients for massless scalars that come from manifestly local interactions, namely interactions of the product of fields and their derivatives at the same spacetime point.In particular, the derivation in Section 3.2 of [12] applies in this case and we conclude that our wavefunction coefficients must satisfy the so-called Manifestly Local Test (MLT) where Bose symmetry allows us to pick the energy k 1 without loss of generality.It will be useful to notice that the different orders in α of the ansatz (2.30) do not need to separately satisfy the MLT, since the derivative ∂/∂k 1 acting on the oscillatory factor mixes them.The MLT thus relates different orders in α of the solution.Notice that, just like the COT, the MLT applies to all finite times as well and so can also be used in the presence of IR divergences.
If we demand our solution (2.30) to satisfy the MLT (2.39), we get . (2.40) Now we can plug the polynomial expansion (2.37) into this expression and solve for the coefficients C r,s .
Since the MLT must be satisfied for general kinematics, we have to demand every monomial in {k 2 , k 3 } of (2.40) to vanish.At order p in derivatives, this results in the following set of relations between coefficients: (2.41) In the next section we will use this to write the general solution.

Bootstrapping resonant bispectra
In Section 2 we have applied the bootstrap rules and we have obtained the most general three-point wavefunction coefficient under our assumptions, which is given by the wavefunction (2.30) with a polynomial (2.37) whose coefficients satisfy (2.38) and (2.41).Nevertheless, this expression contains redundancies, and in the next subsection we will eliminate them in order to present the result in a clearer way.After that, we will perform two consistency checks on our expressions-we will study the flat space limit and count the number of structures.Then we will explicitly write down the lowest-p results for the shapes of resonant non-Gaussianity.Finally, we will comment on parity-odd interactions and IR-divergent terms.

Final result
Our result ψ 3 (k a ), for a given order p in derivatives, is given by (2.30) with a polynomial (2.37) whose coefficients satisfy (2.38) and (2.41).Notice that such a formula contains the new kinematic structures that appear at order p in derivatives (those with k 0 T terms in Poly p+3 ) plus all the structures corresponding to lower p (those with an overall nonzero power of k T in Poly p+3 ).It is desirable to have an expression that just provides the new wavefunction coefficients ψ (p) 3 (k a ) at each order p in derivatives.Before imposing the MLT constraints (2.41), the new structures that appear in Poly p+3 at each order p are where the Kronecker delta ensures that the power of e 2 is integer.Let us now see how the MLT (2.41) ties some of the coefficients of these new structures to other coefficients in Poly p+3 .Notice first that there is an s = 0 term in (3.1)only when p is odd, but in that case (2.41) requires the corresponding coefficient to vanish, C p+3 2 ,0 = 0. We can then just start the sum in (3.1) at s = 1: Among these monomials, those with s > 1 are unconstrained by the MLT (2.41), so we will write them in the final solution with a free coefficient.On the other hand, the structure with s = 1-which is present only when p is even-is tied by the MLT (2.41) to other structures with different powers of k T as follows.
For p = 0 we have while for even p > 0 we have With this information, we are ready to state the main result of this paper.The independent shapes that contribute to the resonant, tree-level, three-point wavefunction coefficient at order p in derivatives are given by e s 3 for even p > 0 . (3.5) In this context, the most general three-point wavefunction coefficient will then be given by a linear combination of these shapes, with coefficients C r,s that must satisfy the COT constraint derived in Subsection 2.5, namely that they must be invariant under the combined action of complex conjugation and α → −α.

Flat space limit
The ansatz for ψ 3 (k a ) generically has a total energy pole of order p equal to the number of derivatives of the operator at hand.This is the same total energy pole encountered in usual scale-invariant wavefunction coefficients [52,53], and its residue is also given by the high-energy limit of a corresponding scattering amplitude.The argument is similar to that of the usual case-the structure of the bulk integral (2.15) implies that its result when k T → 0 is approximated by the integrand times η evaluated at η → −∞, a limit in which we recover the amplitude.Ignoring the phase from the exp(−iα log (k T /k * )) oscillatory factor, we have then The three-point amplitude 3 (k a ) comes from a vertex with p derivatives and hence has dimension p.Since we are considering boost-breaking interactions, we will be interested in the most generic scattering amplitude coming from p-derivative vertices that break Lorentz boosts, which is given by the most general polynomial of {e 2 , e 3 } of degree p [63], 3 (k a ) = Poly p (e 2 , e 3 ) . (3.7) If we look at the k T → 0 limit of the general result (3.5), we indeed find (3.6) with the amplitudes (3.7).This is a nice consistency check of our result, but there is more to it.Notice that our results (3.5) satisfy When the couplings are constants, non-Gaussian correlators are generated around horizon exit (except in the few cases with IR divergent interactions).Around horizon exit, the effect of the expansion are not negligible and the wavefunction coefficients or the corresponding correlators contain both a term corresponding to the flat-space amplitude and other curved spacetime contributions.While these two contributions can be separated using complex momenta and the mathematical limit k T → 0, for general physical configurations they are of the same order and generate a similar amount of non-Gaussianity.The resonant case stands in stark contrast to this expectation, and indeed can be thought of a way to zoom into the flat space limit by taking the large limit of the physical parameter α, namely the frequency of oscillations in units of the Hubble parameter.In the limit α ≫ 1, the non-Gaussianity is produced on scales parametrically shorter that Hubble and so is fully captured by flat space physics, up to higher orders in O(1/α).In other words, to leading order in α, resonant non-Gaussianities are just Minkowski amplitudes dressed with the oscillatory factor e −iα log k T !
Counting structures Here we count the number of possible structures that result from the bootstrap process up to a given order p, and prove that the total number is equal to that of independent cubic vertices with at most p derivatives.The latter can be found by counting the structures in the amplitude (3.7), and is equal to the number of non-negative integer solutions {r, s} to 2r + 3s ≤ p.This number is [64] N amp (p) = The result (3.5) introduces, at every order in p, a number of new structures parameterised by the real constants C r,s .These constants come precisely from the amplitude (3.7), as shown in (3.8), so their number up to p derivatives will be exactly N amp (p).We conclude then that the number of independent structures in ψ 3 up to order p is N amp (p), and hence it coincides with the number of cubic vertices with at most p derivatives.
Notice the difference with the standard bootstrap case, in which the number of structures was equal to the number of scattering amplitudes plus one [12].The extra structure corresponded to the local nongaussianity associated to the only manifestly local field redefinition, ϕ → ϕ + ϕ 2 , which has the form T − 3k T e 2 + 3e 3 .In the resonant case, however, we have11 with b * < 1 a parameter that controls the monotonicity of the potential.The corresponding sum ψ 2 (k 1 ) + ψ 2 (k 2 ) + ψ 2 (k 3 ) cannot be captured by our ansatz (2.10) for ψ 3 and consequently we do not find in our analysis a three-point wavefunction coefficient coming from the above field redefinition.

Comparison with previous results
The literature on resonant non-gaussianities has traditionally presented the bispectrum results as an expansion in α ≫ 1, often computing just the leading order but sometimes providing the NLO and even the NNLO corrections too.The overall phase of the oscillations-which we called k * in our approach-is an undetermined parameter because it depends on the background solution and on the expansion history after inflation.Here we recall that a shift of said phase is degenerate with a change in the functional dependence of the subleading terms in the α expansion.To see this, consider a resonant bispectrum of the form where k a collectively denotes {k 1 , k 2 , k 3 } and f (k a ) and g(k a ) are rational functions independent of α.The overall phase is given by k * and is undetermined, so we are free to shift it.If we shift it by an α-dependent amount, like for example ∆φ = arctan(1/α), we have (3.12) Notice that the kinematic dependence of the NLO term of the expansion has now changed, and a new term has appeared at NNLO.This fact, which was already pointed out in [22] to explain an apparent discrepancy with the result of [9], is general: the form of the subleading terms in the α-expansion can be changed in a specific way by an α-dependence change of the overall phase.Fortunately, this degeneracy is inconsequential for phenomenology since the overall phase of the oscillations is marginalised when these signals are looked for in the data.However, it will be important when comparing our bootstrap results to the previous literature.Recall that our approach does not determine the overall size and phase of the wavefunction coefficients, so comparison must always be done up to overall factors containing powers of iα (as per the COT constraint of Section 2.5) or, equivalently, up to a phase shift (for instance, the phase shift ∆φ = arctan(1/α) in the above example can be understood as coming from an overall factor of (1 − 1/iα) ψ 3 in the wavefunction coefficient).
Finally, let's mention that resonant non-Gaussianity from interactions with a small number of derivatives, as indeed the case in the original example [9], induce IR divergences in the wavefunction coefficient, which had not been previously discussed in the literature.We make up for this omission in Section 3.6.

Lowest-p shapes
We will now explicitly write all the possible resonant bispectra at increasing order p of the total energy pole, up to p = 3. Recall that the relation between the bispectrum and the wavefunction is given by In some cases, we will single out certain linear combinations that correspond to wavefunction coefficients coming from specific vertices in the Effective Field Theory of Inflation (EFToI).In the context of resonant inflation the continuous shift symmetry of the Goldstone boson π becomes a discrete shift symmetry, and the couplings acquire an oscillatory time dependence [19,60].Schematically, the wavefunction coefficients that we will find correspond to interaction terms in the Lagrangian of the form where the phase and overall size of the couplings are immaterial since they are not captured by the bootstrap procedure.Here we focus on the IR finite part of the result, and we discuss IR divergences in Section 3.6 The result (3.5) for p = 0 is simply where we have re-introduced the overall phase k * .The bispectrum is then, up to an overall real factor, (3.16) Figure 1a shows a slice of the corresponding shape A bulk computation confirms that (3.16) corresponds to the bispectrum arising from a vertex π 3 (up to a suitable choice of phase k * as explained in Subsection 3.3): This bispectrum was computed in [19,22] by evaluating the time integrals at the stationary phase point, but the subleading corrections in α ≫ 1 were not correctly taken into account, thus getting an incorrect result at O(1/α).The symmetric squeezed limit of this bispectrum is 12   lim

.19)
12 By symmetric squeezed limit of a bispectrum B(k 1 , k 2 , k 3 ) we mean lim where k l and ks are the long and short modes respectively.The result is generically expressed in terms of k l ≡ |k l |, ks ≡ |ks| and cos θ ≡ kl • ks.
It was argued in [65,66] that the term that diverges as O(1/k 2 l ) must be absent in the symmetric squeezed limit.This was checked explicitly for a resonant model in [65], although the bispectrum was obtained via the stationary phase approximation and the cancellation only worked at leading and next-to-leading orders in 1/α.For our bootstrap solution, however, Equation (3.19) shows that the absence of a O(1/k 2 l ) term is exact.We can see from (3.5) that there is no new wavefunction coefficient for p = 1.In the EFToI there is one vertex with one derivative, π 2 π, and from this point of view the corresponding bispectrum can only be given by the structure that we already had at p = 0, i.e.

.20)
A bulk computation confirms that this is indeed the case, and the bispectra coming from π 3 and π 2 π are the same up to an α-dependent phase shift.The bispectrum (3.20) coincides with the result in [9] at leading order in 1/α but differs starting at next-to-leading order.The reason is again that the bispectrum there was computed via the stationary phase approximation, which gives the correct result only to leading order in 1/α.The difference however can be re-absorbed into an α-dependent phase and so is phenomenologically inconsequential.

p = 2
In this case we have with bispectrum The symmetric squeezed limit is lim where we can see again the absence of the O(1/k 2 l ) term.The EFToI vertex π π2 gives rise to exactly the wavefunction coefficient (3.21), so its bispectrum is The wavefunction corresponding to a π(∂ i π) 2 vertex is obtained by adding up (3.15) and (3.21) with The corresponding bispectrum is This coincides with the expression computed at leading order in α ≫ 1 in [62], but here we provide the full (non-exponentially-suppressed) result.The symmetric squeezed limit of this bispectrum is lim At three derivatives there is a new structure, (3.28) This wavefunction coefficient corresponds to the EFToI vertex π3 [62], and the bispectrum is Taking the sum of (3.15), (3.21), and (3.28) with gives the wavefunction coefficient coming from the EFToI vertex π(∂ i π) 2 , which is The corresponding bispectrum is which coincides with the result computed at leading order in 1/α in [19], but here we find the other power-law-suppressed terms too.The symmetric squeezed limits of both (3.29) and (3.32) satisfy:

Parity-odd interactions
Here we briefly comment on parity-odd correlators in the context of resonant non-Gaussianity.For scalar fields, the power spectrum and bispectrum are necessary parity even and so the first correlator that can be parity-odd is the four-point function or trispectrum.In this sense, this subsection is a bit orthogonal to the main focus of this work.Our main observation is that resonant non-Gaussianity avoids the no-go theorems of [49,67] and can easily induce any parity-odd n-point correlator at tree level for n ≥ 4. Indeed, the oscillating coupling constants invalidate the assumptions of those no-go theorems.As a consequence, for resonant non-Gaussianity the wavefunction is generically complex, with different phases for different values of the kinematics, ψ n ∝ k −iα T and a parity-odd correlator picks up the non-vanishing imaginary part.As a simple example consider the EFT of single-clock inflation.Then any of the following couplings 13(see [49]) leads to a parity-odd trispectrum B P O 4 already for a contact tree-level diagram.Notice that the resulting B P O 4 has a total energy pole and hence is not in factorized form, in contrast to what one finds for constant couplings in de Sitter [50].

Infrared divergences
So far we have neglected the possibility that wavefunction coefficients do not converge as η → 0, a case which we will refer to as IR divergences.These divergences are expected to appear for interactions with few derivatives and more precisely when 2n t + n i < 4, where n t and n i count the number of time and space derivatives, respectively.Here we discuss these terms in detail.The upshot is that there are several contributions to the wavefunction but only one contribution to the bispectrum of a spectator field, which furthermore must vanish for the bispectrum of curvature perturbations in single field inflation.
The main observation is that, in the presence of IR divergences, we have to choose a late time cutoff η 0 to evaluate the wavefunction coefficients.As a result, our discrete scale invariance constraint becomes If we furthermore impose locality in the form of manifestly local interactions, we are guided towards and ansatz containing the possible divergences 1/η 3 0 , 1/η 2 0 , 1/η 1 0 and possibly a logarithmic term, each multiplied by polynomials in the norm of the momenta and possibly a function that breaks continuous scale invariance to the discrete subgroup.Because we are interested in the late time contributions, there is no need to resort to the resonance approximation and instead we proceed considering the full contribution induced by the real coupling constant cos ωt, so that we can use the cosmological optical theorem (COT) in its original form.After carefully studying the bulk integral, we arrive at the following ansatz for the IR divergent terms where Poly 0,1,2 (k a ) are homogeneous polynomials in the momenta of the indicated degree (so that Poly 0 is just a complex number) and φ 1,2,3 are real phases.This ansatz can be justified as follows.Consider the time integral of the form Since we are interested in IR divergences, we focus on the late time contributions coming from the part of the integral at very small η.Then we can expand the exponential and the polynomial function f and keep only the IR divergent terms in the integrand Now the indefinite integral can be easily performed and gives precisely the terms shown in the ansatz.
To proceed, we notice that the COT imposes that the terms even in k a and hence odd in η 0 are purely imaginary.These terms are therefore present in the wavefunction but do not contribute to the correlator, because the latter picks up just the real part of ψ n (for parity-even interactions, which is always the case for scalar bispectra).Then we notice that the manifestly local test demands Poly 1 = 0, because all linear terms in k a must vanish by manifest locality.The expert reader will have noticed that this is verbatim what happens in the case of scale invariant IR divergences in the wavefunction, see [12] for a detailed discussion.Moreover, the MLT imposes that the only possible symmetric polynomial of degree three is the sum of k 3 a and the overall coefficient must be real by virtue of the COT.Hence the only term that can contribute to the correlator is where we have re-absorbed the phase into the constant η * .Notice that this is just the well-known and loved local non-Gaussianity, albeit with an amplitude that oscillates with time Several aspects of this result are novel and interesting.
First, notice that this local contribution is generated by all interactions with 2n t + n i < 4, namely ϕ 3 , φϕ 2 and (∂ i ϕ) 2 ϕ.When these interactions have a coupling that is constant in time, they lead to IR divergences also in the correlator in the form of a term proportional to log(−k T η).These divergences are to be expected because, as we approach the future conformal boundary of dS, the corresponding interaction acts for an infinite amount of time, in contrast to interactions with more derivatives, which redshift away.Often these divergences are interpreted as failures of perturbation theory rather than an actual physical instability, and a non-perturbative semi-classical treatment can sometimes resum these divergences into finite effects, such as for the more studied case of λϕ 4 [68].Here instead the amplitude of the bispectrum never grows because cos(α log(−ηH)) remains bounded at all times.So perturbation theory remains valid but the time integral never converges (indeed η = 0 is an essential singularity).
Second, what we have discussed so far is a generic spectator field and one should ask what happens to curvature perturbations ζ.We know that in (attractor) single field inflation, ζ should become constant on super Hubble scales and this is incompatible with the time dependence in (3.41).Indeed we expect that to compute the bispectrum of ζ one would have to add to the inflation self interaction that generates (3.41) also a second-order gauge transformation from ϕ to ζ that cancels the first contribution exactly.Indeed notice that the field redefinition ϕ → ϕ + cos(α log(−Hη))ϕ 2 , (3.43) generates exactly the contribution in (3.42).This is presumably the reason why this contribution had not been discussed previously in the literature [9,19].
Third, notice that the local contribution in (3.42) oscillates with time, but does not oscillate with the k a .This is in contrast to all the other non-Gaussian shapes we have discussed here and it is a reminder that this term does not come from a resonance, but rather from much slower and much later super Hubble evolution.However, this contribution still breaks scale invariance (dS dilatations involving k a and η as in (3.37)) to a discrete subgroup.For the remainder of this work we set this contribution aside and focus on the IR-convergent resonant shapes.

Corrected mode functions
Throughout this paper we have worked with the assumption that the scalar field mode function is the Bunch-Davies one, cf.Section 2.3.Nevertheless, the background oscillations of the class of models we are considering generically lead to resonant particle production during the inflationary evolution [6].As a result, the mode function of ϕ gets a negative frequency correction: with ϕ + k (η) = (ϕ − k (η)) * the Bunch-Davies solution in (2.16).The function c (−) k (η) quantifies the excitation of the negative frequency mode and was computed in detail in [9,19]; it looks roughly like a step function centered at the particle-production time k η p = −α/2 that interpolates between c (−) k (η) ∼ 0 at early times (before the resonance takes place) and a finite late-time value c (−) k (0).This correction to the mode function gives rise to contributions to the resonant bispectrum that differ from the Bunch-Davies ones-some of them were studied in [19,62].These contributions have two important features.Firstly, they are suppressed with respect to their Bunch-Davies counterparts, since the amplitude of the negative frequency mode is small: with b * < 1 the monotonicity parameter.Secondly, they peak in the folded limit, where the sum of two momenta equals the third.The reason is that, if one substitutes the mode function of one of the interacting fields (let us choose the one with momentum k 3 for concreteness) by the negative frequency one, the bulk integrals (2.15) become This integral has now a stationary phase point at (k in the denominator of ψ 3 .This would naively imply the presence of folded singularities, but a more careful analysis reveals that they are actually not there.Notice that in the limit k 1 + k 2 − k 3 → 0 the time at which the non-gaussianity is generated is pushed to the very far past, η res → −∞, much before the time η p at which the negative frequency mode function is excited.But the function c (−) k3 (η res ) at times η res ≪ η p vanishes, so the integral (4.3) and hence the contribution to ψ 3 will be zero.In other words, the folded limit probes the bulk times before the negative frequency component was excited, and hence it cannot receive a contribution from such component.Despite the suppression (4.2) and the absence of singular behaviour, the analyses of [19,62] revealed that the enhancement near the folded configuration makes the negative frequency mode contributions comparable to the Bunch-Davies bispectra.
One could, in principle, try to bootstrap these negative frequency bispectrum contributions by applying an algorithm similar to the one developed in this paper, but one would find several obstacles.The most important one is that the arguments of Section 2.3 would no longer be valid-the function c (−) k (η) is not rational and hence the result of the integral (4.3) cannot be approximated by the simple ansatz (2.14).In fact, no exact solution has been written for these time integrals and the previous literature has always worked with some approximation of the function c (−) k (η).A possible strategy for the bootstrap would be to follow [62], which modelled the step-like time dependence of c After evaluating (4.3) in the stationary phase approximation, this would lead to an overall factor of which accounts for the ordering of events in the bulk-the negative frequency mode excitation and the nongaussianity production-and ensures that the behaviour described in the previous paragraph is correctly reproduced.The other function in the integrand, f (k a , η), would still be rational in this approximation, so the ansatz for ψ 3 would take the form of a rational function multiplied by the oscillatory exponential and (4.5).The Bose symmetry requirement of Section 2.4 would now be less constraining, but straightforward to implement.The unitarity and locality constraints (sections 2.5 and 2.6 respectively), however, would require a careful analysis of the properties of c (−) k (η). 14We think that this continuous need of invoking the properties of the bulk propagators and integrals partially defeats the purpose of the bootstrap, so we have limited ourselves to outlining the strategy without actually pursuing it.We hope that future developments on our understanding of cosmological correlators will allow to bootstrap these negative frequency bispectra form a more "boundary-centered" point of view.

Constraints from the CMB
In this section we constrain some of the primordial bispectrum shapes bootstrapped in the previous pages using Planck 2018 temperature and polarization data.For that, we make use of the recently-developed CMB-BEST pipeline [55].
We put constraints directly on the bispectra given by the structures that we found with up to three derivatives (cf.Subsection 3.4).By doing this, we do not take into account the fact that in single-field inflation the leading term ∝ 1/k 3 l in the squeezed limit does not contribute to local observables [69,70].Hence, if we want to interpret our bispectra as correlators of the primordial curvature perturbations ζ, they can only arise in a multi-field inflationary theory.We have used the high-resolution "Legendre" basis of CMB-BEST, with mode number p max = 30.This allows for an accurate analysis of our shapes in a reasonable range of frequencies without the need for developing a targeted resonant basis.Nevertheless, the accuracy of the constraints decreases as the frequency α of oscillations increases, mostly due to the difficulty for the chosen basis to capture the rapid oscillations of our shapes at small k T ; we have found in all cases that the "convergence epsilon", which the code uses to estimate the offset in f NL , exceeds ϵ = 0.1 when the frequencies reach α ∼ 35 − 40.At higher frequencies the basis we use does not accurately reproduce the resonant signal.For this reason, and since our results are valid in the α ≫ 1 regime, we restrict the analysis to the range 10 ≤ α ≤ 34.We have maximized our results over the phase of the oscillations, showing constraints for the best-fit value of k * at each frequency α.
We follow the conventions of [55,71], in which f NL is defined by with Φ the Bardeen potential, which on superhorizon scales is related to the curvature perturbation ζ by Φ = 3ζ/5.The bispectrum is normalized following with k p = 0.05 Mpc −1 the pivot scale and A the amplitude of the potential power spectrum,

.3)
Results There are only three independent structures satisfying p ≤ 3: the shape in (3.15) at p = 0, that in (3.21) at p = 2, and finally the shape in (3.28) for p = 3.The results of the analysis are plotted in Figure 2, which shows the constraints on f NL for the three shapes at hand.We plot the detection significance |f NL |/σ(f NL ), i.e. the ratio between the absolute value of the estimated f NL and the sample variance σ(f NL ).We observe that for most shapes this value is below 1σ, with some peaks around 1.75σ.We are scanning over 25 different frequencies α for each shape, so it is natural to find some large detection values as a result of the look-elsewhere effect.A more refined analysis would require quantifying this effect and adjusting the significance accordingly, for example using the treatment in [72], but that goes beyond 14 Consider for example unitarity.It was shown in [14] that the corrected mode function (4.1) satisfies hermitian analyticity if the full solution for c (−) k (η) is used.However, the approximation (4.4), needed to make the integrals attainable, spoils this property, and does not seem to be compatible with the α → −α mapping described in Section 2.5 either.
the scope of this work.The label "raw significance" in our plots emphasizes that we are not adjusting for the look-elsewhere effect.All things considered, we conclude that there is no significant evidence for detection of the three p ≤ 3 resonant shapes.To a large extent this could have been anticipated given the analysis performed by the Planck collaboration [71].Certain linear combinations of the structures found at p ≤ 3 correspond to bispectra generated by vertices present in the effective field theory of inflation (EFToI ) with up to three derivatives, cf.Section 3.4.There are six such vertices.For four of them, π 3 , π 2 π, π π2 , and π3 , the bispectra are degenerate with the p = 0, p = 2 and p = 3 structures that we have already constrained and whose results are plotted in Figure 2. 15 For the other two vertices, π(∂π) 2 and π(∂π) 2 , the bispectra are respectively (3.26) and (3.32), and the bounds on f NL are shown in Figure 3. Just as before, we see low overall significance with a few peaks above 1.5σ, the largest of which is found for π(∂π) 2 at α = 20, with ∼ 2σ.Taking into account the look-elsewhere effect, we do not find significant evidence for rejecting the null hypothesis f NL = 0.

Single field inflation
We conclude with a comment on the non-Gaussian signal in single field inflation.Recall from Subsection 3.4 that the bootstrapped bispectra with p = 0 and p = 2 diverge as 1/k 3 l in the squeezed limit |k l | → 0. References [69,73] showed that, in single-field inflation, the projection to late-time local observables removes that leading term of the squeezed limit.A fully correct comparison of primordial reasonant non-Gaussianity with late time data would involve accurately taking into account the evolution after inflation and would probably leverage the simplification obtained by using conformal Fermi normal  coordinates [74] as shown in [70].However, one could find a plausible estimate of the constraints on singlefield inflation in a much simpler way as follows.One can remove by hand the leading terms of (3.16), (3.22) and (3.26) in the squeezed limits.We have done this and then compared to Planck data.We found significance below 1σ for most models, with peaks around ∼ 1.7σ for a few of them.We conclude that there is also no compelling evidence for primordial non-Gaussianity of this type coming from single-field inflation.

Conclusions
In this work we have demonstrated how to generalize the (boostless) cosmological bootstrap beyond exact scale invariance and we have computed an infinite class of so-called resonant bispectra, characterized by periodic oscillations in the sum of the norms of the momenta.We have re-derived the resonant bispectra that were known in the literature and we have computed in closed form additional shapes corresponding to tree-level interactions with any number of derivatives.We searched Planck data for these signals finding no evidence.In addition to the phenomenological implications of our work, our results are interesting from other points of view.First, we have demonstrated that, to leading order in the large frequency of oscillations, resonant non-Gaussianity is fully specified by a flat-space amplitude for all physical configurations, and not just close to the unphysical total energy pole.In this sense, this a concrete mechanism to "uplift" Minkowski amplitudes to de Sitter wavefunction coefficients without deforming them.Second, our final result can be thought of as giving a meaning to the apparently nonsensical question "what is the bispectrum for an interaction with a complex number of derivatives?"Indeed, we have shown that resonant non-Gaussianity can be obtained by adding an imaginary part to the integer parameter p that counts the order of the total energy singularity.Third, we have demonstrated the power of general constraints such as locality in the form of the manifestly local test, by discovering an (inconsequential) inconsistency at subleading order in the results published in the literature.We have also shown how to modify the constraints of unitarity in the form of the cosmological optical theorem to leverage the analytic dependence of the result on a parameter of the model, namely the frequency of oscillations.
Our results can be extended in a series of promising directions: • As we discussed in Section 4, there are additional contributions coming from the excitation of negative frequency modes away from the Bunch-Davies vacuum.It would be interesting to see if these can be further constrained or even computed explicitly using bootstrap techniques.
• While we relaxed exact scale invariance, we still assumed a discrete form of scale invariance to land on a sufficient constrained ansatz.In so doing, we have ignored the physical effect of frequency drifting, which can be caused for example by slow-roll corrections [75].The power spectrum templates in Section 6.1 of [75] suggest that a simple deformation of our ansatz (that is still compatible with discrete scale invariance (2.8)) might capture this drift.It would be intersting to investigate if helps to provide a better fit to the data.
• We hope that this work is a first step towards extending the cosmological bootstrap into a wider range of inflationary setups.In particular, it would be interesting to apply the lessons learned here to other patterns of scale-invariance breaking, like sharp potential steps or bumps that lead to sinusoidal running of the bispectrum [2,76].Other interesting targets would be the non-Gaussianities induced by heavy particle production via a discrete-shift-symmetric coupling to the inflaton [77,78] or via features in the potential [79,80].
As we consider a "boundary" quantity such as the η → 0 limit of a wavefunction coefficient, the organising principle that reflects (and hence substitutes) the bulk perturbative expansion is the analytic structure of the wavefunction.A boundary strategy then consist in imposing the simplest possible analytic structure in the hope to capture the leading tree level contributions in the bulk.For simplicity, we restrict our discussion here to the two-point wavefunction coefficient ψ 2 (z), which we now take to depend on a single complex variable.This avoids the complications of multi-dimensional complex analysis that arise for higher-point objects.
Let us now rank the possible analytic structures of ψ 2 (z) in increasing level of complexity, until we find a sufficiently rich new ansatz: where z m ≡ λm z and C m is obtained by multiplying every z ∈ C by λm .Hence, we see that for every n, m ∈ Z and some λ ∈ (1, ∞) we must have c n = λ m(n−3) c n .But this is only possible if c n = 0 for all n ̸ = 3, i.e. if the wavefunction coefficient takes the scale-invariant form ψ 2 (z) ∝ z 3 in 0 < |z| < ∞.We conclude that for this simplest analytic structure, a function ψ 2 (z) different from the standard scale-invariant one does not exist.
2. Isolated singularity/ies at z 0 ̸ = {0, ∞} If an isolated singularity is present at a location z = z 0 with 0 < |z 0 | < ∞, then by the property (A.2) there are also infinitely many isolated singularities at z = λm z 0 for every m ∈ Z, see Figure 4.In other words, if isolated singularities are present at locations different from the origin or infinity, then there are infinitely many of them.This infinite sequence of isolated singularities is such that any arbitrarily small circle around z = 0 contains infinitely many of them.This means that the singularities accumulate at z = 0, which is then a cluster point.We don't know of any argument why such a functional form would not be consistent with (A.2), but notice that it would single out a (complex but finite) momentum scale z 0 (or several of them if we consider isolated singularities at several points not related by an integer power of λ).We will then dismiss this analytic structure as unphysical, since neither the setup of the problem nor the symmetry property (A.2) select a particular scale.
3. Branch points at z = 0 and z = ∞ The next level of complexity of the analytic structure comes when we consider multivalued functions and the associated branch points.If a function is multivalued, it must have at least two branch points in the extended complex plane C ≡ C ∪ {∞}.By the previous argument, if we had a branch point at some finite 0 < z < ∞ then we must also have its infinitely many images.The only possibility with a finite number of branch points is hence a function with two branch points at {0, ∞} and a branch (a) Shape of the p = 0 bispectrum B (0) (3.16).(b) Shape of the p = 2 bispectrum B (2) (3.22).(c) Shape of the p = 3 bispectrum B (3) (3.29).

Figure 1 :
Figure 1: Shape S(k 1 , k 2 , k 3 ) ≡ (k 1 k 2 k 3 ) 2 B(k 1 , k 2 , k 3 ) of three resonant bispectra plotted as a function of x a ≡ k a /k 1 for k * = k 1 and α = 70.Notice that due to the absence of scale invariance these plots change under a rescaling k a → λk a or, equivalently, under a change of phase k * .

Figure 2 :
Figure2: Constraints on f NL for the independent resonant shapes found by our bootstrap procedure at zero, two, and three derivatives, whose corresponding bispectra are given in (3.16),(3.22),and (3.29) respectively.We plot the significance |f NL |/σ(f NL ) as a function of the frequency α of oscillations, maximized over phase k * at each frequency.The constraints are obtained from Planck 2018 data using the CMB-BEST pipeline[55].

Figure 3 :
Figure 3: Constraints on f NL for the resonant shapes coming from EFToI vertices π(∂π) 2 and π(∂π) 2 , whose corresponding bispectra are given in (3.26) and (3.32).We plot the significance |f NL |/σ(f NL ) as a function of the frequency α of oscillations, maximized over phase k * at each frequency.The constraints are obtained from Planck 2018 data using the CMB-BEST pipeline [55].

ReFigure 4 :
Figure 4: Possible analytic structure of ψ 2 (z) in the C-plane.An isolated singularity at z 0 ̸ = {0, ∞} repeats itself at z = λm z 0 for every m ∈ Z.This means that the infinite sequence of singularities extends towards z = ∞ and accumulates towards z = 0, which becomes a cluster point.