On primordial black holes and secondary gravitational waves generated from inflation with solo/multi-bumpy potential

Since the beginning of the 21st century, the science of astronomy and astrophysics have been developing rapidly, and researchers have been paying increasing attention to a type of evenly mysterious objects: black holes. From 2015 to date, there have been dozens of gravitational wave events, where gravitational wave signals are released from either binary black holes or other compact objects like white dwarfs and neutron stars [1, 2]. In 2019, the working group of the Event Horizon Telescope (EHT) claimed to obtain the first images for the shadow of the supermassive black hole in the center of the M87 galaxy [3]. Moreover, the recent (2020) Nobel Prize in Physics was also awarded to studies on the black hole. Along with the development of modern science and technology, it is expected that more information about this object will be discovered with time.

Black holes are classified into several categories, based on how they are formed. One of these black-hole categories is of particular interest to astronomers and cosmologists and are called primordial black holes (PBHs) [4-6]. Unlike astrophysical black holes, PBHs are not formed from the collapse of stars. Instead, they are formed in the primordial universe, which exists significantly earlier than the star formation, owing to the gravitational collapse by the overdensity of cosmic spacetime in local patches. For this reason, the formed PBHs can have broad mass ranges, which are not constrained by Chandrasekhar and Oppenheimer limits.

To form PBHs, small-scale inhomogeneities are required, and these can be provided by the mechanism of primordial-perturbation production during the inflation process in the early universe. In this inflation scenario [7-9], where the universe expands to very large extensions within a short time, the quantum fluctuations of the inflation field's vacuum will be stretched out of the horizon and become classical perturbations (see, e.g., [10] for a comprehensive review). If the inflation has some feature that makes these perturbations exceed certain threshold values at small scales, then when the perturbations re-enters the horizon, it will cause large inhomogeneities of the universe, and PBHs will be formed. There have been long discussions on the formation of PBHs from inflation models [11-53]. Among these models, the simplest is a single field inflation in the framework of General Relativity; however, as has been pointed out in [19, 20], the slow-roll condition has to be violated to form PBHs that comprise all the dark matter. Consequently, slow-roll violating models have become an interesting alternative, including ultra-slow-roll inflation [21], inflation with inflection points or bumps [19, 26, 34, 36, 41, 45, 54], etc.

However, although it is well-known that PBHs can function as dark matter, they have been constrained by several experiments, such as Subaru Hyper Suprime-Cam (Subaru-HSC) [55], Experience de Recherche d’Objets Sombres (EROS) [56], Optical Gravitational Lensing Experiment (OGLE) [57], Cosmic Microwave Background (CMB) [58], Femtolensing of Gamma-ray Bursts (FL) [59], White Dwarf Explosions (WD) [60] (also refer to controversies in [61, 62]), 511 keV gamma-ray line [63-66], BH evaporation [67], Neutron Stars (NS) [68], NANOGrav [69, 70], LIGO [71], Leo-I dwarf galaxy [72], Gravitational-Wave Lensing (GW-Lensing) [73], and so on [74]. Owing to Hawking radiation, PBHs smaller than $ 10^{15} $ g can now be evaporated; however, there are still a few observational constraints, such as the big bang nucleosynthesis (BBN) [67], Extragalactic γ-rays [67, 75], CMB spectral distortions, and anisotropies (CMB sda) [76]. For the constraints of the small-mass PBHs, refer to [77, 78]. Currently, the upper limits on the PBH fraction of dark matter $ f_{\rm PBH} $ provided by these experiments have covered almost all the mass ranges of PBHs; however, most of the limits can only reach up to about $ 10^{-5} $. Therefore there are still large blanks for the PBHs in other mass ranges to have larger fractions. In this study, we investigate the generalized bumpy inflation model, whose potential contains multiple bumps. For the multi-bumpy potential, it is expected that the PBHs can be formed in various mass ranges; hence, by comparing it to a single mass range, the detectability of the PBHs can get improved. If the accuracies of these experiments get further improved and the upper limits are put even lower in the future, we can expect that the PBHs generated in our model are more likely to be detected. For example, recently, researchers have studied the electromagnetic signals from the bubbles of a single PBH star and the possibility of their detection via high-energy gamma-ray experiments, especially for small-mass primordial black holes [79]. Note that similar multi-range PBH formations can be realized in [80, 81], uusing a multi-field inflation model.

However, the enhancement of the primordial curvature perturbation will also induce the generation of secondary gravitational waves [82, 83]. Therefore, the observation of SIGWs will further constrain the enhancement of the amplitude of the primordial curvature perturbation during the inflation period.

The remainder of this paper is organized as follows: In Sec. II, we present the general formulation of the inflation process at both background and perturbation levels. For analyticity, we impose slow-roll and ultra-slow-roll conditions. In Sec. III, we introduce our inflation models, where the inflation potential has one or multiple bumps (we take the example of 3 bumps), and then calculate the slow-roll parameter power spectrum numerically. In Sec. IV, we calculate the fraction of generated PBHs from our models in dark matter, using the Press-Schechter formalism. We also consider the formation of PBHs in the case of spherical and ellipsoidal collapses. In Sec. V, we discuss the generation of SIGWs, which is constrained by the latest observations. In Sec. VI, we present our conclusions and discussions. Moreover, in this study, we work in mostly-plus signatures for the metric, (−,+,+,+), and the unit where $ c = \hbar = 1 $ and $ M_{pl} = 1/\sqrt{8\pi G} $ in this paper.

First, we start with the inflation scenario driven by a single scalar field with an arbitrary potential. The action of the inflation model is expressed as

From the above action, we can obtain the Friedmann equation:

where $ H = \dot a/a $ represents the Hubble parameter; a is the scale factor and dot denotes a derivative with respect to cosmic time t. We can also obtain the equation of motion for the inflation field ϕ:

where $ V'(\phi) = {\rm d}V(\phi)/{\rm d}\phi $. To maintain the inflation process, it is usually required that the slow-roll condition should be satisfied. We define the slow-roll parameters as:

where the subscript H denotes that these parameters are defined for the Hubble parameter, or in other words, the "whole universe." The slow-roll condition states that during inflation, all these parameters have values that are significantly less than 1, namely

However, another set of slow roll parameters can also be defined using the inflation potential itself:

and when the slow-roll condition (7) holds, $ \epsilon_H\simeq \epsilon_V $, $ \eta_H\simeq\eta_V-\epsilon_V $ is obtained. Moreover, the total number of e-foldings is usually adopted to describe the duration of inflation:

where the subscripts i and e denote the parameters at the beginning and ending time of the inflation, respectively. To solve the notorious big-bang problems, it is usually required that $ N\geq 60 $ [84].

Now, we turn to the perturbations generated during inflation. The perturbed FRW metric in its 3+1 decomposed form is [85]:

where α, β, ζ denote the scalar-type perturbations, and $ h_{ij} $ represents the tensor-type perturbation. The comoving time τ is defined by the relationship $ {\rm d}\tau\equiv a^{-1}(t){\rm d}t $. In uniform-ϕ gauge, α and β become constraint degrees of freedom only, and the only dynamical scalar perturbation is ζ [86, 87]. To obtain an analytical solution of ζ, we define the well-known Mukhanov-Sasaki (MS) variable u as:

and this variable satisfies Mukhanov-Sasaki equation in the Fourier space as [88, 89]:

where prime denotes the derivative with respect to comoving time τ, and the effective potential term is given by the following expression:

For a given mode k, and for a sub-Hubble region with $ k\gg aH $ at sufficiently early times, we can assume $ u_{k} $ to be quantum fluctuations in the Bunch-Davies vacuum [90] satisfying

The above solution can be viewed as the initial condition of the MS variable $ u_k $. However, because the inflation process will lead to the decline of the comoving Hubble radius, the k-mode will exit the horizon and enter the super-Hubble region with $ k\ll aH $, while the quantum fluctuations will be decoherent and become classical perturbations. By solving the MS equation and adopting the relationship between the MS variable and its dimensionless primordial power-spectrum [16]

we obtain the power spectrum expression under the slow roll approximation:

On the large cosmological scale accessible to CMB observations, the power spectrum usually takes the power-law form:

where $ A_S = P_S(k_*) $ is the amplitude of the scalar power spectrum at the pivot scale. In the slow-roll approximation, the scalar spectral tilts $ n_{s} $ are given by [84]

Recent CMB observations suggest $ P_S \simeq 2.1\times 10^{-9} $ at $ k_* = 0.05\;{\rm{Mpc}}^{-1} $ [91]; however, to increase the abundance of PBH to an appropriate order of magnitude, $ P_S $ must reach the amount of $ O(10^{-2}) $. Note that when the slow-roll condition holds, $ P_S $ will be conserved in time and will not deviate significantly from its pivot value; hence, it is impossible to generate sufficient PBHs, and we must break the condition. This is also pointed out in recent papers such as [19, 20]. A possible approach is to have $ \epsilon_H $ decrease rapidly at small scales, while results at large scales are kept unchanged. One such model is the ultra-slow-roll (USR) inflation model [92-94]. This model has an absolutely flat potential; therefore,

Therefore, $ \epsilon_H $ decreases with the explosion of a, triggering an abrupt increase in $ P_{S} $ according to Eq. (17) [21]. Moreover, in this case, the slow-roll condition gets violated because $ \eta_H = 3+\epsilon_H>1 $, and the above relationship between two sets of the slow-roll parameters does not hold any longer. Another example is the featured potiential like the one with a bump [41, 45]. When the inflation passes through the bump, it will enter a state similar to USR, which will achieve a local enhancement of the power spectrum.

In this subsection, we construct an inflation field with one bump in its potential. We consider the potential to be in the following form:

Here, the basic potential is of the power-law form while there is a Gaussian function-like bump on the potential, with the parameters b, c, and d denoting the height, central value in ϕ, and width of the bump. In general, the power spectrum of the slow-roll inflation potential does not contribute to the peak in the power spectrum, which is required for the formation of the primordial black holes. In this case, it is common to add a local Gaussian bump in the inflation potential, which has been studied in [41, 54, 95]. In the following part, we will see that by choosing the parameters practically, the addition of local Gaussian bumps can effectively slow down the inflation field and enhance the scalar power spectrum associated with the formation of the primordial black hole. From this potential, the following can be obtained:

In Fig. 1, we plot various relationships between variables involved in this case. In the numerical calculation, we select the parameters as follows:

Figure 1.  (color online) Upper left figure: the inflation potential of $ p = 2 $, in which the bump is located at $ \phi = 11M_{pl} $. Upper right figure: the evolution of ϕ with the parameter N, which enters the USR-like stage near $ N = 40 $. Lower left figure: the evolution of the slow roll parameters $ |\epsilon_{H}| $ and $ |\eta_{H}| $, in which the blue line corresponds to the value of 3 on the vertical axis. Lower right figure: the relationship between the power spectrum $ P_{S} $ (in Logarithmic scale) and the wave-number k. Our results are consistent with the observational constraints.

From the plot, it can be easily observed that there is a bump in the potential at $ \phi = 11M_{pl} $. Owing to the bump, the inflation will stay around $ \phi = 11M_{pl} $ for a short period. However, if the initial velocity of ϕ is insufficient, inflation may stop at the bump forever and will not continue to roll, leading to an eternal inflation. Conversely, if the initial velocity of the inflation is too large, the inflation will pass through the bump very quickly, which will make the peak value of the power spectrum too low to have sufficient abundance of PBHs. Therefore, both the initial velocity and the parameters of the bump need to be carefully settled. In the numerical calculation, we set the initial value of ϕ to be $ \phi_i = 16.2M_{pl} $, from which ϕ rolls down, and the initial velocity of ϕ to be $ \dot\phi_i = 4.5\times 10^{-6} $, to satisfy the observational constraint from CMB.

We solve the equation of motion (3) for the scalar field ϕ to obtain its evolution with respect to the e-folding number N. From the plot, it can be obseved that the inflation meets the bump near $ N = 40 $. We also plot the evolution of $ \epsilon_H $ and $ \eta_H $ with respect to N. To increase the power spectrum amplitude to $ 10^{-2} $, the inflation will stay at the bump for approximately 7 e-folding numbers. Because of the very slow motion of the inflation, the slow roll parameter $ \epsilon_H $ gets a negligible value, while $ \eta_H $ becomes significantly large ($ \sim O(1) $), just like the case of USR inflation models. Moreover, the inflation can stop when $ \epsilon = 1 $ at approximately $ N = 70 $.

We also solve the MS Eq. (13) numerically and obtain the variation of the power spectrum $ P_{S} $ with respect to the wave number k. It can be observed that the power spectrum fits the constraints from CMB on large scales with a peak value of $ O(10^{-2}) $ (the numerical value is 0.02284) at $ k = 1\times 10^{14} $ Mpc⁻¹.

For a further verification, we choose a second set of parameters as

and plot the same figures in Fig. 2. In this case, it can be observed that the decrease in $ \epsilon_H $ occurs around $ N = 20 $ during inflation, and a peak appears in the power spectrum at around $ k = 1\times 10^{6}\;{\rm{Mpc}}^{-1} $. Therefore, we find that the Gaussian-function-like bump on the potential does lead to the peak on the spectrum at small scales, provided we correctly choose the parameters. Such a feature in the spectrum can in principle generate PBHs, which we will discuss in the next section.

Figure 2.  (color online) Upper left figure: the inflation potential of $ p = 2/3 $, where the bump is located at $ \phi = 6M_{pl} $. Upper right figure: the evolution of ϕ with the parameter N, which enters the USR-like stage near $ N = 20 $. Lower left figure: the evolution of the slow roll parameters $ |\epsilon_{H}| $ and $ |\eta_{H}| $, where the blue line corresponds to the value of 3 on the vertical axis. Lower right figure: the relationship between power spectrum $ P_{S} $ (in Logarithmic scale) and the wave-number k. Our results are consistent with the observational constraints.

In this subsection, we extend the model discussed above to the case where the potential have more than one bump. In this case, there will be multiple peaks in the power spectrum, leading to the formation of PBHs in multiple mass ranges. It is interesting in the sense that, in the future, as the accuracy of the experiments improve continuously, it is expected that PBHs can be caught in various mass ranges, which can thus be explained by our model. In the following, we take the potential to be of the following form:

In principle, there could be arbitrary numbers of bumps in the potential; however, in this paper, we take an example of three Gaussian function-like bumps for simplicity and analycity. The specific parameters of Eq. (26) are given as

In this case, the inflation rolls down from the value of $ \phi = 16.2M_{pl} $ and satisfies the observation constraint of CMB from Fig. 3. We set $ V_{0} = 1.45\times 10^{-11}M_{pl}^{4} $.

Figure 3.  (color online) Upper left figure: the inflation potential with three bumps locating at $ \phi = 11M_{pl} $, $ \phi = 12M_{pl} $ and $ \phi = 13M_{pl} $. Upper right figure: the evolution of ϕ with the parameter N, which enters the USR-like stage three times near $ N = 20 $, $ N = 35 $ and $ N = 50 $. Lower left figure: the evolution of the slow roll parameters $ |\epsilon_{H}| $ and $ |\eta_{H}| $, where the blue line corresponds to the value of 3 on the vertical axis. Lower right figure: the relationship between power spectrum $ P_{S} $ (in Logarithmic scale) and the wave-number k. Our results are consistent with the observational constraints.

It is easy to observe that there are three USR stages in Fig. 3, which means that there will be three peaks in the power spectrum. For initial conditions, we still set $ \phi_i = 16.2M_{pl} $ and $ \dot\phi_i = 4.5\times 10^{-6} $, which is the same as for the solo-bumpy model. However, it differs from the solo-bumpy ones, where we have to take care of not only the shape of each bump (height, width, etc), as well as the relative distance of the bumps. This is important because, when the inflation field passes through one bump, it will lose kinetic energy, and if the bumps are far from each other, it may not have enough energy to pass through the next ones. Therefore, in this study, we set the bumps close to each other, namely at $ \phi = 11M_{pl} $, $ \phi = 12M_{pl} $ , and $ \phi = 13M_{pl} $. It can be observed from the plot that we will have three times $ \eta_H $, which breaks the slow-roll condition at approximately $ N = 20 $, $ N = 35 $ and $ N = 50 $, leading to three peaks on the power spectrum in small scales, while the inflation finally ends at approximately $ N = 87 $, and both $ \epsilon_H $ and $ \eta_H $ exceed unity. As will be demonstrated in the next section, this will form PBHs in three different mass ranges.

First, we briefly review how the PBHs are formed from inflation and how to relate the mass of PBHs to the inflation power spectrum. Generally, if the fluctuations re-entering the horizon become sufficiently large, the primordial black holes with extensive mass will form due to gravitational collapse. The mass of PBHs formed at a certain epoch in the radiation dominated era with fluctuation mode $ k_{\rm PBH} $ is given by the Hubble mass at that epoch up to an efficiency factor γ [96, 97]:

such that when the PBHs are formed, the Hubble radius and the energy density of the universe are denoted by $ H_{\rm form}^{-1} $ and $ \rho_{\rm form} $, respectively, and $ M_H $ is the total mass inside the Hubble radius. Moreover, in the radiation dominated era, we obtain $ \gamma\simeq0.2 $. With the Friedmann equation: $ 3H_{\rm form}^{2} = 8\pi G\rho _{\rm form} $, Eq. (28) can be written in another form:

Note that in the radiation dominated era, $ \rho = \dfrac{\pi ^2}{30}g_*T^4 $, and in the adiabatic environment of the early universe, entropy density: $ s\sim g_*T^3\sim a^{-3} $, therefore $ T\sim g_{*}^{-1/3}a^{-1} $ [98]. Then, we can obtain [99]

where $ g_* $ is the total effective degree of freedom of the universe in the radiation dominated era. The subscript "form" denotes variables that are evaluated when PBHs are formed, and the subscript "eq" denotes variables that are evaluated when radiation and matter are equal to each other. Note that in the radiation/matter equality, $ g_{*{\rm eq}}\simeq3.38 $ and $ k_{\rm eq}\simeq0.07\varOmega_{m0}h^2\;{\rm{Mpc}}^{-1} $. Hence, it can be observed that the mass of PBHs are determined by the scale of their formation denoted by $ k_{\rm form} $. For a given $ k_{\rm form} $, we can obtain the corresponding $ M_{\rm PBH} $.

We are interested in whether the PBHs formed can act as dark matter, or more precisely, how much fraction of dark matter can be contributed from the PBHs. To evaluate this, we need to define the fraction of PBHs in dark matter, namely

and this can be connected to the fraction of PBHs in the entire universe

via the equation:

where subscript "0" denotes variables evaluated today, and $ \Omega_{\rm DM0} $ is today's density fraction of dark matter that can be obtained from the observations. Combining Eqs. (30) and (33) and doing some manipulations (the details of which could be found in [67]), we can obtain the relationship between $ \beta(M) $ and $ f_{\rm PBH} $:

The standard treatment of $ \beta(M) $ is based on the Press-Schechter formalism [100] of the gravitational collapse that is widely adopted in large-scale structure studies [101]. In this formalism, $ \beta(M) $ is given by the probability that the fractional overdensity $ \delta\equiv\delta \rho/\rho $ is above a certain threshold $ \delta_{c} $ for PBH formation [102, 103]. For Gaussian primordial fluctuations, $ \beta(M) $ is given by [104, 105]:

As can be deduced from the above formula, the value of $ \beta(M) $ is uniquely determined by the variance $ \sigma^2\left(M\left(k\right)\right) $, which is assumed to be coarse-grained variance smoothed on a scale of $ R = k^{-1} $. During the radiation dominated era, it is given by the following expression [106]:

where $ P_{S}(q) $ is the power spectrum of curvature perturbation and $ W(x) $ is the smoothing window function [107, 108]. In this study, we use the Gaussian form: $ W(x) = $$ \exp(-x^2/2) $. Moreover, in general, in the radiation dominated period, the threshold density $ \delta _c = c_{s}^{2} = 1/3\approx 0.33 $. However, some studies have shown that the threshold density $ \delta _c $ can be between 0.33 and 0.66 [109, 110].

In Fig. 4, we numerically plot the abundances of $ f_{\rm PBH} $ with different threshold energies $ \delta_c $, taking a solo-bumpy potential with $ p = 2 $ as an example. Values of $ \delta_c $ are as presented in Table 1. We find that the larger the $ \delta_c $ value, the smaller the abundance, and this is easy to understand: the larger the threshold energy, the more difficult it is to form black holes. For a small value of $ \delta_c = 0.41 $, the abundance of PBHs can reach approximately 10% of dark matter. The mass range of the PBHs formed is approximately $ 10^{-15}M_\odot $, namely the asteroid mass.

Figure 4.  (color online) We plot $ f_{\rm PBH} $ for the potential (21) with $ p = 2 $ using different threshold densities $ \delta _c $, where the yellow line corresponds to $ \delta _c = 0.41 $, the purple line corresponds to $ \delta _c = 0.46 $, and the blue line corresponds to $ \delta _c = 0.486 $. Our results are consistent with the constraints from current observations [111]. The specific parameters are presentd in Table 1. Constraints are obtained from the publicly available Python code PBHbounds [112].

$k_{\rm form}/{\rm Mpc}^{-1}$	$M_{\rm PBH}$	$ \delta _c $	$f_{\rm PBH}$
$ 1.1539\times 10^{14} $	$ 3.68\times 10^{-16}M_{\odot} $	0.41	$ 8.0455\times 10^{-2} $
$ 1.1539\times 10^{14} $	$ 3.68\times 10^{-16}M_{\odot} $	0.46	$ 4.6573\times 10^{-6} $
$ 1.1539\times 10^{14} $	$ 3.68\times 10^{-16}M_{\odot} $	0.486	$ 1.8905\times 10^{-8} $

Table 1.  $ f_{\rm PBH} $ for potential (21) with $ p = 2 $ for different threshold densities $ \delta _c. $

Moreover, we also considered the case where the PBHs formed in the early universe were in an ellipsoidal collapse instead of spherical collapse [113, 114]. The difference between these two collapse models is that the threshold density for forming PBHs is different. Compared with the spherical collapse, the PBHs formed by the ellipsoidal collapse will increase the ellipticity of the formed PBHs, which will lead to the correction of the threshold density. In the following, we uniformly use PBHs/e-PBHs to represent the primordial black hole formed by the spherical/ellipsoidal collapse and $ \delta _c $/$ \delta_{ec} $ to describe the threshold density of PBHs/e-PBHs, respectively. According to [115, 116], we have the following relationship:

where $ \kappa = 9/\sqrt{10\pi} $ and $ \nu = 1/2 $. From Eq. (34), we can easily obtain the abundance of PBHs and e-PBHs with a given $ \delta _c $. The calculation shows that the abundance of e-PBHs is lower than that of PBHs, owing to the difference in their threshold densities.

In Fig. 5, we plot the abundances of $ f_{\rm PBH} $ and $ f_{\rm e-PBH} $ from solo-bumpy potential (21), with parameter sets {$ p = 2 $, $ \delta _c = 0.41 $} and {$ p = 2/3 $, $ \delta _c = 0.33 $}, respectively, as well as the constraints from various observations. For the $ p = 2 $ case, the fraction of PBHs could reach O(10%) at the mass range of $ \sim 10^{-15}M_\odot $ (asteroid mass), while for $ p = 2/3 $ case, the fraction is also nearly O(0.1%) at the mass range of $ \sim 1M_\odot $ (solar mass). The specific parameters are presented in Table 2 and Table 3.

Figure 5.  (color online) The figures above show the constraints on primordial black holes acting as dark matter, in which the colored region is excluded by various observations. The blue lines correspond to $ f_{\rm PBH} $, and the red lines correspond to $ f_{\rm e-PBH} $. The left figure is for potential (21) with $ p = 2 $ and $ \delta _c = 0.41 $, while the right figure is for potential (21) with $ p = 2/3 $ and $ \delta _c = 0.33 $. Constraints are obtained from the publicly available Python code PBHbounds [112].

$k_{\rm form}/{\rm Mpc}^{-1}$	$M_{\rm PBH}$	$ \delta _c $	$f_{\rm PBH}$	$f_{\rm e-PBH}$
$ 1.1539\times 10^{14} $	$ 3.68\times 10^{-16}M_{\odot} $	0.41	$ 8.0455\times 10^{-2} $	$ 1.1781\times 10^{-8} $

Table 2.  $ f_{\rm PBH} $ and $ f_{\rm e-PBH} $ for potential (21) with $ p=2 $ and $\delta _c=0.41.$

$k_{\rm form}/{\rm Mpc}^{-1}$	$M_{\rm PBH}$	$ \delta _c $	$f_{\rm PBH}$	$f_{\rm e-PBH}$
$ 6.0256\times 10^{5} $	$ 9.2437M_{\odot} $	0.33	$ 7.2256\times 10^{-4} $	$ 2.9602\times 10^{-8} $

Table 3.  $ f_{\rm PBH} $ and $ f_{\rm e-PBH} $ for potential (21) with $ p=2/3 $ and $\delta _c=0.33.$

Similarly, we plot the abundances of $ f_{\rm PBH} $ and $ f_{\rm e-PBH} $ from the multi-bumpy potential (26) (with parameter sets (27)) in Fig. 6. As can be deduced from the figure, PBHs can be formed in three mass ranges. The leftmost peak in Fig. 6 can have a large amount of PBHs, the abundance of which can reach 27% of dark matter; however, because the mass range of PBHs is negligible (around $ 3.6975\times $$ 10^{-27}M_{\odot} $), if one considers the effects of Hawking radiation, such PBHs may not have any chance to live till today. Nonetheless, it may have some non-trivial effects in the evolution of the Universe, which we will mention in the next section. However, the other two peaks refer to PBHs of $ 10^{-15}M_\odot $ (asteroid mass) and $ 10^{-3}M_\odot $ (planet mass), which can live till today and may act as dark matter. The specific values are presented in Table 4.

Figure 6.  (color online) The figures above present the constraints on primordial black holes acting as dark matter, in which the colored region is excluded by various observations. The plot is for potential (26) and $ \delta _c = 0.465 $. Potential (26) can form three mass ranges of PBHs. Because the PBHs in one of the mass ranges is negligible, it cannot be constrained by existing observations; hence, we plot it separately and use $ initial-f_{\rm PBH} $ to represent it (left figure). The other two mass ranges of PBHs can still be constrained by existing observations (right figure). The blue lines correspond to $ f_{\rm PBH} $ while the red lines correspond to $ f_{\rm e-PBH} $. From left to right, the masses of PBHs are $ 3.6975\times 10^{-27}M_{\odot} $, $ 5.8601\times 10^{-16}M_{\odot} $ and $ 3.6975\times 10^{-3}M_{\odot} $ respectively. Constraints are obtained from the publicly available Python code PBHbounds [112].

$k_{\rm form}/{\rm Mpc}^{-1}$	$M_{\rm PBH}$	$ \delta _c $	$f_{\rm PBH}$	$f_{\rm e-PBH}$
$ 3.1548\times 10^{19} $	$ 3.6975\times 10^{-27}M_{\odot} $	0.465	$ 2.7181\times 10^{-1} $	$ 4.0947\times 10^{-10} $
$ 7.9245\times 10^{13} $	$ 5.8601\times 10^{-16}M_{\odot} $	0.465	$ 9.0930\times 10^{-7} $	$ 1.4448\times 10^{-15} $
$ 3.1548\times 10^{7} $	$ 3.6975\times 10^{-3}M_{\odot} $	0.465	$ 4.9587\times 10^{-13} $	$ 8.3559\times 10^{-22} $

Table 4.  $ f_{\rm PBH} $ and $ f_{\rm e-PBH} $ for potential (26) and $ \delta _c=0.465. $

Scalar-induced gravitational waves (SIGWs) are generated by linear scalar perturbations of the second-order [82, 83]. After inflation, the universe should go through the reheating epoch and will be dominated by radiation. The energy density of the SIGWs in the radiation-dominated (RD) era can be estimated by [27, 117-121]

where the kernel function in the RD era takes the form [122]

where Θ denotes the Heaviside theta function.

The current energy density of the induced GWs can be related to their value in the RD era [123]

where $ \Omega_{r,0}h^{2}\approx4.2\times10^{-5} $ is the current energy density parameter of radiation and $ g_{c}\approx106.75 $ is the effective degrees of freedom of the total radiation energy density at time $ \tau_{c} $ [124]. We convert to frequency f from the comoving wavenumber k via the relationship $ k = 2\pi f $, namely

thus, we can plot the energy spectra of the SIGWs against its frequency f in Fig. 7. We consider the SIGWs sourced from the multi-bumpy potential. Figure 7 presents that three peaks of SIGWs that be generated correspondingly by the three bumps of the inflation potential. Interestingly, the curve of the first peak at $ f\sim [10^{-9},10^{-8}] $ Hz lies in the allowed region of the observational constraints by the NANOGrav results [125].

Figure 7.  (color online) The energy spectra of scalar-induced GWs from multi-bumpy potential. The dashed region represents observational constraints by the NANOGrav data [125]. The colored lines are designed sensitivity curves of current and future gravitational wave detectors including FAST [126], SKA [127], EPTA [128], TianQin [129], LISA [130], Taiji [131], NEMO [132], Einstein Telescope(ET) [133], Cosmic Explorer (CE) [134], LIGO Design [135].

In this study, we investigated the formation of PBHs from an inflation model with bumpy features on its potential. Although the case of a potential with one bump has been discussed in the literature, we discussed the extended case of multiple bumps, which can generate PBHs at different mass ranges. This will enlarge the detectability of PBHs from the coming development of astronomical observations. Specifically, we considered the power-law potential as the basic potential and added one or several bumps of Gaussian type, which makes the inflation roll from the slow-roll stage to the USR-like stage. With the potential, we constructed the power spectrum with single or multiple peaks in small scales, while keeping the large-scale power spectrum consistent with CMB data. We numerically calculated the abundances of PBHs (fraction to dark matter) at the mass range given by the solo-bumpy potential, as well as three mass ranges given by the multi-bumpy potential. Owing to these peaks, PBHs can be formed at different mass ranges, including the asteroid mass range ($ 10^{-16}-10^{-14}M_\odot $), planet mass range ($ 10^{-6}-10^{-3}M_\odot $), and solar mass range (around $ 1M_\odot $), some of which can reach significant abundance. We also demonstrated that the abundance has a negative correlation with the threshold energy density $ \delta_c $. We confronted the results to the current observational data, and all the abundances were determined to be consistent with the data constraints. Although we take three bumps for simplicity, in principle, from such kind of potential, abundances at other mass ranges are also expected.

Moreover, PBHs formed in the early universe can also be formed by ellipsoidal collapse. Therefore, in addition to considering the formation of spherical collapsed PBHs, we also considered the formation of PBHs under ellipsoidal collapse, owing to the correction brought by the ellipticity of the latter; hence, according to the Press-Schecther formalism, the abundance will be lower. We obtained the abundances of PBHs in the two types of collapse at all the mass ranges discussed above, and determined the numerical results to be consistent with the analytical analysis.

Considering the generation of scalar-induced gravitational waves, we have drawn the SIGWs generated by three bumps of the inflation potential. Interestingly, the curve of the first peak lies within the area allowed by the NANOGrav observation constraint.

Before ending, we briefly discuss the evaporation of the formed PBHs. According to Hawking's theory of black-hole radiation, these black holes will inevitably emit virtual particles and continuously reduce their masses until they vanish, which is known as Hawking evaporation [136]. It has been pointed out that, considering the age of the universe, PBHs with initial mass less than $ 10^{15}\text{g} $ $ (\sim 10^{-18}\text{M}_\odot) $ have been completely evaporated today. Therefore, for the multi-bumpy inflation model in our analysis, the PBHs standing for the leftmost peak in Fig. 6 will actually vanish, and this cannot explain the dark matter today. However, it does not imply that they are not important at all. Actually, they may still have a significant impact on the early universe, such as the process of Big-Bang nucleosynthesis, reheating, baryogenesis, etc. [136-139]. With the further development of observational techniques, we may be able to detect the traces left by the PBHs of this type, to find more evidence of their existence. Meanwhile, for other mass ranges, the PBHs are hardly evaporated to date, and can therefore act as dark matter. More details about the influence of the PBHs in our model on these other aspects of our universe, although very interesting, are postponed to an upcoming paper.

T. Q. thanks Yun-Gui Gong and Shi Pi for useful discussions.