Primordial Black Holes and Gravitational Waves in Hybrid Inflation with Chaotic Potentials

We study the formation of primordial black hole (PBH) dark matter and the generation of scalar induced secondary gravitational waves (SIGWs) in a non-supersymmetric model of hybrid inflation with chaotic (polynomial-like) potential, including one-loop radiative corrections. A radiatively corrected version of these models is entirely consistent with Planck's data. By adding non-canonical kinetic energy term in the lagrangian, the inflaton experiences a period of ultra-slow-roll, and the amplitude of primordial power spectrum is enhanced to $O(10^{-2})$. The enhanced power spectra of primordial curvature perturbations can have both sharp and broad peaks. %In the enhancement mechanism we explore two possible extensions by employing a Guassian and a step size kinetic function. A wide mass range of PBHs is realized in our model with the frequencies of scalar induced gravitational waves ranged from nHz to kHz. We present several benchmark points where the PBH mass generated during inflation is around $(1 - 100) \, M_{\odot}$, $(10^{-9} - 10^{-7}) \, M_{\odot}$ and $(10^{-16} - 10^{-11}) \, M_{\odot}$. The PBHs can make up most of the dark matter with masses around $(10^{-16} - 10^{-11}) \, M_{\odot}$ and $(1 - 100) \, M_{\odot}$, and their associated SIGWs can be probed by the upcoming ground and space-based gravitational wave (GW) observatories. The evidence of stochastic process recently reported by NANOGrav may be interpreted as SIGWs associated with the formation of PBHs. These SIGWs may also be tested by future interferometer-type GW observations of SKA, DECIGO, LISA, BBO, TaiJi, TianQin, CE and ET.


Introduction
From astronomical and cosmological observations, there is convincing evidence that 85% of the matter in the Universe is in the form of cold, non-baryonic dark matter (DM) [1]. The study of Primordial Black Holes (PBHs) dates back to the 1960's and 70's [2,3,4], and shows that the PBHs may form due the collapse of large overdensities in the early universe. The early universe may contain regions with high densities at small scales that can trigger gravitational collapse to form PBHs. This PBH production can be tested through their effects on a variety of cosmological and astronomical processes, and therefore, can serve as an inspiring tool to probe physics in the very early Universe [5,6]. In particular, PBH could be a potential candidate for (a fraction of) dark matter (DM), which has drawn a lot of attention [7,8].
PBHs are non-baryonic, as they form before matter-radiation equality. The PBHs with masses 10 15 g would have evaporated by now emitting Hawking radiation [9]. The emitted particles may impact the gamma-ray background [10] and the abundance of light elements produced by the big bang nucleosynthesis (BBN) [11]. The PBHs with masses greater than 10 15 g, on the other hand may survive up to the present epoch and are expected to be constrained by their gravitational effects, such as gravitational lensing [12], dynamical effects on baryonic matter [13], or the fast radio bursts created by mergers of charged PBHs [14].
The Laser Interferometer Gravitational-Wave Observatory (LIGO), the Scientific Collaboration and the Virgo Collaboration have detected several events of GWs coming from the merger of black holes (BHs) [15,16,17,18,19]. Recently, the North American Nano hertz Observatory for Gravitational Wave (NANOGrav) Collaboration [36] has published an analysis of the 12.5 yrs pulsar timing array (PTA) data, where strong evidence of a stochastic process with a common amplitude and a common spectral slope across pulsars was found. These two observations moved the physicists attention toward the gravitational waves (GWs) generated by PBH-PBH mergers [20,21,22], as well as the scalar-induced GWs from the enhanced primordial density perturbations associated with PBH formation [23,24,25,26,27,28]. The GWs survey shall be a promising window to reveal the physical processes of PBH formations.
The scalar induced gravitational waves (SIGWs) associated with the formation of PBHs may be the source of NANOGrav signal [29], or the GWs detected by LIGO/Virgo. In order to produce PBHs in the radiation era from the gravitational collapse of overdense regions, it is required that the density of overdense regions exceed the threshold value at the horizon re-entry. The initial conditions for these overdense regions are produced during the inflationary era. To produce the desired abundance of PBHs one needs the primordial scalar power spectrum at small scales to be enhanced to P ζ (k > 1) ∼ O(0.01). This condition is also required to explain the NANOGrav signal if it is regarded as a SIGW. On the other hand, the constraint on the amplitude of power spectrum at large scales from the cosmic microwave background (CMB) anisotropy measurements from Planck [30] is P ζ (0.05) ∼ O(10 −9 ).
To produce enough abundance of PBH dark matter (DM) and SIGWs measurable by NANOGrav, the amplitude of the power spectrum at small scales should be enhanced at least seven orders of magnitude to reach the threshold value.
In this paper, we study the non-supersymmetric model of hybrid inflation and the formation of primordial black hole (PBH) dark matter with the generation of their associated scalar induced secondary gravitational waves (SIGWs). In order to produce the required abundance of PBHs, we enhance the power spectrum using the mechanism discussed in [31] and extend the mechanism to study two more possibilities. This mechanism relies on the incorporation of a non-canonical kinetic term with a function G(φ) which exhibits a peak at some value of the field φ p . Beyond this point, G(φ) falls exponentially, suppressing the scalar perturbations to the value observed today. With the extended mechanism, the predictions of our model are in much better agreement with the experimental bounds. Several other enhancement mechanisms are also discussed in the literature as well. The enhancement by ultraslow-roll inflation with an inflection point is discussed in [32]. The other possibility is fine-tuning the model parameters while keeping the total number of e-folds around 50-60 [33].
The layout of the paper is as follows. In Sec. 2 we present the non-supersymmetric hybrid inflation model. The inflation setup and PBH production is described in Sec. 3. PBH abundance is computed in Sec. 4 while the generation of SIGWs is studied in Sec. 5. Our conclusions are summarized in Sec. 6.

Description of the Model
The scalar potential of hybrid inflation (HI) can be expressed as a combination of Higgs potential V (χ) and inflaton potential δV (φ) with an additional term, g 2 χ 2 φ 2 , which represents the interaction between the Higgs field χ and inflaton φ. The treelevel hybrid inflation potential, therefore, can be written as where δV (φ) is the inflaton potential and is taken to be a chaotic polynomial-like potential, i.e., δV (φ) = λ p φ p with p > 0. Here, the role of the interaction term is to generate an effective (squared) mass, for the χ field in the χ = 0 direction. This direction is a local minimum for φ > φ c = √ 2κM g and can be used for inflation with effective single field potential given by where V 0 = κ 2 M 4 . The chaotic potential, here, provides the necessary slope for the slow-roll inflation in the otherwise flat-valley. We consider suitable initial conditions for inflation to occur only in the χ = 0 valley until φ = φ c is reached where inflation is terminated abruptly, followed by a waterfall phase transition. We now tend to include one-loop radiative corrections, as the tree level predictions are not consistent with the Planck 2018 results [47]. The radiative corrections arise from the possible coupling of inflaton with other fields. These couplings can contribute to the reheating process in order to recover the hot big bang initial conditions. The corrections arising from the coupling of inflaton to fermions or bosons may be termed as fermionic or bosonic radiative corrections. The one-loop radiative corrections to V (φ) in the inflationary valley can be found from the following form of Coleman-Weinberg formula [48], where A < 0 (A > 0) for fermonic (bosonic) radiative corrections. The fermionic radiative corrections have already been seen to play an important role for the chaotic inflation driven by the quadratic and the quartic potentials [49]. The fermionic radiative corrections generally reduce both r and n s in the chaotic inflation. In the following, we study the effect of fermionic radiative corrections on the tree-level predictions and compare them with the Planck's latest bounds on r and n s . Using Eqs. (3) and (4), the one-loop radiatively corrected hybrid inflationary (RCHI) potential can be written as, In order to discuss the predictions of the model, some discussion of the effective number of independent parameters is in order. Apart from the parameter λ p of the chaotic potential, the fundamental parameters of the potential in Eq. (1) are κ, g and M , which can be reduced to V 0 and φ c for the effective potential in Eq. (3). We, however, take V 0 and κ c ≡ g 2 /κ as the effective independent parameters with φ c = 2V 1/2 0 /κ c . With this choice we can develop a simple correspondence for the supersymmetric hybrid inflation for which κ c = g = κ [50].  Table 1: Hybrid inflation with chaotic potential parameters V 0 , λ p , A and φ c for p = 1, p = 2 and p = 2/3. φ * corresponds to the value of φ at the pivot scale k * = 0.05 Mpc −1 . The scalar spectral index n s and power spectra P ζ (k * ) = 2.15×10 −9 are evaluated at the pivot scale k * for the three chaotic potentials.

Inflation and PBH Production
PBHs are formed from the gravitational collapse of over-dense regions when their density contrasts at the horizon re-entry during radiation domination exceeding the threshold value. The overdense regions may be seeded from the primordial curvature perturbations generated during inflation. The feasible way to produce enough abundance of primordial black hole (PBH) dark matter (DM) is by enhancing the amplitude of the power spectrum at least seven orders of magnitude to reach the threshold at small scales. We employ the enhancement mechanism of the power spectrum at small scales proposed in Refs [31] using kinetic or K/G inflation. The kinetic inflation is defined when inflaton field's kinetic part is coupled to inflaton field function where m p = 1/ √ 8πG = 1. The background equations of motion are where K ,φ = dK(φ)/dφ. These equations can be written in terms of derivatives with respect to number of efolds n = ln(a) as follows The slow-roll parameters are defined as where the last slow roll parameter is specific for Kinetic Inflation. The slow-roll inflation is realized when | i | 1, where i = 1, 2, K. The second order action S (2) in perturbation theory of the curvature perturbation ζ, for Kinetic Inflation is [52], wherez = a(t)φ/H and ζ represents derivative with respect to the conformal time τ = dt/a(t). Using z = √ Kz in the quadratic action (15) and varying with respect to curvature perturbation ζ k = v k /z and v k in the Fourier space, we obtain the famous Mukhanov equations for scalar perturbations The above mode equation for ζ k can be written in terms of derivatives with respect to number of e-folds n as follows where we have used z = K(φ)aφ ,n , Numerically it is more convenient to solve the equation for ζ k rather than standard Mukhanov variable v k , as the former yields stable results. Solving the mode Eq. (17) with Bunch Davies vacuum [53], we obtain the scalar power spectrum on super horizon scales |kτ | 1, where * marks the horizon crossing values for each mode k. In the slow-roll approximation 4 , the scalar power spectrum can be estimated as, The scalar spectral index n s , tensor to scalar ratio r and the total number of e-folds N for KG inflation at the pivot scale are given by The CMB power spectrum as reported by Planck 2018 [30] is P ζ (k * ) = 2.15 × 10 −9 at the pivot scale k * = 0.05 Mpc −1 . For the PBH production, the power spectrum needs to be enhanced to P ζ (k) 10 −2 at small scales k > 10 2 Mpc −1 . This enhancement in the power spectrum can be achieved using K(φ) = 1 + G(φ), as some kind of peak function. We employ a polynomial peak function K q (φ), from [31] and introduce two functions; a Gaussian peak function K g (φ) and a step function K s (φ), where w is the width and h is the height of the peaks. These peak functions are plotted in Fig. 1 as a function of φ. It can be seen that the polynomial peak functions K q=1,2 (φ) have wider base and falls off slowly. The Gaussian peak function 4 5 6 7 8 9 10 0  The curves are drawn for Gaussian peak function K g (φ). The numerical values of the relevant parameters are given in Table 2.
K g (φ) falls off exponentially with a narrow base whereas, the step peak function K s (φ) has a flat top with the width exactly equal to w. Due to these peak functions, the first slow roll parameter becomes very small, O(∼ 10 −10 ). This leads to ultra slow-roll inflation; the inflaton field gets trapped in this local minimum for 15 to 40 e-folds as shown in Fig. 2, where , |η| and φ are plotted against the number of e-folds n = log(a) for Gaussian peak function K g (φ). Because of inverse relation between the power spectrum P ζ and , the sudden fall in the value of enhances the scalar power spectrum P ζ by seven order of magnitude. This enhancement in power spectra  Table 2: Benchmark points and predictions of gravitational wave spectrum Ω GW and PBH abundance Y PBH for models p = 1, p = 2 and p = 2/3 generated using Gaussian peak function K g (φ) and the step function K s (φ).  Table 3: Benchmark points and predictions of gravitational wave spectrum Ω GW and PBH abundance Y PBH for models p = 1, p = 2 and p = 2/3 generated using polynomial peak function K q=2 (φ). The labels NG, TL, ET and SE corresponds to peaks location of GW fraction Ω GW (f ) curve in regions of experiments; NANOGrav (NG), TaiJi/Lisa (TL), Einstein Telescope (ET) and SKA/EPTA (SE), respectively.

Model
can be realized in string theory inspired inflationary models due to large parametric space and the possibility of colliding branes [54]. We will evaluate power spectrum P ζ (k), energy spectrum of induced gravitational wave Ω GW (f ) and PBH abundance Y PBH using peak functions (25) - (27). We ensure that the total number of e-folds vary between 56 to 64, enough to solve the horizon  Figure 3: The scalar power spectrum P ζ (k) generated using the Gaussian K g (φ) and step function K s (φ). The panels are drawn using parameter sets for the three models; p = 1 (upper left), p = 2 (upper right) and p = 2/3 (bottom) as listed in Table 2. These functions generate broad power spectrum which appears extended and flat.
problem. The benchmark points and predictions of gravitational waves spectrum Ω GW and PBH abundance Y PBH for the hybrid inflation model with p = 1, p = 2 and p = 2/3 are listed in Table 2 and 3 for the peak functions K g (φ), K s (φ) and K q=2 (φ) along with the parameters h, w and φ p . These parameters are chosen to enhance the power spectrum O(10 −2 ) at small scales. The results for scalar power spectrum P ζ (k) are shown in Figs. 3 and 4 for our hybrid inflation model described by the scalar potential in Eq. (5). The curves in Fig. 3 are generated using the Gaussian peak function K g (φ) and step function K s (φ), defined in Eqs. (25) and (27). The curves in Fig. 4 are generated using the polynomial peak function K q (φ) (defined in Eq. (26)) with q = 2, for the models p = 1, p = 2 and p = 2/3. The four peaks in each panel correspond to the parameter sets; NG, TL, ET and SE, listed in Table 3. It can be seen that at large scales (0.05 Mpc −1 ), the power spectrum is of the order of O(10 −9 ), compatible with CMB constraints. At small scales (k > 10 2 Mpc −1 ) however, the power spectrum is enhanced to the order O(10 −2 ), large enough to produce PBHs after the horizon re-entry as discussed in the next section. Moreover, the Gaussian and step functions generate broad power spectrum which appears extended and flat, whereas the polynomial function generates peaked and narrow power  Table 3.
spectrum. The broad power spectrum has important implications for GWs induced by the curvature perturbations responsible for PBH formation, as discussed in Sec. 5. Note that the models satisfy the constraints from CMB µ-distortion, big bang nucleosynthesis (BBN) and pulsar timing array (PTA) observations, except for the parameter set SE for polynomial peak function.

PBH Abundance
The gravitational collapse of primordial curvature perturbation upon horizon re-entry during radiation dominated era may yield PBHs. The PBH mass is equal to γM hor , where M hor is the horizon mass and we take γ = 0.2 [55]. The current fractional energy density of PBHs with mass M to DM is [7] Y PBH (M ) = β(M ) 3.94 × 10 −9 where M is the solar mass, g * is the effective degrees of freedom at the formation time g * = 10.75 0.5 MeV < T < 300 GeV 107. 5 T > 300 GeV , (29) and Ω DM is the current energy density of the dark matter. The fractional energy density of PBHs at the formation is given by [56,57] where δ c denotes the critical density perturbation for the PBH formation and σ(k) being the mass variance associated with the PBH mass M (k) smoothing on the comoving horizon length k −1 = 1/(aH), given by [56] σ 2 (k) = 4 9 with the Gaussian window function W (x) = exp(−x 2 /2). To calculate PBH abundance, we take the observational values; Ω DM h 2 = 0.12 [58] and δ c = 0.4 [57,59]. The relation between PBH mass M and the scale k is given by [56] M (k) = 3.68 Using the approximation that the power spectrum is scale invariant, we obtain where µ c = 9δ c /4. Substituting the power spectrum obtained above into Eqs. (28), (30), (31) and (32), we obtain the PBH abundances as displayed in Fig. 5 for polynomial peak function K q (φ) with q = 2 and in Fig. 6 for Gaussian K g (φ) and step peak function K s (φ). The three panels in Fig. 5 are drawn for the models p = 1, p = 2 and p = 2/3 whereas, the peaks correspond to the parameter sets ET, TL and SE from Table 3. The shaded regions in the background represent the observational constraints on the PBH abundance from different experiments such as; accretion constraints by CMB, extragalactic gamma-rays by PBH evaporation (EG γ ), galactic center 511 keV gamma-ray line (INTEGRAL), white dwarf explosion (WD), microlensing events with Subaru HSC, the Kepler satellite and EROS/MACHO. For the polynomial peak function K q (φ) with q = 2 the models p = 1, p = 2 and p = 2/3 produce PBHs with masses M PBH (0.03, 0.04, 0.1), respectively. Finally, for the step peak function K s (φ), the parameter sets for models p = (1, 2, 2/3) produce PBHs with masses M PBH (10 −8 , 10 −16 , 10 −5 ) M and abundances Y peak PBH (0.004, 0.14, 0.098), respectively. Note that there are no observational constraints on PBH abundances in these mass ranges and therefore, PBHs can constitute all of the dark matter (DM).

Production of Secondary Gravitational Waves
The production of primordial black holes (PBHs) due to large curvature or density perturbations can induce secondary GWs due to second order mode coupling. These Secondary Induced GW (SIGW) are gauge invariant [29] stochastic background waves that could be observed by future GW experiments.
The tensor mode h k for SIGW is sourced by quadratic scalar perturbations function S k (Φ k ) [24,23] where H = aH, w = p/ρ, e ij (k) is polarization tensor and Φ k is the gauge invariant Bardeen potential [60]. During the radiation dominated era, these SIGWs decouple from their scalar part and plateau after horizon crossing. The energy density of these SIGWs today is given by [29,61] Ω GW (k) = 0.387 Ω r 6 g 4 * ,s g −3 * 106.75 where In the above expression, Ω r = 5.38 × 10 −5 is the radiation energy density and g * , g * ,s are effective degrees of freedom at horizon crossing for each mode. Because of second order mode coupling, Ω GW has quadratic dependence on scalar power spectrum P ζ (k).  Table 2. The shaded regions in the background represent the sensitivity of current and future gravitational waves observatories.
The GW spectrum for our hybrid inflation model is shown in Fig. 7 with Gaussian K g (φ) and step peak function K s (φ) and in Fig. 8 with the polynomial peak function K q (φ), q = 2. The peaks in these figures correspond to the parameter sets listed in Table 2 for functions K g (φ) and K s (φ) and Table 3 for the polynomial function K q (φ).
LIGO O1/LIGO O2 are also included, although the peaks lie outside these bounds. Note that all the peaks in Fig. 7 lie in the sensitivity bounds of NANOGrav and these SIGWs signals, associated with the formation of primordial black holes (PBHs), may be the source of stochastic process recently reported by NANOGrav analysis of 12.5 yrs of data [36]. The polynomial peak function K q (φ), q = 2 generates narrow Ω GW (k)h 2 peaks in different frequency ranges as shown in Fig. 8. The three models p = 1, p = 2 and p = 2/3, with parameter sets listed in Table 3, generate similar gravitational wave spectrum which can be seen in various experiments. For example, the Ω GW (k)h 2 peak generated by parameter sets NG and SE can be seen in NANOGrav and SKA and may explain the stochastic process recently reported by NANOGrav collaboration. Similarly, the peak generated by parameter set TL can bee seen in future gravitational wave detectors such as, TaiJi, TianQin, LISA, AEDEG, DECIGO and BBO. Finally, the peak corresponding to parameter set ET lies in the sensitivity region of future experiments such as, AEDEG, DECIGO, BBO, CE and ET. : Gravitational wave signals from hybrid inflation for the models p = 2, p = 1 and p = 2/3 compared to the NANOGrav observations. The dark and light shaded regions represent 1-and 2-σ bounds reported by NANOGrav [36]. The left panel corresponds to the Gaussian peak function K g (φ), whereas the right panel corresponds to the step function K s (φ).
We can also compare the SIGW signal associated with the formation of primordial black holes (PBHs) to the recent NANOGrav results [36], which constrain the amplitude and slope of a stochastic process. The amplitude of the SIGW from [62] is given as, which allows direct comparison of our results to the NANOGrav bounds in the Ω yr gw − n gw plane as shown by the dark (1-σ) and light (2-σ) shaded regions in Fig. 9. We extract the amplitude and slope by comparing the amplitude at pivot scale f * = 5.6 × 10 −9 Hz and taking the logarithmic derivative of Ω GW (f ) at the pivot scale, Ω yr gw = Ω GW (f * ) f yr f * ngw .
(41) Fig. 9 shows comparison of SIGW predictions, associated with the primordial black hole (PBH) formation, from hybrid inflation model with the constraints on the amplitude and tilt from NANOGrav [36]. The curves are drawn for the models with p = 1, p = 2 and p = 2/3 using Gaussian peak function K g (φ) (left panel) and step function K s (φ) (right panel). It is evident that the predictions from all three models lie within the 1-σ bound of NANOGrav. It is also important to emphasis that the broader peaks generated by Gaussian and step functions provide much better fitting to the NANOGrav bounds as compared to the sharp peaks of polynomial peak function.

Summary
To summarize, we have investigated production of primordial black holes (PBHs) and their associated induced secondary gravitational waves (SIGWs) using a background of hybrid inflation. The cosmological observables, such as the tensor-to-scalar ratio r and the scalar spectral index n s , are computed for various effective potentials. To produce the required abundance Y pbh for PBHs as a dark matter (DM) and their associated SIGWs, the curvature power spectrum is enhanced by seven order of magnitude, P ζ ∼ 0.01 at the scale k > Mpc −1 by employing a non-canonical kinetic energy term. We have utilized a well known polynomial peak function K q (φ) and employed two new functions; a Gaussian peak K g (φ) and a step function K s (φ), with a peak at φ p in order to enhance the power spectrum at small scales. The peak functions K(φ) induce an inflection point (with flat plateau) in the potential and effectively lead to ultra slow-roll inflation. However, the functions have a minor role away from the peak value where the usual slow-roll inflation is recovered, constrained by the CMB observations at large scales. Based on our analysis, it is quite possible that PBHs constitute most of the DM if the mass of PBH lies within the range (