Non-Gaussianity in inflationary scenarios for primordial black holes

Working in an idealised framework in which a series of phases of evolution defined by the second slow-roll parameter $\eta$ are matched together, we calculate the reduced bispectrum, $f_{\rm NL}$, for models of inflation with a large peak in their primordial power spectra. We find $f_{\rm NL}$ is typically approximately constant over scales at which the peak is located, and provide an analytic approximation for this value. This allows us to identify the conditions under which $f_{\rm NL}$ is large enough to have a significant impact on the resulting production of primordial black holes (PBHs) and scalar induced gravitational waves (SIGWs). Together with analytic formulae for the gradient of the rise and fall in the power spectrum, this provides a toolkit for designing or quickly analysing inflationary models that produce PBHs and SIGWs.


Introduction
Models of inflation capable of producing a large peak in the primordial power spectrum at small scales are the subject of much active research. At large scales, the power spectrum is tightly constrained by observations of the Cosmic Microwave Background (CMB) to be approximately scale-invariant, with a slight red tilt, and with an amplitude A R ∼ O(10 −9 ) [1,2]. These constraints do not apply, however, to the power at much smaller scales.
One potential consequence of an enhanced power spectrum is the formation of primordial black holes (PBHs). These form from the gravitational collapse of perturbations at enhanced scales upon horizon re-entry [3]. Recently, there has been a renewed interest in PBHs, especially owing to their potential as a dark matter candidate [4][5][6][7][8][9]. PBHs have also been suggested as a source for gravitational wave (GW) events witnessed by LIGO and Virgo collaborations [10][11][12][13][14] and as contributors to the stochastic GW background which may have been detected by NANOGrav [15,16].
If PBHs are to be produced in the numbers necessary to be responsible for these effects, then the power spectrum at small scales likely needs to be enhanced by ∼ O(10 7 ) relative to CMB measurements of A R [17][18][19][20][21], the precise enhancement required depending on the finer details of both the collapse and the shape of the power spectrum. For power spectra with a large peak at short scales, scalar induced gravitational waves (SIGWs) will also be produced at second order in perturbation theory [22][23][24][25][26][27][28][29][30][31][32].
As well as the necessary growth in the power spectrum, it is well understood that non-Gaussianity of the scalar perturbations can have important effects on both PBH production [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48] and the spectrum of SIGWs [32,[49][50][51][52][53][54]. Although in reality the full probability distribution of perturbations is needed to accurately determine the production of PBHs, the reduced bipsectrum, f NL , can be used as a guide for when non-Gaussianity becomes important to the calculation. An f NL of O(1) will significantly affect both the abundance of PBHs, and their mass distribution [37,38,41,42,46]. In the case of SIGWs, if non-Gaussianity is local, then f NL together with the power spectrum is sufficient to fully determine the spectrum of GWs. In this context, for example, an f NL ∼ O (1) to O(10) will have a significant effect when A R ∼ 10 −3 [49,53]. Since PBHs and SIGWs form at enhanced scales, the most relevant value of f NL is its value around the peak of the power spectrum.
There are many models for generating enhancements of the primordial power spectrum that have already been explored in the literature (see Refs. [46,48,[55][56][57] for recent examples). In general, these are based on specific single-field inflationary potentials. A framework more amenable to analytic studies, however, was introduced by Byrnes et al. in Ref. [58]. Instead of working from a pre-determined potential, inflation is modelled as a series of distinct "phases" each characterised by a constant value of the second slow-roll parameter η; transitions between phases occur instantly [55,59]. This approach can be used to approximate a wide variety of models, and we refer to models in this framework as "idealised".
At present, works studying idealised models have focused on the power spectrum [56-58, 60, 61]. Our primary aim in this paper, therefore, is to extend this approach to compute the bispectrum as well, focusing in particular on the value of f NL at scales around the peak of the power spectrum. We find that f NL takes a similar form for all idealised models with a large peak in their power spectrum, and find a simple expression for the value of f NL at scales around this peak. We also generalise the idealised formalism further to allow η to have a linear dependence on conformal time, using this to study the effect of relaxing the assumption of instantaneous transitions. This paper is structured as follows: in Section 2 we present the general idealised framework for generating models of inflation with peaked primordial power spectra. In Section 3 we derive the form of the bispectrum using standard cosmological perturbation theory and the in-in formalism, f NL is also defined. In Section 4 we present plots of f NL in a selection of idealised models and derive a simple expression for its value around peak scales. In Section 5 we consider models with more realistic transitions which occur linearly rather than instantly in conformal time and discuss the effect of this on f NL . Finally we summarise our results and discuss their implications for PBH formation and the spectrum of SIGWs.

Inflationary Action
In this paper we consider single-field, minimally-coupled models of inflation defined by the action where g is the determinant of the metric, R is the Ricci scalar, φ is the inflaton field and V is its self-interaction potential. We define the first and second slow-roll parameters as = −Ḣ H 2 and η =˙ H respectively. A dot indicates a derivative with respect to cosmic time and H is the Hubble parameter. It is assumed in all of our models that 1, but that η can be much larger (at least of O(1), or higher). Using cosmological perturbation theory the action (2.1) can be expanded in terms of the gauge-invariant comoving curvature perturbation R about a flat Friedmann-Robertson-Walker background. The power spectrum and bispectrum of R is what we wish to compute.

Mode Functions
After expanding the action to second order in R we can obtain an equation of motion for the Fourier modes R k . Defining z 2 ≡ 2a 2 , where a is the scale factor and v k ≡ zR k , the auxiliary field v k satisfies the Mukhanov-Sasaki equation [62,63] where a dash indicates a derivative with respect to conformal time τ . This equation can be rewritten in terms of the quantity where we have used the de-Sitter space approximation aH = −1/τ . For constant η, the general solution to this equation is written in terms of Hankel functions of the first and second kind as Note also that depends on conformal time as ∝ (−τ ) −η for constant η.

Idealised Approach
Following the approach of Byrnes et al. in Ref. [58], we can now use this solution to construct models of inflation with peaked power spectra in our idealised framework. We model inflation as a succession of phases defined by a constant η and assume that transitions between these phases occur instantaneously. This means that the solutions for the mode functions in each phase are given by Eq. (2.5), with constants C 1 and C 2 taking different values in each phase.
In the first phase the constants are fixed by assuming the Bunch-Davies initial conditions at some initial time, and in subsequent phases by requiring continuity of R and its first derivative at each transition time, according to the Israel junction conditions [64,65].

Linear Transitions
If we want to consider more realistic transitions that are not instantaneous, we must allow for η, and hence ν, to be time dependent. If η varies linearly in conformal time, an analytical solution to Eq. (2.4) exists, allowing us to continue with a fully analytic approach. The solution is given in terms of Whittaker functions as where the values of p, q, and r depend on the details of the transition. If we take a general linear transition between two values of η, η 1 and η 2 , then we have where the transition between η 1 and η 2 takes place between the conformal times τ 1 and τ 2 , τ 1 < τ 2 . In this case the values of p, q, r are: (2.10) C 1 and C 2 are once again fixed by continuity.
In this paper, in addition to instantaneous transitions, we consider models where η varies linearly in conformal time during a transition with non-zero duration. We refer to these as "linear transition" models.

Power Spectrum
The power spectrum P R (k) for a given inflationary model can be calculated from the mode functions. It is defined by its relation to the two-point correlator (2.11) The dimensionless power spectrum P R (k) is then defined in terms of the power spectrum as In this paper, "power spectrum" will always refer to the dimensionless power spectrum, unless otherwise specified. This can be expressed neatly in terms of the mode functions as (2.13)

Scenarios with a Large Peak
As mentioned, we are interested in models of inflation which produce a large peak in the power spectrum. Sharp peaks can be generated in a variety of ways, but in general they require the first slow-roll parameter to decay rapidly (see for example the discussion in Ref. [60]). The simplest case is therefore to begin in a phase of slow-roll, with negligible η, and transition rapidly to a phase with η 0. When η < −3, rather than being frozen, modes grow on superhorizon scales. This not only affects modes which exit the horizon during the the large negative η phases, but also modes which exited during some period before it began. The number of modes affected is determined by the length of the negative η phase. As described in Ref. [60] it is these modes that gain a large blue spectral index, and lead to a steep power spectrum. Modes that exit during slow-roll initially freeze, but can then grow again later during a large negative η phase and gain a P R (k) ∝ k 4 dependence [58]. This is almost the steepest sustained k-dependence the power-spectrum can achieve (so long as there is only one large negative η phase -see Appendix A for a full discussion); a marginally steeper growth of P R (k) ∝ k 5 (log k) 2 is possible if an η = −1 phase precedes the large negative η phase [60]. In general, modes that exit the horizon before a large η phase which they are affected by gain a k-dependence of the form [60] P R (k) ∝ k 5−|η+1| , (2.14) where η in this formula is the η-value of the phase in which modes exit. The scale dependence is as described by Eq. (2.14) so long as the growing phase lasts a sufficiently long time for the initially sub-leading term of the asymptotic expansion of the mode functions to become dominant [60]. This compares to the standard expression for modes whose spectral index is not affected by future phases of evolution, which follows from the leading term in an asymptotic expansion of the mode functions. We review and refine the discussion of Ref. [60] in Appendix A. We now highlight a selection of scenarios in which a large power spectrum peak can be generated.
The standard scenario The simplest way to get a large peak in the power spectrum follows from the discussion above. We begin in a phase of slow-roll η ∼ 0 inflation in order to be consistent with CMB observations of approximate scale-invariance at large scales. During this period of slow-roll inflation the inflaton rolls down its self-interaction potential. This is followed by a period with a substantial decrease in the kinetic energy of the inflaton. In terms of η values this corresponds to a phase with a large, negative value of η. The simplest case of η = −6 is referred to as "ultra-slow-roll" (USR) and is one of the examples moststudied in the literature. This phase is achieved by introducing an extremely flat region in the inflaton's potential. Modes which exit during USR have a scale invariant spectrum, and so, in order bring the power spectrum down again after the peak, a final phase with a large positive η-value is required. This corresponds to a steep fall off in the inflationary potential and to the inflaton gaining kinetic energy again. Modes which exit during this phase freeze after horizon crossing in the usual way. It isn't necessary for the power spectrum's amplitude to be reduced back to CMB levels, but a model must end in a phase of η > −3 otherwise super-horizon modes will continue to evolve. This phase could be followed by other phases, such as a further slow-roll phase, but as long as there are no more phases with η < −3 the modes which exited in previous phases are unaffected. We can also change the value of η in the large negative η phase so that we get something other than USR. Values with η < −6 correspond to the inflaton slowing down more quickly than it would on a flat potential, and require the potential to start to rise. This option is appealing because modes which exit during this phase have a red spectral index, and the spectrum begins to decay from its peak due to these modes alone (i.e. without necessarily needing modes to exit in a subsequent large positive η phase). A further phase of η > 0 is still required, however, so that the modes freeze.
In summary, in the "standard scenario" we have three phases with η ∼ 0, η ≤ −3, and η > 0. Modes which exit towards the end of the slow-roll phase ultimately gain a k 4 dependence, as predicted by Eq. (2.14); modes which exit during the η < −3 phase gain a scale dependence as predicted by Eq. (2.15); and modes which exit in the large positive η phase are red tilted also as predicted by Eq. (2.15). The overall appearance is a sharp rise followed by a sharp fall, although possibly with an extended peak in between depending on the value of the large negative η. In the simplest case of USR, there is an extended, flat peak.
Variants on the standard scenario It is, of course, possible to generalise the standard scenario. We can consider adding a further phase after the slow-roll phase but before the large negative η phase. Modes exiting during these phases will be affected by the large negative η phase and lead to a different growth than the k 4 behaviour. If, for example, an η = −1 phase is added, the slightly steeper growth of P R (k) ∝ k 5 (log k) 2 is found as mentioned above.
Inflaton Falls One particular variant of the standard scenario is noteworthy, and we refer to it as the "inflaton falls" model after the study of Inomata et al. in Ref. [57]. In this case an initial phase of slow-roll inflation is followed immediately by a phase with a large positive value of η, after which there is a large negative η stage, in that order, before the model returns back to slow-roll (or indeed any other η > −3 value). This corresponds to the inflaton "falling" off its slow-roll plateau and its kinetic energy initially increasing, before reaching a sufficiently flat (or growing) region of the potential in which the inflaton decelerates again. If the large negative η phase is large enough to affect not only the modes which exit immediately before it during the large positive η phase, but also some of the modes that exit at the end of the slow-roll phase, the result is that the modes which exit during slow-roll gain a k 4 dependence as usual, and the modes which exit during the large η phase gain a dependence according to Eq. (2.14). With a large positive η, this implies a strongly red, decaying spectrum (but with a different tilt than if these modes simply freeze and the tilt was determined by Eq. (2.15)).
In Ref. [57] the discussion implies that the peak in the power spectrum is produced by a rapid increase in the kinetic energy of the inflaton thanks to a sharp dip in its potential, but this hides the origin of the peak somewhat. Really it is modes that exit during slow-roll which ultimately account for the steep growth to a peak, and the modes which exit during the phase in which the kinetic energy increases which account for a steep decay. As in the standard scenario the large negative η phase is required to generate the steep rise (by changing the tilt of the slow-roll modes).
Repeating Model In Ref. [60] it is argued that only a marginally steeper sustained growth than k 4 could possibly be generated during canonical inflation. This statement remains true if there is only one large negative η phase, but in later work Tasinato [56] showed that a faster growth can be achieved by successive phases of SR and USR. We attempt to give a full explanation of how this works in Appendix A. This model begins in a phase of SR, before switching twice between phases of USR and SR. It can be generalised to have other values of η in these alternating phases, and will also require a large positive η-value added in one of the ways described in the simpler scenarios above in order for the spectrum to decay.
We will now compute the bispectra for each of the kinds of model just introduced. The remarkable thing we find is that, despite the numerous possibilities for generating a large peak in the power-spectrum and the differences in the precise form of the spectrum between them, in all of them the reduced bispectrum, f NL , looks remarkably similar and possesses the same structure. As discussed, of particular interest is the non-Gaussianity at scales around the peak of these power spectra. We will show that the level of non-Gaussianity around the peak, in all of these models, can be approximately calculated from the same extremely simple formula.

Third Order Action
To compute the bispectrum we need to expand the action (2.1) to third order in the curvature perturbation [66,67] where (3. 2) The last term in the cubic action is usually removed by making a field redefinition [66] of the form Once the final term has been removed, the remaining terms in Eq. (3.1) are at least of O( 2 ) or higher with the exception of fourth term, which is proportional to η. In our models η varies rapidly between phases of inflation (either instantaneously, or via rapid linear transitions as discussed in §2.4), meaning that contributions from this term dominate over all others. Calculation of the bispectrum can now proceed using the redefined quantity.
In the final answer, however, we must remember to convert back to the original physical curvature perturbation, and hence f (R n ) will also contribute to the bispectrum. In the models we consider in this work we will only need to retain the first term in f (R n ). This is because evaluation of the final bispectrum will always take place at late times, once the scales being considered are super-horizon, and during a phase in which R is conserved on super-horizon scales (a phase where η > −3). This means the third term can be neglected due to spatial derivatives, while the second and fourth terms can be neglected as both are proportional tȯ R n . After these considerations, we end up with just two contributions to the bispectrum: one coming from the only surviving term in the cubic order action and the other coming from the first term in the field redefinition (3.5)

Form of the Bispectrum
The bispectrum B(k 1 , k 2 , k 3 ) can be expressed in terms of the three-point correlator of the mode functions This can be computed by using the in-in formalism to calculate the three point function of the redefined field and then accounting for the field redefinition using Wick's theorem [66]. The contribution from Eq. (3.4) (dropping the subscript n for simplicity) is given by the integral [68] (see also [69]) 7) where R k i (τ ) are the mode functions to be calculated in the idealised approach of §2.3.
Since we are interested in the bispectrum evaluated at late times, τ 0 will be some conformal time long after the transitions have taken place and all modes of interest have left the horizon. In the case of instantaneous transitions the integral Eq. (3.7) becomes a sum [68] 8) where τ i labels the time of the ith transition, and ∆η i is the change in η at this transition.
Next we account for the field redefinition. For a field redefinition of the form: R → R n + qR 2 n , Wick's Theorem tells us that In our case the first term here is the contribution from the cubic order action of our redefined field given above, and q = η/4. Returning to Fourier space, the contribution of the field redefinition to the bispectrum can be written in terms of mode functions R k as The sum of Eq. (3.7) and Eq. (3.10) gives the full bispectrum for any model of the form we're considering. Finally, we define the reduced bispectrum f NL (k 1 , k 2 , k 3 ), which is useful a measure of the relative size of non-Gaussianity . (3.11)

Non-Gaussianity in Idealised Models of Inflation
Using the expressions given above we now calculate the value of f NL for a selection of idealised inflationary models, as described in §2.6. We consider first the case of instantaneous transitions between phases, and focus on the equilateral configuration 1 (k 1 = k 2 = k 3 = k).  Remarkably, we find that in all models investigated that give a very large peak, the non-Gaussianity consistently exhibits the same structure. Moreover, at scales around the peak of the power spectrum, f NL is quasi-local (we observe this by investigating f NL away from the equilateral limit, finding only weak dependence on configurations for these scales).
To present the results 2 , in Fig. 1 we plot both the power spectrum and the value of f NL for a selection of models which produce power spectrum peaks in the standard scenario and variations. In Fig. 2 we do the same for an inflaton falls model, and a repeating model. In all of these cases the structure of f NL looks very similar; the features common to all these plots include: • A large peak for modes exiting during the initial phase -this corresponds to the dip in the power spectrum where the leading and sub-leading terms in the mode functions cancel each other out; • A long, approximately flat region that begins for modes exiting during the initial phase and ends for modes exiting during the final phase. We refer to this as the "plateau" region, and the value of f NL within it is what we call the "plateau value". It extends over the scales at which the peak in the power spectrum occurs; • A growing, highly oscillatory region at short scales.
Similar features, such as a scale-invariant f NL over peak scales, have also been noted in previous studies e.g. [45]. It is already known that the presence of sharp features in the inflaton's self-interaction potential will lead to oscillations in f NL at very small scales [71], in our models this is due to the instantaneous nature of the transitions between phases. These oscillations can clearly be seen at short scales in our plots in Fig. 1, for example. They begin to appear at scales k ∼ −1/τ n where τ n is the conformal time at which the final transition between phases occurs. This means that they appear for modes exiting the the horizon during the final phase of a given model, which in our case will always be a phase with η > −3 where the power spectrum is decaying. Consequently, these oscillations will always appear at scales smaller than those at which the peak of the power spectrum is located. These are scales smaller than those we are interested in.
One can also see in Fig. 1 that once the oscillations in f NL begin at short scales, the amplitude of these oscillations appears to grow without bound. This is a result of dealing with models in which the transitions between phases are instantaneous; such a growth is not a feature of realistic physical models where oscillations will be damped on short scales.
Of most interest is the value of f NL in the plateau region. Since the plateau extends over the scales at which the power spectrum reaches its peak, it is this value of f NL which is important for the formation of PBHs and SIGWs. It is natural, therefore, to ask what determines its value. We find that, as long as the large negative η phase lasts sufficiently long to form a peak of the size needed to produce PBHs, a plateau region in f NL always forms. Moreover, we see that the value of f NL in this region is sensitive only to the large negative η-value, and the value of η in the phase immediately following it. The plateau value is almost completely independent of the η-values of any previous phases, or their duration.

Plateau Value
An important feature of the plateau region of f NL is that it extends across modes which exited the horizon in different phases. For example, in standard scenarios it includes modes which exited during the large negative η phase, and the phase immediately prior.
To obtain an approximate expression for f NL in the plateau phase, therefore, we use this observation and study k-modes which exit the horizon during a phase of η 1 < −3, and are matched at the transition to a phase of η 2 > −3 evolution. Since R is ultimately conserved in this second phase, as long as this phases last sufficiently long, any further phases that also conserve R cannot affect the answer. We also assume that such modes are unaffected by transitions while deep inside the horizon and so are normalised to the Bunch-Davis vacuum at early times. Under these assumptions, the answer for the two phase model gives us an approximate value for the f NL in the plateau of a model where η 1 is the η-value in its most recent negative η phase, and η 2 the η-value of the phase immediately following it.
To find an analytic expression for f NL using the the two-phase model, the general mode functions (2.5) for both phases are expanded in the superhorizon limit (kτ → 0). Using these expansions of the mode functions the power spectrum and bispectrum are computed as described above. The bispectrum receives contributions from two terms: one from the single instantaneous transition in this model, and the other from the field redefinition term. It is clear that the redefinition term provides a scale-invariant contribution f NL, redef = 5 12 η 2 to the non-Gaussianity [45,46], while the contribution from the transition is harder to estimate 3 .
We proceed by considering the full expression obtained from both contributions using Mathematica. The expression is rather long and unwieldy, but it can be simplified drastically by taking suitable limits. First, we expand the reduced bispectrum in terms of τ 0 , the time at which it is evaluated. Since we are interested in the non-Gaussianity at times long after the transitions have taken place, only the leading order O(τ 0 0 ) term of this expansion is retained. A barrier to further simplification is the presence of gamma functions in the resulting expression. This can be overcome by initially restricting to integer values of η and dealing with odd and even values separately. Subsequently, the k → 0 limit is applied to reveal a very neat expression. Finally, we verify that the result holds for non-integer η. The expression is given by Recall that this estimate is for η 1 < −3 and η 2 > −3, and hence that it breaks down as η 1 → η 2 where the expression exhibits a singularity.
A few further comments are in order. As discussed above if additional phases are present after the η 2 phase, they do not affect the power spectrum or bispectrum as long as they conserve R. We note that in such cases the scale invariant contribution to f NL from the field redefinition term alters in value, but this change is precisely accounted for by contributions from the additional transition terms.
It only remains then to understand why this approximation works for modes which exit the horizon in the phase immediately before the large negative η phase, and why all other transitions before this phase do not contribute to the bispectrum. We believe the former is true because the only dependence on the previous phases comes from the previous transition term in Eq. (3.8), which is negligible when compared to the final transition term as well as to the field redefinition terms, for scales that are affected by the rapid growth. The reason for this can be seen from the time dependence of the transition term in Eq. (3.8). The three copies of R are evaluated long after the transition and, have thus grown far beyond the amplitudes of the other copies of R evaluated at the transition time. However, when computing f NL , we divide by the final power spectrum, which effectively includes four copies of R. This ratio is thus proportional to 1/R meaning it decays as R grows. The same does not happen for the final transition term because all copies of R have similar sizes, since the phase after the transition conserves R. It is also clear that earlier transitions provide even smaller contributions to the bispectrum, since the amplitude of R which enters Eq. (3.8) is so much smaller before modes enter the η < −3 growing phase.
Finally, we note that the virtue of the expression in Eq. (4.1) is that it can be used to rapidly estimate the value of f NL around the peak of a power spectrum in a wide range of Figure 3: Dependency of f NL on η-values on the relevant phases of inflation. In general, η 1 is the most recent large negative η-value of the model and η 2 the η-value of the phase immediately following. In all cases, the value of f NL is larger when the magnitude of any η-values are increased. When one η is kept fixed, as in a (top left) and b (top right), f NL tends to a certain value as the other η's magnitude is increased. In c (bottom left) one can see a section of the 3D parameter space. In d (bottom right) we see the dependence of f NL on η in a model of the form 0 → −η → η. When η is increased f NL grows with no bound on its value. inflationary models. Below, we discuss the implications of this expression for f NL for a few different classes of models present in the literature.

Standard Scenario
In Fig. 1, in addition to the full f NL , we plot the approximate plateau value obtained by substituting the values of η during the final two phases into Eq. (4.1). Note that we have also verified that this value does not depend on the initial phase being SR (any other phase, such as η = −1 will lead to the same plateau value of f NL ), and that extra η > −3 phases can be glued to the end of the model also without affecting the plateau value.
Next, in Fig. 3, we use Eq. (4.1) to see what effect varying the values of η in the final two phases will have on the plateau f NL . If we fix one of these two values, and vary the other, then the plateau f NL grows with the magnitude of the varying η. The greater the magnitude of η, the greater the plateau f NL , but only up to a limit. When one of the η values is kept constant there is a maximum f NL that can be achieved. If we fix η in the middle phase to be η m , then the maximum non-Gaussianity that can be obtained is − 5(4+ηm) 4 . This is the value of f NL that would be obtained if the perturbations in the growing mode phase froze immediately at the end of this phase [67,68]. In other words, the limiting value of f NL in standard scenario models is the value of f NL if one only had the growing phase.
In addition, if we fix η in the final phase to be η f , then the maximum value of f NL is . There is, however, no limit to the value of f NL if both η values are allowed to vary. The value of f NL can be made arbitrarily large by increasing the magnitude of both.

Inflaton Falls Models
In Fig. 2 (top left and top right), we plot the value of f NL for two inflaton falls models. In these cases the approximate plateau value is still obtained by using Eq. (4.1), except now substituting η 2 = 0 (since the growing phase is followed by slow-roll) and whatever the value of η is in the third phase of the model (the one with a large negative η value) for η 1 . We have verified that the plateau value of f NL does not depend at all on the large positive η value, or on the fact that the model begins in a phase of SR inflation.
This means that, so long as the slow-roll phase immediately after the growing phase is retained, there is a limit on the f NL that can be achieved in an inflaton falls model. Since η 2 is fixed to 0 in Eq. (4.1), the maximum f NL is obtained in the limit as η 1 → −∞. In this limit, f NL = 5 2 . If this slow roll phase is instead replaced by a different η > −3 phase then this value should be used for η 2 in Eq. (4.1) and values of f NL larger than 5 2 can, in principle, be obtained.

Repeating Models
In repeating models a faster than k 4 growth in the power spectrum is achieved thanks to the alternating SR and USR phases. However, once again, the non-Gaussianity has the same structure in repeating models as for the inflaton falls and standard scenario models. A long plateau in f NL extends over scales at which the peak in the power spectrum occurs, and the value of f NL here is given by substituting η 1 = −6 and η 2 = 0 into Eq. (4.1). This is the same non-Gaussianity that one would expect for a standard scenario model with a final transition from USR to SR.
Note that the value of f NL around the peak of such a model still only depends on the final large negative η-value and the η-value of the phase immediately following. A 0 → −6 → 0 → −6 → 0 model will have the same f NL around its peak as a 0 → −10 → 0 → −6 → 0, for example. The significance of the first negative η-value is for the scale dependence of the rise in the power spectrum, not for the value of f NL .

Linear Transitions
Up until this point we have only been considering models of inflation with instant transitions between phases. While this simplified our analysis, these kinds of transition are unphysical approximations to realistic, continuous transitions. To investigate how our results are affected if we don't assume instantaneous transitions, we performed bispectrum calculations for models with transitions of finite duration between phases where the value of η varied linearly in conformal time. This way the value of η is continuous throughout the evolution.
We say that these models have "linear" transitions, since the value of η evolves according to Eq. (2.7) during the transition. The setup and mode functions during the transitions were described in §2.4, while the bispectrum calculation proceeds by integrating Eq. (3.7) with the mode functions having a dependence on τ . What we notice is that, although the reduced bispectrum in the linear case has a similar form to that in the instant case, the value of f NL in the plateau region is always reduced in the linear case. The reduced value of f NL depends only on the duration of the transition between the two phases whose η-values would be used in Eq. (4.1) to calculate f NL in the instant case. The details of any other linear transitions in the model do not affect f NL . The longer the relevant transition lasts, the smaller the plateau value of f NL compared to its value in the instant case.
In Figure 4 we plot the dependence of f NL on the duration of the transition for a selection of models in the standard scenario. In most cases studied f NL is reduced to around 50% of its value for transitions lasting longer than an e-fold. In most cases f NL is reduced at most to 90% of its value for transitions lasting less than a fifth of an e-fold. For models with larger η-values however, the decline in f NL is steeper.
Our findings agree with earlier work that indicates that the value of f NL is very sensitive to the details of the transitions between phases [68]. Slow, smooth transitions are liable to wipe out large amounts of non-Gaussianity. So long as the final transition is sufficiently short (much less than an e-fold), however, our formula gives a good estimate for f NL around the peak of the power spectrum.

Conclusions
In this work we have shown that the level of non-Gaussianity around the peak of a power spectrum, as measured by f NL , can be predicted in general idealised models of inflation using a simple expression. This expression can be applied to estimate f NL in any model which has a correspondent in our idealised framework.
An f NL of O(1), or higher, can have a significant impact on PBH formation and the spectrum of SIGWs produced. In our framework there is, in principle, no limit to the magnitude of f NL that can be obtained. The larger the magnitude of the relevant η-values, the larger f NL . O(1) values of f NL can be obtained for fairly reasonable η-values; a model with η 1 = −8 and η 2 = 10 has f NL = 3.26, a model with η 1 = −30 and η 2 = 30 has f NL = 17.3.
Our idealised models are clearly not fully realistic, since transitions between phases occur instantaneously. When we smooth out transitions by giving them a non-zero duration, we find invariably that f NL is reduced. If the transitions are long and smooth, this is liable to wipe out a significant fraction of f NL , as predicted by our simple expression. So long as the transitions are short (much less than an e-fold), then our expression is still a good estimate for f NL . It is worth asking whether short transitions between phases are realistic when constructing models directly from the inflationary potential. The answer requires investigation beyond the scope of this work, but it is clear that for some transitions this likely can be realised, while for others it would require much more fine tuning. For example transitioning rapidly into a slow-roll phase requires the field velocity to match precisely that of the of the slow-roll predicted value. More likely is that the field enters a regime of the potential that supports slow-roll, but with incorrect velocity, and there is a decay time-period over which the field velocity decays or grows to reach the slow-roll attractor. On the other hand, moving from to a phase of very large positive or negative η can likely happen very rapidly. Such phases are defined by a large acceleration or deceleration of the field, and this only requires the potential to change suddenly such that it applies the required accelerating (or decelerating) force on the field. We defer more detailed investigation to future work, and conclude by reiterating we expect that our formula (2.14) for the tilt of the power spectrum (and the standard expression (2.15)), together with our formula for f NL (4.1), will help identify the important features of the spectra produced by models of inflation with a large peak, and aid in the designing of such models.

A Detailed explanation of scaling of power spectrum
In this appendix, we describe the scaling of the power spectrum induced by a fast growing stage.
We begin by re-writing Eq. (2.5) as which becomes, in the super-horizon limit kτ → 0 where the variables with tildes are combinations of A η and B η depending only on η. It is clear here that the constant mode is always present in theB η term, but time-dependent contributions can arise in both terms. In particular, for −3 < η < −1 the leading decaying mode comes from theÃ η term, whereas for η > −1, the slowest decaying mode comes from theB η contribution. However, when a growing mode is present (η < −3), it only arises in theÃ η , making it the crucial component for stages like USR. Matching a stage with constant leading behaviour (e.g. SR) to one with growing behaviour (e.g. USR) at τ = τ 1 , one finds the following solution for theÃ η term of the second stage where the barred variables are proportional to the tilde variables. It can be easily checked that, at leading order in the limit kτ 1 → 0, this expression always has the same scale dependence as the leading decaying term in the previous stage, as shown in Ref. [60]. Squaring the expression in Eq. (A.3) and multiplying by k 3 gives the scaling of the power spectrum P ∼ k α as A hidden assumption used in deriving the result shown above is that the coefficients A η 1 and B η 1 are scale-independent. This is true for general values of η in the case of Bunch-Davies initial conditions, but may not be correct if other transitions have occurred before the first stage considered here. This is particularly important if there are several stages with growing behaviour, as has been considered by Ref. [56], finding a scaling that does not obey Eq. (A.4). We will now proceed to explain the origin of this effect.
After an USR-like stage it is often useful to match to an SR-like stage, in which the scale dependence of the growing term is naturally passed to the constant mode in the SR-like stage. What is crucial is the fact that the coefficient of the constant mode in the SR-like stage (B η ) is the same as the leading decaying term (for η > −1). And if that coefficient gains an additional scale dependence due to a previous USR-like stage, then so does the leading decaying term, contrary to what happens with Bunch-Davies initial conditions. For this reason, when there is a second USR-like stage following the second SR-like stage, the scaling of the power spectrum becomes modified.
To demonstrate this in more detail, we perform these two additional matching calculations. We are assuming a scenario with three transitions in which the values of η obey η 1 , η 3 > −3 and η 2 , η 4 < −3. It is clear from Eq. (A.3), that one would need to know the scale dependence of both coefficients in the previous stage. Under the assumption that η 1 > −3, we findB η 2 ≈B η 1 k (η 2 −η 1 )/2 , (A.5) where once moreB η 1 is a place-holder variable that is proportional toB η 1 . It is not equal to the one shown in Eq. (A.3), but we are interested here only on how these coefficients depend on each other, not their specific expressions. After the second transition, the coefficients arẽ By multiplying by the scale dependence of the constant mode and neglecting theB η 2 term, we find k −3/2−η 3 /2B η 3 ≈Ā η 2 k 3/2+η 2 /2 , (A.8) where we can see that this mode inherits the same scale dependence of the previous growing mode, as expected (the right-hand side above is the same as the left-hand side of Eq. (A.3)).
To finalize, we put everything together into Eq. (A.3), sinceÃ η 4 obeys the same equation. For η 3 > −1, we find k 3/2+η 4 /2Ã η 4 ≈ k 2 k 3/2+η 2 /2Ā η 2 +B η 2 k −3/2−η 2 /2 , (A.9) where we can see that the contribution of the first growing stage is multiplied by an additional k 2 , thus potentially increasing the steepness of the final power spectrum. Performing the full substitution and multiplying by an additional k 3/2 results in k 3+η 4 /2Ã η 4 ≈B η 1 k 4−η 1 /2 +B η 1 k 2−η 1 /2 , (A.10) where we have also assumed η 1 > −1. Squaring this expression we see that the scaling of the power spectrum can now reach α = 8 − η 1 . In practice, the second term in Eq. (A.10) can dominate because these are super-horizon scales, unless the duration of the different stages is tuned. However, this tuning is possible and a power spectrum with a k 8 scaling can be found for η 1 = η 3 = 0 and η 2 = η 4 = −6 as shown in the main text, as well as first found by Ref. [56]. For that case, the conditions are where ∆N i denotes the duration in e-folds of the stage i. Therefore, to make this new scaling appear, one requires a very short first USR stage followed by an SR stage that is at least twice as long. There is a also a requirement that the final USR stage is sufficiently long, so that a sufficient range of scales is affected by both USR stages. In those scenarios, both scalings are typically present, as seen in Fig. 2.