Primordial black holes and secondary gravitational waves from natural inflation

The production of primordial black hole (PBH) dark matter (DM) and the generation of scalar induced secondary gravitational waves by using the enhancement mechanism with a peak function in the non-canonical kinetic term in natural inflation is discussed. We show explicitly that the power spectrum for the primordial curvature perturbation can be enhanced at $10^{12}$ Mpc$^{-1}$, $10^{8}$ Mpc$^{-1}$ and $10^{5}$ Mpc$^{-1}$ by adjusting the model parameters. With the enhanced primordial curvature perturbations, we show the production of PBH DM with peak masses around $10^{-13}\ M_{\odot}$, the Earth's mass and the stellar mass, and the generation of scalar induced gravitational waves (SIGWs) with peak frequencies around mHz, $10^{-6}$ Hz and nHz, respectively. The PBHs with the mass scale $10^{-13}\ M_{\odot}$ can make up almost all the DM and the associated SIGWs is testable by spaced based gravitational wave observatory.

To produce abundant PBH DM, the amplitude of primordial curvature perturbations at small scales needs to be in the order of 0.01 [28,29] while at large scales the cosmic microwave background (CMB) constraint on the amplitude of the power spectrum is A s = 2.1 × 10 −9 [30]. Inflationary models with inflection point or non-minimal coupling are usually considered to enhance the power spectrum at small scales [31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55]. Usually it is very difficult to enhance the power spectrum to the order of 0.01 at small scales while keeping the total number of e-folds to be 50-60 [41,46]. For canonical single field inflationary models with inflection point, the power spectrum either cannot be enhanced to the order of 0.01 [34,35,36,46] or the model parameters need to be fine tuned by more than six decimal digits [31,40,43,46]. By introducing a non-canonical kinetic term like k inflation [56,57] and G inflation [58,59,60,61], a new mechanism with a peak function in the non-canonical kinetic term was proposed to enhance the primordial power spectrum at small scales [62,63,64]. The forms of the peak function and the inflationary potential in the mechanism are not restricted, and both sharp and broad peaks are possible [63,64,65]. The mechanism works for both Higgs inflation and T-model. To the second order of perturbations, the first order perturbations are the sources of the second order tensor perturbation, so associated with the formation of PBHs the large curvature perturbations at small scales induce secondary GWs after the horizon reentry during the radiation dominated epoch [66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93]. These scalar induced GWs (SIGWs) consist of the stochastic background and they can detected by Pulsar Timing Arrays (PTA) [94,95,96,97,98] and the space based GW observatories like Laser Interferometer Space Antenna (LISA) [99,100], Taiji [101] and TianQin [102].
Nambu-Goldstone bosons are ubiquitous because they arise whenever a global symmetry is spontaneously broken. These particles become pseudo-Nambu-Goldstone bosons (PNGBs) if additional explicit symmetry is broken. For axions, the mass of PNGBs arises from nonperturbative instantons through the chiral anomaly. Instanton effects produce a periodic potential of height Λ 4 for PNGBs when the associated gauge group becomes strong at a mass scale Λ. Taking the symmetry breaking scale f a as the Planck scale M pl and the mass scale Λ as the scale of grand unification 10 15 GeV, we get a naturally flat potential without any fine tuning. The inflationary model with this flat potential is called natural inflation [103].
In this paper, we use the enhancement mechanism [63,64] to discuss the production of PBHs and SIGWs in natural inflation. The paper is organized as follows. In Sec. 2, we discuss the enhancement of the power spectrum, the production of PBHs and the generation of SIGWs in natural inflation. The conclusions are drawn in Sec. 3.

Natural Inflation
In this section, we consider natural inflation with a non-canonical kinetic term. The action is where X = −g µν ∇ µ φ∇ ν φ/2, the reduced Planck mass is taken to be M −2 Pl = 8πG = 1, and we introduce noncanonical kinetic terms like those in scalartensor theory of gravity, k inflation [56,57] or G inflation [58]. The potential for natural inflation is [103] where Λ ∼ 10 15 GeV is the energy scale of the potential, and f a ∼ M Pl is the symmetry breaking scale. The noncanonical kinetic function G(φ) = G p (φ) + f (φ) [63,64] with the function and the peak function G p (φ). Motivated by the noncanonical kinetic term X/φ in Brans-Dicke theory, we can take G p (φ) = hw/(φ − φ p + w) with hw ∼ O(1) and the dimensionless parameter w 1 to avoid singularity at φ = φ p . In this paper, we take more general peak function [63,64] where the dimensionless parameter h determines the height of the peak, the dimensionless parameter w controls the width of the peak, φ p normalized by the reduced Planck mass M Pl controls the peak position of the enhanced power spectrum, the power index q controls the shape of the enhanced power spectrum, and f 0 has the dimension of mass. Note that at low energies after inflation, the non-canonical kinetic function G(φ) is negligible and the standard canonical kinetic term is recovered. The background equations are where G φ = dG(φ)/dφ, V φ = dV /dφ.

The enhancement of the scalar power spectrum
To understand how the mechanism can enhance the power spectrum qualitatively, first we use slow-roll formulae to explain it. For the calculation of the power spectrum, we do not assume slow-roll and instead invoke numerical method to solve both the background and perturbation equations. The equation for the curvature perturbation ζ is where η is the conformal time η = dt/a(t), u k = zζ k , z = aφ(1 + G) 1/2 /H. In the slow-roll approximation, the power spectrum for the curvature perturbation is The scalar spectral index n s = 1 + d ln P ζ /d ln k. The tensor perturbations are not affected by the non-canonical kinetic term and the power spectrum for tensor perturbations is The tensor to scalar ratio is r = P T /P ζ . From Eq. (9), we see that the non-canonical function G(φ) can enhance the scalar power spectrum. Around the peak φ p , the major contribution to G(φ) is from G p (φ) and G(φ) ≈ h, so the scalar power spectrum can be enhanced by h. To achieve seven orders of magnitude enhancement, h should be at least the order of 10 7 . On the other hand, the number of e-folds before the end of inflation when the pivotal scale exits the horizon is so the enhancement also increases N and the contribution from the peak function is about 20 e-folds, here φ * and φ e are the values of the scalar field at the horizon exit and the end of inflation, respectively. This means that when the peak function G p (φ) enhances the scalar power spectrum at small scales, it also effectively moves φ * closer to φ e and the remaining number of e-folds N ef f ∼ 40 for the slow-roll inflation since the usual slow-roll inflation is recovered away from the peak. This changes the predictions of the scalar spectral tilt and the tensor to scalar ratio for natural inflation at large scales. Away from the peak, the peak function is negligible and the non-canonical function f (φ) dominates, we can take the transformation [63,64] dΦ = 1 + f (φ)dφ (12) to change the non-canonical scalar field φ to be canonical scalar field Φ with the potential for the function f (φ) is derived from the above transformation (12). For the power law potential U (Φ) = U 0 Φ 1/3 and N ef f ∼ 40, we get n s = 0.971 and r = 0.033 which are consistent with Planck 2018 results. If we use the full peak function G(φ) in Eq. (12) to change the non-canonical scalar field φ to be canonical scalar field Φ, then we can get the effective potential U ef f (Φ) of the canonical field Φ. In general, there is no analytical relation between φ and Φ and we cannot obtain the expression for the effective potential. However, we can use numerical method to get the effective potential U ef f (Φ) as shown in Fig. 1 for the model N1 defined below. The effective potential U ef f (Φ) of the canonical field Φ has an inflection point. The existence of inflection point explains the enhancement of the power spectrum. Therefore, we see that for natural inflation it is possible to satisfy the large scale constraints and enhance the scalar power spectrum at small scales. To show this, we numerically solve the background Eqs. (5)-(7) and the perturbation Eq. (8) to get both the scalar and tensor power spectra P ζ and P T . Then we numerically calculate the scalar spectral index n s and the tensor to scalar ratio r from P ζ and P T . We take f 0 = 7500, φ * = 14.5, the symmetry breaking scale f a = 7, and the values of the parameters Λ, h and φ p as shown in Table 1. In table 1, we use the label "N" to represent the model with the parameter q = 1 which produces sharp peaks at small scales in the scalar power spectrum, and the label "WN" to represent the model with the parameter q = 6/5 which produces broad peaks at small scales in the scalar power spectrum. We also use the labels 1, 2, 3 to distinguish the models with different peak scales in the scalar power spectrum, the peak scale for the model labelling as 1 is around 10 12 Mpc −1 , the peak scale for the model labelling as 2 is around 10 8 Mpc −1 , and the peak scale for the model labelling as 3 is around 10 5 Mpc −1 . The numerical results for N , n s = 1 + d ln P ζ /d ln k, the tensor to scalar ratio r = P ζ /P T , and the peak scale k peak are summarized in Table 1, and the scalar power spectra for these models are shown in Fig. 2. We also show the peak values of the power spectra in Table 2. From these results, we see that the results n s ≈ 0.968 and r ≈ 0.04 at the pivotal scale k * = 0.05 Mpc −1 are consistent with the observational constraints [30] n s = 0.9649 ± 0.0042 (68%CL), r 0.05 < 0.06 (95%CL).
As discussed above and in Ref. [63], the parameter q controls the shape of the peak and the parameter φ p determines the peak scales. By choosing different values of φ p and q, we can produce either sharp peaks or broad peaks with different peak scales.  Figure 2: The results for the scalar power spectrum. The solid lines denote the models with the parameter q = 1 and the dashed lines denote the models with the parameter q = 6/5. The blue, red, black lines denote the models with the peaks around 10 12 Mpc −1 , 10 8 Mpc −1 and 10 5 Mpc −1 , respectively. The parameters for the models and the peak scales k peak are shown in Table 1. The peak values of the power spectra are shown in Table 2. The lightgreen shaded region is excluded by the CMB observations [30]. The pink, cyan and orange regions show the constraints from the PTA observations [104], the effect on the ratio between neutron and proton during the big bang nucleosynthesis (BBN) [105] and µ-distortion of CMB [106], respectively.
We also calculate the three-point correlation function numerically to get the bispectrum B ζ [107,108,109], whereζ k is the corresponding quantum operator of the curvature perturbation ζ k . The non-Gaussianity parameter f NL is [107,110] f NL (k 1 , k 2 , k 3 ) = 5 6 where P ζ (k) = |ζ k | 2 = 2π 2 P ζ /k 3 . The non-Gaussianity parameter f NL in the equilateral and squeezed limits are shown in Figs. 3 and 4, respectively.

Primordial black hole dark matter
During radiation dominated era, PBHs may form by gravitational collapse from the enhanced primordial curvature perturbations at small scales when they reenter the horizon. Ignoring the mass accretion and evaporation in the Press-Schechter formalism, the fractional energy density of PBHs at the formation is [111,112,113]  where the critical density perturbation for the PBH formation is δ c = 0.4 [113,114,115,116] and the mass variance σ(M ) associated with the PBH mass is W (kR) is the window function to smooth the density contrast, the smoothing scale R = (aH) −1 , and the power spectrum P δ of the matter perturbation is related with the primordial curvature as Although there exist different window functions [117], we take the Gaussian window function Substituting the Gaussian window function (19) and Eq. (18) into Eq. (17), the mass variance becomes where x = kR. Note that the main contribution to the integral (20) is from the scale x ∼ 1 and the result is not affected much by other scales, so to a good approximation, the primordial curvature power spectrum can be treated as scale invariant even though it may change rapidly [29]. With this assumption, the mass variance (20) becomes Substituting Eq. (21) into the definition (16), we obtain where µ c = 9δ c /2 √ 2. The current fractional energy density of PBHs with mass M to DM [22,31] where M is the solar mass, γ = 0.2 [118], the current energy density parameter of DM is taken to be Ω DM h 2 = 0.12 [119], the effective degrees of freedom g * = 107.5 for T > 300 GeV and g * = 10.75 for 0.5 MeV < T < 300 GeV. The peak mass scale M for PBHs and the peak scale k peak is estimated to be [31] M (k) = 3.68 γ 0.2 g * 10.75 Note that there are subtleties for the simple relationship between β(M ) and P ζ . When we consider the non linearities between the Gaussian curvature perturbation and the density contrast and the non-linear effects arising at horizon crossing, the value of δ c may become larger [120]. On the other hand, the abundance of PBHs also depends on the shape and non-Gaussianity of P ζ and non-linear statistics [121,122,123,124]. If we consider non-Gaussianities of the primordial power spectrum, then the PBH abundance at the peak is [125] whereζ = ζ/ P ζ (k peak ) and f NL (k peak , k peak , k peak ) P ζ (k peak ).
For a good estimation, we can use the parameter J peak to characterize the effect of the non-Gaussianities. If J peak 1, then we can ignore the effect of non-Gaussianities on the PBH abundance. From Figs. 3 and 4, we see that non-Gaussianities around the peak scale are small and henceforth J peak 1, so the effect of non-Gaussianities on the PBH abundance is negligible.
Substituting the numerical results of the power spectra obtained in the previous subsection into Eqs. (22), (23) and (24), we obtain the abundance and the peak mass of PBHs and the results are shown in Fig. 5 and Table  2. In Fig. 5, the constraint on the abundance of PBHs from white dwarf explosion [126] is not included because it is not robust [127]. From Fig. 5 and Table 2, we see that different peak scales in the scalar power spectrum correspond to different masses of PBHs. In particular, the models generate PBHs with masses around 10 −13 M , 10 −6 M and 10 M , respectively. The stellar mass PBHs may explain BHs observed by LIGO/Virgo collaboration. The peak abundance of PBHs with the mass scale 10 −13 M is Y peak PBH ≈ 1, so they can make up almost all the DM. PBH DM with the mass around O(1)M ⊕ may explain the planet 9. The results for the peak value of the primordial scalar power spectra, the peak mass and peak abundance of PBH and the peak frequency of SIGWs. The parameters for the models are shown in Table 1

Scalar induced secondary gravitational waves
For the perturbation to the second order, the first order and second order perturbations are mixed and the first order scalar perturbations become the   Table 1, the peak abundance and the peak mass of PBHs are shown in Table 2. The shaded regions show the observational constraints on the PBH abundance: the yellow region from accretion constraints by CMB [128,129], the red region from extragalactic gamma-rays by PBH evaporation [130] (EGγ), the cyan region from galactic center 511 keV gamma-ray line (INTEGRAL) [131,132], the green region from microlensing events with Subaru HSC [133], the blue region from the Kepler satellite [134], the gray region from the EROS/MACHO [135]. The magenta solid line shows the constraints on the stochastic gravitational wave background by LIGO/Virgo [9,136] and the black dots shows the limits from the LIGO/Virgo merger rate [9,137].
source of the second order tensor perturbations. Therefore, the large primordial scalar perturbations at small scales induce secondary GWs associated with the production of PBHs when they reenter the horizon. The equation for the Fourier components of the second order tensor perturbations h k is [68,69] h where the source from the first order scalar perturbations is H = aH, w = p/ρ = 1/3 during radiation domination, e ij (k) is the polarization tensor, Φ is the gauge invariant Bardeen potential. The energy density of SIGWs in the radiation domination is [78,83,84] where the kernel function I 2 RD is [84,28] T c (u, v, x) andT s (u, v, x) are given in Ref. [28]. The current energy density of SIGWs is where Ω r is the fraction energy density of radiation. Plugging the power spectra in Fig. 2 into Eqs. (29) and (31) and using Eq. (30) we obtain current energy densities of SIGWs and the results are shown in Fig. 6. The peak frequencies are shown in Table 2. From Table 2, we see that the peak frequencies of the SIGWs are around mHz, µHz and nHz respectively. For the models WN1, WN2 and WN3, Ω GW has a broad shape which spans a wide frequency bands because the enhanced power spectrum has a broad peak. The broad shape for the model WN3 which produces the stellar mass PBHs with the abundance in the order of 10 −3 leads to its exclusion by the EPTA data [94,95,96,97], but the broad shape for the model WN2 which produces the earth-mass PBHs to explain the planet 9 makes the model testable by LISA and Taiji. The models N1 and WN1 can explain DM in terms of PBHs and they can be tested by LISA/Taiji/TianQin. The model N3 may explain the possible stochastic GW background detected in the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) 12.5-year data [138]. If the abundance of PBHs produced in the model WN3 is in the order of 10 −15 , then the model WN3 can also explain the possible stochastic GW background detected in the NANOGrav 12.5-year data. For comparison, we also show the energy density of the primordial GWs from the model N1 in Fig. 6. The energy density of the primordial GWs is in the order of 10 −16 and the peak values of the energy densities of the SIGWs generated in our models are in the order of 10 −8 .  Figure 6: The energy densities of SIGWs. The parameters for the models are shown in Table 1 and the peak frequencies are shown in Table 2. The dashed pink curve denotes the EPTA limit [94,95,96,97] , the dotted cyan curve denotes the SKA limit [98], the shaded region is the observational result from NANOGrav 12.5-year data [138], the dashed green curve in the middle denotes the TianQin limit [102], the dot-dashed magenta curve shows the TaiJi limit [101], the dashed brown curve shows the LISA limit [100], and the dashed gray curve denotes the aLIGO limit [139,140]. For comparison, the energy density of primordial GWs from the model N1 is also shown with the gray solid line.

Conclusion
The amplitude of primordial curvature perturbations at small scales can become large by the enhancement mechanism with a peak function in the non-canonical kinetic term [62,63,64]. We apply the enhancement mechanism to natural inflation to produce abundant PBHs and observable SIGWs. The power spectrum at large scales is consistent with observational constraint and the power spectrum at small scales is enhanced to the order of 0.01. Either sharp peak or broad peak is possible with different peak shapes for the peak function by choosing different values of q. By adjusting the peak position φ p in the peak function, the power spectrum is enhanced at different scales, henceforth associated with the generation of SIGWs with different peak frequencies, PBHs with different masses are produced. We choose three different values of φ p to get enhanced power spectrum at 10 12 Mpc −1 , 10 8 Mpc −1 and 10 5 Mpc −1 , respectively. The enhanced curvature perturbations produce PBH DM with peak masses around 10 −13 M , the Earth's mass and the stellar mass, and SIGWs with peak frequencies around mHz, µHz and nHz. The stellar mass PBHs may explain BHs observed by LIGO/Virgo collaboration, and the earth-mass PBHs may explain the planet 9. The PBHs with the mass scale 10 −13 M can make up almost all the DM. The peak energy densities of SIGWs for the models discussed in this paper are around 10 −8 while the energy density of primordial GWs is around 10 −16 . The SIGWs with the peak frequency around nHz is testable by PTA observations, and SIGWs with the peak frequency around mHz is testable by space based GW observatory. In particular, the SIGWs produced in the model N3 can explain the stochastic GW background observed by NANOGrav 12.5year experiment. The broad shape for the model WN3 which produces the stellar mass PBHs with the abundance in the order of 10 −3 leads to its exclusion by the EPTA data and the broad shape for the model WN2 makes the model testable by LISA and Taiji. If the abundance of PBHs produced in the model WN3 is in the order of 10 −15 , then the model WN3 can also explain the possible stochastic GW background detected in the NANOGrav 12.5-year data. These results are similar to those obtained in [62,63,64,65]. The detailed shape of the primordial power spectrum and non-Gaussianity in different models are not exactly the same, more accurate measurements may be able to distinguish different models.
In conclusion, the enhancement mechanism with a peak function in the non-canonical kinetic term works for natural inflation.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.