Primordial black holes and associated gravitational waves in axion monodromy inflation
Shu-Lin Cheng Wolung Lee Kin-Wang Ng
astro-ph.CO hep-ph hep-th









^1Department of Physics, National Taiwan Normal University,
Taipei 11677, Taiwan

^2Institute of Physics, Academia Sinica, Taipei 11529, Taiwan

^3Institute of Astronomy and Astrophysics, Academia Sinica, Taipei 11529, Taiwan




In the axion monodromy inflation, the inflation is driven by the axion with super-Planckian field values in a monomial potential with superimposed sinusoidal modulations. The coupling of the axion to massless gauge fields can induce copious particle production during inflation, resulting in large non-Gaussian curvature perturbation that leads to the formation of primordial black holes. In this paper, we explore the parameter space in the axion monodromy inflation model that favors the formation of primordial black holes with masses ranging from 10^8 grams to 20 solar masses. We also study the associated gravitational waves and their detection in pulsar timing arrays and interferometry experiments.


98.80.Cq, 97.60.Lf, 04.30.-w
Primordial black holes and associated gravitational waves in axion monodromy inflation
    Kin-Wang Ng^2,3
    December 30, 2023
======================================================================================





§ INTRODUCTION


The inflation scenario is generally accepted for explaining the observed spatial flatness and homogeneity of the Universe.
A simple model of the scenario such as the slow-roll inflation driven by a flat inflaton potential predicts quasi de Sitter
vacuum fluctuations during the inflation which could give rise to Gaussian and
nearly scale-free metric perturbations containing both matter density fluctuations (scalar modes)
and gravitational waves (GWs or tensor modes) <cit.>. 

The simplest model that provides the flat potential is perhaps the large-field inflation with a monomial potential. The drawback is that the inflaton needs 
super-Planckian field values to fulfill the slow-roll conditions. However, these large field values can be realized in the string theory and hence are utilized in the axion monodromy inflation <cit.>. The inflaton potential in a single-field axion monodromy inflation has a monomial form with superimposed sinusoidal modulations whose size is model-dependent, given by <cit.>


    V(φ)=V_0+μ^4-pφ^p +Λ(φ)^4 cos[φ/f(φ)+γ_0],

where V_0, μ, and γ_0 are constants, and p=3,2,4/3,1,2/3. The modulation contains the energy scale Λ and the axion decay constant f, both of which are in general functions of φ. The cosmological phenomena derived from this kind of inflation models, such as the tensor-to-scalar ratio, and the resonantly enhanced modulations of the scalar power spectrum and bispectrum due to the sinusoidal modulations of the potential, have been studied <cit.>.

The axion is expected to be coupled to some gauge field via a pseudoscalar-vector coupling,


    ℒ_int = - α/4 f φ F̃^μν F_μν,

where α is a dimensionless coupling constant, F_μν = ∂_μ A_ν - ∂_ν A_μ is the vector field strength tensor, and F̃^μν = ϵ^μναβ F_αβ/2 is its dual. 
Recently, there has been a lot of studies on the phenomenological effects of this coupling in the axion inflation model. The coupling leads to the particle production with a rate proportional to the rolling speed of axions, which induces a new source of metric perturbations while backacting to the axion dynamics. 
This results in very interesting effects such as the generation of non-Gaussian and non-scale-invariant scalar power spectrum <cit.>, 
the production of GWs that could be detected at the ground-based gravity-wave interferometer <cit.>,
and the formation of primordial black holes (PBHs) with masses ≲ 10^8g near the end of inflation <cit.>.
More recently, the authors in Refs. <cit.> have made the axion potential steeper locally, or involved non-minimal coupling to gravity <cit.> to boost the particle production rate at certain wavenumbers in order to seed PBHs with much higher masses. Furthermore, the strength of the GWs associated with the formation of these PBH seeds may reach the sensitivity of future pulsar timing arrays and interferometry experiments.

The cosmic inflation is the most efficient way to seed the formation of PBHs. There have been many inflationary models that produce PBHs with various masses and associated GWs. In a single-field slow-rolling inflation model, matter density perturbations are generally well below the threshold to form PBHs, though they can be formed at rare density peaks. Modifications of the inflation potential to achieve blue-tilted matter power spectra or running spectral indices may lead to large enough density perturbations at the end of inflation <cit.>. However, the resulting PBH masses in most of these models are many orders of magnitude below M_⊙ <cit.>. To boost the PBH mass into the astrophysical and even the cosmological mass scales, several scenarios involving multi-field inflations have been proposed, such as the hybrid inflation <cit.>, the double inflation <cit.>, and the curvaton models <cit.>, in which small-scale density perturbations can be inflated to a scale ranging from the size of a stellar-mass PBH to that of a supermassive PBH. 

In this paper, we study the formation of PBHs in the axion monodromy inflation governed by the potential (<ref>) and the interaction (<ref>). We will not consider a specific axion model; instead, we treat the model parameters as free as possible. Then, we will probe the parameter space favorable to the formation of PBHs. At the end, we will discuss briefly the physically motivated model with the favorable model parameter. The paper is organized as follows. In the next section, we lay out our base axion inflation model. In Sec. <ref>, we calculate the curvature perturbation induced by particle productions and then consider the formation of PBHs. In Sec. <ref>, GWs sourced by the particle production are computed. Section <ref> is our conclusion.



§ THE MODEL


We consider the axion monodromy inflation in which the axion couples to a U(1) gauge field via the interaction (<ref>).
The action is given by


    𝒮 =  ∫ d^4 x √(-g)[ M_p^2/2  R  -1/2∂_μφ∂^μφ- V(φ) - 1/4F^μνF_μν - 1/√(-g)α/4 fφ F̃^μν  F_μν],

where R is the curvature scalar and M_p=2.435× 10^18 GeV is the reduced Planck mass. For the inflaton potential (<ref>), we choose a linear or concave potential (p=1 or p=2/3) because convex potentials (p>1) are disfavored by cosmic microwave background (CMB) experiments <cit.>, assuming that the ground state is at φ=0, and taking the approximation that Λ and f are constants. Hence, we have


    V(φ)= M_p^pμ^4-p c^-p[√(1+(cφ/M_p)^2p) - 1] +Λ^4 [1+cos(φ/f+π)],

where we have introduced a numerical factor c to adjust the smoothness of the potential at φ=0.
The second term in the above equation is a standard axion potential that has been widely used in axion inflation models as mentioned in the Introduction. The drifts in Λ and f during the course of inflation occur in a broad range of axion monodromy scenarios in string theory <cit.>. We have neglected the drifts and assumed constant Λ and f. This is a good approximation in the context of the production of PBHs as long as the drifts ΔΛ,Δ f ≪ M_p.  It would be interesting to include any possible drifts in the future work.

Here we assume a spatially flat Friedmann-Robertson-Walker metric:


    ds^2=-g_μν dx^μ dx^ν= a^2(η) (dη^2- d x⃗^2),

where a(η) is the cosmic scale factor and η is the conformal time related to the cosmic time by dt=a(η)dη.
The Hubble parameter is defined by H ≡ (da/dt) / a or H≡ (da/dη) / a.

From the action (<ref>), we can write down the Friedmann equation, the equation of motion for the inflaton, and the Maxwell equations, respectively:


    H^2 = 1/3 M_p^2[ 1/2(∂φ/∂η)^2 + 1/2( ∇⃗φ)^2   + a^2   V(φ) + a^2/2( E⃗^2 + B⃗^2 ) ],
          ∂ ^2φ/∂η^2 + 2  H∂φ/∂η - ∇⃗^2 φ  + a^2  d V/d φ  = a^2 α/fE⃗·B⃗,
          ∂^2 A⃗/∂η^2 - ∇⃗^2 A⃗
    + ∇⃗ (∇⃗·A⃗) = α/f∂φ/∂η∇⃗×A⃗ - α/f∇⃗φ×∂A⃗/∂η,
          ∂/∂η(∇⃗·A⃗) = α/f∇⃗φ·(∇⃗×A⃗),

where for the Maxwell equations we have chosen the temporal gauge, i.e. A_μ = (0, A⃗), and we have introduced the physical “electric” and “magnetic” fields,


    E⃗ = -1 a^2∂A⃗/∂η,  B⃗ = 1 a^2∇⃗×A⃗.


In Ref. <cit.>, we have calculated the production of gauge quanta by the rolling inflaton via the interaction during a slow-roll inflation, taking into account self-consistently the backreaction of the gauge quanta production on inflation.  To calculate the production of gauge quanta, we separate the inflaton into a mean field and its fluctuations:


    φ = ϕ (η) + δφ (η,x⃗).

Under the linear approximation, we decompose A⃗(η,x⃗) into two circularly polarized Fourier modes, A_± (η,k⃗), which satisfy
the equation of motion as


    [ d^2/dη^2 + k^2 ∓ 2aH kξ] A_±(η,k) = 0,   ξ≡α/2 f Hdϕ/dt .

This implies that either one of the two modes, when satisfying the condition k/(aH)< 2|ξ| for the spinoidal instability, grows exponentially fast.
The energy density and the interaction term of the produced gauge quanta are given by the vacuum expectation values of the electric and magnetic fields, respectively,


    1/2⟨E⃗^2+B⃗^2 ⟩   =   ∫dk k^2/4 π^2 a^4∑_λ=±( |dA_λ/dη|^2 + k^2 | A_λ|^2 ), 
    ⟨E⃗·B⃗⟩   =    -∫dk k^3/4 π^2 a^4d/d η(| A_+ |^2-| A_- |^2).

As a consequence, the production of gauge quanta gives rise to a backreaction on the background, whose evolution is then governed by


    d^2ϕ/dt^2 + 3 H dϕ/dt +dV/dϕ = α/f⟨E⃗·B⃗⟩, 
    
           3 H^2 = 1/M_p^2[ 1/2(dϕ/dt)^2 + V(ϕ) + 1/2⟨E⃗^2 + B⃗^2 ⟩].




§ CURVATURE PERTURBATION AND PRIMORDIAL BLACK HOLES



It is well known that the vacuum quantum fluctuation of an inflaton field during the inflation gives rise to the curvature perturbation whose power spectrum is governed by <cit.>


    Δ_ζ  vac^2(k)≡⟨ζ(x)^2⟩=1/4π^2H^4/(dϕ/dt)^2.

In this work, fluctuations of the gauge quanta production during the inflation would lead to a new source for inflaton fluctuations. 
The induced inflaton perturbation satisfies an inhomogeneous equation:


    [ ∂^2/∂ t^2 +3 β H ∂/∂ t -∇⃗^2/a^2.    +   . d^2V/dϕ^2] δφ(t,x⃗) 
       =   α f(  E⃗·B⃗ - ⟨E⃗·B⃗⟩),

where the frictional term can be derived as <cit.>


    β≡ 1-2πξα/f⟨E⃗·B⃗⟩/3H(dϕ/dt).

The particular solution to this equation can be well approximated by <cit.>


    δφ= α3β f H^2(  E⃗·B⃗ - ⟨E⃗·B⃗⟩),

which contributes to the power spectrum of the curvature perturbation an amount given by


    Δ_ζ^2(k)=H^2 ⟨δφ^2⟩/(dϕ/dt)^2
    =[α⟨E⃗·B⃗⟩/3β f H(dϕ/dt)]^2.


PBHs can be arisen from the collapse of over-dense regions originating in the curvature perturbation spawned by the inflaton as they re-entering the horizon in the subsequent expanding universe. If a PBH seed is created during the inflation, the mass of the PBH can be estimated as follows. 
The energy contained within the comoving seed volume that leaves the horizon N e-foldings before the end of inflation is given by


    4 3π H^-3 e^3Nρ_e    with   ρ_e=3H_e^2 M_p^2,

where H=H(N) is the Hubble parameter at N e-foldings before the end of inflation and H_e at the end of inflation. Let a_0=1 be the scale factor at the onset of the inflation, which will eventually span N_0 e-foldings over the entire inflationary course. After the inflation has terminated, the universe reheats and is becoming radiation-dominated. Whenever the conformal time η>η_e, the radiation dominant universe would expand to a size of 


    a(η)=H_e e^2N_0 (η-2η_e),   where  η_e=-(H_e e^N_0)^-1.
 
The comoving volume re-enters the horizon when its scale k=H e^N_0-N satisfies the condition kη∼ 1, i.e., 
when a=e^N_0+N (H_e/H) or the temperature of the thermal bath is red-shifted by a factor of e^N (H_e/H) <cit.>. 
Therefore, the mass of the PBH presumably formed at this time is 


    M_ BH=4 π M_p^2 H_e/H^2 e^2N=2.74× 10^-38 e^2N(M_p H_e/H^2) M_⊙.


For the purposes of forming PBHs in certain mass ranges associated with GWs of observational interest, we have worked out three specific cases for the potential V(ϕ). Figure 1 shows the first case. In this figure and hereafter, we rescale all dynamical variables in terms of the reduced Planck mass. With this inflaton potential, the background solutions for ϕ and ξ in Eqs. (<ref>) and (<ref>) are plotted in Fig. 2, for the given initial position and speed of the inflaton, ϕ_0 and (dϕ/dt)_0, and the coupling constant α.
The number of e-foldings after the onset of the inflation is defined by ∫_0^t H(t') dt'. Note that a_0=1 and we find that N_0≃ 61. 
Also, in Fig. 3, we have evaluated the scalar spectral index n_s = 1 - 2ϵ_1 - ϵ_2, the tensor-to-scalar ratio r = 16ϵ_1, and the running of the scalar spectral index dn_s/dln k, where the slow-roll parameters ϵ_1 and ϵ_2 are derived from


    ϵ_1= - 1/H^2dH/dt,  ϵ_2= 1/Hϵ_1dϵ_1/dt.

Figure 4 presents the total power spectrum of the curvature perturbation, which is the sum of the vacuum contribution in Eq. (<ref>) and the induced power in Eq. (<ref>). 
Note that the amplitude of the total power spectrum of the curvature perturbation at the largest scales as well as the levels of the scalar spectral index, the tensor-to-scalar ratio, and the scalar spectral index running for cosmologically interesting scales that correspond to about the first 7 e-foldings are all consistent with the Planck measurements: 
Δ_ζ^2≃ 2× 10^-9, n_s≃ 0.97, r<0.1, and |dn_s/dln k|<0.013 <cit.>. 

The perturbation ζ in Eq. (<ref>) is obviously non-Gaussian. In Ref. <cit.>, since A⃗ is a Gaussian vector field, it was argued that 
ζ can be written as ζ=g^2-⟨ g^2⟩, with g a Gaussian field. Let P(g) be its Gaussian distribution function. Then, the probability distribution function of ζ is given by P(ζ)dζ=P(g)dg, with g_c^2=ζ_c+⟨ g^2⟩, where ζ_c∼ 1 denotes the critical value leading to the PBH formation. This value of g_c for a Gaussian field determines the energy fraction that can collapse to form horizon-sized PBHs with mass M at the re-entry:


    β(M)=∫_ζ_c^∞ P(ζ)dζ = ∫_g_c^∞ P(g)dg.
 
The astrophysical and cosmological constraints on the PBH energy fraction β(M) can then be translated into an upper bound on the power spectrum of the curvature perturbation <cit.>, as shown by the short-dashed line in Fig. 4.













In Fig. 2, the first 10 e-foldings of the inflation are in standard slow-roll stage. When the inflaton enters the steeper slope, it speeds up and hence ξ reaches a maximum value of 5.8 at N≃ 40. This produces a huge amount of gauge quanta and induces a high peak of inflaton fluctuations at N≃ 38.5. Consequently, the strong backreaction due to the gauge field production almost stops the inflaton motion at ϕ≃ -2.6 for about half of the inflationary duration. After then, the inflaton speeds up again and the inflation ends up with a rather complicated dissipation-fluctuation processes. This is a typical example of the so-called trapped inflation in which the production of gauge quanta can sustain a nearly steady thermal bath during the inflationary epoch and exhaust the vacuum energy to terminate the inflation gracefully without undergoing a large-scale preheating or perturbative reheating <cit.>.

In Ref. <cit.>, the authors have shown that in the in-in formalism the validity of the perturbative treatment of axion and gauge modes restricts ξ≲ 5. However, this result does not apply to our present work. Instead of using the perturbation approach, we have taken only the linear approximation for the axion and gauge modes in the equations of motion to carry out all calculations including the strong backreaction. The validity of this linear approximation was recently discussed in Ref. <cit.> for the same gauge-field axion inflation, in which the authors have proposed a consistency condition for the linearity. This condition was later analyzed in a detailed numerical calculation by including the backreaction self-consistently, which shows that the linearity condition can be maintained for strong couplings as large as ξ≲ 10 <cit.>. In Fig. 2, our main concern is the peak at N≃ 40, where ξ≃ 5.8 is well within the linearity condition. Near the end of the inflation ξ gets larger values that drive the scalar power spectrum to high spikes at which PBHs of much smaller masses are likely to be copiously formed. However, these high-density spikes should be damped to a lower level because of the local strong gravity. This gravity effect has been ignored in the present work; indeed, it can be properly taken into account by including a higher-order gravitational term in the inflaton perturbation Eq. (<ref>) <cit.>. From Eq. (<ref>), a black hole with M_ BH≃ 19.7 M_⊙ can be produced by taking H_e=1.48× 10^-6 and H=2.38× 10^-6 when N=38.5. Hence, the peak of the scalar power spectrum in Fig. 4 will seed the formation of higher solar-mass PBHs,
with the black hole fraction β(M)∼ 10^-9, which gives the fraction of dark matter in PBHs f(M)∼ 0.1 <cit.>. In Appendix, we will present two other inflaton potentials that can produce PBHs with smaller masses.



§ ASSOCIATED GRAVITATIONAL WAVES



In addition to de Sitter vacuum fluctuations, GWs can be directly sourced by the gauge field production <cit.>. The GW equation reads 


    [∂^2/∂η^2 + 2 ada/dη∂/∂η -∇⃗^2] h_ij
    = 2a^2/M_p^2( -E_i E_j - B_i B_j )^TT,

where TT denotes the transverse and traceless projection of spatial components of the energy-momentum tensor of the gauge field. 

While the curvature perturbation is sufficiently large for seeding the formation of PBHs, GWs simultaneously produced at the horizon crossing as the second-order effect in the metric perturbation theory cannot be neglected <cit.>. The second-order effect is mainly contributed by the transverse-traceless part of a source term S_ij(δφ) involving quadratic terms of the curvature perturbation that should appear on the right-hand side of the GW Eq. (<ref>). In the present consideration, we will show that the amount of GWs sourced by S_ij(δφ) is subdominant to that by the gauge field production. From naive power counting, the GW amplitude induced by the gauge field can be estimated as h_A∼ A^2∼ζ, whereas that induced by S_ij(δφ) is h_δφ∼δφ^2∼ζ^2. As long as Δ_ζ^2∼ζ^2 is beneath the PBH bound, ζ≪ 1 and thus we have h_δφ≪ h_A. In fact, this simple estimation is supported by detailed calculations <cit.>. In Ref. <cit.>, it was shown that the energy densities of inflaton perturbation ρ_δφ and of gauge quanta ρ_A during inflation scale as ρ_δφ/ρ_A ∼ζ≪ 1. This implies that the energy-momentum tensor that sources the generation of GWs mainly comes from the gauge quanta rather than the inflaton perturbation. More recently, the authors in Ref. <cit.> have confirmed the sub-dominance of the second-order effect by explicitly computing the GW power spectra induced by the gauge field production as well as by the second-order curvature perturbation.

Rather than numerically solving Eq. (<ref>) in conjunction with the gauge mode equation (<ref>), we use the approximate analytic gauge mode solutions to estimate the present relative GW energy density per logarithmic k interval, given by <cit.>


    Ω_ GW h^2≃ 3.5× 10^-7H^2/M_p^2(1+4.3× 10^-7H^2/M_p^2e^4πξ/ξ^6).

To evaluate Ω_ GW h^2, we treat ξ as a function of N given by the numerical results plotted as ξ(N) in Figs. 2, 7, and 11. The Hubble scale H=H(N) is calculated by Eq. (<ref>). This spectral energy density is at k=H e^N_0-N with H_0 corresponding to the present horizon of size 0.002  Mpc^-1.  In Fig. 5, we plot Ω_ GW h^2 against the frequency f=k/(2π)=3× 10^-18 (H/H_0) e^N_0-N Hz, where N_0 is the total number of e-foldings for the inflationary duration with inflaton potentials plotted in Figs. 1, 6, and 10, respectively. This gives rise to three prominent peaks, which from left to right are associated with the production of PBHs with masses of 19.7 M_⊙, 2.4× 10^-13 M_⊙, and  1.2× 10^8  g, respectively. In the figure, we list the current upper limits on GW background inferred from pulsar timing array data <cit.> and aLIGO O1 data <cit.>. Also shown are the projected sensitivities of on-going and future GW experiments such as aLIGO O3, O5 <cit.>, LISA <cit.>, and SKA radio telescope <cit.>.






§ CONCLUSION



The axion monodromy inflation is a well-motivated inflationary scenario in string theory. The inflation is driven by the axion with super-Planckian field values in a monomial potential with superimposed sinusoidal modulations. While the monomial potential provides a slow-roll inflation consistent with CMB data and large-scale-structure surveys, the modulations may give rise to interesting observational signatures. 

Flauger et al. in Ref. <cit.> have found that the periodic modulations can drive resonant enhancements of inflaton perturbations, with characteristic scale-dependent modulated amplitude, giving rise to oscillations in the CMB anisotropy power spectrum. For an axion decay constant f≪ M_p, they have calculated the modulated curvature power spectrum, and determined the limits that CMB data places on the amplitude and frequency of modulations. In some favorable axion monodromy models, resonant contributions to the CMB power spectrum and bispectrum could be detected in near-future CMB experiments.

In this paper, we have studied the consequences from the periodic modulations of the axion potential with f∼ M_p, when the axion is coupled to a massless gauge field with a coupling constant α∼ 20. The modulations with Planckian frequencies have little effect on the slow-roll inflation at CMB scales. After the inflaton leaves the slow-roll regime and slides down a steeper slope of the modulated potential, a huge amount of gauge quanta is produced. As a consequence, the backreaction of the particle production to the inflaton motion induces large inflaton fluctuations, creating a peak in the curvature power spectrum that seeds the formation of primordial black holes. We have given permissible model parameters that produce PBHs with masses ranging from 10^8 grams to 20M_⊙. Interestingly, the amplitudes and frequencies of gravitational waves sourced by the stress-energy tensor of the generated gauge quanta could lie within the sensitivities of on-going and future gravity-wave detectors such as aLIGO/VIRGO, LISA, and pulsar timing arrays.

The axion monodromy provides a sound theoretical framework for a large-field inflation model to work with. In addition to the monomial potential that realizes the slow-roll regime, inherent nonperturbative effects generate small sinusoidal modulations of the potential with model-dependent axion decay constant f and energy scale Λ. In the present consideration, we have introduced an axion-photon interaction with a coupling constant α. This interaction strength can be rewritten as α/f=α_e/(2π f_a), where α_e≃ 1/137 is the fine structure constant and the rescaled axion decay constant f_a∼ 10^-4 f. It is expected that the axion decay constant f is in the order of the grand unification scale, which is generally smaller than the Planck scale. However, there are some examples in string theory that allow f to be near the Planck scale <cit.>. The preferred values for the model parameters that can produce primordial black holes and associated gravitational waves of astrophysical interest are f∼ M_p, Λ∼ 10^-3 M_p, and f_a∼ 10^-4 M_p. These disparate energy scales may be constructed, for example, in the clockwork mechanism <cit.>. Overall, the axion monodromy inflation has very rich astrophysical and cosmological implications that may be tested in on-going and future CMB and gravity-wave experiments.


The authors would like to thank K. Choi, M. Peloso, and E. Silverstein for useful conservations.
This work was supported in part by the Ministry of Science and Technology, Taiwan, ROC under Grants No. MOST104-2112-M-001-039-MY3 (K.W.N.) and No. MOST104-2112-M-003-013 (W.L.).


*


§ TWO MORE CASES WITH SMALLER PBH MASSES


Figure 6 and Fig. 10 are the second and third potentials, respectively.  Similar to the first case, we have calculated the background and the perturbation in these two potentials. The results are summarized in Figs. 7-9 and Figs. 11-13. Figure 8 and Fig. 12 show n_s, r, and dn_s/dln k, all of which satisfy the Planck constraints on CMB scales. In Fig. 7, the ξ bump at N≃ 25 induces a Δ_ζ peak at N≃ 23 with a peak value saturating the PBH bound as seen in Fig. 9, which corresponds to M_ BH≃ 2.4× 10^-13 M_⊙, by taking
H_e=3.99× 10^-7 and H=2.15× 10^-6. 
This contributes to the fraction of dark matter in PBHs of amount f(M)∼ 0.06 <cit.>.
In Fig. 11, the ξ peak at N≃ 9 produces a Δ_ζ peak just outside the PBH bound as shown in Fig. 13. This density peak at N≃ 8 seeds PBHs with M_ BH≃ 1.2× 10^8  g when H_e=8.54× 10^-7 and H=1.78× 10^-6. These small PBHs would have evaporated through the
emission of Hawking radiation.




















olive
For reviews see:  K. A. Olive, Phys. Rep. 190, 307 (1990);
D. H. Lyth and A. Riotto, Phys. Rep. 314, 1 (1999).

silverstein
For a review, see: E. Silverstein, in Proceedings of the 2015 Theoretical Advanced Study Institute in Elementary Particle Physics,
(World Scientific, Singapore, 2017).

mcallister
L. McAllister, E. Silverstein, A. Westphal, and T. Wrase, J. High Energy Phys. 09 (2014) 123;
F. Marchesano, G. Shiu, and A. M. Uranga, J. High Energy Phys. 09 (2014) 184.

flauger
R. Flauger, L. McAllister, E. Silverstein, and A. Westphal, J. Cosmol. Astropart. Phys. 10 (2017) 055.

cosmology
E. Silverstein and A. Westphal, Phys. Rev. D 78, 106003 (2008);
L. McAllister, E. Silverstein, and A. Westphal, Phys. Rev. D 82, 046003 (2010);
R. Flauger, L. McAllister, E. Pajer, A. Westphal, and G. Xu, J. Cosmol. Astropart. Phys. 06 (2010) 009;
R. Flauger and E. Pajer, J. Cosmol. Astropart. Phys. 01 (2011) 017;
N. Kaloper, A. Lawrence, and L. Sorbo, J. Cosmol. Astropart. Phys. 03 (2011) 023;
T. Higaki, T. Kobayashi, O. Seto, and Y. Yamaguchi, J. Cosmol. Astropart. Phys. 10 (2014) 025;
E. Palti and T. Weigand, J. High Energy Phys. 04 (2014) 155;
Y. Wan, S. Li, M. Li, T. Qiu, Y. Cai, and X. Zhang, Phys. Rev. D 90, 023537 (2014);
Q. E. Minor and M. Kaplinghat, Phys. Rev. D 91, 063504 (2015);
L. C. Price, Phys. Rev. D 92, 103507 (2015);
K. Choi and H. Kim, Phys. Lett. B 759, 520 (2016).

barnaby
N. Barnaby and M. Peloso, Phys. Rev. Lett. 106, 181301 (2011);
N. Barnaby, R. Namba, and M. Peloso, J. Cosmol. Astropart. Phys. 04 (2011) 009.

meer
P. D. Meerburg and E. Pajer, J. Cosmol. Astropart. Phys. 02 (2013) 017.

barnaby2
N. Barnaby, E. Pajer, and M. Peloso, Phys. Rev. D 85, 023525 (2012).

cook
J. L. Cook and L. Sorbo, Phys. Rev. D 85, 023534 (2012);
86, 069901(E) (2012);
M. M. Anber and L. Sorbo, Phys. Rev. D 85, 123537 (2012);
N. Barnaby, J. Moxon, R. Namba, M. Peloso, G. Shiu, and P. Zhou, Phys. Rev. D 86, 103508 (2012);
D. Jiménez, K. Kamada, K. Schmitz, and X.-J. Xu,  J. Cosmol. Astropart. Phys. 12 (2017) 011.

linde
A. Linde, S. Mooij, and E. Pajer, Phys. Rev. D 87, 103506 (2013).

lin
C.-M. Lin and K.-W. Ng, Phys. Lett. B 718, 1181 (2013);
E. Bugaev and P. Klimai, Phys. Rev. D 90, 103501 (2014);
E. Erfani, J. Cosmol. Astropart. Phys. 04 (2016) 020.

cheng
S.-L. Cheng, W. Lee, and K.-W. Ng, Phys. Rev. D 93, 063510 (2016).

cheng2
S.-L. Cheng, W. Lee, and K.-W. Ng, J. High Energy Phys. 02 (2017) 008;
J. García-Bellido, M. Peloso, and C. Unal, J. Cosmol. Astropart. Phys. 12 (2016) 031.

juan
J. García-Bellido, M. Peloso, and C. Unal, J. Cosmol. Astropart. Phys. 09 (2017) 013.

domcke
V. Domcke, F. Muia, M. Pieroni, and L. T. Witkowski, J. Cosmol. Astropart. Phys. 07 (2017) 048.

drees
See, for example, M. Drees and E. Erfani,
J. Cosmol. Astropart. Phys. 01 (2012) 035.

hybrid
See, for examples, S. Clesse and J. García-Bellido, Phys. Rev. D 92, 023524 (2015);
M. Kawasaki and Y. Tada, J. Cosmol. Astropart. Phys. 08 (2016) 041.

double
M. Kawasaki, N. Sugiyama, and T. Yanagida, Phys. Rev. D 57, 6050 (1998);
J. Yokoyama, Phys. Rev. D 58, 083510 (1998);
T. Kawaguchi et al., Mon. Not. Roy. Astron. Soc. 388, 1426 (2008);
P. H. Frampton, M. Kawasaki, F. Takahashi, and T. T. Yanagida, J. Cosmol. Astropart. Phys. 04 (2010) 023;
M. Kawasaki, A. Kusenko, Y. Tada, and T. Yanagida, Phys. Rev. D 94, 083523 (2016).

curvaton
M. Kawasaki, N. Kitajima, and T. T. Yanagida, Phys. Rev. D 87, 063519 (2013);
K. Kohri, C.-M. Lin, and T. Matsuda, Phys. Rev. D 87, 103527 (2013).

planckinflation
Planck Collaboration, P. A. R. Ade, et al., Astron. Astrophys. 594, A20 (2016).

anber2
M. M. Anber and L. Sorbo, Phys. Rev. D 81, 043534 (2010).

ng
For example, see K.-W. Ng, Int. J. Mod. Phys. A 11, 3175 (1996).

planckparameter
Planck Collaboration, P. A. R. Ade, et al., Astron. Astrophys. 594, A13 (2016).

perturbativity
R. Z. Ferreira, J. Ganc, J. Noreña, and M. S. Sloth, J. Cosmol. Astropart. Phys. 04 (2016) 039;
M. Peloso, L. Sorbo, and C. Unal, J. Cosmol. Astropart. Phys. 09 (2016) 001.

fugita
T. Fujita, R. Namba, Y. Tada, N. Takeda, and H. Tashiro,
J. Cosmol. Astropart. Phys. 05 (2015) 054.

carr16
B. Carr, F.  Kühnel, and M. Sandstad, Phys. Rev. D 94, 083504 (2016);
B. Carr, M. Raidal, T. Tenkanen, V. Vaskonen, and H. Veermäe, Phys. Rev. D 96, 023514 (2017).


sorboGW
L. Sorbo,  J. Cosmol. Astropart. Phys. 06 (2011) 003.

inomata
K. Inomata, M. Kawasaki, K. Mukaida, Y. Tada, T. Yanagida, Phys. Rev. D 95, 123510 (2017);
N. Orlofsky, A. Pierce, and J. D. Wells, Phys. Rev. D 95, 063518 (2017);
T. Nakama, J. Silk, and M. Kamionkowski, Phys. Rev. D 95, 043511 (2017).

PTA
L. Lentati et al., Mon. Not. Roy. Astron. Soc. 453, 2576 (2015);
Z. Arzoumanian et al. (NANOGrav), Astrophys. J. 821, 13 (2016);
R. M. Shannon et al., Science 349, 1522 (2015).

ligoback17
LIGO Scientific Collaboration and Virgo Collaboration: B. P. Abbott et al.,
Phys. Rev. Lett.  118, 121101 (2017).

ligoback16
LIGO Scientific Collaboration and Virgo Collaboration: B. P. Abbott et al.,
Phys. Rev. Lett.  116, 131102 (2016).

lisa
 N. Bartolo et al., J. Cosmol. Astropart. Phys. 12 (2016) 026.

ska
G. Janssen et al., Advancing Astrophysics with the Square Kilometre Array, PoS AASKA14 (2015) 037.

banks
T. Banks, M. Dine, P. J. Fox, and E. Gorbatov, J. Cosmol. Astropart. Phys. 06 (2003) 001.

clockwork
K. Choi, H. Kim, and S. Yun, Phys. Rev. D 90, 023545 (2014); 
K. Choi and S. H. Im, J. High Energy Phys. 01 (2016) 149;
D. E. Kaplan and R. Rattazzi, Phys. Rev. D 93, 085007 (2016).



