{'fraction_non_alphanumeric': 0.04286493501608141, 'fraction_numerical': 0.016993623864427182, 'mean_word_length': 4.313296357367711, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Large slow-roll corrections to the bispectrum of noncanonical inflation

Nongaussian statistics are a powerful discriminant between inflationary models, particularly those with noncanonical kinetic terms. Focusing on theories where the Lagrangian is an arbitrary Lorentz-invariant function of a scalar field and its first derivatives, we review and extend the calculation of the observable three-point function. We compute the"next-order"slow-roll corrections to the bispectrum in closed form, and obtain quantitative estimates of their magnitude in DBI and power-law k-inflation. In the DBI case our results enable us to estimate corrections from the shape of the potential and the warp factor: these can be of order several tens of percent. We track the possible sources of large logarithms which can spoil ordinary perturbation theory, and use them to obtain a general formula for the scale dependence of the bispectrum. Our result satisfies the next-order version of Maldacena's consistency condition and an equivalent consistency condition for the scale dependence. We identify a new bispectrum shape available at next-order, which is similar to a shape encountered in Galileon models. If fNL is sufficiently large this shape may be independently detectable.


Introduction
In the last decade, our picture of the early universe has become numerically precise-a quantitative revolution made possible by analysis of the cosmic microwave background (CMB) radiation. The CMB is usually interpreted as a relic of an earlier hot, dense, primordial era [1]. During this era the universe is believed to have been filled with an almost-smooth plasma, in which tiny perturbations were seeded by an unknown mechanism. Eventually, photons decoupled from the cooling plasma to form the CMB, while matter particles collapsed under gravity to generate bound structures. The origin of these perturbations is a matter of dispute. A possible candidate is the inflationary scenario, in which an era of quasi-de Sitter expansion smoothed the universe while allowing quantum fluctuations to grow and become classical on superhorizon scales. In its simplest implementation, inflation forecasts an approximately scale-invariant, Gaussian distribution of perturbations [2,3], characterized entirely by their two-point correlations. Although falsifiable, these predictions are in agreement with present-day observations [4][5][6].
The increasing sophistication of CMB experiments has suggested that nontrivial threeand higher n-point correlations could be detected [7,8], and in the medium-term future it is possible that competitive constraints will emerge from the nongaussian statistics of collapsed structures [9]. In principle, valuable information is encoded in each n-point function. In practice, extracting information from the four-point function is already computationally challenging [10,11], and it is unclear whether useful constraints can be obtained for the fiveand higher n-point functions. For this reason attention has focused on the amplitude of three-point correlations, f NL [7,8]. For a Gaussian field these correlations are absent and f NL = 0.
In simple models both the de Sitter expansion and quantum fluctuations originate from a single scalar field, φ, although more complicated possibilities exist. 1 But whatever combination of ingredients we choose, the perturbations inherit their statistical properties from a microphysical Lagrangian by some mapping of the correlation functions of a field such as φ [13][14][15]. The two-point function for a scalar field with canonical kinetic term was calculated during the early days of the inflationary paradigm [2,[16][17][18][19][20][21], followed by partial results for the three-point function [22][23][24][25]. The complete three-point function was eventually computed by Maldacena [26], and later extended to multiple fields [27] and four-point correlations [28,29].
Maldacena's calculation showed that f NL would be unobservably small in the simplest model of inflation-of order the tensor-to-scalar ratio, r, which is constrained by observation to satisfy r 0.2 [6]. However, it was clear that information about the Lagrangian operators responsible for generating fluctuations in more exotic theories could be encoded in the threeand higher n-point functions. Which operators could be relevant? Creminelli argued that if the dominant kinetic operator for the background field was the canonical term (∇φ) 2 , where ∇ is the spacetime covariant derivative, then the three-point correlations of its perturbations were effectively indistinguishable from Maldacena's simplest model [30]. However, in theories where other operators make significant contributions, the amplitude of three-point correla-tions can be radically different. This gave valuable diagnostics for more complicated scenarios such as "ghost inflation" [31,32] or models based on the Dirac-Born-Infeld (DBI) action [33].
In these examples, and many others, the scalar field action controlling the dynamics of both background and perturbations can be written where X ≡ −g ab ∇ a φ∇ b φ, the spacetime Ricci scalar is R, and P (X, φ) is an arbitrary function. More recently, diagnostic signatures for related models based on "Galileon" actions have been obtained [34][35][36][37][38]. We describe models in which the background kinetic structure is dominated by (∇φ) 2 as canonical, and more general examples of the form (1.1) as noncanonical. Eq. (1.1) was suggested by Armendáriz-Picón et al. [39], who described it as 'k-inflation.' The corresponding two-point function was obtained by Garriga & Mukhanov [40]. Three-point correlations were studied by Gruzinov in a decoupling limit where mixing with gravity could be neglected [41]. Gravitational interactions were included in Ref. [42].
In this paper we return to actions of the form (1.1) and reconsider the three-point function. The analyses discussed above controlled their calculations by invoking some form of slow-roll approximation, typically restricting attention to the lowest power of ε ≡ −Ḣ/H 2 (or other quantities of similar magnitude) where H is the Hubble parameter and an overdot denotes a time derivative. One may expect this procedure to yield estimates accurate to a fractional error of order ε, which in some models may be as large as 10 −1 to 10 −2 . Where f NL is small, as in canonical scenarios, this is an excellent approximation. However, in models where f NL is numerically large, a fractional error of order ε may be comparable to the resolution of forthcoming data from the Planck survey satellite, especially if the O(ε) term enters with a relatively large coefficient. In the equilateral mode, Planck is expected to measure the amplitude of three-point correlations with an error bar in the range 25 -30, and a future CMB satellite such as CMBPol or CoRE may even achieve ∆f NL ≈ 10 [53,54]. Ideally, we would like the theoretical uncertainty in our predictions to fall below this threshold.
Corrections to the power spectrum at subleading order in ε are comparatively wellunderstood. For the two-point function of a canonical field, the necessary expressions were obtained by Stewart & Lyth [55], who worked to second order in the slow-roll expansion. 2 Even more precise results, accurate to third order, were given by Gong & Stewart [57,58]. Stewart & Lyth's calculation was generalized to the noncanonical case by Chen et al. 3 [43], although third-order results analogous to Gong & Stewart's are not yet available. The secondorder calculation was systematized by Lidsey et al. [60], who organized their expansion into "lowest-order" terms containing the fewest slow-roll parameters, followed by "next-order" 2 Stewart & Lyth were obliged to assume that the slow-roll parameter ε = −Ḣ/H 2 was small. Grivell & Liddle [56] gave a more general construction in which this assumption was not mandatory, but were unable to obtain analytic solutions. Their numerical results confirmed that the Stewart-Lyth formulae were valid within a small fractional error. 3 Earlier results had been obtained in certain special cases, by Wei et al. [59].
terms containing a single extra parameter, "next-next-order" terms containing two extra parameters, and so on. In what follows we adopt this organizational scheme. Corrections to f NL have received less attention, but are more likely to be required. The prospect of a sizeable "theory error" in a lowest-order slow-roll calculation was identified by Chen et al. [43], who computed next-order corrections for P (X, φ) models. But their final expressions were given as quadratures which were not easy to evaluate in closed form. How large a correction should we expect? The third-order action derived in Ref. [42] is perturbative in the amplitude of fluctuations, but exact in slow-roll quantities: it does not rely on an expansion in slow-roll parameters, although it mixes lowest-order and next-order terms. Therefore a subset of O(ε) corrections to the vertices of the theory may be obtained without calculation, and can be used to estimate the size of the unknown remainder. 4 In the case of DBI inflation [33,61], it can be checked that these terms generate a fractional correction 101ε/7 ≈ 14ε in the equilateral limit. For the value ε ≈ 1/20 suggested by Alishahiha, Silverstein & Tong [33] this is of order 70%. 5 In the absence of further information one should expect the unknown terms to be of comparable magnitude but unknown sign, making the uncertainty in a purely lowest-order prediction of order 10%-20% even if ε is substantially smaller than that required for the model of Ref. [33]. In models where f NL is observable, but ε is not negligible, the lowest-order prediction is likely to survive unaltered only if an accidental cancellation removes most of the O(ε) terms.
There are other reasons to pursue next-order corrections. In models based on (1.1), only two "shapes" of bispectrum can be produced with amplitude enhanced by a low sound speed. At subleading order, we may expect more shapes to appear-the question is only with what amplitude they are produced. If these shapes are distinctive, and the amplitude is sufficiently large, it may be possible to detect them directly. Were such a signal to be found, with amplitude appropriately correlated to the amplitude generated in each lowestorder shape, it would be a striking indication that an action of the form (1.1) was responsible for controlling the inflationary fluctuations.
The calculation presented in this paper enables a precise evaluation of the slow-roll corrections, resolving the large theoretical uncertainties and providing detailed information about their shapes. For canonical models where |f NL | 1 we verify that the effect is small. For noncanonical models we find that the next-order corrections may be enhanced by a numerical prefactor, and can be several tens of percent in interesting cases. Although nextorder calculations are likely to be sufficient for Planck, this suggests that next-next-order 4 These are the next-order components of the gi in Eqs. (3.10)-(3.11), for which it is not necessary to compute corrections from the propagator or account for time-dependence of the interaction vertices. If the subset is representative of the next-order terms in sign and magnitude then it gives a rough estimate of the possible effect. In fact, we will see in § §6-7 that the signs in this subset are probably systematically biased: when computed in standard scenarios the remaining terms are of almost exactly the same magnitude but opposite sign. But this could not be predicted in advance. 5 Baumann & McAllister later suggested that the Lyth bound [62] placed a limit on ε [63]. Lidsey & Huston [64] argued that in combination with the large value of ε implied by Alishahiha, Silverstein & Tong [33] (see also Ref. [65]) this made the "UV" version of the DBI model microscopically unviable. The UV model has other difficulties. Bean et al. [66] noted that backreaction could invalidate the probe brane approximation, spoiling inflation. For the purpose of making an estimate we are ignoring these details.
results could be required to bring the theoretical uncertainty below the observational resolution of a fourth-generation satellite such as CMBPol or CoRE. If desired, these could be computed using the Green's function formalism of Gong & Stewart [57,58], although we do not attempt this here. We find that P (X, φ) models can produce a distinctive-but not unique-bispectrum shape at next-order. A similar shape arises in certain Galileon models.
Outline.-This paper is organized as follows. In §2 we discuss the background model and specify our version of the slow-roll approximation. In §2.2 we give a short account of the twopoint function, which is required to study the squeezed limit of the bispectrum. We review the derivation of the third-order action in §3, and give an extremely brief introduction to the "in-in" formulation for expectation values in §3.1.1. We give a careful discussion of boundary terms in §3.1.2.
Next-order corrections to the bispectrum are calculated in §3.2. Readers familiar with existing calculations of the three-point function may wish to skip directly to this section, and refer to § §2-3.1.2 for our notation and definitions. We translate the next-order bispectrum to a prediction for f NL in §3.4, and study its shape and scale dependence in § §3.5-3.6. We identify a new bispectrum shape orthogonal to the lowest-order possibilities, and give cosines for the next-order shapes with standard templates. In §3.6 we obtain a general formula for the overall running of f NL with scale.
In §4 we briefly recapitulate the next-order computation for tensor fluctuations. This extends the range of observables to include the tensor-to-scalar ratio, r, and the tensor spectral index n t , which are helpful in assessing observational signatures. In § §5-6 we apply our results to a selection of concrete models. §5 concentrates on canonical single-field inflation with an arbitrary potential. The next-order corrections are small, but we are able to verify Maldacena's consistency condition to second-order in slow-roll parameters. In §6 we discuss models which can imprint a significant nongaussian signature, focusing on DBI inflation and k-inflation. For DBI inflation one can choose a late-time power-law attractor solution, or a generic quasi-de Sitter background which depends on the shape of the potential and warp factor. In §6.1.1 we calculate the next-order corrections for both cases, obtaining corrections to f NL arising from the shape of the potential and warp factor for the first time. We compare our results to known exact solutions. In §6.1.2 we apply our calculation to the power-law solution of k-inflation. Finally, we summarize our findings and conclusions in §7.
A number of appendices collect extra material. We briefly recount the derivation of next-order corrections to the propagator in Appendix A. Appendices B-C give mathematical details of certain integrals and discuss the role of analytic continuation in obtaining expressions valid for arbitrary momentum shapes. In Appendix D we provide explicit formulae for the quadratures used in Ref. [43].

Single-field inflation: an overview
Consider the theory (1.1), with a homogeneous background solution given by the Robertson-Walker metric, The corresponding Friedmann equations are where H ≡ȧ/a is the Hubble parameter, an overdot denotes differentiation with respect to cosmological time, t, and a comma denotes a partial derivative. The nontrivial kinetic structure of P causes fluctuations of the scalar field φ to propagate with phase velocitywhich we will refer to as the sound speed, c s -different to unity [67], In the special case of a canonical field, we have P = X/2 − V (φ) and the fluctuations in φ propagate at the speed of light.

Fluctuations
During inflation the universe rapidly isotropizes [68], making φ spatially homogeneous to a good approximation. At the same time, quantum fluctuations generate small perturbations whose statistics it is our intention to calculate. Since φ dominates the energy density of the universe by assumption, its fluctuations must be communicated to the metric. Our freedom to make coordinate redefinitions allows these metric and field fluctuations to be studied in a variety of gauges [69,70]. Without commitment to any particular gauge, a perturbed metric can be written in terms of the Arnowitt-Deser-Misner (ADM) lapse and shift functions [71], respectively denoted N and N i , In the absence of perturbations, Eq. (2.4) reduces to the background metric (2.1). The utility of this formulation is that N and N i do not support propagating modes, which are restricted to φ and the spatial metric h ij . There are three scalar and two tensor modes, of which one scalar can be gauged away by choice of spatial coordinates and another by choice of time. At quadratic order the tensor modes decouple from the remaining scalar. The tensors contribute at tree-level for n-point functions with n 4 [29,48], but only at loop-level 6 to the threepoint function. In what follows we work to tree-level and discard tensor modes. Whether or not they are retained, however, one uses constraint equations (discussed below) to express the lapse and shift algebraically in terms of the propagating modes. 7 Where only a single scalar field is present, it is convenient to work on spatial slices of uniform field value. In this gauge the propagating scalar mode is carried by h ij , and expresses local modulations in the expansion a(t), Writing ln a(t) = t dN , where dN = H dt and N (t) − N (t 0 ) expresses the number of e-folds of expansion between times t and t 0 , it follows that ζ = δN . 8 Slow-variation parameters.-We work perturbatively in ζ. During inflation the universe is typically well-approximated by a de Sitter epoch, in which the Hubble parameter is constant. To express the action (1.1) in terms of ζ, it is convenient to introduce dimensionless parameters which measure the variation of background quantities per Hubble time. We define Note that ε 0, provided H decreases with time. There is no requirement of principle that any of ε, η or s are small, although whenever background quantities are slowly varying we may expect ε, |η|, |s| 1.
Action and constraints.-In ADM variables, the action (1.1) can be written where E ij satisfies where | denotes a covariant derivative compatible with h ij , and indices enclosed in brackets (· · · ) are to be symmetrized. The extrinsic curvature of spatial slices, K ij , can be written Varying the action with respect to N and N i yields constraint equations [42], Eqs. (2.9a)-(2.9b) can be solved order-by-order. We write N = 1 + α. Likewise, the shift vector can be decomposed into its irrotational and solenoidal components, We expand α, β and θ perturbatively in powers of ζ, writing the terms of n th order as α n , β n and θ n . To study the three-point correlations it is only necessary to solve the constraints to first order [26,43]. One finds The quantity Σ measures the second X-derivative of P . In what follows it will be helpful to define an analogous quantity for the third derivative, denoted λ [42],

Two-point correlations
To obtain the two-point statistics of ζ one must compute the action (2.7) to second order, which was first accomplished by Garriga & Mukhanov [40]. Defining a conformal time variable by τ = t ∞ dt/a(t), one finds where z ≡ ε/c 2 s , a prime denotes differentiation with respect to τ , and ∂ represents a spatial derivative. Although our primary interest lies with P (X, φ)-type models of the form (1.1), we leave intermediate expressions in terms of z and its corresponding slow-variation parameter, The last equality applies in a P (X, φ) model. When calculating three-point functions in §3, we will do so for an arbitrary z. In order that the fluctuations are healthy and not ghostlike we must choose z positive. To simplify certain intermediate expressions, it will be necessary to have an expression for the variation δS 2 /δζ,χ where we have introduced a variable χ, which satisfies The equation of motion follows by setting δS 2 /δζ = 0 in (2.14).
Slow-variation approximation.-Although Eqs. (2.12) and (2.14)-(2.15) are exact at linear order in ζ, it is not known how to solve the equation of motion (2.14) for arbitrary backgrounds. Lidsey et al. [60] noted that the time derivative of each slow-variation parameter is proportional to a sum of products of slow-variation parameters, and therefore if we assume 0 < ε 1, |η| 1, |s| 1 and |v| 1, (2.16) and work to first order in these quantities, we may formally treat them as constants. Expanding a general background quantity such as H(t) around a reference time t , one finds where ∆N (t) = N (t) − N (t ). Similar formulae apply for the slow-variation parameters ε, η, s and v themselves, and their derivatives. Physical quantities do not depend on the arbitrary scale t , which plays a similar role to the arbitrary renormalization scale in quantum field theory. Eq. (2.17) yields an approximation to the full time evolution whenever |ε ∆N (t)| 1 [57,58,82], but fails no later than ∼ 1/ε e-folds after the reference time t . In order that solutions of (2.14) are sufficiently accurate for their intended use-to compute correlation functions for t ∼ t -we must demand that O(ε) quantities are sufficiently small that (2.17) applies for at least a few e-folds around the time of horizon crossing, but it is not usually necessary to impose more stringent restrictions. This approach will fail if any slow-variation quantity becomes temporarily large around the time of horizon exit, which may happen in "feature" models [83][84][85][86][87].
Working in an arbitrary model, results valid many e-folds after horizon exit typically require an improved formulation of perturbation theory obtained by resumming powers of ∆N [88][89][90], for which various formalisms are in use [13,[91][92][93]. However, as is well-known (and we shall see below), this difficulty does not arise for single-field inflation.
Two-point function.-The time-ordered two-point correlation function is the Feynman propagator, T ζ(τ, (2.18) The δ-function enforces conservation of three-momentum, and k = |k 1 | = |k 2 |. One finds The mode function u k is a positive frequency solution of (2.14) with δS 2 /δζ = 0. Invoking (2.16) and working to first-order in each slow-variation parameter, we find 20) in which H (2) ν is the Hankel function of the second kind of order ν, and ≡ ε + v/2 + 3s/2. At sufficiently early times, for which |kc s τ | 1, an oscillator of comoving wavenumber k cannot explore the curvature of spacetime and feels itself to be in Minkowski space. In this limit (2.20) approaches the corresponding Minkowski wavefunction [94].
Power spectrum.-When evaluated at equal times, the two-point function defines a power spectrum P (k, τ ) by the rule P (k, τ ) = G k (τ, τ ). (2.21) In principle, the power spectrum depends on time. It is conventionally denoted P (k) and should not be confused with the Lagrangian P (X, φ). Using (2.19) for τ → τ , working in the limit |kc s τ | → 0 and expanding uniformly around a reference time τ , one finds The Euler-Mascheroni constant is γ E ≈ 0.577. We have introduced a quantity k satisfying |k c s τ | = 1, where we recall that τ is arbitrary and need not be determined by k. Since |k c s τ | ≈ |k c s /a H |, we describe τ as the horizon-crossing time associated with the wavenumber k . Inclusion of next-order effects marginally shifts the time of horizon exit, leading to a small mismatch between |k c s τ | and |k c s /a H |. 9 In the terminology of Lidsey et al. [60], the "lowest-order" result is the coefficient of the square bracket [· · · ], and the "next-order" correction arises from the term it contains of first-order in slow-variation parameters. Therefore the lowest-order result can be recovered by setting the square bracket to unity. This convention was introduced in Ref. [60], and when writing explicit expressions we adopt it in the remainder of this paper. Although expanded around some reference time τ , Eq. (2.22) does not depend on ∆N and therefore becomes time-independent once the scale of wavenumber k has passed outside the horizon. This is a special property of single-field inflation. Working in classical perturbation theory it is known that ζ becomes constant in the superhorizon limit provided the fluctuations are adiabatic and the background solution is an attractor [95][96][97]. 10 We are not aware of a corresponding theorem for the correlation functions of ζ, computed according to the rules of quantum field theory. However, it appears that in all examples compatible with the classical conservation laws discussed in Refs. [95][96][97], a time-independent limit is reached. We will return to this issue in §3.3 below.
Scale dependence.-The logarithmic term in P (k) indicates that the power spectrum varies weakly with scale k, making (2.22) quantitatively reliable only if the reference time k is chosen sufficiently close to k that | ln(2k/k )| 1. Defining a "dimensionless" power spectrum P by the rule P = k 3 P (k)/2π 2 , the variation of P(k) with scale is conventionally described in terms of a spectral index, which is valid to lowest-order provided k ≈ k. Eq. (2.23) is a renormalization group equation describing the flow of P with k, where β P ≡ (n s − 1)P plays the role of the β-function. An expression for n s valid to next-order can be obtained by setting k = k , making ' ' the time of horizon exit of wavenumber k. Having made this choice, the k-dependence in (2.22) appears only through the time of evaluation ' ' and is accurate to next-order [55]. We define 9 Chen et al. [43]   and quote results in Table 1 for n s − 1 at lowest-order and next-order.

Third-order action
Three-point statistics can be obtained from the third-order action. This calculation was first given in Ref. [42], where the three-point function was obtained under certain hypotheses. Chen et al. later computed the full three-point function [43]. After integration by parts, using both the background equations of motion (2.2) and the solutions of the constraints given in (2.10), we find where ∂ denotes an integral over a formal boundary, whose role we will discuss in more detail below. We have temporarily reverted to cosmic time t, rather than the conformal time τ . The parameters Σ and λ were defined in Eqs. (2.11a) and (2.11b). After further integration by parts, and combining (2.14) and (2.10) with (2.15), one finds 11 where f is defined by  [26,42,43]. Adopting a procedure initially used by Maldacena, these calculations discarded all boundary terms, retaining only contributions proportional to δS 2 /δζ in the bulk component of (3.2). The δS 2 /δζ terms were subtracted by making a field redefinition. This procedure can be misleading. The terms proportional to δS 2 /δζ give no contribution to any Feynman graph at any order in perturbation theory, because δS 2 /δζ is zero by construction when evaluated on a propagator. Therefore these terms give nothing whether they are subtracted or not. On the other hand, a field redefinition may certainly shift the three-point correlation function. Therefore, in general, the subtraction procedure will yield correct answers only if this nonzero shift reproduces the contribution of the boundary component in (3.2), which need not be related to f . This argument was first given in Ref. [27], and later in more detail in Refs. [90,98], but was applied to the third-order action for field 11 In Refs. [42,43], a further transformation was made to rewrite the terms proportional to η. Using the field equation (2.14) and integrating by parts, these can be consolidated into the coefficient of a new operator ζ 2ζ -which does not appear in (3.2)-together with corresponding new contributions to f and the boundary term. In this paper, we will leave the action as in (3.2) for the following reasons. First, when computing the three-point correlation function, the contribution of each operator must be obtained separately. Therefore nothing is gained by introducing an extra operator, whose contribution we can avoid calculating by working with (3.2). Second, after making the transformation, the contribution from the boundary term is nonzero and must be accommodated by making a field redefinition. This redefinition must eventually be reversed to obtain the correlation functions of the physical field ζ. If we leave the action as in (3.2) then it transpires that no field redefinition is necessary. fluctuations in the spatially flat gauge. In this gauge only a few integrations by parts are required. The boundary term is not complicated and the subtraction procedure works as intended. In the present case, however, it appears impossible that the subtraction procedure could be correct, because the boundary term contains operators such as ζ 3 which are not present in f . Indeed, because ζ approaches a constant at late times, the ζ 3 term apparently leads to a catastrophic divergence which should manifest itself as a rapidly evolving contribution to the three-point function outside the horizon.
This potential problem can be seen most clearly after making the redefinition ζ → π − f under which the quadratic action transforms according to The bulk term proportional to δS 2 /δζ disappears by construction. After the transformation, the boundary term becomes in which χ is to be interpreted as a function of π [cf. Eq. (2.15)]. Eq. (3.5) is not zero. To satisfy ourselves that it does not spoil the conclusions of Refs. [26,42,43], we must determine how it contributes to the three-point correlation function. Before doing so, we briefly describe the in-in formalism which is required. Readers familiar with this technique may wish to skip to §3.1.2.

Schwinger's in-in formalism
The correlation functions of interest are equal time expectation values taken in the state corresponding to the vacuum at past infinity. At later times, persistent nontrivial correlations exist owing to gravitational effects associated with the time-dependent background of de Sitter.
Feynman's path integral computes the overlap between two states separated by a finite time interval, which is taken to infinity in scattering calculations. Schwinger obtained expectation values at a finite time t * by inserting a complete set of states |i, t f at an arbitrary time where |in is the "in" vacuum in which one wishes to compute the expectation value and O is an arbitrary local functional. Our choice of t f is irrelevant. Choosing a basis of energy eigenstates, |i, This phase cancels in (3.6). 12 Expressing each overlap as a Feynman path integral, we obtain [99][100][101][102] 12 One could just as well insert a complete set of states at an arbitrary time t < t * , but the resulting overlap would not be expressible in terms of a path integral. Had we retained t f > t * , the resulting contributions would have yielded only cancelling phases in (3.7).
The δ-function restricts the domain of integration to fields φ + and φ − which agree at time t * , but are unrestricted at past infinity. The requisite overlap of an arbitrary field configuration with the vacuum is obtained by deforming the contour of time integration. Cosmological applications of Schwinger's formulation were considered by Jordan [103] and by Calzetta & Hu [104], to which we refer for further details. Applications to inflationary correlation functions were discussed by Weinberg [73,74] and have been reviewed elsewhere [105,106].

Removal of boundary terms
If boundary terms are present they appear in (3.7) as part of the action S with support at past infinity and t = t * . The deformed contour of integration kills any contribution from past infinity, leaving a boundary term evaluated precisely at t * , where the δ-function constrains the fields to agree. Therefore, at least for O containing fields but not time derivatives of fields, any boundary operators which do not involve time derivatives produce only a phase which cancels between the + and − contours. 13 This cancellation is a special property of the in-in formulation: it would not occur when calculating in-out amplitudes, for which the uncancelled boundary term would diverge near future infinity. Rapid oscillations of e iS induced by this divergence would damp the path integral, yielding an amplitude for any scattering process which is formally zero. This can be understood as a consequence of the lack of an S-matrix in de Sitter space [107].
In virtue of this cancellation we may disregard the first three operators in the boundary part of (3.5). However, the δ-function in (3.7) in no way requires that time derivatives of the + and − fields are related at t = t * . Therefore operators involving time derivatives need not reduce to cancelling phases. To understand their significance we subtract them using a further field redefinition.
Inspection of Eq. (3.5) shows that the time-derivative terms are of the schematic form ππ 2 , and therefore lead to a field redefinition of the form ζ → π + ππ. We now argue that boundary operators with two or more time derivatives are irrelevant on superhorizon scales. Using the schematic field redefinition, the three-point correlation functions of ζ and π are related by ζ 3 = π 3 + 3 π 2 ππ plus higher-order contributions. However, Eq. (2.22) implies ππ → 0 on superhorizon scales, and therefore ζ 3 = π 3 up to a decaying mode. This field redefinition will inevitably produce bulk terms proportional to δS 2 /δζ, but we have already seen that these do not contribute to Feynman diagrams at any order in perturbation theory. Therefore, on superhorizon scales, the correlation functions of the original and redefined fields agree. It follows that after subtraction by a field redefinition, the unwanted boundary terms in (3.5) are irrelevant and can be ignored. Equivalently, one may confirm this conclusion by checking that operators with two or more time derivatives give convergent contributions to the boundary action at late times. Similar arguments apply for any higher-derivative combination.
We conclude that the only non-negligible field redefinitions are of the schematic form ζ → π+π 2 , which arise from boundary operators containing a single time derivative. Eq. (3.5) contains no such operators. However, had we made a further transformation to consolidate the η dependence, a contribution of this form would have appeared in f and the boundary action. It can be checked that this single-derivative term would be correctly subtracted by Maldacena's procedure, and in this case the subtraction method applied in Refs. [26,42,43] yields the correct answer. However, in theories where single-derivative terms already appear in the boundary component of (3.1) there seems no guarantee it will continue to do so.

The bispectrum beyond lowest-order
We define the bispectrum, B, in terms of the three-point function, Observational constraints are typically quoted in terms of the reduced bispectrum, f NL , which satisfies [8,13] f NL ≡ 5 6 .
Current constraints on f NL in the simplest inflationary models have been discussed by Senatore et al. [108].
Next-order corrections.-We are now in a position to compute the three-point function of (3.2) to next-order. The arguments of the previous section show that the boundary action and δS 2 /δζ contributions can be discarded. Next-order terms arise in the remaining operators from a variety of sources. First, the coefficients of each vertex in (3.2) contain a mixture of lowest-order and next-order contributions. Second, the lowest-order part of each vertex is a time-dependent quantity which must be expanded around a reference time, as in (2.17), producing next-order terms. Third, there are next-order corrections to the propagator, obtained by expanding (2.19)- (2.20) in the neighbourhood of the chosen reference time. Propagator corrections appear on both the external and internal legs of the diagram.
Reference time, factorization scale.-To proceed, we must choose a reference point τ around which to expand time-dependent quantities. Consider an arbitrary correlation function of fields ζ(k i ). Whatever our choice of τ , the result (2.22) for the power spectrum shows that we must expect logarithms of the form ln k i /k which account for the difference in time of horizon exit between the mode k i and the reference wavenumber k . To obtain a reliable answer we should attempt to minimize these logarithms. If all fields participating in the correlation function carry momenta of approximately common magnitude k i ∼ k-described as the "equilateral limit"-the logarithm will be small when k ∼ k. In this case, naïve perturbation theory is not spoiled by the appearance of large logarithms. In the opposite limit, one or more fields have "soft" momenta of order k IR which are much smaller than the remaining "hard" momenta of order k UV . When k IR /k UV → 0 it will not be possible to find a choice of k which keeps all logarithms small and the calculation passes outside the validity of ordinary perturbation theory. We have encountered the problem of large logarithms which led to the renormalization group of Gell-Mann & Low [109].
In the study of inflationary correlation functions, configurations mixing hard and soft momenta with k IR k UV are referred to as "squeezed," and are of significance because they dominate the bispectrum for canonical inflation [26]. In principle one could study the behaviour of a correlation function as its momenta are squeezed by setting up an appropriate renormalization group analysis [110]. But this is more complicated than necessary. Maldacena argued that, as the momentum carried by one operator becomes soft, the three-point function would factorize: it can be written as a "hard subprocess," described by the two-point correlation between the remaining hard operators on a background created by the soft operator [26]. Factorization of this kind is typical in the infrared dynamics of gauge theories such as QCD, where it plays an important role in extracting observational predictions. The various factorization theorems for QCD correlation functions have been comprehensively reviewed by Collins, Soper & Sterman [111]. 14 Maldacena's argument was later generalized by Creminelli et al. [112]. The factorization property can be exhibited by an explicit decomposition of the field into hard and soft modes [113,114].
Because the squeezed limit can be described by Maldacena's method, the outcome of this discussion is that the reference scale should usually be chosen to minimize the logarithms when all momenta are comparable. In the remainder of this paper we quote results for arbitrary k , but frequently adopt the symmetric choice k = k 1 + k 2 + k 3 where numerical results are required. 15 Having done so, we will be formally unable to describe the squeezed limit. Nevertheless, because there is no other scale in the problem, our results must be compatible with the onset of factorization in appropriate circumstances-a property usually referred to as Maldacena's consistency relation. We will see below that this constitutes a nontrivial check on the correctness of our calculation; see also Renaux-Petel [114] for a recent discussion of Maldacena's condition in the case of P (X, φ) models.
Operators.-To simplify our notation, we rewrite the cubic action (3.2) as (3.10) In a P (X, φ) model the interaction vertices are but our calculation will apply for arbitrary g i . Although ζ is dimensionless, it is helpful for power-counting purposes to think of it as a field of engineering dimension [mass], obtained after division by the Hubble rate H. In this counting scheme, the ζ 3 operator is dimension-6, 14 The background created by soft modes is typically described by some version of the DGLAP (or Altarelli-Parisi) equation. A similar phenomenon seems to occur in the inflationary case [110]. Equally, the separate universe method can be thought of as a factorization theorem for secular time-dependent logarithms ∼ ln |kcsτ |. The δN rules which translate correlation functions of the field perturbations into correlation functions of ζ are an important special case. In this sense, factorization is as important in extracting observable quantities for inflation as it is for QCD. 15 In the gauge theory language discussed above, the scale k can be thought of as the factorization scale.
Operators carrying momentum k k should not be included as part of the hard subprocess, but factorized into the background.
whereas the remaining four operators are dimension-5. At low energies one would naïvely expect the dimension-6 operator to be irrelevant in comparison to those of dimension-5. However, the dimension-5 operators are suppressed by the scale H making all contributions equally relevant. This manifests itself as an extra power of H in the denominator of g 1 .
The vertex factors g i are themselves time-dependent background quantities. We define slow-variation parameters h i which measure their rate of change per e-fold, and take these to be O(ε) in the slow-variation approximation.

Three-point correlations
We use these conventions to compute the next-order bispectra for each operator in (3.10).
The resulting three-point functions are complicated objects, and when quoting their values it is helpful to adopt an organizing principle. We divide the possible contributions into broadly similar classes. In the first class, labelled 'a,' we collect (i) the lowest-order bispectrum; (ii) effects arising from corrections to the wavefunctions associated with external lines; and (iii) effects arising from the vertex corrections. In the second class, labelled 'b,' we restrict attention to effects arising from wavefunctions associated with internal lines. These are qualitatively different in character because wavefunctions associated with the internal lines are integrated over time. Adapting terminology from particle physics, we occasionally refer to the lowest-order bispectrum as the "LO" part, and the next-order piece as the "NLO" part.
Large logarithms, infrared singularities.-The computation of inflationary n-point functions has been reviewed by Chen [115] and Koyama [106]. At least three species of large logarithms appear, disrupting ordinary perturbation theory. We carefully track the contribution from each species. The most familiar types-already encountered in the two-point functionmeasure time-and scale-dependence. A third type of large logarithm is associated with the far infrared limit k IR /k UV → 0 discussed in §3.2. This is Maldacena's "squeezed" limit, discussed in §3.2, in which the behaviour of the three-point function obeys a factorization principle. We will show that the various large logarithms arrange themselves in such a way that they can be absorbed into the scale-dependence of background quantities. Time-dependent logarithms appear after expanding background quantities near a fixed reference scale, as in (2.17), where at conformal time τ we have N = ln |k c s τ |. In §2.2 we explained that the correlation functions of ζ are expected to become time-independent outside the horizon. Therefore one should expect all ln τ dependence to disappear. Some Ntype logarithms cancel among themselves but others cancel with time-dependent logarithms arising from wavefunction corrections associated with internal lines. The internal lines are aware only of the intrinsic geometrical scale k t ≡ k 1 + k 2 + k 3 and cannot depend on the arbitrary reference scale k , so the outcome of such a cancellation leaves a residue of the form ln k t /k . These are scaling logarithms, entirely analogous to the logarithm of (2.22), describing variation of the three-point function with the geometrical scale k t . Scale logarithms can also occur in the form ln k i /k . The third species of logarithm takes the form ln k i /k t . Each side of the triangle must scale linearly with the perimeter, so despite appearances these have no dependence on k tthey are unaffected by rigid rescalings of the momentum triangle (cf. Eqs. (3.23a)-(3.23c) below). We describe them as 'purely' shape dependent. The 'pure' shape logarithms become large in the squeezed limit k i /k t → 0.
a-type bispectrum.-Collecting the a-type contributions to the bispectrum, we find and Z a (k 1 ) are functions of all three momenta k i and are symmetric under the exchange k 2 ↔ k 3 . We adopt the convention, used through the remainder of this paper, of writing only the asymmetric momentum explicitly. The notation 'cyclic permutations' denotes addition of the preceding term under cyclic permutations of the k i . The result is symmetric under interchange of any two momenta. We give explicit expressions for the coefficient functions in Table 2. The quantity E is a combination of slow-variation parameters, E = (2 − γ E − ln 2) − ε − s, and also appears in the power spectrum (2.22). The term proportional to X a (k 1 ) includes the entire lowest-order bispectrum.
b-type bispectrum.-The b-type bispectrum must be added to the a-type terms. It has no lowest-order contributions, and can be written (3.14) The same convention applies to the arguments of the coefficient functions We give explicit expressions in Table 3. Eq. (3.14) depends on three logarithmic functions J i (which are not Bessel functions) defined by Table 2. Coefficients of the leading order bispectrum.
where ϑ i = 1 − 2k i /k t . These exhaust the 'pure' shape logarithms of the form ln k i /k t , discussed above, which appear only in the J i . There is an obvious logarithmic divergence in the squeezed limit k i → 0, which we will show to be responsible for factorization of the correlation function. There is potentially a power-law divergence in the limit k t → 2k i . This is also a squeezed limit-in which the i th side stays fixed while a different momentum goes to zero. In this limit ϑ i → 0, making the J i naïvely divergent. If present, such power-law divergences would be puzzling. However, it can be checked that-in combination with the logarithm-each J i is finite. This infrared-safe behaviour relies on a resummation procedure which is discussed in more detail in Appendix B.
The function J satisfies This function is discussed in Appendix C. The coefficients γ 0 , γ 1 , γ 2 , γ 3 , δ 0 , δ 1 , δ 2 and δ 3 depend on the operator under consideration. We quote values for each operator in Table 4. Table 3. Coefficients of the subleading corrections to the bispectrum.
The operators ζζ 2 , ζ ∂ j ζ∂ −2 ∂ j ζ and ∂ 2 ζ(∂ −2 ∂ j ζ ) 2 are all dimension-5, and differ only in the arrangement of spatial gradients. For arbitrary shapes their three-point functions will not coincide, but for equilateral triangles the arrangement of gradients is irrelevant and the resulting f NL should agree. This will represent a minimal check of our expressions. We will carry out further checks in §3.6 and § §5-6.
A subset of these terms were calculated by Chen et al. [43]. Our calculations exhibit two principal differences. First, Chen et al. worked to fixed order in slow-roll quantities, keeping terms of O( ) only. In a model where c s 1 this gives the next-order corrections. However, in a model where c s ∼ 1 the leading terms are themselves O( ) and the formulae of Chen et al. reduce to these leading contributions. In our calculation, we work uniformly to next-order rather than a fixed order in powers of . When c s 1 our next-order corrections are O( ), and we have verified that they agree with those computed in Ref. [43]. (We give more details of the relation between our calculations in Appendix D.) When c s ∼ 1 the next-order corrections are O( 2 ). These were not included in the formulae of Ref. [43].
Second, we retain a floating reference scale k . In Appendix B 16 of Ref. [43] this was cho- Table 4. Coefficients appearing in the function J for each operators. Note that the γ i contain complex numbers. The imaginary part is cancelled on addition of the + and − Feynman diagrams, and only the real part of these coefficients contribute. In an intermediate step for the three-point function of ζ(∂ζ) 2 , the cancellation of power-law divergences in the conformal time τ (which is required by Weinberg's theorem [73]) depends on a real contribution generated from the product of two imaginary terms.
sen to be k = k t . Retaining this scale allows us to extract the scale-and shape-dependence of f NL ( § §3.5-3.6).

Formulae for f NL
The individual bispectra, with their detailed shape-dependence, are the principal observable objects. However, for simple model comparisons it is helpful to have an explicit expression for the nonlinearity parameter f NL defined in Eq. (3.9). Accounting for scale-dependent logarithms present in the power spectrum, one finds This expression is to be expanded uniformly to O(ε) in slow-variation parameters.
There is another reason to study f NL . We have explained that large logarithms of the form ln k i /k or ln k i /k t are to be expected in the squeezed limit k i → 0, describing variation of the bispectrum with shape. The power spectrum P (k) contains similar large logarithms. Since copies of the power spectrum must be factored out to obtain f NL , one may expect it to be more regular in the squeezed limit. Indeed, a stronger statement is possible. Partitioning the momenta into a single soft mode of order k IR and two hard modes of order k UV , Maldacena's consistency condition requires [26] as k IR → 0, where the right-hand side is to be evaluated at horizon exit for the mode of wavenumber k UV . Eq. (3.18) is finite and independent of any logarithms associated with the limit k IR → 0, which is why this behaviour is described as factorization. It imposes the nontrivial requirement that all large logarithms can be absorbed into P (k IR ). Such logarithms are subtracted by the denominator of (3.17), making f NL finite.
For each operator i, we write the corresponding f NL as f NLi and quote it in the form In Tables 5 and 6 we give explicit expressions for the coefficient functions f NL | i0 and κ i in the case of equilateral 17 and squeezed triangles. Table 6 confirms that (3.19) is finite in the squeezed limit, as required. In the equilateral case, we find that the operators ζζ 2 , ζ ∂ j ζ∂ −2 ∂ j ζ and ∂ 2 ζ(∂ −2 ∂ j ζ ) 2 agree, for the reasons explained above.

Shape dependence
Stewart & Lyth's interest in next-order corrections to the power spectrum lay in an accurate estimate of its amplitude. In comparison, next-order corrections to the bispectrum could be relevant in at least two ways. First, they could change the amplitude of three-point correlations, as for the power spectrum. Second, they could lead to the appearance of new "shapes," by which is meant the momentum dependence of B(k 1 , k 2 , k 3 ) [118], defined in (3.8).
In principle, both these effects are measurable.

Inner product and cosine
Babich et al. introduced a formal "cosine" which may be used as a measure of similarity in shape between different bispectra [118]. Adopting Eq. (2.21) for the power spectrum P (k), one defines an inner product between two bispectra B 1 and B 2 as

20)
17 The quoted quantity is fNL(k, k, k), which is not the same object as f equi NL for which constraints are typically derived from data [6,108]. To obtain f equi NL , one should take an appropriately normalized inner product (see §3.5 for a simple example) between the full next-order bispectrum and the equilateral template [116,117]. At this level of precision, it may even be desirable to include experiment-dependent information in the inner product.
where the sum is to be taken over all triangular configurations of the k i . The cosine between B 1 and B 2 is These expressions require some care. In certain cases the result may be infinite, requiring the summation to be regulated.
Inner product.-We define the sum over triangles as an integral over triangular configurations in a flat measure, so . It is sometimes useful to introduce a more complicated measure, perhaps to model observational effects [119]. In this paper we retain the flat measure for simplicity. The δ-function can be integrated out immediately, leaving a space parametrized by two vectors forming a planar triangle which we choose to be k 1 and k 2 . The triangle is invariant under a group SO(2) × U(1), representing arbitrary rotations of k 1 combined with azimuthal rotations of k 2 ; these change our representation of the triangle but not its intrinsic geometry. The volume of this group may be factored out of the measure and discarded. Reintroducing k 3 in favour of the remaining angular integration, we conclude The k i can be parametrized geometrically in terms of the perimeter, k t , and two dimensionless ratios. We adopt the parametrization of Fergusson & Shellard [119], where 0 β 1 and β − 1 α 1 − β. The measure dk 1 dk 2 dk 3 is proportional to k 2 t dk t dα dβ. Also, on dimensional grounds, each bispectrum B i scales likeB i k −6 t , wherẽ B i is dimensionless, and each power spectrum P scales likeP k −3 t whereP is dimensionless. In the special case of scale-invariance,P is constant and theB i depend only on α and β. Therefore where N is a harmless infinite normalization which can be divided out. With this understanding we use (3.24) to determine the cosine of Eq. (3.21). In practice, our bispectra are not scale invariant and therefore (3.24) does not strictly apply. However, the violations of scale invariance (to be studied in §3.6 below) are small.
Divergences.-Eq. (3.24) may be infinite. For example, the well-studied local bispectrum diverges like (1 + α + β) −2 or (1 − α + β) −2 in the limit β → 0, α → ±1, or like (1 − β) −2 in the limit β → 1, α → 0 [118]. These correspond to the squeezed limits discussed in §3.2. Eq. (3.24) therefore exhibits power-law divergences on the boundaries of the region of integration, and in such cases the integral must be regulated to obtain a finite answer. For simplicity, we adopt a sharp cutoff which requires k i /k t > δ min . As δ min → 0 the cosine (3.21) may converge to a nonzero limit if B 1 · B 2 , B 1 · B 1 and B 2 · B 2 diverge at the same rate. Otherwise, except in finely-tuned cases, it converges to zero. For this reason, where divergences exist, the value assigned to cos(B 1 , B 2 ) is largely a matter of convention. However, to resolve the practical question of whether two shapes can be distinguished by observation it should be remembered that experiments cannot measure arbitrarily small wavenumbers. Therefore their ability to distinguish shapes peaking in the squeezed limit is limited. In this case, to obtain accurate forecasts of what can be distinguished, one should restore the k t -dependence in (3.24) and restrict the integration to observable wavenumbers, yielding a manifestly finite answer [119,120].

Bispectrum shapes from slow-variation parameters
In a model with arbitrary g i , the bispectrum is a linear combination of the shapes produced by the five operators in (3.10). Of these, ζζ 2 and ζ(∂ζ) 2 are predominantly correlated with the local template and the remainder correlate strongly with the equilateral template. The ζ 3 operator has some overlap with the enfolded template, yielding a cosine of order 0.75.
With generic values of the slow-variation parameters the situation at next-order is similar, and each next-order shape is largely correlated with its parent lowest-order shape.
Lowest-order shapes.-A P (X, φ) is not generic in this sense, but imposes strong correlations among the g i . At lowest order g 4 and g 5 do not contribute. We focus on a model with small sound speed, in which next-order corrections are most likely to be observable, and retain only contributions enhanced by c −2 s . The remaining three operators organize themselves into a family of shapes of the form S 1 + αS 2 , where S 2 arises only from ζ 3 but S 1 is a linear combination of the shapes produced by ζ 3 , ζζ 2 and ζ(∂ζ) 2 . The parameter α is the enhanced part of λ/Σ, that is In the DBI model α = 1/2. We plot the shapes S 1 and S 2 in Table 9. Note that although S 1 involves a linear combination of the local-shape operators ζζ 2 and ζ(∂ζ) 2 , the P (X, φ) Lagrangian correlates their amplitudes in such a way that there is no divergence in the squeezed limit. Both S 1 and S 2 are strongly correlated with the equilateral template. They are similar to the M 1 -and M 2 -shapes studied in a Galileon theory by Creminelli et al. [37].
Next-order shapes. 18 -At next-order, more shapes are available. Naïvely, the family of enhanced bispectra is labelled by ε, η, s and also (following Chen et al. we define =λ/Hλ [43]). In practice there is some degeneracy, because the shapes corresponding to these independent parameters may be strongly correlated. We will see these degeneracies emerge naturally from our analysis.
The c −2 s -enhanced next-order shape can be written as a linear-combination of shapes proportional to the ε, η, s and parameters, modulated by α, εS ε + ηS η + sS s + α εS ε + ηS η + sS s + S . (3.26) We give overlap cosines of the S i , S i with the standard templates in Table 7 and plot their shapes in Table 10. Because the cutoff dependence complicates comparison between different analyses we list the cosines between templates in Table 8, computed using the same conventions. Generally speaking, these shapes have strong overlaps with the equilateral template. However, two are quite different in appearance and have a slightly smaller cosine ∼ 0.85 with this mode: these are S ε and S s . We fix two coefficients in (3.26) by choosing a linear combination orthogonal to both S 1 and S 2 . Without loss of generality we can choose these to be η and s. We find the required combination to be approximately It is possible this procedure is stronger than necessary. Both S 1 and S 2 are correlated with the equilateral template, and it may be sufficient to find a linear combination orthogonal to that. In what follows, however, we insist on orthogonality with S 1 and S 2 and defer generalizations to future work. For certain values of α the denominator of both η and s may simultaneously vanish, making the required η and s very large. This implies that, near these values of α, no shape orthogonal to both S 1 and S 2 can be found within the validity of next-order perturbation theory. Therefore we restrict attention to those α which allow acceptably small η and s. This process leaves two linear combinations proportional to ε and . In principle these can be diagonalized, yielding a pair of shapes orthogonal to each other and {S 1 , S 2 }. However, the 2×2 matrix of inner products between these linear combinations is degenerate. Therefore, only one member of this pair is physical and can be realized in a P (X, φ) model. The other is not: it has zero inner product with (3.26), and is impossible to realize because of enforced correlations between coefficients. We denote the physical orthogonal combination O. It has a vanishing component proportional to . This was expected, because the shape S is the same as S 2 . For this reason, Chen et al. [43] absorbed into a redefined λ/Σ. It is indistinguishable from the lowest-order prediction and could never be observed separately, which is the origin of the degeneracy. We could have arrived at the same O by excluding S from (3.26). Demanding the inner product with S 1 and S 2 be zero reproduces the physical linear combination obtained from diagonalization.
We plot the shape of O in Table 11. Its dependence on α is modest. As a function of the k i there are multiple peaks, and therefore O is not maximized on a unique type of triangle. In Table 13 we give the overlap cosine with common templates. The lowest-order shapes S 1 and S 2 are strongly correlated with the equilateral template, and since O is orthogonal to these by construction it also has small cosines with the equilateral template, of order 10 −2 . There is a moderate cosine with the local template of order ∼ 0.3 -0.4. The precise value depends on our choice of δ min , but the dependence is not dramatic. In Table 13 we have used our convention δ min = 10 −3 . For δ min = 10 −5 the local cosines change by roughly 25%. Overlaps with the remaining templates are stable under changes of δ min . There is a cosine of order 0.35 -0.40 with the orthogonal template, and of order 0.30 -0.35 with the enfolded template. We conclude that O is not strongly correlated with any of the standard templates used in CMB analysis. To find subleading effects in the data, it will probably be necessary to develop a dedicated template for the purpose.
The O-shape is very similar to a highly orthogonal shape constructed by Creminelli et al. [37] in a Galileon model, although O contains marginally more fine structure. For comparison, we plot the Creminelli et al. shape in Table 12 and include its cosine with O in Table 13. For varying α we find a cosine in the range 0.8 -0.9, which indicates it would be difficult to distinguish these shapes observationally. In particular, even if a bispectrum with this shape were to be detected, further information would be required to distinguish between candidate P (X, φ) or Galileon models for its origin.

Scale dependence
In this section we use the logarithms ln k i /k and ln k t /k to study the scale-dependence of the three-point function. In the squeezed limit this is determined by (3.18). The only scale which survives is the common hard momentum k UV , and the variation of f NL with this scale is determined by the variation of n s − 1. This is typically called the running of the scalar spectral index [122], and leads to a further consistency relation inherited from Maldacena'sand, in general, a hierarchy of such consistency equations generated by taking an arbitrary number of derivatives. In the case of single-field canonical inflation, discussed in §5 below, we are able to verify this explicitly.
Away from the squeezed limit, deformations of the momentum triangle may change either its shape or scale. Scale dependence occurs in even the simplest models for the same reason that the spectrum P and spectral index n s depend on scale [13,26,27,42,43]. Chen introduced a 'tilt,' n f NL , defined by 19 [123] For a fixed triangular shape, this measures changes in f NL as the perimeter varies. Scale dependence of this type was subsequently studied by several authors [124,125]. Observational constraints have been determined by Sefusatti et al. [126]. Byrnes et al. performed a similar analysis in the special case of multiple-field models producing a local bispectrum [127,128]. They allowed for deformations of the momentum triangle including a change of shape, but found these to be less important than rescalings of k t . Shape dependence is often substantially more complicated than scale dependence. Eq. (3.9) makes f NL dimensionless, but contains both powers and logarithms of the k i . The powers occur as dimensionless ratios in which k t divides out, but the shape dependence remains. The argument of each logarithm is also a dimensionless ratio, but an extra scale is available: 19 Chen implicitly worked in the equilateral limit ki = k, where kt = 3k and d ln kt = d ln k. We are defining n f NL to be the variation of fNL with perimeter for an arbitrary triangle if the shape is kept fixed. the reference scale k . When present, this gives rise to the scaling logarithms ln k i /k and ln k t /k described above, which depend on k t as well as the shape. It follows that a simple way to track the k t -dependence of f NL is to study the k -logarithms, yielding the identity For a general model, n f NL can be written where the f i (k i ) functions are dimensionless ratios of polynomials in the k i which are listed in Table 14.
Squeezed and equilateral limits.-In the equilateral limit, we find The squeezed limit gives a simple result, We define the running of the spectral index, α s , by [122] α s = d(n s − 1) d ln k . (3.33) Compatibility with (3.18) in the squeezed limit requires In §5 we will verify this relation in the special case of canonical single-field inflation.

Tensor modes
Inflation will inevitably produce tensor fluctuations to accompany the scalar fluctuation ζ.
Detection of the B-mode polarization signal produced by these fluctuations is a major aim of the Planck satellite and future CMB experiments. If it can be measured, this signal will provide important constraints on the energy scale of inflation.
In certain models the tensor sector provides sufficient observables to allow one or more quantities, such as f NL , to be written in terms of other observables. In the inflationary literature such relationships are typically known as consistency relations, and were introduced by Copeland et al. [129,130]. In the language of particle physics they are "observables in terms of observables"-predictions which are independent of how we parametrize the theory, and which the renormalization programme has taught us represent the physical content of any quantum field theory. Such consistency relations represent important tests of entire classes of models. To be used effectively with the next-order results of this paper we will require next-order predictions for the tensor modes. These were obtained by Stewart & Lyth [55], and are unchanged by the noncanonical action (1.1).
In this section our aim is to obtain the next-order consistency relation, Eq. (4.6). The tensor fluctuation is a propagating spin-2 mode which belongs to the ADM field h ij of (2.4). We write h ij = a 2 e 2ζ (e γ ) ij , where tr γ ij = 0. At quadratic order, the action is [131] There are two polarizations, traditionally denoted '+' and '×,' making γ ij transverse in the sense ∂ i γ ij = 0. Introducing a reference scale k and adding the power in each polarization incoherently, the resulting dimensionless spectrum can be written [55] This is the sum of two copies of the power spectrum for a massless scalar field with c s = 1, and is conserved on superhorizon scales. Including next-order corrections, the scale dependence of P g is measured by the tilt n t , It is conventional to measure the amplitude of tensor fluctuations relative to ζ. One defines the tensor-to-scalar ratio, r, by the rule [130] r ≡ P g P , (4.4) where P is the dimensionless version of Eq. (2.22). We find r 16ε c s 1 − 2η + (s + η ) γ E + ln 2k k . (4.5) In canonical models, r can be written purely in terms of observable quantities. In the noncanonical case this is not automatically possible without the addition of new observables. In general, One may use the lowest-order result for n t to eliminate ε. To eliminate η would require the scalar spectral index, n s . It is possible to use f NL to rewrite c s in the prefactor [65], but in doing so one introduces dependence on the parameter . Therefore at least two extra observables would be required to eliminate the dependence on s and . If these depend on t, ξ or similar parameters, then further observables could be required. We conclude that at next-order, for a general P (X, φ) Lagrangian, the observables {r, n s , n t , f NL } do not form a closed set.

Canonical single-field inflation
The simplest model of inflation comprises a single scalar field with canonical kinetic terms. Maldacena showed that the fluctuations in this model are almost Gaussian, with f NL of order r [26]. This is unobservably small. In a canonical model, c s = 1.
-To calculate f NL we require the flow parameters h i , which measure time dependence in the vertex factors g i . These are In this model the time-dependence of z is described by v = η. Collecting contributions from where we have used the numerical constant ω = coth −1 5, and k should be regarded as the common momentum scale, k i = k. The squeezed limit may be recovered from Table 6. We find where k should now be regarded as the scale of the hard momenta in the correlation function. In §3.4 and §3. 6 we emphasized that f NL should be finite in this limit, containing no large logarithms, because these factorize into the power spectrum and are subtracted. The remaining logarithms [the ln 2k/k terms in (5. 3)] track the dependence of f NL on the hard scale, and will be studied below. Using (3.18) and comparing with the spectral indices quoted in Table 1, it is easy to check that our formula correctly reproduces the Maldacena limit. We note that, strictly, one should regard agreement in this limit as an accident which happens because the simple slow-roll model contains no other scale which could interfere with factorization of the correlation function.
Scale dependence of f NL .-Specializing to the equilateral limit of n f NL , we find In the squeezed limit one obtains which correctly describes the running of the scalar spectral index, α s , in agreement with (3.34). The consistency conditions (5.3) and (5.5) represent a nontrivial check on the correctness of our calculation. In particular, throughout the calculation we have cleanly separated the conceptually different scales k t and k . Therefore the correct formula (5.5) is not simply a consequence of obtaining the correct lowest-order terms in (5.3)-although, as described in §2.2, it can be obtained from these by differentiating with respect to ln k after setting k = k t .

Non-canonical single-field inflation
For the noncanonical action (1.1), the h i can be written where we have defined ξ =η/Hη, t =ṡ/Hs.
In the canonical case, it was possible to verify Maldacena's consistency condition to nextorder. In the noncanonical case this is not possible without a next-next-order calculation, because for c s = 1 the leading contribution to f NL is O(1) in the slow-variation expansion. Therefore our calculation of subleading corrections produces a result valid to O(ε), which is short of the O(ε 2 ) accuracy required to verify the consistency condition at next-order. Chen et al. gave the subleading corrections in terms of undetermined integrals [43]. Expanding these asymptotically, they argued that the consistency relation would be satisfied at lowestorder. More recently, Renaux-Petel [114] gave an equivalent demonstration. Here, we have knowledge of the full bispectrum to subleading order. Using (6.1), it can be verified that in the squeezed limit, and expanding around a reference scale k , One may check that this agrees with Eq. (3.18) and Table 1. We expect n f NL = O(ε 2 ), and therefore a next-next-order calculation is required to estimate the running of f NL in noncanonical models.

Asymptotically power-law models
Power-law inflationary models were introduced by Lucchin & Matarrese [132,133], who studied potentials producing an expansion history of the form a(t) ∝ t 1/ε . The exponent 1/ε is the usual parameter ε = −Ḣ/H 2 . It need not be small, but should be taken as constant which makes η = ξ = 0. The solution is inflating provided ε < 1. Exact solutions can be found in the canonical case, which form the basis of the next-order calculation [60]. In this section we study two examples which are asymptotically described by noncanonical power-law inflation at late times. The first is Dirac-Born-Infeld ("DBI") inflation, which produces a scale-invariant power spectrum at lowest-order. Departures from scale invariance appear at next-order. These properties imply that we can compare our results to a formula of Khoury & Piazza which was obtained without invoking the slow-roll approximation [134]. Our second example is k-inflation, for which the power spectrum is not scale invariant at lowest-order, and to which Khoury & Piazza's result does not apply.

Dirac-Born-Infeld inflation
The DBI action is a low-energy effective theory which describes a D3-brane moving in a warped throat. It was proposed as a model of inflation by Silverstein & Tong [61], and subsequently developed by Alishahiha et al. [33]. The action is of the form (1.1), with P (X, φ) satisfying where f is an arbitrary function of φ known as the warp factor, and V (φ) is a potential arising from couplings between the brane and other degrees of freedom. The DBI Lagrangian is algebraically special [65,125,135] and enjoys a number of remarkable properties, including a form of nonrenormalization theorem [136][137][138]. In principle non-minimal curvature couplings can be present, of the form Rφ 2 , which spoil inflation [139]. This gives a form of the ηproblem, and we assume such terms to be negligible. Eq. (6.3) makes 2λ/Σ = (1 − c 2 s )/c 2 s , which requires g 1 → 0 but causes the denominator of h 1 in (6.1) to diverge. Only the finite combination g 1 h 1 appears in physical quantities, and it can be checked that g 1 h 1 → 0 as required. The square root in (6.3) must be real, giving a dynamical speed limit for φ. It is conventional to define a Lorentz factor When γ ∼ 1 the motion is nonrelativistic. When γ 1, the brane is moving close to the speed limit. The Lorentz factor is related to the speed of sound by c s = γ −1 .
Silverstein & Tong [61] argued that (6.3) supported attractor solutions described at late times by power-law inflation. In this limit, the slow-variation parameters ε and s are constant, with η = ξ = t = 0 but not zero. Variation of the sound speed gives s = −2ε, making = 0 and yielding scale-invariant fluctuations at lowest-order [cf. (2.23)]. In the equilateral limit, 20 we find that f NL satisfies In §1 we estimated the relative uncertainty in f NL to be ∼ 14ε, working in the limit γ 1, based on O(ε) terms from the vertices only. Eq. (6.5) shows that, due to an apparently fortuitous cancellation, this large contribution is almost completely subtracted to leave a small fractional correction ∼ 0.69ε.
In the squeezed limit we find which is O(ε 2 ), as predicted by the Maldacena condition (3.18) and the property = 0 [33].
Comparison with previous results.-Khoury & Piazza estimated the bispectrum in a powerlaw model satisfying = 0 without invoking an expansion in slow-variation parameters [134]. 21 They quoted their results in terms of a quantity f X which replaces λ, For the DBI model, f X = 1 − c 2 s . They assumed constant f X , making their result valid to all orders in ε but only lowest-order in the time dependence of f X . Working in the equilateral limit for arbitrary constant f X we find Adopting the evaluation point k = 3k, this precisely reproduces (8.4) of Khoury & Piazza [134] when expanded to order ε. Although (6.8) does not strictly apply to DBI, where c s is changing, it can be checked that effects due to the time dependence of f X do not appear at next-order. Indeed, Eq. (6.8) yields (6.5) when f X = 1 − γ −2 .
Generalized DBI inflation.-The foregoing analysis was restricted to the asymptotic powerlaw regime, but this is not required. Using an arbitrary potential V (φ) in (6.3) one can construct a generic quasi-de Sitter background. Many of their properties, including the attractor behaviour, were studied by Franche et al. [141]. However, in the absence of a controlled calculation of next-order terms it has not been possible to estimate corrections from the shape of V (φ) or f (φ). Analogous effects have been computed for Galileon inflation [34], but our computation enables us to determine them in the DBI scenario for the first time.
For γ 1 the noncanonical structure suppresses background dependence on details of the potential. But small fluctuations around the background cannot be shielded from these details, which induce three-body interactions whether or not they are relevant for supporting the quasi-de Sitter epoch. These interactions generate relatively unsuppressed contributions to the three-point function.
We adapt the notation of Franche et al., who defined quantities measuring the shape of the potential V and the warp factor f , 22 The same branch of f 1/2 should be chosen in computing f 3/2 and sgn(φf 1/2 ). Note that these shape parameters do not coincide with the global slow-variation parameters ε and η. 21 See also Baumann et al. [140]. 22 The factor sgn(φf 1/2 ) was not used by Franche et al. [141], but is necessary here because the relativistic background solution requires f X = 1 + O(γ −2 ) up to corrections suppressed by O(ε) which are higher-order than those we retain. Depending on the direction of motion, this yieldsφ = ±f −1/2 + O(γ −2 ).

Franche et al. argued that ∆
1 was required to obtain attractor solutions, which we will assume to be satisfied in what follows. In addition, we work in the equilateral limit and take γ 1, which is the regime of principal interest for observably large f NL . We find The leading term of (6.5) is unchanged, but the subleading terms are dominated by shapedependent corrections, .
These subleading terms are more important than those of (6.5), which began at relative order γ −2 and are therefore strongly suppressed for γ 1. Moreover, inflation can occur even for relatively large values of ε v and η v -roughly, whenever ε v /γ < 1-so these corrections need not be extremely small. For large |f NL |, we estimate the relative correction to be For negative equilateral f NL , current constraints approximately require |f NL | 1/2 12 [6]. Therefore, these corrections can be rather important unless the potential is tuned to be flat, although some cancellation occurs because ε v and η v enter with opposite signs.
To obtain an estimate, suppose that ∆ is negligible. Taking the extreme 95%-confidence value f NL = −151 [6] and ε v ≈ |η v | ∼ 1 to obtain an estimate for a "generic" potential, the correction is of order 14% if η v > 0 and 20% if η v < 0. To reduce these shifts one might be prepared to tolerate a small tuning, giving perhaps ε v ≈ |η v | ∼ 0.1 and suppressing the correction to the percent level. However, the corrections grow with decreasing |f NL |. Keeping the generic estimate ε v ≈ |η v | ∼ 1, and using |f NL | ≈ 50, for which γ ≈ 10 and the approximation γ 1 used to derive (6.12) is at the limit of its applicability, we find the corrections to be of order 24% for η v > 0 and 36% for η v < 0.
Although Franche et al. argued that ∆ must be small to obtain attractor behaviour, it need not be entirely negligible. In such cases it introduces a dependence on the shape of the warp factor in addition to the shape of the potential. This may be positive or negative. If the ∆ and ε terms add constructively, the next-order correction can become rather large.
Infrared model.-The DBI scenario can be realized in several ways. The original "ultraviolet" model is now disfavoured by microscopic considerations [63,65]. Chen introduced [123, 142,143] an alternative "infrared" implementation which evades these constraints and remains compatible with observation [66,[144][145][146]. In this model the warp factor f (φ) is λ/φ 4 , where λ is a dimensionless parameter. The potential is V (φ) = V 0 − β 2 H 2 φ 2 /2, in which the mass is expressed as a fraction β 1/2 of the Hubble scale. The constant term V 0 is taken to dominate, making ε negligible. However the remaining slow-variation parameters need not be small. It is convenient to express our results in terms of the number, N e , of e-folds to the end of inflation. Background quantities evaluated at this time are denoted by a subscript 'e'. Computing (6.10) with these choices of V and f we find η e ≈ 3/N e and s e ≈ 1/N e . The infrared model is an example where ∆ e is not negligible, being also of order 1/N e . Specializing (6.12) to this case, we find

k-inflation
The k-inflation model of Armendáriz-Picón et al. [39] also admits power-law solutions. The action is where γ is a constant, no longer related to the speed of sound by the formula c s = γ −1 which applied for DBI. Unlike the DBI Lagrangian, Eq. (6.14) is unlikely to be radiatively stable and its microscopic motivation is uncertain. Nevertheless, nongaussian properties of the inflationary fluctuations in this model were studied by Chen et al. [43]. There is a solution with X = 2 − γ 4 − 3γ , (6.15) making ε = 3γ/2 and c s constant. Therefore this model has s = 0 but ε = 0, and is not scale-invariant even at lowest-order. Inflation occurs if 0 < γ < 2/3. The lowest-order contribution to f NL is of order 1/γ, making the next-order term of order unity. A next-nextorder calculation would be required to accurately estimate the term of order γ.
In the equilateral limit, Chen et al. quoted the lowest-order result f NL ≈ −170/81γ. Proceeding as in §1, one can estimate the fractional theoretical uncertainty in this prediction to be ∼ 9γ, or roughly ±20. This is comparable to the Planck error bar, and is likely to exceed the error bar achieved by a subsequent CMB satellite. Still working in the equilateral limit, we find As for DBI inflation, a fortuitous cancellation brings the fractional correction down from our estimate ∼ 9γ to ∼ 0.5γ. It was not necessary to choose a reference scale k in order to evaluate (6.16). In DBI inflation, to the accuracy of our calculation, the power spectrum is scale invariant but f NL is not. For the power-law solution of k-inflation, with the same proviso, it transpires that f NL is scale invariant even though the the power spectrum is not.

Because
= 0 in this model, the analysis of Khoury & Piazza does not apply. In a recent preprint, Noller & Magueijo gave a generalization which was intended to be valid for small and constant but otherwise arbitrary ε and s [147]. Their analysis also assumes constant f X . We set the reference scale to be k = 3k and work in the equilateral limit for arbitrary constant f X . One finds This occurs partially because they approximate the propagator (2.20) in the superhorizon limit |kc s τ | 1 where details of the interference between growing and decaying modes around the time of horizon exit are lost. For example, their approximation discards the Ei-contributions of (A.2) although these are O( ) and as large as other contributions which are retained. But were these terms kept, the superhorizon limit |kc s τ | 1 could not be used to estimate them. Infrared safety of the J i integrals in (3.15) is spoiled if truncated at any finite order, causing divergences in the squeezed limit ϑ i → 0 and a spurious contribution to the bispectrum with a local shape. As we explain in Appendices A and B, it appears that-as a point of principle-if = 0 corrections are kept then the shape of the bispectrum can be accurately determined only if the full time-dependence of each wavefunction around the time of horizon exit is retained.

Conclusions
In the near future, we can expect key cosmological observables to be determined to high precision. For example, the Planck satellite may determine the scalar spectral index n s to an accuracy of roughly one part in 10 3 [148]. In a formerly data-starved science, such precision is startling. But it cannot be exploited effectively unless our theoretical predictions keep pace.
Almost twenty years ago, Stewart & Lyth developed analytic formulae for the twopoint function accurate to next-order in the small-parameter ε = −Ḣ/H 2 . Subsequent observational developments have restricted attention to a region of parameter space where ε 1 is a good approximation, making the lowest-order prediction for the power spectrum an accurate match for experiment. The same need not be true for three-and higher n-point correlations, where the imminent arrival of data is expected to improve the observational situation. The results of §6 show that Planck 's observational precision in the equilateral mode may be comparable to next-order corrections. For a future CMB satellite it is conceivable that the data will be more precise than a lowest-order estimate. In this paper we have reported a next-order calculation of the bispectrum in a fairly general class of single-field inflationary models: those which can be described by a Lagrangian of the form P (X, φ), where X = g ab ∇ a φ∇ b φ is twice the field's kinetic energy. Our calculation can be translated to f NL , and in many models it provides a much more precise estimate than the lowest-order result.
Our major results can be categorized into three groups. The first group comprises tests of the accuracy, or improvements in the precision, of existing lowest-order calculations. The second generate new shapes for the bispectrum at next-order. The third involves technical refinements in calculating scalar n-point functions for n 3.

Accuracy and precision
Except in special cases where exact results are possible-such as the result of Khoury & Piazza for constant f X discussed in §6-predictions for observable quantities come with a "theory error" encapsulating uncertainty due to small contributions which were not calculated. For inflationary observables the typical scale of the theory error is set by the accuracy of the slow-variation approximation, where the dimensionless quantities ε = −Ḣ/H 2 , η =ε/Hε, s =ċ s /Hc s (and others) are taken to be small.
Power-law DBI inflation and k-inflation.-In §1 and § §6.1.1-6.1.2 we estimated the precision which could be ascribed to the lowest-order formula for f NL in the absence of a complete calculation of next-order effects. To do so one may use any convenient-but hopefully representative-subset of next-order terms, estimating the remainder to be of comparable magnitude but uncertain sign. Using the next-order contributions from the vertex coefficients g i , which can be obtained without detailed calculation, we estimated the fractional uncertainty to be of order 14ε for DBI and 9γ for k-inflation. The next-order terms neglected in this estimate come from corrections to the propagator, and from the time-dependence of each vertex. The prospect of such large uncertainties implies one has no option but to carry out the full computation of all next-order terms.
In the power-law DBI and k-inflation scenarios, we find that the terms omitted from these estimates generate large cancellations, in each case reducing the next-order contribution by roughly 95%. Similar large cancellations were observed by Gong & Stewart in their calculation of next-next-order corrections to the power spectrum [57]. After the fact, it seems reasonable to infer that the contributions from g i systematically overpredict the nextorder terms. But this could not have been deduced without a calculation of all next-order effects. Therefore, for power-law DBI and k-inflation models we conclude that the lowestorder calculation is surprisingly accurate.
For DBI inflation the status of the power-law solution is unclear, being under pressure from both observational and theoretical considerations. More interest is attached to the generalized case to be discussed presently. For k-inflation, taking present-day constraints on the spectral index into account, the next-order correction is of order 1%. Estimating the contribution of next-next-order terms using all available contributions from our calculation gives a fractional uncertainty-measured with respect to the lowest-order term-of order 22γ 2 . If similar large cancellations occur with terms not included in this estimate, the nextorder result could be rather more precise than this would suggest. 23 Without knowledge of such cancellations, however, we conclude that the uncertainty in the prediction for f NL has diminished to ∼ 250γ% 10% of the uncertainty before calculating next-order terms. In the case of power-law DBI inflation, the same method yields an estimate of next-next-order corrections at ∼ 40ε 2 , also measured from the dominant lowest-order term in the limit γ 1. This reduces the uncertainty to ∼ 300ε% 15% of its prior value.
Generalized DBI inflation.-The situation is different for a generalized DBI model with arbitrary potential V (φ) and warp factor f (φ). The largest next-order corrections measure a qualitatively new effect, not included in the power-law solution, arising from the shape of V and f . The fractional shift was quoted in (6.12) and can be large, because the DBI action supports inflation on relatively steep potentials. Indeed, if one were to tune the potential to be flat in the sense ε v ∼ |η v | 10 −2 then much of the motivation for a higher-derivative model would have been lost. Even for rather large values of |f NL | the correction can be several tens of percent for an "untuned" potential with ε v ∼ |η v | ∼ 1. For slightly smaller |f NL | the correction is increasingly significant, perhaps growing to ∼ 35%. The formulas quoted in §6.1.1 assume γ 1 and would require modification for very small f NL where O(γ −1 ) corrections need not be negligible. If desired, these can be obtained from our full formulae tabulated in §3. 2.
In a concrete model-the infrared DBI scenario proposed by Chen [123,142,143]-we find the correction to be ∼ 12% for parameter values currently favoured by observation, which translates to reasonably large shifts in f NL . For the maximum-likelihood mass suggested by the analysis of Bean et al., we find next-order corrections increase the magnitude of f NL by a shift |∆f NL | ≈ 19. This is a little smaller than the error bar which Planck is expected to achieve, but nevertheless of comparable magnitude. We conclude that a nextorder calculation will be adequate for Planck, but if the model is not subsequently ruled out a next-next-order calculation may be desirable for a CMBPol -or CoRE -type satellite.
For models producing small f NL , such as the "powerlike" Lagrangian discussed in §VI.B of Franche et al. [141], we find similar conclusions. However, in such models f NL is unlikely to be observable and the subleading corrections are of less interest.

New bispectrum shapes
Because our final bispectra capture the shape dependence in the squeezed limit, we are able to determine the relationship between the lowest-order and next-order shapes, discussed in §3.5.
Working in a P (X, φ) model, the enhanced part of the lowest-order bispectrum is wellknown to correlate with the equilateral template. Only two shapes are available, plotted in Table 9, and the bispectrum is a linear combination of these. The next-order bispectrum is a linear combination of seven different shapes, although these cannot be varied independently: strong correlations among their coefficients are imposed by the P (X, φ) Lagrangian. Many of these shapes also correlate with the equilateral mode, but two of them are different: in the language of §3.5, these are S ε and S s . For typical values of α (which is the c −2 s -enhanced corrections provided by Gong & Stewart [57]. part of λ/Σ) one can obtain a linear combination of the next-order shapes orthogonal to both lowest-order shapes. This is the orthogonal shape O appearing in Table 11. This shape represents a distinctive prediction of the next-order theory, which cannot be reproduced at lowest-order. It is very similar to a new shape constructed by Creminelli et al. [37], working in an entirely different model-a Galilean-shift invariant Lagrangian with at least two derivatives applied to each field. We conclude that even a clear detection of a bispectrum with this shape will not be sufficient, on its own, to distinguish between models of P (X, φ)-and Galileon-type. In practice, a network of interlocking predictions is likely to be required. For example, this shape cannot be dominant in a P (X, φ) model. It must be accompanied by a mixture of S 1 -and S 2 -shapes of larger amplitude.
The orthogonal shape is not guaranteed to be present at an observable level in every model, even in models where the next-order corrections are detectable. This may happen because α takes a value for which no bispectrum orthogonal to S 1 and S 2 can be generated perturbatively, or because the slow-variation parameters accidentally conspire to suppress its amplitude. Nevertheless, it will be present in many models. If detected, it would play a role similar to the well-known consistency relation between r and n t . Searching for the O-shape in real data is likely to require a dedicated template, and until this is constructed it is not possible to estimate the signal-to-noise and therefore the amplitude detectable by Planck or a subsequent experiment. It would be interesting to determine the precise amplitude to be expected in motivated models, but we have not attempted such an estimate in this paper.

Technical results
Our calculation includes a number of more technical results.
Treatment of boundary terms.-In §3.1 we gave a systematic treatment of boundary terms in the third-order action. Although these terms were properly accounted for in previous results [26,42,43], these calculations used a field redefinition which was not guaranteed to remove all terms in the boundary action.
Pure shape logarithms.-The subleading correction to the propagator contains an exponential integral contribution [of the form Ei(x)] whose time dependence cannot be described by elementary functions [34,43,57]. This term must be handled carefully to avoid unphysical infrared divergences in the squeezed limit where one momentum goes to zero. In Appendix B we describe how this contribution yields the J i functions given in Eq. (3.15). These are obtained using a resummation and analytic continuation technique introduced in Ref. [34]. The possibility of spurious divergences in the squeezed limit shows that, as a matter of principle, one should be cautious when determining the shape of the bispectrum generated by an approximation to the elementary wavefunctions.
In the present case, the J i contain 'pure' shape logarithms which are important in describing factorization of the three-point function in the limit k i → 0. Obtaining the quantitatively correct momentum behaviour requires all details of the interference in time between growing and decaying modes near horizon exit. The possibility of such interference effects, absent in classical mechanics, is a typical feature of quantum mechanical processes.
This interference correctly resolves the unwanted divergences in the opposing infrared limit k j → 0 with j = i. We discuss these issues in more detail in Appendices A and B.
Comparison with known results.-We have verified Maldacena's consistency condition (3.18) to next-order in canonical models by explicit calculation of the full bispectrum. This agrees with a recent calculation by Renaux-Petel [114]. In the case of power-law inflation with = 0 and constant ε and s, we reproduce a known result due to Khoury & Piazza [134]. A subset of our corrections were computed by Chen et al. [43]. Up to first-order terms in powers of slow-roll quantities-where our results can be compared-we find exact agreement.

Note added
Immediately prior to completion of this paper, a preprint by Arroja & Tanaka appeared [149] which appears to present arguments regarding the role of boundary terms which are equivalent to those of §3.1.2.  Evaluated at the conventional reference scale k = k t Table 5. Equilateral limit of f NL at lowest-order and next-order. The numerical constant ω satisfies ω = 1 a Evaluated at the reference scale k = 2k UV , where k UV is the common hard momentum Table 6. Squeezed limit of f NL at lowest-order and next-order. The numerical constant ω satisfies ω = 1 2 ln 3 2 = coth −1 5.  [108]. e The local template, and the operators ζζ 2 and ζ(∂ζ) 2 , are strongly peaked in the "squeezed" limit where one momentum becomes much softer than the other two. For these shapes the inner product which defines the cosine is divergent, and must be regulated. The resulting cosines are almost entirely regulator-dependent. See the discussion in §3. 5.
The values we quote are meaningful only for our choice of regulator. For the values quoted above we have used δ min = k/k t = 10 −3 , where δ min was defined in the main text.  [118]. The plotted quantity is x 2 y 2 B(x, y, 1), where x = k 1 /k 3 , y = k 2 /k 3 and B is the bispectrum, and normalized to unity at the equilateral point [119]. The plotted quantity is k 2 1 k 2 2 k 2 3 B(k 1 , k 2 , k 3 ) as a function of the α and β parameters defined in (3.23a)-(3.23c). Table 9. Lowest-order bispectrum shapes enhanced by c −2 s in P (X, φ) models. Table 10. Bispectrum shapes enhanced by c −2 s at next-order in a P (X, φ) model. The Babich et al. and Fergusson-Shellard plots are defined in Table 9.   Table 13. Overlap cosines between the orthogonal shape O and common templates, defined in Table 7. Sign information has been discarded.

A Propagator corrections
At lowest-order, the wavefunctions can be obtained from Eq. (2.20) by setting all slowvariation parameters to zero. Choosing a reference scale k , this yields the standard result (A.1) Next-order corrections to the propagator were discussed by Stewart & Lyth [55], who quoted their result in terms of special functions and expanded uniformly to next-order after taking a late-time limit. The uniform next-order expansion at a generic time was given by Gong & Stewart [57] for canonical models and by Chen et al. [43] in the noncanonical case. Their result was cast in a more convenient form in Ref. [34] whose argument we briefly review.
The next-order correction is obtained after systematic expansion of each quantity in (2.20) to linear order in the slow-variation parameters. Contributions arise from each timedependent factor and from the order of the Hankel function. Collecting the formulae quoted in Ref. [43], one finds where Ei(x) is the exponential integral, This is well-defined for x < 0. For x 0 it should be understood as a Cauchy principal value. For complex argument-required by (A.2)-one promotes Eq. (A.3) to a contour integral, taken on a path running between t = x ∈ C and |t| → ∞ in the half-plane Re(t) < 0. Using Cauchy's theorem to rotate this contour onto the negative real t-axis, one finds The next-order correction to (A.1), expanded uniformly to O(ε) in slow-variation parameters but including the exact time-dependence, is and quantities labelled ' ' are evaluated at the horizon-crossing time for the reference scale k . We have used ∆N to denote the number of e-folds which have elapsed since this time, so ∆N = ln |k c s τ |. The limit → 0 should be understood after the integration has been performed and merely guarantees convergence as ξ → −∞.
Setting the integral aside, Eq. (A.5) has the appearance of an expansion in powers of kτ . In fact, the time dependence is exact. Although terms of high-order in kτ become increasingly irrelevant as |kτ | → ∞, all such terms make comparable contributions to a generic n-point function at the time of horizon crossing, where |kτ | ∼ 1, and should not be discarded. This reflects interference effects between the growing and decaying modes at horizon exit. Indeed, a high-order term such as (kτ ) n typically generates a contribution ∼ n!(k/k t ) n . Such terms are suppressed by the factor (k/k t ) n which typically varies between 0 and 1/2, but are enhanced by the rapidly growing factorial. Therefore their contribution to the shape dependence must usually be retained. An infinite series of such terms may converge only for certain ratios k/k t , requiring a continuation technique to obtain the momentum dependence-and hence the shape of the bispectrum-for arbitrary k i . We will see an explicit example in Appendix B below.
(A. 7) In order that the power spectrum remains conserved on superhorizon scales, one should verify that δu approaches a constant and δu → 0 as |kτ | → 0. The explicitly time-dependent term ∆N apparently spoils the required behaviour, but is compensated by a logarithmic divergence from the integral. Indeed, for τ → 0 one finds τ −∞(1+i ) dξ ξ e 2ikcs ξ = ln |2ikc s τ | + O(kc s τ ), (A. 8) which precisely cancels the time dependence arising from ∆N . Note the incomplete cancellation of the logarithm, which leaves a residual of the form ln 2 and is the origin of the ln 2 term in the Stewart-Lyth constant C = −2 + ln 2 + γ E [55,60]. This is a primitive form of the incomplete cancellation which leads to residual ln k i /k t and ln k t /k terms after cancellation of ln τ logarithms in the three-point function.

B Integrals involving the exponential integral Ei(z)
The principal obstruction to evaluation of the next-order corrections using standard methods is the Ei-term in (A.5). This can not be expressed directly in terms of elementary functions whose integrals can be computed in closed form. It was explained in Appendix A that although Ei(2ikc s τ ) contains terms of high orders in kτ which become increasingly irrelevant as |kτ | → 0, these cannot usually be neglected when computing n-point functions. A term of order (kτ ) n generates a contribution ∼ n!(k/k t ) n and as part of an infinite series will converge at best slowly for certain values of k/k t -typically, when the scale k itself is squeezed to zero. The ratio k/k t is maximized when k is left fixed but one of the other momenta is squeezed. To obtain a quantitative description of this limit may require continuation. Hence, in order to understand the shape generated by arbitrary k i , it follows that the entire kτ -dependence must be retained. It was for this reason that we were careful to expand uniformly only in slow-roll parameters in Appendix A, while keeping all time-dependent terms.
The integrals we require are of the form [34] I m (k 3 ) = We are again using the convention that, although I m depends on all three k i , only the asymmetric momentum is written as an explicit argument. In the following discussion we specialize to I 0 . Comparable results for arbitrary m can be found by straightforward modifications and are given in Ref. [34]. We introduce the dimensionless combinations After contour rotation and some algebraic simplification, we can express I 0 in terms of ϑ 3 and θ 3 [34] 3) The v-integral has a Puiseaux series representation where the sum converges uniformly for all complex θ 3 u. The right-hand side has a singularity at θ 3 u = 0, and the logarithm generates a branch cut along | arg(θ 3 u)| = π. The u-independent term γ E is the Euler-Mascheroni constant γ E ≡ 0 ∞ e −x ln x dx and is obtained by expanding e −2v in series, integrating term-by-term, and matching the undetermined constant of integration with the left-hand side of (B.4) in the limit u → 0.
Because (B.4) is uniformly convergent we may exchange integration and summation, evaluating I 0 using term-by-term integration. We find [34] This is singular when ϑ 3 → 0, which corresponds to the squeezed limit where either k 1 or k 2 → 0. There is nothing unphysical about this arrangement of momenta, for which θ 3 → ∞, but we will encounter a power-law divergence unless the bracket [· · · ] vanishes sufficiently rapidly in the same limit. The sum converges absolutely if −1/2 < Re(θ 3 ) < 1/2, corresponding to the narrow physical region 0 < k 3 < k t /4 where k 3 is being squeezed. For θ 3 satisfying this condition the summation can be performed explicitly, yielding ln(1 + 2θ 3 ) −1 = ln ϑ 3 . It can be verified that I 0 (k 3 ) is analytic in the half-plane Re(θ 3 ) > 0. Therefore it is possible to analytically continue to θ 3 1/2, allowing the behaviour as θ 3 → ∞ and ϑ 3 → 0 to be studied. We find where we have introduced a function J 0 , which is related to I 0 by a simple change of normalization. It is more convenient to express our final bispectra in terms of J 0 rather than I 0 . In the limit ϑ 3 → 0, Eq. (B.6) is finite with J 0 → −1. The θ 3 → ∞ behaviour of the sum has subtracted both the power-law divergence of the prefactor and the logarithmic divergence of ln(2θ 3 ). To obtain this conclusion required information about the entire momentum dependence of (B.1) in order to perform a meaningful analytic continuation. 24 Had we truncated the k 3 ζ dependence of (B.1), or equivalently the θ 3 u dependence of (B.4), we would have encountered a spurious divergence in the squeezed limit and a misleading prediction of a strong signal in the local mode. Several of our predictions, including recovery of the Maldacena limit described by (5.3), depend on the precise numerical value of J 0 as θ 3 → ∞ and therefore constitute tests of this procedure.
We can now evaluate I m (k 3 ) for arbitrary m, although only the cases m 2 are required for the calculation presented in the main text. The necessary expressions are (B.7b) In the limit ϑ 3 → 0 one can verify that J 1 → −1/2 and J 2 → −2/3.

C Useful integrals
To simplify evaluation of the various integrals which arise in computing the bispectrum, it is helpful to have available a formula for a master integral, J , for which the integrals of interest are special cases [34]. We define dξ e iktcs ξ γ 0 + iγ 1 c s ξ + γ 2 c 2 s ξ 2 + iγ 3 c 3 s ξ 3 + γ 4 c 4 s ξ 4 + N δ 0 + iδ 1 c s ξ + δ 2 c 2 s ξ 2 + iδ 3 c 3 s ξ 3 + δ 4 c 4 s ξ 4 (C.1) 24 A precisely analogous situation occurs at finite temperature, where one first computes the imaginary time Green's function G(iωn) on a discrete spectrum of Matsubara frequencies {ωn}. To obtain the realtime Green's function one must extend G into the complex plane (in the sense that one finds a function G(z), where z ∈ C, which coincides with the imaginary time Green's function when z = iωn), and analytically continue to the real axis. Unless one has knowledge of G on the complete set of Matsubara frequencies-which may be impossible in an effective field theory-it is difficult to extract meaningful results for the real-space correlation function. In the present case, the requirement to know the series expansion in (B.4) or (B.5) for all n is analogous to the requirement to know G(iωn) for all Matsubara frequencies.
where N represents the number of e-folds elapsed since horizon exit of the reference mode k , giving N = ln |k c s ξ|. Eq. (C.1) is an oscillatory integral, for which asymptotic techniques are well-developed, and which can be evaluated using repeated integration by parts. To evaluate the N -dependent terms one requires the standard integral which can be obtained by contour rotation and use of the Euler-Mascheroni constant defined below Eq. (B.4). We conclude We use this result repeatedly in §3.3 to evaluate integrals in closed form.