Stochastic dynamics of multi-waterfall hybrid inflation and formation of primordial black holes

We show that a hybrid inflation model with multiple waterfall fields can result in the formation of primordial black holes (PBHs) with an astrophysical size, by using an advanced algorithm to follow the stochastic dynamics of the waterfall fields. This is in contrast to the case with a single waterfall field, where the wavelength of density perturbations is usually too short to form PBHs of the astrophysical scale (or otherwise PBHs are overproduced and the model is ruled out) unless the inflaton potential is tuned. In particular, we demonstrate that PBHs with masses of order 1020 g can form after hybrid inflation consistently with other cosmological observations if the number of waterfall fields is about 5 for the case of instantaneous reheating. Observable gravitational waves are produced from the second-order effect of large curvature perturbations as well as from the dynamics of texture or global defects that form after the waterfall phase transition.


Introduction
The seeds of cosmological large-scale structures are generated by curvature perturbations during inflation, and their signals are imprinted and observed as the temperature anisotropies of the cosmic microwave background (CMB).The spectrum of curvature perturbations is not exactly flat, which is consistent with the prediction of slow-roll inflationary scenario [1].However, there are several puzzles in the formation of small-scale structures, such as the formation of supermassive black holes (SMBHs) [2,3] and black holes (BHs) observed by the LIGO-Virgo-KAGRA colllaboration [4,5].Although their evolution is highly non-linear, it is said that their number densities are in tension with the standard scenario of BH formation in astrophysics.One therefore needs to consider a primordial origin of such BHs called primordial black holes (PBHs) [6][7][8] (see also the recent review article on PBH [9]).It is also discussed that astrophysical size PBHs of order 10 17-23 g can explain the dark matter (DM) in the Universe [10] (see, e.g., Refs.[11][12][13][14]).Among several scenarios of the formation of PBHs, it would be minimal to consider the formation of those small-scale structures by generating large curvature perturbations in the small scales during inflation (see, e.g., Refs.[15][16][17][18][19][20] and also references in Ref. [9]). 1 The curvature perturbations are converted to density perturbations after inflation, and the dense regions collapse to form PBHs when the corresponding wavelength enters the horizon [6][7][8].
To generate large curvature perturbations, a bump or at least a very flat region for inflaton potential is required [16-18, 20, 31-49].This can be naturally realised by hybrid inflation models because a waterfall field has a vanishing mass at the waterfall phase transition [50].Remarkably, it can generate large density perturbations at small scales that result in the formation of PBHs [51][52][53].Since large density perturbations are generated during the waterfall phase transition, the waterfall field experiences the stochastic dynamics that requires detailed numerical simulations (see Refs. [54][55][56][57][58][59][60][61][62][63] for the first papers on the subject) and was omitted in some literature.Several studies on hybrid inflation models reveal that the stochastic effect on the waterfall fields significantly alters the predictions of hybrid inflation models (see, e.g., Ref. [64]).Ref. [65] provided an analytic method to calculate the spectrum of curvature perturbations in a model with a single waterfall field, considering the stochastic effects under some approximations.Their analytic formula is partially consistent with the numerical simulations for the stochastic dynamics performed in Ref. [66].These studies are performed in models with a single waterfall fields, and we thereby concluded that the peak amplitude of curvature perturbations is too large to provide consistent astrophysical consequences even though the stochastic effect is considered. 2 Our recent work [67] shows that this problem can be evaded by considering a higher-order term for the inflaton potential, though it requires a tuning in a parameter.As another aspect, the hybrid inflation model with a single waterfall field suffers from the so-called domain-wall problem because the Z 2 symmetry is spontaneously broken after the waterfall phase transition.This problem can be evaded by introducing a small explicit Z 2 -symmetry breaking term, which however may affect the dynamics of waterfall phase transition (see, e.g., Ref. [68]).
In this paper, we analyze a hybrid inflation model with multiple waterfall fields, considering the stochastic effect on the waterfall fields.They are assumed to be symmetric under O(N ) symmetry.The idea to introduce multiple waterfall fields was considered in Ref. [69], although they did not consider the stochastic dynamics of waterfall fields.We show that the spectrum of curvature perturbations can be reduced with a fixed peak frequency for this case.Alternatively, the peak frequency can be reduced for a fixed peak amplitude of curvature perturbations.This can result in the formation of PBHs at astrophysically interesting mass scales when the number of waterfall fields is sufficiently large.We employ a sophisticated numerical method proposed and developed in Refs.[70][71][72][73][74] to incorporate the stochastic dynamics of multiple waterfall fields.We also confirm that the numerical results can be qualitatively understood by the analytical calculation of Ref. [65], extending it for the case with multiple waterfall fields.In particular, the peak amplitude of curvature perturbations is proportional to the inverse of the number of waterfall fields.We provide a semi-analytic formula for the spectrum of curvature perturbations, which is useful to consider the formation of PBHs.From the formula, we particularly demonstrate that PBHs with mass of order 10 20 g can form after hybrid inflation consistently with other cosmological observations if the number of waterfall fields is about 4 or 5 for the case of instantaneous reheating.Moreover, the domain wall problem is absent for the case with O(N ) symmetric waterfall fields.The spontaneous symmetry breaking of O(N ) symmetry after the waterfall phase transition results in the formation of texture (for N = 4) or global defects (for N ≥ 5).Gravitational waves (GWs) are generated via the dynamics of those structure as well as the second-order effect of large curvature perturbations.We demonstrate that the GW signals can be observed by future GW experiments, including LISA.
The rest of the paper is organised as follows.In Sec. 2, we first summarise the formalism for the stochastic equations for O(N ) symmetric fields and write down the Fokker-Planck and Langevin equations.These equations are solved analytically and numerically in Sec. 3, in hybrid inflation models with multiple waterfall fields.We provide a semi-analytic formula for the spectrum of curvature perturbations, where numerical coefficients are fitted by the results of numerical calculations.Furthermore, in Sec. 4 we use it to consider the formation of PBHs and the prediction of GW spectrum.Sec. 5 is devoted to discussion and conclusions.

Stochastic equations for symmetric fields
Let us first see the Fokker-Planck (FP) and Langevin equations that describe the stochastic dynamics of O(N ) symmetric fields during inflation.One will find that the system reduces to the one only for the radial mode in the slow-roll limit with an effective centrifugal force proportional to the number N of the fields due to the diffusion term.If the N fields correspond to the waterfall fields in the hybrid inflation, such a centrifugal force prevent them from staying around the origin of the potential and accordingly suppresses the curvature perturbation as we will see.
We employ the slow-roll approximation throughout this paper.The Langevin equation for N (canonical) inflaton fields ψ i (i = 1, 2, • • • , N ) (which we will identify as waterfall fields in the context of hybrid inflation model) is expressed as (see, e.g., Ref. [72]) where N is the e-folding number as the time variable, M Pl is the reduced Planck mass, V is the inflatons' potential, is its first derivative, and ξ i (N ) is an independent stochastic noise: 2) The corresponding FP equation for the inflatons' probability distribution This is equivalent to the following adjoint FP equation for the distribution P FPT (N | ψ) of the first passage time N from a certain field space point ψ to the end-of-inflation surface: Note that the forward e-folding number N is a deterministic time variable while the backward one N is stochastic as it depends on the realisation of the inflatons' stochastic dynamics.
Let us then suppose that the inflatons are symmetric and the potential (and hence all physical variables in the slow-roll limit) depends only on the radius ψ r = ψ 2 i .In this case, P FPT should solely be a function of N and ψ r .Therefore, the adjoint FP equation reduces to (refer also to Ref. [75]) (2.5) One finds a new term ∝ V N −1 ψr ∂ ψr from the Laplacian V ∂ 2 i .It can be understood as an effective single field system with a new centrifugal force dictated by the FP equation or the corresponding Langevin equation with a normalised noise We below see the effect of this centrifugal force in the hybrid inflation model.

Multi-waterfall hybrid inflation
In order for a general discussion, we consider the (1 + N )-field hybrid inflation with a phenomenological potential [65,66] characterised by the five dimensionful parameters Λ, M , ϕ c , µ 1 , and µ 2 .Here and hereafter, the summation in terms of i is implicit for ψ 2 i (≡ i ψ 2 i ).The inflaton potential is expanded around the critical point ϕ c and the higher-order terms than the quadratic one are neglected as we are mainly interested in the dynamics around the waterfall phase transition.Among the inflatons ϕ and ψ i (i = 1, 2, • • • , N ), the latter is often called waterfall fields as they terminate inflation by their tachyonic instability in the standard scenario.For a large ϕ, the waterfall fields stay around the origin and inflation occurs by the potential energy (≃ Λ 4 ).The Hubble parameter during inflation is given by H inf ≃ Λ 2 /( √ 3M Pl ).When the field ϕ reaches the critical point ϕ c , ψ i becomes tachyonic (called the waterfall phase) and starts to roll down toward the potential minimum ψ 2 i = M 2 .Depending on the parameters, the waterfall phase does not necessarily end quickly but can last as a second phase of inflation.The inflation ends when the slow-roll condition is violated.It is often controlled by the second slow-roll parameter along the waterfall direction In particular, the potential of interest is symmetric with respect to ψ i and therefore the model is effectively a two-field system of ϕ and ψ r = ψ 2 i with the centrifugal force disucssed around Eq. (2.7) in the slow-roll limit.We below study this two-field system.

Analytic estimations
The spectrum of curvature perturbations can be estimated by an analytic method proposed in Ref. [65] (refer also to Refs.[66,67]) for the hybrid inflation model.We extend their analysis to the case with multiple waterfall fields, including the leading-order effect of quadratic term for the inflaton potential.
We first estimate the distribution of the waterfall fields at the critical point, which determines the dynamics in the waterfall phase and then the final curvature perturbation.It follows the Langevin equations, with the independent noise The stochastic noise for the inflaton ϕ is usually negligible for the waterfall dynamics, whereas that for the waterfall fields is important around the waterfall phase transition.We denote the e-folding number at the time of the waterfall phase transition as N c .Neglecting the stochastic noise for ϕ, it is solved as for N ≈ N c .Eq. (3.2) derives the evolution equation for ⟨ψ 2 r ⟩ as where we assume ⟨ψ 2 r ⟩ ≪ M 2 and define With an approximation V ≃ Λ 4 , it can be solved as where we adopt the asymptotic condition x 0 e −t 2 dt is the error function.ψ r,c is defined by hence represents the amplitude of ⟨ψ 2 r ⟩ at the phase transition (N = N c ).This is enhanced by a factor of N for multiple waterfall fields.
The solution Eq. (3.7) approaches to 2ψ 2 r,c exp 8 2) ≡ N cl .The time scale N cl represents the time after which the classical dynamics (namely the effect of the first term in the right-hand side of Eq. (3.2)) begins to dominate.We will observe later that Π should be 10-40, so that N c = 7-28 for the case we are interested in.The stochastic effect is therefore important even well after the waterfall phase transition.As will be observed below, large curvature perturbations are generated during this regime, so that the numerical simulation for the FP equation is necessary to calculate the spectrum of curvature perturbations.This will be done in Sec.3.2.In this subsection, we nevertheless derive an analytic rough estimation on the curvature perturbation by omitting the stochastic effect and solving the classical equation of motion in the waterfall phase with the initial condition of ψ r = ψ r,c at N = N c , following Refs.[65][66][67].We will see that it indeed provides some information of the spectrum of the curvature perturbation.
Here we quote some results presented in Ref. [67] with some corrections by a factor of N including the leading-order correction from the quadratic potential for the inflaton.The e-folding number from the time of the waterfall phase transition N c to the end of inflation N end is calculated as where c (< 1) and χ 2 are given by c ≡ The δN approach gives the spectrum of curvature perturbations as where the backward e-folds N k is defined by the time when k = aH and ∆N 1 is the solution to the following equation: The peak amplitude reads where we neglect O(c/χ 2 ) terms in the second line.We will see in Sec. 4 that the parameter χ 2 is ∼ 10 almost constantly for the parameters of interest and also µ 2 is fixed to ≃ 10M Pl by the CMB observation.Therefore, the peak amplitude P (peak) R is uniquely determined by the corresponding peak wavelength (through N PT ) except for the suppression due to N .If N is equal to unity or two, the peak amplitude is too large for a desirable amount of PBHs with astrophysical scales [65,66].This challenge can be evaded when we consider multi-waterfall fields because of the factor of 1/N in Eq. (3.15).
As we commented above, the stochastic effect is not negligible for several e-foldings after the waterfall phase transition.Nevertheless, the above analytic form works reasonably well as we confirm with the detailed numerical calculations in the stochastic formalism in the next subsection.We will use the analytic formula (3.12) as a fitting function of the numerical result, regarding N PT and P (peak) R as fitting parameters rather than giving them by Eqs.(3.9) and (3.14).

Numerical implementation of the stochastic-δN formalism
In the full stochastic approach to inflation, the power spectrum of the curvature perturbation can be calculated in the so-called stochastic-δN formalism [70][71][72][73][74].The essential idea is that the probability distribution of the stochastic backward e-folds N is nothing but that of the curvature perturbation thanks to the δN formalism [76][77][78][79].The power spectrum is obtained by extracting the variance of N corresponding to the scale of interest.Specifically, it is given by the field-space integration (see Eq. (3.19) of Ref. [74]) While the backward probability P bw is formulated by a combination of the probability density function of ϕ (the solution of the FP equation) and that of N (the solution of the adjoint FP equation) in a rigorous way (see Ref. [73]), it is easy to sample the solution of the Langevin equation in a practical, numerical simulation.There, the backward probability is approximated as where S is the sampling number, ϕ i (N ) is the ith solution of the Langevin equation, and N i is the first passage time from the initial field value ϕ 0 for that solution.Its derivative can be further approximated by the finite difference: with a certain small step ∆N .Substituting these approximations back into the formula (3.16), the power spectrum is approximated as where Note that the average δN 2 ϕ 0 → ϕ ± i (k) is defined only for sample paths passing the point ϕ ± i (k), which is a non-trivial condition in a two or more than two fields case.We further try to approximate it by the average without a condition.Let us first define the total perturbation for paths passing ϕ ± i (k) by δN (ϕ 0 ) Mean squares of both sides read where we used the Markovian property of the Langevin equation for the last line.We then expect that the difference between δN 2 (ϕ 0 ) ϕ + i (k) and δN 2 (ϕ 0 ) ϕ − i (k) will be suppressed as O(∆N 2 ).This would be statistically justified because for typical points ϕ ± i (k), δN 2 (ϕ 0 ) ϕ ± i (k) asyptotes to a function only of ϕ 0 and becomes stationary against the change between ϕ + i (k) and ϕ − i .Therefore, one finds at the leading order in ∆N .
The variances on the right-hand side are found as a solution of the adjoint FP equation.In fact, recalling that the average ⟨f (N (ϕ))⟩ of an arbitrary function f (N (ϕ)) is defined by the probability density P FPT as the adjoint FP equation with the boundary condition P FPT (N | ϕ) = δ(N ) on the end-of-inflation surface ∂Ω leads to the following recursive partial differential equations: with the boundary condition ⟨N (ϕ)⟩ ϕ∈∂Ω = ⟨δN 2 (ϕ)⟩ ϕ∈∂Ω = 0.In this paper, we numerically solve these partial differential equations in the Jacobi method with use of the C ++ package, StocDeltaN [80] (see also the author's GitHub page [81]).

Results
We show the numerical results of the stochastic-δN formalism in this subsubsection.We fixed M and ϕ c by M = ϕ c / √ 2 = 10 16 GeV while µ 1 is determined by Π parameter through its definition (3.6).Λ and µ 2 are fixed by Eqs.(4.4) and (4.7) due to the CMB constraint.In the δN formalism, the end surface should be given by a uniform-density slice.We first solve the classical slow-roll equations of motion (Eq.(3.2) without the noise terms) from the initial condition ϕ = ϕ c and ψ = ψ r,c until η ψψ ≡ M 2 Pl V ψr ψr V reaches −2 (for sufficient violation of the slow-roll condition) for the first time and define the end-of-inflation energy density for the stochastc-δN scheme by the potential value at that time.
Fig. 1 shows the contours of ⟨N (ϕ, ψ r )⟩ and ⟨δN 2 (ϕ, ψ r )⟩ for Π 2 = 100 and N = 1 or 100 as solutions of the recursive partial differential equations (3.25) obtained by the StocDeltaN package.One finds that inflatons' sample trajectories represented by red lines pass the larger-ψ r region for N = 100 than N = 1 due to the noise-induced centrifugal force.The interesting feature is that though ⟨N ⟩ is not so sensitive to the number of waterfall fields N , the perturbation variance ⟨δN 2 ⟩ indeed decreases in the large N case.Therefore, one can expect that a large N can reduce P R , keeping its peak scale large enough.
Fig. 2 shows the resultant power spectra of curvature perturbations for several Π 2 and N .We adopt the sampling number S = 10 5 and the e-folds step ∆N = 0.5.Each error bar represents the propagated error in the finite difference equation (3.22) from the estimated standard error in the sampling averages 1 S S i=1 ⟨δN 2 (ϕ ± i (k))⟩.We also show the fitting formula (3.12) by the dotted lines with the corresponding colours.The fitting parameters N PT and P (peak) R are chosen by the numerically obtained peak wavenumber and peak amplitude.One finds that the fitting formula reasonably works well, particularly for a large Π 2 and a large N where the slow-roll approximation and the classical approximation (to neglect the noise terms) are better respectively.Table 1 shows the numerical results of the fitting parameters N PT and P (peak) R for Π 2 = 50, 100, 200, 500, 1000, and 2000 with N = 1, 3, 10, 100, and 1000.Those results are also shown in Fig. 3.The shaded region on the top-left panel represents the analytical result Eq. (3.14) multiplied by a factor of 1/3 for Π 2 ∈ (50, 2000).One finds that the analytic formula reasonably works again.Furthermore, the behaviour P (peak) R ∝ 1/N is confirmed by the numerical result.The thin dashed (solid) curves on the bottom panel represent the analytic results Eq. (3.9) in the leading ∼ Π (next-to-leading ∼ Π 2 ) order approximations with c = 0. Though the leading order approximation is not enough, the next-to-leading order approximation is well consistent with the numerical result.One can also see that N PT is almost independent of N and hence the degeneracy between N PT and P  chosen by the numerically obtained peak wavenumber and peak amplitude.

Parameters of interest for astrophysical PBH formation
In this section, we discuss which parameters should be chosen to generate the desired mass of PBHs consistently with other observed quantities.We then apply our result of the previous section to calculate the PBH mass function.We also show a GW spectrum predicted in our model.for given Π 2 and N .50,2000).The thin dashed (solid) curve on the bottom panel represents the analytic formula Eq. (3.9) in the leading ∼ Π (next-to-leading ∼ Π 2 ) order approximations with c = 0.

CMB constraints
Since perturbations at the CMB scale exits the horizon before the waterfall phase transition, they are generated by the quantum fluctuation of the inflaton ϕ.The amplitude of curvature perturbations is therefore given by the standard textbook formula where is a slow-roll parameter.This should be consistent with the observed one: where k * (= 0.05 Mpc −1 ) represents the wavenumber at the pivot scale [1].This implies or In addition, using Eq.(4.3) and Eq.(3.11), we obtain Thus we expect χ 2 ∼ 10.
The spectral index in our model is provided by We can consider µ 2 ≃ 10M Pl to explain the observed spectral index of n s = 0.9649±0.0042[1].

Constraint from self-ordering scalar fields
Our model has the global O(N ) symmetry, which is spontaneously broken to O(N − 1) at the waterfall phase transition.For the case with N ≤ 4, topologically non-trivial field configurations form after the spontaneous symmetry breaking (SSB).For higher N , no stable topological defect exist, but the Nambu-Goldstone (NG) modes are randomly distributed over the horizon scale and their gradient energy is proportional to the SSB scale squared.When those modes come into horizon at a late time, they tend to move in order to reduce the gradient energy [82].The dynamics of NG modes as well as topological defects can affect the CMB temperature anisotropy.Especially for the case with N = 2 and 3, where cosmic strings (for N = 2) and monopoles (for N = 3) form at the phase transition, Ref. [83] discussed a constraint on the SSB scale based on the Planck results such as M < 2.9 × 10 15 GeV for N = 2, M < 6.4 × 10 15 GeV for N = 3, ( We expect that the constraint for the case with large N is weaker by a factor of a few.Therefore we are interested in the case with M ≲ 10 16 GeV.Note that this constraint can be evaded if one introduces a tiny but nonzero explicit O(N ) breaking term in the Lagrangian because all NG modes are aligned by that term.The dynamics of self-ordering scalar fields also emits gravitational waves [84][85][86] (see also Ref. [87] in a different context).It has a flat spectrum in the frequency range of interest.Its amplitude is given by where Σ N is a numerical factor that approaches to unity for a large N limit.For example, Σ N ≃ 4.1 (1.7) for N = 4 (8) [88].This should be added to the gravitational wave signals from second-order effect from large curvature perturbations, which we will discuss shortly.

PBH formation
PBHs form when an overdense region with an O(1) density perturbation enters the horizon.
A sizable amount of PBHs can form when [8] 4 From the top-right panel of Fig. 3, the number of waterfall fields should be O(1-10).
The perturbation of comoving wave length λ enters the horizon when a(t)λ(t) = 1/H(t) ≡ 1/H HC .Assuming that this happens during the radiation dominated epoch, we obtain the corresponding e-folding number N PT at which the relevant perturbation exits the horizon: where we use t end ∼ 1/H inf and t RH ∼ 1/H RH .The PBH mass is given by the total energy enclosed within the Hubble horizon at the horizon crossing: where γ (= O(1)) is a numerical constant [8].We thus obtain where we assume γ = 0.2.From the bottom panel of Fig. 3, we are interested in Π ∼ 10-20.
We should take care of the tail of the spectrum for the curvature perturbations so that the spectrum around the CMB scale is not affected.Since the width of the spectrum ∆N 1 is of the same order with the peak frequency N PT , the latter one cannot be arbitrary large.For the case of instantaneous reheating (i.e., H RH = H inf ), we obtain Π

Summary of parameters of interest
The parameters we are interested in are roughly given by Λ ≃ 5.4 × 10 15 GeV Π 10 ) M ≲ 10 16 GeV, ( This is the number of waterfall fields that is requied to generate PBHs with a mass of interest.

Results
We show an example that provides PBHs with a mass of order 10 20 g, where we choose N = 5, Π 2 = 100, and µ 2 = 10M Pl with ϕ c / √ 2 = M = 10 16 GeV.Interporating the numerical results in Table 1, we obtain N PT ≃ 14 and P GeV.The peak amplitude and frequency are determined by the interpolation for numerical results of stochastic-δN formalism, whereas the shape of the spectrum is given by the analytic result Eq. (3.12).results in Table 1, whereas the shape of the spectrum is given by the analytic result Eq. (3.12) to show a smooth curve.The shaded region above the dashed curve on the top of the figure is excluded by the overproduction of PBHs (see Ref. [9] and references therein).The redshaded region in the left corner of the figure is excluded by the CMB observation [112][113][114][115].The blue and green shaded region in the upper left corner is excluded by constraint for µ and y distortions [116,117]. 5The constraint and future sensitivity curves for gravitational wave experiments are shown as the other dense and light-shaded regions, respectively (see Ref. [119] and references therein).
In the right panel of Fig. 4, we show the PBH mass function generated by the collapse of overdense regions.The peak mass for the PBH is about 5.3 × 10 19 g, which is within the window for the PBH DM.Note that, in order to keep the model parameters simple, we chose the threshold value δ c (∼ 0.3) so that the total PBH energy density is equal to the observed DM density rather than fixing δ c and fine-tuning the parameters.Hence this mass function should be understood as a schematic illustration.In fact, the threshold value with respect to the density contrast is known to be non-universal but depend on the profile of the overdensity and anyway one has to go beyond the Press-Schechter approach for a precise evaluation of the PBH abundance.The non-Gaussianity of the curvature perturbation will also affect the PBH abundance.We leave all these possible corrections for future works (see also the comments in footnote 4).The shaded regions are excluded by overproduction of PBHs (see, e.g., Ref. [9] and references therein for detail).
Fig. 5 shows the spectrum of stochastic gravitational waves predicted in our model.As we discussed around Eq. (4.9), gravitational waves are generated from the dynamics of re-alignment of NG modes when the modes enter into the horizon.In the figure, we plot the signal as the red dashed line for the case with M = 10 16 GeV and N = 5.We assume Σ N = 4. Stochastic gravitational waves are also generated by the second-order effect from large curvature perturbations.The solid blue curve represents the prediction for the case with the above parameter sets.The dense shaded regions are excluded by current experiments and the light-dense regions are future sensitivity curves for planned experiments.

Discussion and conclusions
We have analytically and numerically calculated the spectrum of curvature perturbations in the hybrid inflation model with O(N )-symmetric waterfall fields.The numerical simulation is based on the stochastic-δN algorithm proposed and developed in Refs.[70][71][72][73][74].We formulate the stochastic dynamics under the O(N ) symmetry and demonstrate that the resulting amplitude of curvature perturbations is reduced by a factor of order √ N for a fixed peak frequency.This factor is important to generate PBHs within a desirable mass range from a collapse of the overdense region at the horizon entry.It has been observed that the PBHs with a mass of order 10 20 g can be generated consistently with other cosmological constraints if the number of waterfall fields is about 5 for the case of instantaneous reheating.Such PBHs can explain all energy density of dark matter in the Universe.We have also calculated the spectrum of stochastic gravitational waves generated by the second-order effect from large curvature perturbations.The resulting GW signals can be observed in future GW experiments, including LISA.
The extended sector for the waterfall field is also motivated by avoiding the domain wall problem that arises in a model for a single waterfall field with Z 2 symmetry.The O(N ) symmetry is spontaneously broken by the vacuum expectation value of waterfall fields after the waterfall phase transition.This implies that stochastic NG modes are produced, which results in a continuous re-ordering of NG modes after the phase transition for the case of N ≥ 5 [82].The cosmological effect of those dynamics is qualitatively similar to the case with N = 2, 3, and 4, where topologically non-trivial field configurations form after the spontaneous symmetry breaking.Those dynamics particularly affect the CMB spectrum and give an upper bound on the symmetry-breaking scale.Moreover, it generates GWs with a flat spectrum [84][85][86], which would be observed by future GW experiments as well as the pulsar-timing array experiments.This is an interesting smoking-gun signal of the present model.
One may consider that any global symmetries should not be exact to ensure consistency with quantum gravity.The NG boson then has a small mass from an explicit symmetrybreaking term, and the stochastic modes are ordered by the mass term at a late time.In this case, the constraints from the re-ordering of NG modes may be evaded and the GW spectrum from the dynamics of self-ordering NG modes is suppressed at a small frequency.
Although we have specifically considered the case with global O(N ) symmetry for waterfall fields, the qualitative results, such as the scaling of the amplitude of curvature perturbations as ∝ N −1 , are expected to hold for other symmetries such as SU(N ).Moreover, one can consider the gauged O(N ) or other symmetries instead of the global ones.The stochastic effect on the NG modes as well as the radial mode should be the same in a covariant gauge at least in the regime where the mass of the gauge boson is much smaller than the Hubble parameter during inflation (see, e.g., Refs.[120,121]).Since the vacuum expectation value of the waterfall fields is comparable to the Hubble parameter at the waterfall phase transition, our calculation can be justified at least for the case with a sufficiently small gauge coupling constant.The gauge field becomes heavier after the waterfall phase transition than the Hubble parameter.The stochastic structure of NG modes in the covariant gauge can be understood as the massive gauge fields in a unitary gauge.The heavy gauge fields are then supposed to decay into light fields.Therefore the constraint from CMB anisotropy and prediction of GW signals from NG modes are absent for N ≥ 4 in this case.Still, the GW signals from the second-order effect of large curvature perturbations are the prediction of the model.
)weighted by the backward probability distribution P bw (ϕ | N ), the probability such that the inflatons take certain field values ϕ = {ϕ, ψ r } at the time N e-folds before the end of inflation.Ω denotes the field-space region where the slow-roll inflation can be realised.ϕ 0 is certain initial field values of the inflatons and δN (ϕ 0 → ϕ * ) is a fluctuation in the e-folds N (ϕ 0 → ϕ * ) from ϕ 0 to ϕ * , i.e., δN (ϕ 0 → ϕ * ) ≡ N (ϕ 0 → ϕ * ) + ⟨N (ϕ * )⟩ − ⟨N (ϕ 0 )⟩, where N (ϕ) is a first passage time from ϕ to the end of inflation.Note that this formula requires the leading-order slow-roll approximation, N k ≃ − ln(k/k f ), where k f is the comoving Hubble scale aH at the end of inflation.
by introducing N .

Figure 3 .
Figure 3.The numerical result for the fitting parameters P (peak) R (top) and N PT (bottom) as functions of Π and N .On the top-left panel, the shaded region represents the analytical result Eq. (3.14) multiplied by a factor of 1/3 for Π 2 ∈ (50, 2000).The thin dashed (solid) curve on the bottom panel represents the analytic formula Eq. (3.9) in the leading ∼ Π (next-to-leading ∼ Π 2 ) order approximations with c = 0.

2 ≲ 250 for ϕ c / √ 2 =
M = 10 16 GeV, which leads to N PT ≲ 20 and M PBH ≲ 2 × 10 25 g.For the case of the lowest possible reheating temperature of H RH ∼ 1 MeV 2 /M Pl , we obtain Π 2 ≲ 200 for ϕ c / √ 2 = M = 10 16 GeV, which leads to N PT ≲ 18 and M PBH ≲ 1 × 10 27 g.The upper bounds on Π or N PT do not change much for different values of ϕ c and M (with Π fixed).Because of the relations of Eq. (4.13) and Π ≡ M √ µ 1 ϕ c /M 2 Pl , a larger PBH mass can be generated for a smaller M and ϕ c under the same upper bound on N PT .

Figure 4 .
Figure 4. Spectra of curvature perturbations (left panel) and PBH mass function (right panel) for the case with N = 5, Π 2 = 100, µ 2 = 10M Pl , and ϕ c / √ 2 = M = 10 16GeV.The peak amplitude and frequency are determined by the interpolation for numerical results of stochastic-δN formalism, whereas the shape of the spectrum is given by the analytic result Eq. (3.12).

Figure 5 .
Figure5.Spectra of stochastic gravitational waves generated by the dynamics of NG modes (red dashed line) and the second-order effect from large curvature perturbations (solid blue curve).

Table 1 .
Numerical results of the fitting parameters N PT and P PT given by Eq.(4.13).The first two conditions, Eqs.(4.14) and (4.15) are very rough estimation and are shown for illustrative purposes.More accurate conditions should be read by interpolating the numerical results in Table1.We added a sub-Planckian condition for the critical point: ϕ c ≲ M Pl .The parameters we have are µ 1 , µ 2 , M , ϕ c , Λ, and N .Four of them, µ 1 , µ 2 , Λ, and N , are (approximately) determined by the first four conditions.We have two free parameters, M and ϕ c , which have the upper bounds.