A Density Spike on Astrophysical Scales from an N-Field Waterfall Transition

Hybrid inflation models are especially interesting as they lead to a spike in the density power spectrum on small scales, compared to the CMB, while also satisfying current bounds on tensor modes. Here we study hybrid inflation with $N$ waterfall fields sharing a global $SO(N)$ symmetry. The inclusion of many waterfall fields has the obvious advantage of avoiding topologically stable defects for $N>3$. We find that it also has another advantage: it is easier to engineer models that can simultaneously (i) be compatible with constraints on the primordial spectral index, which tends to otherwise disfavor hybrid models, and (ii) produce a spike on astrophysically large length scales. The latter may have significant consequences, possibly seeding the formation of astrophysically large black holes. We calculate correlation functions of the time-delay, a measure of density perturbations, produced by the waterfall fields, as a convergent power series in both $1/N$ and the field's correlation function $\Delta(x)$. We show that for large $N$, the two-point function is $<\delta t({\bf x})\,\delta t({\bf 0})>\,\propto\Delta^2(|{\bf x}|)/N$ and the three-point function is $<\delta t({\bf x})\,\delta t({\bf y})\,\delta t({\bf 0})>\,\propto\Delta(|{\bf x}-{\bf y}|)\Delta(|{\bf x}|)\Delta(|{\bf y}|)/N^2$. In accordance with the central limit theorem, the density perturbations on the scale of the spike are Gaussian for large $N$ and non-Gaussian for small $N$.

Inflation, a phase of accelerated expansion in the very early universe thought to be driven by one or several scalar fields, is our paradigm of early universe cosmology [1][2][3][4]. It naturally explains the large scale homogeneity, isotropy, and flatness of the universe. Moreover, its basic predictions of approximate scale invariance and small non-Gaussianity in the ∼ 10 −5 level departures from homogeneity and isotropy are in excellent agreement with recent CMB data [5,6] and large scale structure.
While the basic paradigm of inflation is in excellent shape, no single model stands clearly preferred. Instead the literature abounds with various models motivated by different considerations, such as string moduli, supergravity, branes, ghosts, Standard Model, etc [7][8][9][10][11][12][13][14][15][16][17]. While the incoming data is at such an impressive level that it can discriminate between various models and rule out many, such as models that overpredict non-Gaussianity, it is not clear if the data will ever reveal one model alone. An important way to make progress is to disfavor models based on theoretical grounds (such as issues of unitarity violation, acausality, etc) and to find a model that is able to account for phenomena in the universe lacking an alternate explanation. It is conceivable that some version of the so-called "hybrid inflation" model may account for astrophysical phenomena, for reasons we shall come to. The hybrid inflation model, originally proposed by Linde [18], requires at least two fields. One of the fields is light and another of the fields is heavy (in Hubble units). The light field, called the "timer", is at early times slowly rolling down a potential hill and generates the almost scale invariant spectrum of fluctuations ob-served in the CMB and in large scale structure. The heavy field, called the "waterfall" field, has an effective mass that is time-dependent and controlled by the value of the timer field. The waterfall field is originally trapped at a minimum of its potential, but as its effective masssquared becomes negative, a tachyonic instability follows, leading to the end of inflation; an illustration is given in Figure 1. The name "hybrid inflation" comes from the fact that this model is a sort of hybrid between a chaotic inflation model and a symmetry breaking inflation model.
As originally discussed in Refs. [19,20], one of the most fascinating features of hybrid models is that the tachyonic behavior of the waterfall field leads to a sharp "spike" in the density power spectrum. This could seed primordial black holes [21][22][23][24][25]. For generic parameters, the length scale associated with this spike is typically very small. However, if one could find a parameter regime where the waterfall phase were to be prolonged, lasting for many e-foldings, say N w ∼ 30 − 40, then this would lead to a spike in the density perturbations on astrophysically large scales (but smaller than CMB scales). This may help to account for phenomena such as supermassive black holes or dark matter, etc. Of course the details of all this requires a very careful examination of the spectrum of density perturbations, including observational constraints.
Ordinarily the spectrum of density perturbations in a given model of inflation is obtained by decomposing the inflaton field into a homogeneous part plus a small inhomogeneous perturbation. However, for the waterfall fields of hybrid inflation, this approach fails since classically the waterfall field would stay forever at the top of a ridge in its potential. It is the quantum perturbations themselves that lead to a non-trivial evolution of the waterfall field, and therefore the quantum perturbations cannot be treated as small. Several approximations have been used to deal with this problem [19,20,[26][27][28][29][30][31][32][33][34][35].
Here we follow the approach presented by one of us recently in Ref. [36], where a free field time-delay method was used, providing accurate numerical results.
In this paper, we generalize the method of Ref.
[36] to a model with N waterfall fields sharing a global SO(N ) symmetry. A model of many fields may be natural in various microscopic constructions, such as grand unified models, string models, etc. But apart from generalizing Ref. [36] to N fields, we also go much further in our analysis: we derive explicit analytical results for several correlation functions of the so-called time-delay. We formulate a convergent series expansion in powers of 1/N and the field's correlation function ∆(x). We find all terms in the series to obtain the two-point correlation function of the time-delay for any N . We also obtain the leading order behavior at large N for the three-point function timedelay, which provides a measure of non-Gaussianity. We find that the non-Gaussianity is appreciable for small N and suppressed for large N .
We also analyze in detail constraints on hybrid inflation models. We comment on how multiple fields avoids topological defects, which is a serious problem for low N models. However, the most severe constraint on hybrid models comes from the requirement to obtain the observed spectral index n s . We show that at large N , it is easier to engineer models that can fit the observed n s , while also allowing for a prolonged waterfall phase. This means that large N models provide the most plausible way for the spike to appear on astrophysically large scales and be compatible with other constraints.
With regards to tensor modes, the confirmation of the recent detection and amplitude of primordial gravitational waves [37] will depend largely on understanding dust foregrounds [38]. So we will only briefly discuss possibilities of obtaining observable tensor modes from hybrid inflation in Section V.
The organization of our paper is as follows: In Section II we present our hybrid inflation model and discuss our approximations. In Section III we present the timedelay formalism, adapting the method of Ref. [36] to N fields. In section IV we derive a series expansion for the two-point function, we derive the leading order behavior of the three-point function, and we derive results in k-space. In Section V we present constraints on hybrid models, emphasizing the role that N plays. In Section VI we discuss and conclude. Finally, in the Appendices we present further analytical results.

II. N FIELD MODEL
The model consists of two types of fields: The timer field ψ that drives the first slow-roll inflation phase, and the waterfall field φ that becomes tachyonic during the second phase causing inflation to end. In many hybrid models, φ is comprised of two components, a complex field, but here we allow for N real components φ i . We assume the components share a global SO(N ) symmetry, and so it is convenient to organize them into a vector For the special case N = 2, this can be organized into a complex field by writing φ complex = (φ 1 + i φ 2 )/ √ 2. The dynamics is governed by the standard two derivative action for the scalar fields ψ, φ minimally coupled to gravity as follows (signature + − −−) The potential V is given by a sum of terms: V 0 providing false vacuum energy, V ψ (ψ) governing the timer field, V φ (φ) governing the waterfall field, and V int (ψ, φ) governing their mutual interaction, i.e., During inflation we assume that the constant V 0 dominates all other terms. The timer field potential V ψ and the waterfall field potential V φ can in general be complicated. In general, they are allowed to be non-polynomial functions as part of some low energy effective field theory, possibly from supergravity or other microscopic theories. For our purposes, it is enough to assume an extremum at ψ = φ = 0 and expand the potentials around this extremum as follows: The timer field is assumed to be light m ψ < H and the waterfall field is assumed to be heavy m 0 > H, where H is the Hubble parameter. In the original hybrid model, this quadratic term for V ψ was assumed to be the entire potential. This model leads to a spectrum with a spectral index n s > 1 and is observationally ruled out. Instead we need higher order terms, indicated by the dots, to fix this problem. This also places constraints on the mass of the timer field m ψ , which has important consequences. We discuss these issues in detail in Section V. As indicated by the negative mass-squared, the waterfall field is tachyonic around φ = ψ = 0. This obviously cannot be the entire potential because then the potential would be unbounded from below. Instead there must be higher order terms that stabilize the potential with a global minimum near V ≈ 0 (effectively setting the late-time cosmological constant).
The key to hybrid inflation is the interaction between the two fields. For simplicity, we assume a standard dimension 4 coupling of the form This term allows the waterfall field to be stabilized at φ = 0 at early times during slow-roll inflation when ψ is displaced away from zero, and then to become tachyonic once the timer field approaches the origin; this is illustrated in Figure 1.

A. Approximations
We assume that the constant V 0 is dominant during inflation, leading to an approximate de Sitter phase with constant Hubble parameter H. By assuming a flat FRW background, the scale factor is approximated as a = exp(Ht).
At early times, ψ is displaced from its origin, so φ = 0. This means that we can approximate the ψ dynamics by ignoring the back reaction from φ. The fluctuations in ψ establish nearly scale invariant fluctuations on large FIG. 1. An illustration of the evolution of the effective potential for the waterfall field φ as the timer field ψ evolves from "high" values at early times, to ψ = ψc, and finally to "low" values at late times. In the process, the effective mass-squared of φ evolves from positive, to zero, to negative (tachyonic).
scales, which we shall return to in Section V. However, for the present purposes it is enough to treat ψ as a classical, homogeneous field ψ(t). We make the approximation that we can neglect the higher order terms in the potential V ψ in the transition era, leading to the equation of motionψ Solving this equation for ψ(t), we insert this into the equation for φ i . We allow spatial dependence in φ i , and ignore, for simplicity, the higher order terms in V φ , leading to the equation of motion Here we have identified an effective mass-squared for the waterfall field of where the dimension 4 coupling g has been traded for the parameter ψ c as g 2 = m 2 0 /ψ 2 c . The quantity ψ c has the physical interpretation as the "critical" value of ψ such that the effective mass of the waterfall field passes through zero. So at early times for ψ > ψ c , then m 2 φ > 0 and φ i is trapped at φ i = 0, while at late times for ψ < ψ c , then m 2 φ < 0 and φ i is tachyonic and can grow in amplitude, depending on the mode of interest; Figure 1 illustrates these features.

B. Mode Functions
Since we ignore the back-reaction of φ onto ψ and since we treat ψ as homogeneous in the equation of motion for φ (eq. (9)), then by passing to k-space, all modes are decoupled. Each waterfall field φ i can be quantized and expanded in modes in momentum space as follows where c † k (c k ) are the creation (annihilation) operators, acting on the φ i = 0 vacuum. By assuming an initial Bunch-Davies vacuum for each φ i , the mode functions u k are the same for all components due to the SO(N ) symmetry. To be precise, we assume that at asymptotically early times the mode functions are the ordinary Minkowski space mode functions, with the caveat that we need to insert factors of the scale factor a to convert from physical wavenumbers to comoving wavenumbers, i.e., , at early times.
At late times the mode functions behave very differently.
Since the field becomes tachyonic, the mode functions grow exponentially at late times. The transition depends on the wavenumber of interest. The full details of the mode functions were explained very carefully in Ref. [36], and the interested reader is directed to that paper for more information.
Since we are approximating φ i as a free field here, its fluctuations are entirely Gaussian and characterized entirely by its equal time two-point correlation function φ i (x) φ j (y) . Passing to k-space, and using statistical isotropy and homogeneity of the Bunch-Davies vacuum, the fluctuations are equally well characterized by the socalled power spectrum P φ (k), defined through This means the power spectrum is It is simple to show that the power spectrum is related to the mode functions by In Figure 2 we plot a rescaled version of P φ , where we divide out by the root-mean-square (rms) of φ as defined in the next Section.

III. THE TIME-DELAY
We would now like to relate the fluctuations in the waterfall field φ to a fluctuation in a physical observable, namely the density perturbation. An important step in this direction is to compute the so-called "time-delay" δt(x); the time offset for the end of inflation for different parts of the universe. This causes different regions of the universe to have inflated more than others, creating a difference in their densities (though we will not explicitly compute δρ/ρ here). This basic formalism was first introduced by Hawking [39] and by Guth and Pi [40], and has recently been reviewed in Ref. [41]. In the context of hybrid inflation, it was recently used by one of us in Ref. [36]. It provides an intuitive and straightforward way to calculate primordial perturbations and we now use this to study perturbations established by the N waterfall fields.
In its original formulation, the time-delay formalism starts by considering a classical homogeneous trajectory φ 0 = φ 0 (t), and then considers a first order perturbation around this. At first order, one is able to prove that the fluctuating inhomogeneous field φ(x, t) is related to the classical field φ 0 , up to an overall time offset δt(x), In the present case of hybrid inflation, the waterfall field is initially trapped at φ = 0 and then once it becomes tachyonic, it eventually falls off the hill-top due to quantum fluctuations. This means that there is no classical trajectory about which to expand. Nevertheless, we will show that, to a good approximation, the field φ(x, t) is well described by an equation of the form (16), for a suitably defined φ 0 . The key is that all modes of interest grow at the same rate at late times. Further information of the time evolution of the mode functions can be found in Ref. [36].
To show this, we need to compute the evolution of the field φ according to eq. (9). This requires knowing m φ (t), which in turn requires knowing ψ(t). By solving eq. (8) for the timer field and dispensing with transient behavior, we have where (note p > 0). We have set the origin of time t = 0 to be when ψ = ψ c and assume ψ > ψ c at early times. Substituting this solution into m φ (t), we can, in principle, solve eq. (9). In general, the solution is somewhat complicated with a non-trivial dependence on wavenumber. However, at late times the behavior simplifies. Our modes of interest are super-horizon at late times. For these modes, the gradient term is negligible and the equation of motion reduces tö So each mode evolves in the same way at late times.
Treating m φ (t) as slowly varying (which is justified because the timer field mass m ψ < H and so p is small), we can solve for φ i at late times t in the adiabatic approximation. We obtain where Here t 0 is some reference time. For t > 0, we have λ(t) > 0, so the modes grow exponentially in time. Later in Section V we explain that in fact λ is roughly constant in the latter stage of the waterfall phase, i.e., the exp(−2 p t) piece becomes small. We now discuss fluctuations in the time at which inflation ends. For convenience, we define the reference time t 0 to be the time at which the rms value of the field reaches φ end ; the end of inflation where we have included a factor of N to account for all fields, allowing φ rms to refer to fluctuations in a single component φ i , i.e., φ 2 rms = φ 2 i (0) . In terms of the power spectrum, it is If we were to include arbitrarily high k, this would diverge quadratically, which is the usual Minkowski space divergence. However, our present analysis only applies to modes that are in the growing regime. For these modes, P φ (k) falls off exponentially with k and there is no problem. So in this integral, we only include modes that are in the asymptotic regime, or, roughly speaking, only super-horizon modes. Using eq. (20), we can express the field φ i (x, t) at time If t is chosen to be the time t end (x) at which inflation ends at each point in space, then φ · φ x, t end (x) = φ 2 end = N φ 2 rms (t 0 ), and the above equation becomes which can be solved for the time-delay field δt( This finalizes the N component analysis of the timedelay, generalizing the two component (complex) analysis of Ref. [36].

IV. CORRELATION FUNCTIONS
We now derive expressions for the two-point and threepoint correlation functions of the time-delay field. To do so, it is convenient to introduce a rescaled version of the correlation function ∆(x) defined through By definition ∆(0) = 1, and as we vary x, ∆ covers the interval ∆(x) ∈ (0, 1]. (Ref. [36] used a different convention where ∆ covers the interval ∆(x) ∈ (0, 1/2]).
In Figure 3, we show a plot of ∆ at t end as a function of x, measured in Hubble lengths (H −1 ), for various combinations of masses.

A. Two-Point Function
We now express the time-delay correlation functions as a power series in ∆ and 1/N . An alternative derivation of the power spectra of the time-delay field in terms of an integral, which is closer to the language of Ref. [36], can be found in Appendix C.
Using the above approximation for δt in eq. (26), the two-point correlation function for the time-delay is where we have used the abbreviated notation φ x ≡ φ(x).
The two-point function will include a constant (independent of x) for a non-zero δt . This can be reabsorbed into a shift in t 0 , whose dependence we have suppressed here, and so we will ignore the constant in the following computation. This means that we will compute the connected part of the correlation functions. We would like to form an expansion, but we do not have a classical trajectory about which to expand. Instead we use the following idea: we recognize that φ · φ should be centralized around its mean value of N φ 2 rms , plus relatively small fluctuations at large N . This means that it is convenient to write and treat the term in parenthesis on the right as small, as it represents the fluctuations from the mean. This allows us to Taylor expand the logarithm in powers of with Φ = 0. Now recall that the series expansion of the logarithm for small Φ is allowing us to expand eq. (28) to any desired order in Φ.
The leading non-zero order is quadratic ∼ Φ 2 . It is where we are implicitly summing over indices in the second line (for simplicity, we will place all component indices (i, j) downstairs). Using Wick's theorem to perform the four-point contraction, we find the result where x ≡ |x|. This provides the leading approximation for large N , or small ∆. This should be contrasted to the RSG approximation used in Refs. [19,20] and summarized in [36], in which the correlation function is approximated as ∼ ∆, rather than ∼ ∆ 2 , at leading order. For brevity, we shall not go through the result at each higher order here, but we report on results at higher order in Appendix A. By summing those results to different orders, we find (where ∆ = ∆(x) here). We find that various cancellations have occurred. For example, the −8∆ 2 /N 2 term that enters at cubic order (see Appendix A) has canceled. Note that at a given order in 1/N there are various powers of ∆ 2 . However, by collecting all terms at a given power in ∆ 2 , we can identify a pattern in the value of its coefficients as functions of N . We find We then identify the entire series as where the coefficients are with a b the binomial coefficient. This series is convergent for any N and ∆ ∈ (0, 1], and is one of our central results. For the case of a single scalar or a doublet (complex field), the full series organizes itself into known functions. For N = 1 we find For N = 2 we find (2λ) 2 δt(x)δt(0) = ∆ 2 (x) + ∆ 4 (x) 2 2 + ∆ 6 (x) 3 2 + . . . = Li 2 (∆ 2 (x)), where Li s (z) is the polylogarithm function. We also find that for any even value of N , the series is given by the polylogarithm function plus a polynomial in ∆; this is described in Appendix B.
Using the power series, we can easily obtain plots of the two-point function for any N . For convenience, we plot the re-scaled quantity N (2λ) 2 δt(x)δt(0) as a function of ∆ in Figure 4 (top) for different N . We see convergence of all curves as we increase N , which confirms that the leading behavior of the (un-scaled) two-point function is ∼ 1/N . In Figure 4 (bottom) we plot N δt(x)δt(0) as a function of x for different masses and two different N .

B. Three-Point Function
The three-point function is given by a simple modification of eq. (28), namely Here we will work only to leading non-zero order, which in this case is cubic. We expand the logarithms as before to obtain where "perms" is short for permutations under interchanging x, y, 0. Using Wick's theorem to perform the various contractions, we find the result We would now like to use the three-point function as a measure of non-Gaussianity. For a single random variable, a measure of non-Gaussianity is to compute a dimensionless ratio of the skewness to the 3/2 power of the variance. For a field theory, we symmetrize over variables, and define the following measure of non-Gaussianity in position space S ≡ δt(x)δt(y)δt(0) δt(x)δt(y) δt(x)δt(0) δt(y)δt(0) .
(Recall that we are ignoring δt , as it can be just reabsorbed into t 0 , so the three-point and two-point functions are the connected pieces). Computing this at the leading order approximation using eqs. (33, 42) (valid for large N , or small ∆), we find Curiously, the dependence on x, y has dropped out at this order. We see that for small N there is significant non-Gaussianity, while for large N the theory becomes Gaussian, as expected from the central limit theorem.

C. Momentum Space
Let us now present our results in k-space. We shall continue to analyze the results at high N , or small ∆, which allows us to just include the leading order results. For the two-point function, we define the power spectrum through We use eq. (33) and Fourier transform to k-space using the convolution theorem. To do so it is convenient to introduce a dimensionless fieldφ i ≡ φ i /φ rms and introduce the corresponding power spectrum Pφ(k) = P φ (k)/φ 2 rms = |u k | 2 /φ 2 rms , which is the Fourier transform of ∆(x). We then find the result A dimensionless measure of fluctuations is the following The factor of k 3 /(2π 2 ) is appropriate as this gives the variance per log interval: (Hδt) 2 = d ln k P δt (k). By studying eq. (46), one can show P δt ≈ const for small k and falls off for large k. This creates a spike in P δt (k) at a finite k * and its amplitude is rather large. (This is to be contrasted with the usual fluctuations in de Sitter space, P δt (k) ∝ 1/k 3 , making P δt (k) approximately scale invariant.) We see that the amplitude of the spike scales as ∼ 1/N , and so it is reduced for large N . A plot of P δt (k) is given in Figure 5.
For the three-point function, we define the bispectrum through where we have indicated that the bispectrum only depends on the magnitude of the 3 k-vectors, with the constraint that the vectors sum to zero. We use eq. (42) and Fourier transform to k-space, again using the convolution theorem. We find (49) To measure non-Gaussianity in k-space, it is conventional to introduce the dimensionless f N L parameter, defined as By substituting the above expressions for P δt and B δt , we see that f N L is independent of N at this leading order. However, this belies the true dependence of non-Gaussianity on the number of fields. This is because f N L is a quantity that can be large even if the non-Gaussianity is relatively small (for example, on CMB scales, any f N L smaller than O(10 5 ) is a small level of non-Gaussianity). Instead a more appropriate measure of non-Gaussianity in k-space is to compute some ratio of the bispectrum to the 3/2 power of the power spectrum, analogous to the position space definition in eq. (43). For the simple equilateral case, k 1 = k 2 = k 3 , we define where we have inserted a factor of k 3/2 /( √ 2 π) from measuring the fluctuations per log interval. Using eqs. (46,49), we see that F N L ∝ 1/ √ N , as we found in position space. In Figure 6, we plot this function. We note that although the non-Gaussianity can be large, the peak is on a length scale that is small compared to the CMB and so it evades recent bounds [6].

V. CONSTRAINTS ON HYBRID MODELS
Hybrid inflation models must satisfy several observational constraints. Here we discuss these constraints, including the role that N plays, and discuss the implications for the scale of the density spike.

A. Topological Defects
The first constraint on hybrid models concerns the possible formation of topological defects. Since the waterfall A factor of 6/5 is often included when studying the gauge invariant quantity ζ that appears in cosmological perturbation theory, but it does not concern us here. field starts at φ = 0 and then falls to some vacuum, it spontaneously breaks a symmetry. For a single field N = 1, this breaks a Z 2 symmetry; see Figure 1. For multiple fields N > 1, this breaks an SO(N ) symmetry. N = 1 leads to the formation of domain walls, which are clearly ruled out observationally, so these models are strongly disfavored; N = 2 leads to the formation of cosmic strings, which have not been observed and if they exist are constrained to be small in number. This would require a very large number of e-foldings of the waterfall phase to make compatible with observations, and seems unrealistic; N = 3 leads to the formation of monopoles, which are somewhat less constrained; N = 4 leads to the formation of textures, which are relatively harmless; N > 4 avoids topological defects altogether. So choosing several waterfall fields is preferred by current constraints on topological defects.

B. Inflationary Perturbations
Inflation generates fluctuations on large scales which are being increasingly constrained by data. An important constraint on any inflation model is the bound on the tensor-to-scalar ratio r. CMB measurements from Planck place an upper bound on tensor modes of r < 0.11 (95% confidence) [6]. The amplitude of tensor modes is directly set by the energy density during inflation. Typical hybrid models are at relatively low energy scales, without the need for extreme fine tuning, and so they immediately satisfy this bound. Recent data by the BICEP2 experiment [37] claim a detection of gravitational waves with r ∼ 0.2 at a high confidence level. Going back to the simplest potential of hybrid inflation V (ψ) = V 0 + 1 2 m 2 ψ ψ 2 . .., we have so far considered the regime V 0 ≫ 1 2 m 2 ψ ψ 2 , which generates a very small amount of tensor modes. However one can operate in a regime with V 0 1 2 m 2 ψ ψ 2 , where the tensor to scalar ratio can be pushed to be closer to O(0.1). Furthermore, new realizations of hybrid inflation [42,43] can generate appreciable tensor modes with appropriate parameter choices. Since it is quite conceivable that the tensor mode signal will be lowered after subtraction of the dust foreground [38], we will leave a complete discussion of producing the correct value of r in hybrid inflation models for future work.
Although the detection of tensor modes is still a debatable issue, scalar modes are pinned down with great accuracy. The tilt of the scalar mode spectrum is characterized by the primordial spectral index n s . WMAP [5] and Planck measurements [6] place the spectral index near n s,obs ≈ 0.96, giving a red spectrum. Here we examine the constraints imposed on hybrid models in order to obtain this value of n s . The tilt on large scales is determined by the timer field ψ. For low scale models of inflation, such as hybrid inflation, the prediction for the spectral index is where This is to be evaluated N e e-foldings from the end of inflation, where N e = 50−60 in typical models. Combining the above equations, we need to satisfy V ′′ ψ ≈ −0.06 H 2 . If we take V ψ = m 2 ψ ψ 2 /2, then V ′′ ψ > 0, and n s > 1, which is ruled out. So we need higher order terms in the potential to cause it to become concave down at large values of ψ where η is evaluated, while leaving the quadratic approximation for V ψ valid at small ψ. For most reasonable potential functions, such as potentials that flatten at large field values, we expect |V ′′ ψ | m 2 ψ . So this suggests a bound which can only be avoided by significant fine tuning of the potential. Hence although the timer field is assumed light (m ψ < H), it cannot be extremely light.

C. Implications for Scale of Spike
The length scale associated with the spike in the spectrum is set by the Hubble length during inflation H −1 , red-shifted by the number of e-foldings of the waterfall phase N w . Since the Hubble scale during inflation is typically microscopic, we need the duration of the waterfall phase N w to be significant (e.g., N w ∼ 30 − 40) to obtain a spike on astrophysically large scales. Here we examine if this is possible.
Since we have defined t = 0 to be when the transition occurs (ψ = ψ c ), then N w = Ht with final value at t = t end . To determine the final value, we note that modes grow at the rate λ, derived earlier in eq. (21). For m ψ < H, we can approximate the parameter p (eq. (18)) as p ≈ m 2 ψ /3H. Using the above spectral index bound in eq. (55), we see that for significantly large N w (e.g., N w ∼ 30 − 40) the exponential factor is somewhat small and we will ignore it here. In this regime, the dimensionless growth rate λ/H can be approximated as a constant The typical starting value for φ is roughly of order H (de Sitter fluctuations) and the typical end value for φ is roughly of order M P l (Planck scale). For self consistency, φ must pass from its starting value to its end value in N w e-foldings with rate set by λ/H. This gives the approximate value for N w as This has a clear consequence: If we choose m 0 ≫ H, as is done in some models of hybrid inflation, then H/λ ≪ 1. So unless we push H to be many, many orders of magnitude smaller then M P l , then N w will be rather small. This will lead to a spike in the spectrum on rather small scales and possibly ignorable to astrophysics. Note that if we had ignored the spectral index bound that leads to eq. (55), then we could have taken m ψ arbitrarily small, leading to an arbitrarily small p value. In this (unrealistic) limit, it is simple to show So by taking m ψ arbitrarily small, λ could be made small, and N w could easily be made large. However, the existence of the spectral index bound essentially forbids this, requiring us to go in a different direction.
The only way to increase N w and satisfy the spectral index bound on m ψ is to take m 0 somewhat close to H. This allows H/λ to be appreciable from eq. (57). For instance, if we set m 0 = 1.3 H, then H/λ ≈ 2. If we A better approximation comes from tracking the full time dependence of λ and integrating the argument of the exponential exp( t dt ′ λ(t ′ )), but this approximation suffices for the present discussion.
then take H just a few orders of magnitude below M P l , say H ∼ 10 −6,7 M P l , which is reasonable for inflation models, we can achieve a significant value for N w . This will lead to a spike in the spectrum on astrophysically large scales, which is potentially quite interesting. It is possible that there will be distortions in the spectrum by taking m 0 close to H, but we will not explore those details here. However, there is an important consequence that we explore in the next subsection.

D. Eternal Inflation
Since we are being pushed towards a somewhat low value of m 0 , near H, we need to check if the theory still makes sense. One potential problem is that the theory may enter a regime of eternal inflation. This could occur for the waterfall field at the hill-top. This would wipe out information of the timer field, which established the approximately scale invariant spectrum on cosmological scales.
The boundary for eternal inflation is roughly when the density fluctuations are O(1), and this occurs when the fluctuations in the time delay are (H δt) 2 = O(1). To convert this into a lower bound on m 0 , let us imagine that m 0 is even smaller than H. In this regime, the growth rate λ can be estimated using eq. (57) as λ ∼ m 2 0 /H 2 . Using eq. (33) this gives (H δt) 2 ∼ H 4 /(N m 4 0 ). This implies that eternal inflation occurs when the waterfall mass is below a critical value m c 0 , which is roughly So when N ∼ 1 we cannot have m 0 near H, because we then enter eternal inflation. On the other hand, for large N we are allowed to have m 0 near H and avoid this problem. This makes sense intuitively, because for many fields it is statistically favorable for at least one of the fields to fall off the hill-top, causing inflation to end. Hence large N is more easily compatible with the above set of constraints than low N .

VI. DISCUSSION AND CONCLUSIONS
In this work we studied density perturbations in hybrid inflation caused by N waterfall fields, which contains a spike in the spectrum. We derived expressions for correlation functions of the time-delay and constrained parameters with observations. Density Perturbations: We derived a convergent series expansion in powers of 1/N and ∆(x), the dimensionless correlation function for the field, for the two-point function of the time-delay for any N , and the leading order behavior of the three-point function of the time-delay for large N . These correlation functions are well approximated by the first term in the series for large N (even for N = 2 the leading term is moderately accurate to ∼ 30%, and much more accurate for higher N ). In this regime, the fluctuations are suppressed, with two-point and three-point functions given by Although this reduces the spike in the spectrum, for any moderate value of N , such as N = 3, 4, 5, the amplitude of the spike is still quite large (orders of magnitude larger than the ∼ 10 −5 level fluctuations on larger scales relevant to the CMB), and may have significant astrophysical consequences. Also, the relative size of the three-point function to the 3/2 power of the two-point function scales as ∼ 1/ √ N . In accordance with the central limit theorem, the fluctuations become more Gaussian at large N . This will make the analysis of the subsequent cosmological evolution more manageable, as this provides a simple spectrum for initial conditions. We note that since we are considering small scales compared to the CMB, then this non-Gaussianity evades Planck bounds [6].
Constraints: We mentioned that hybrid models avoid topological defects for large N , while tensor mode constraints can be satisfied in different models. A very serious constraint on hybrid models comes from the observed spectral index n s ≈ 0.96, which requires the potential to flatten at large field values. One consequence of this is that the timer field mass m ψ needs to be only a little smaller than the Hubble parameter during inflation H, or else the model is significantly fine tuned. For a large value of the waterfall field mass m 0 , this would imply a large growth rate of fluctuations, a rapid termination of inflation, and in turn a density spike on very small scales. Otherwise, we need to make the waterfall field mass m 0 somewhat close to H, but this faces problems with eternal inflation. However, by using a large number of waterfall fields N , it is safer to make the waterfall field mass m 0 somewhat close to H. This reduces the growth rate of fluctuations, prolonging the waterfall phase for many e-foldings.
Thus large N presents a plausible setup to establish a spike in the density perturbations on astrophysically large length scales that is consistent with other constraints.
Outlook: It may be possible that these perturbations seed primordial black holes, which may be relevant to seeding supermassive black holes, or an intriguing form of dark matter. An investigation into these topics is underway. It would be important to fully explore the eternal inflation bound and the effects on the spectrum for relatively light waterfall field masses. Finally, it would be of interest to try to embed these large N models into fundamental physics, such as string theory, and to explore reheating [44][45][46][47] and baryogenesis [48][49][50][51][52][53][54] in this framework.
Earlier we computed the two-point correlation function for the time-delay at quadratic order, and then stating the results at all orders. Here we mention the results order by order.

a. Cubic Order
At cubic order ∼ Φ 3 we find When N is even the expression always involves the polylog function that we found for N = 2, plus corrections that depend on N . We find that the form of the answer is For N = 4, we find For N = 6, we find For N = 8, we find Here we write the two-point function as a multidimensional integral. It is convenient now to switch to a vector notation thus making the components of φ explicit. φ(0, t) ≡ φ x ≡ (X 1 , X 2 , X 3 , . . . , X N ), φ(x, t) ≡ φ 0 ≡ (X N +1 , X N +2 , X N +3 , . . . , X 2N ), (C2) X ≡ (X 1 , X 2 , X 3 , . . . , X 2N ).
The average value of a function F of a random variable X with probability distribution function p(X) is given by Since we are using a free field approximation, X follows a joint Gaussian distribution where is the correlation matrix. The components of Σ can be easily calculated using the commutation relations for the creation and annihilation operators in φ(x, t), from Eq. (11). Due to the high degree of symmetry the matrix itself has a very simple structure: Finally, we can write the two-point function as: (2λ) 2 δt(x)δt(0) Then we can Taylor expand the integrand in powers of ∆. We see that all odd powers of ∆ vanish, leaving only even powers in the expansion. Performing these integrals term by term in the Taylor expansion, leads to the results reported earlier in eqs. (35,36,37).