Inflation and Primordial Black Holes

Abstract

We review conceptual aspects of inflationary scenarios able to produce primordial black holes by amplifying the size of curvature fluctuations to the level required to trigger black hole formation. We identify general mechanisms to do so, both for single- and multiple-field inflation. In single-field inflation, the spectrum of curvature fluctuations is enhanced by pronounced gradients of background quantities controlling the cosmological dynamics, which can induce brief phases of non-slow-roll inflationary evolution. In multiple-field inflation, the amplification occurs through appropriate couplings with additional sectors characterized by tachyonic instabilities that enhance the size of their fluctuations. As representative examples, we consider axion inflation and two-field models of inflation with rapid turns in field space. We develop our discussion in a pedagogical manner by including some of the most relevant calculations and by guiding the reader through the existing theoretical literature, emphasizing general themes common to several models.

1. Introduction

• Primordial black holes: history of the concept

Inflation, a short period of accelerated expansion in the very early moments of the universe, has become one of the main pillars of modern cosmology [1,2]. Leaving aside its success in addressing the puzzles of the standard hot Big Bang cosmology, inflation provides an explanation for the quantum mechanical origin of structures such as galaxies (including our own) and the anisotropies in the cosmic microwave background (CMB) radiation [3]. In the last two decades, advances in observational cosmology and, in particular, observations of the CMB and of the large-scale structure (LSS) of our universe have confirmed the predictions of inflation and arguably established its status as the main theoretical framework describing the very early universe [4,5]. These successes notwithstanding, CMB and LSS probes only provide us with information on the early universe at the largest cosmological scales (${10}^{-4}\lesssim k\left[{\mathrm{Mpc}}^{-1}\right]\lesssim {10}^{-1}$) corresponding to a small fraction of the early stages of inflationary dynamics. Hence, while inflation provides us with a consistent, testable framework to understand the initial conditions in the universe at the largest scales, we do not have direct access to most of the inflationary dynamics or to the evolution of the universe in the early post-inflationary era. Importantly, these stages could be host to a number of interesting phenomena, including the production of stable relics such as dark matter (see, e.g., [6] for a historical review of dark matter), which is essential in understanding the world we observe today, as well as for establishing new physics. Indeed, the existence of non-luminous, cold dark matter (CDM) that constitutes a quarter of the total energy budget in the universe [7] is one of the most glaring pieces of evidence beyond Standard Model physics [8]. The absence of signatures from collider experiments, along with unsuccessful direct and indirect detection searches, have all made the DM puzzle particularly compelling [9].

An intriguing and economical explanation that might account for DM density in our universe is a scenario wherein DM is made of compact objects such as primordial black holes (PBHs). Pioneered by the works of Y. Zel’dovic and I. Novikov [10] and S. Hawking [11], the initial ideas in this direction began with the realization that PBHs could form as a result of the gravitational collapse of over-dense inhomogeneities in the early universe. In the mid 1970s, it was realized in the works of B. Carr [12,13] and G. Chapline [14] that PBHs could contribute to DM density and provide the seeds for the supermassive BHs populating our universe [15]. Following these theoretical progressions, the interest of the scientific community in PBHs increased in the mid 1990s with the reported detection of microlensing events from the MACHO collaboration [16]. An immediate interpretation of these results suggested the possibility that a significant fraction of mass density in our galaxy could be composed of subsolar-mass PBHs. However, these considerations were later rendered invalid by the findings of the EROS [17] and OGLE [18,19] collaborations, concluding that only a small fraction of mass in the Milky Way could be in the form of PBHs.

Stimulated both by the absence of signals for well-motivated particle DM candidates and the first detection of gravitational waves (GWs) from merging BHs by the the LIGO/VIRGO collaboration [20], a second surge of interest in PBHs was ignited (see Figure 1). In particular, different groups suggested that merging PBHs could be responsible for the observed GW signals, constituting a significant fraction of DM density in our universe [21,22,23]. Since the first appearance of these articles, a significant amount of effort has been pushed forward by the community to search for and constrain the abundance of PBHs by utilizing their gravitational and electromagnetic effects on the environment at small scales. Various experiments set stringent constraints on PBH abundance for the solar- and subsolar-mass range, leaving a viable window for this scenario of tiny PBH masses (${10}^{-17}\lesssim M_{\mathrm{pbh}}\left[M_{\odot }\right]\lesssim {10}^{-12}$ ($M_{\odot }\simeq 1.98\times {10}^{33}\mathrm{gr}$)) as the totality of DM (see, e.g., [24,25,26]).

It should be noted that some of the constraints derived in the literature make specific assumptions about the formation process and the subsequent evolution of PBHs (such as monochromatic mass functions, clustering and accretion processes, etc.) and about other model-dependent specifics (such as non-Gaussianity) that could relax or tighten these constraints1. Since we mostly focus on the subject of inflationary model building, we will not review these issues and aforementioned constraints, but the interested reader can find more details in excellent reviews published recently, see e.g., [28,29,30,31,32,33,34].

PBHs are likely to form well before the end of the radiation-dominated era (i.e., before the so-called matter–radiation equality) and behave like cold and collision-less matter. Therefore, they constitute an interesting DM candidate if they are massive enough ($M_{\mathrm{pbh}}\gtrsim {10}^{15}\mathrm{g}\simeq {10}^{-18}{M}_{\odot }$) to ensure a lifetime comparable with the age of the universe [35]2. In this context, a particularly appealing aspect of PBH dark matter is its economical and minimal structure in the sense that this scenario does not require any additional beyond-Standard-Model (BSM) physics (such as new particles and interactions), provided that one alters the not-so-well-constrained early universe at small scales by introducing a viable mechanism to account for the production of large density fluctuations required for PBH formation.

Similar to the generation of CMB anisotropies, a compelling and natural source of these perturbations in the early universe could be the quantum fluctuations that are stretched outside the horizon during inflation. However, in order to generate such over-dense regions that can collapse to form PBHs in the post-inflationary universe, one needs to devise a mechanism to enhance, by several orders of magnitude, the inflationary scalar perturbations at small scales ($k\gg k_{\mathrm{cmb}}$) (corresponding to late stages of inflation) far above the values required to match CMB observations. As the observed temperature anisotropies prefer a red-tilted power spectrum at CMB scales, this situation generally requires a blue-tilted power spectrum or some specific features at scales associated with PBH formation.

In the context of canonical single scalar field inflation, the first ideas in this direction appeared in works by P. Ivanov, P. Naselsky and I. Novikov [38] (see also [39]). In particular, these authors showed that if the inflation potential has a very flat plateau-like region for field ranges corresponding to the late stages of accelerated expansion, the inflationary dynamics enters a “non-attractor” regime called ultra-slow-roll (USR) inflation [40]. This leads to super-horizon growth of scalar perturbations [41,42,43] that can eventually trigger PBH formation in the post-inflationary universe3. Many explicit single-field inflationary models that exhibit similar local features were subsequently studied in the literature; for a partial list of popular works, see, e.g., [46,47,48,49,50,51,52,53,54,55,56] (see also [57,58,59,60,61,62,63] for earlier constructions). In the context of single scalar field inflation, another possibility to generate an enhancement in the scalar power spectrum is to invoke a variation of the sound speed of scalar fluctuations, for example, through a reduction in the speed of sound ($c_s^2$) [52,64,65] or through rapid oscillation ($c_s^2$), which triggers a resonant instability in the scalar sector [66,67].

From a top-down model-building perspective, a rich particle content during inflation is not just an interesting possibility, but appears to be a common outcome of many BSM theories [68]. Since the early days of research on PBHs, multifield inflationary scenarios have also attracted considerable attention as a natural way to realize enhancement in the scalar perturbations at small scales. For instance, large scalar perturbations may arise through instabilities arising in the scalar sector, e.g., during the waterfall phase of hybrid inflation [57,69,70] or due to turning trajectories in the multiscalar inflationary landscape, as reported recently in [71,72,73,74,75,76,77,78]. Another intriguing possibility in this context is to employ axion-gauge field dynamics during inflation [79,80,81,82,83,84,85,86]. In these models, particle production in the gauge field sector act as a source for the scalar fluctuations and can therefore be responsible for PBH formation.

A common feature of all inflationary scenarios capable of producing PBH populations is the inevitable production of a stochastic GW background (SGWB) induced through higher-order gravitational interactions between enhanced scalar and tensor fluctuations of the metric [87,88,89]. Interestingly, this signal may carry crucial information about the properties of its sources, including the amplitude, statistics and spectral shape of scalar perturbations (see, e.g., [90,91,92,93,94,95,96,97,98]) and could provide invaluable information about the underlying inflationary production mechanism. Furthermore, since the resulting GW background interacts very weakly with the intervening matter between the time of its formation and today, it leads to a rather clean probe of the underlying PBH formation scenario. This allows us to access inflationary dynamics on scales much smaller than those currently probed with CMB and LSS experiments through space- and ground-based GW interferometers including the Laser Interferometer Space Antenna (LISA) [99,100], Pulsar Timing Array (PTA) experiments [101,102] and DECIGO [103,104]. For a detailed review of induced SGWB and the dependence of its properties on the post-inflationary expansion history, see [105].

• The structure of this review

If their origin is attributed to the large primordial fluctuations, PBHs may offer us a unique window to probe inflationary dynamics at sub-CMB scales. In this work, focusing mainly on the activity in the literature within the last few years, we aim to revisit and review different inflationary production mechanisms of PBHs4 and their main predictions, in a heuristic and pedagogical manner. The audience we have in mind is graduate students or researchers in related fields who wish to learn about inflation and primordial black holes and to be guided through the large literature on the subject by emphasizing common conceptual themes behind many different realizations.

This review is organized as follows. In Section 2, we present a simplified, intuitive picture of PBH formation in the inflationary universe and provide an approximate estimate of the required conditions to produce PBHs from the perspective of inflationary dynamics. In Section 3, we discuss ideas to enhance the curvature power spectrum within single-field inflation, as required for PBH formation. These mechanisms exploit large gradients in background quantities that are converted into an amplification of fluctuations. Besides reviewing analytic findings, we also develop numerical analysis and provide a link to a code for reproducing our results (see Notes 21). In Section 4, we focus on multifield inflationary scenarios that can generate PBH populations including particle production during axion inflation or sudden turns in the multiscalar inflationary landscape. Finally, we end with a discussion on future directions in Section 5. We supplement this work with several technical appendices, where we provide useful formulas and calculations used in the main body of this work.

• Conventions

Throughout this review, we work with natural units ($\hslash =c=1$). We use the reduced Planck mass defined as $M_{\mathrm{pl}}^2={\left(8\pi G\right)}^{-1}$ and retain it in the equations unless otherwise stated. For time-dependent quantities, over dots and primes denote derivatives with respect to cosmological time (t) and conformal time ($\mathrm{d}\tau \equiv \mathrm{d}t/a\left(t\right)$), respectively, where $a\left(t\right)$ is the scale factor of the background FLRW metric ($g_{\mu \nu }=\mathrm{diag}(-1,a^2,a^2,a^2)$).

2. PBH Formation in the Early Universe

We start by providing a physical description of PBH formation in the early universe, as comprised of an early stage of inflation, followed by radiation and matter domination (for a mini review of background cosmology, see Appendix A). Our aim is to set the stage and relate basic properties of a PBH population—as their mass and abundance—with the features of primordial curvature fluctuations originating from inflation. For this purpose, in Section 2.1, we describe the mechanism of PBH formation in the post-inflationary universe, emphasizing its nature as a causal process controlled by inflationary quantum fluctuations. In Section 2.2, we discuss relevant concepts such as the threshold for collapse into black holes and the corresponding mass and collapse fraction of PBHs relevant for a computation of their abundance. Finally, in Section 2.3, we relate the PBH abundance to primordial physics during inflation with the aim of determining the amplitude of scalar fluctuations required to produce a population of PBHs with interesting consequences for cosmology. All the concepts we discuss form the basis and motivations for our analysis of inflationary mechanisms for PBH production, which we develop in Section 3 and Section 4.

• Main References: In compiling the materials for this section and to set the main framework for our discussion, we benefited from the ideas presented in the reviews by C. Byrnes and P. Cole [126] and M. Sasaki et al. [30], as well as the Ph.D. thesis of G. Franciolini [127].

2.1. PBH Formation as a Causal Process

An important concept in an expanding space–time is the horizon scale, which is crucial for understanding the causal properties of the dynamics of perturbations that are responsible for PBH formation. As an indicator of the rate of expansion of our universe, the Hubble rate ($H\left(t\right)\equiv \dot{a}\left(t\right)/a\left(t\right)$) has dimensions of inverse length (or time${}^{-1}$ in natural units). This makes the quantity $H^{-1}$ (Hubble horizon) a natural candidate for a physical length scale in an expanding universe. Commonly referred to as the Hubble distance, the quantity $1/H$ (or $c/H$ if one wishes to recover physical units) measures the distance that light travels within one unit of Hubble time. Therefore, it can be considered a good proxy for a (time-dependent) length scale controlling the size of a causal patch in our universe. Bearing in mind that we relate physical quantities to comoving quantities by the scale factor ($a\left(t\right)$), a useful quantity that guides us in this direction is the comoving Hubble horizon (${\left(aH\right)}^{-1}$) and, in particular, its time evolution. When expressed in terms of the second derivative of the scale factor, the time derivative of the comoving horizon can be written as

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{1}{aH}\right)=-\frac{\ddot{a}}{a^2{H}^2}.$$

(1)

Notice that during inflation, $\ddot{a}>0$; hence, the comoving horizon scale is a decreasing function of time, whereas in a decelerating universe with $\ddot{a}<0$ (i.e., in the post-inflationary universe before dark energy domination), this quantity is an increasing function of time. The property that the comoving horizon decreases during an accelerated expansion is perhaps the most important element to understand inflation as a solution of the horizon problem of hot Big Bang cosmology and a framework for the quantum mechanical origin of structures in our universe5. The time dependence of the comoving Hubble horizon is controlled by the value of the background equation of state (w) (EoS) as (see Appendix A)

$${\left(aH\right)}^{-1}\propto a^{(1+3w)/2}.$$

(2)

Therefore, during inflation, $w\simeq -1$ and ${\left(aH\right)}^{-1}\propto a^{-1}$, while during the subsequent phases of the radiation-dominated (RDU) and matter-dominated universe (MDU), the comoving horizon evolves as ${\left(aH\right)}^{-1}\propto a^1$ and ${\left(aH\right)}^{-1}\propto a^{1/2}$, respectively. The evolution of the comoving horizon with respect to the logarithm-of-scale factor ($\ln \left(a\right)$) is illustrated in Figure 2.

When we study the statistical properties of fluctuations in Fourier space, we often label a given perturbation mode with a comoving length scale ($k^{-1}$) measured in units of megaparsecs ($\mathrm{Mpc}=3.26\times {10}^6\mathrm{light}\mathrm{years}\simeq 3.1\times {10}^{19}\mathrm{km}$). Therefore, a crucial quantity to conceptualize the behavior of perturbations in the inflationary universe is the ratio of the wavelength of a given mode with respect to the size of the comoving Hubble horizon, i.e., $\left(aH\right)/k$. Fluctuations with wavelengths larger than the comoving horizon are referred to as superhorizon modes ($k<aH$), while subhorizon perturbations satisfy $k>aH$. Each mode crosses the horizon at $k=aH$. As shown in Figure 2, a typical fluctuation with a comoving size of $k^{-1}$ (horizontal lines) begins its life deep inside the horizon (typically as a quantum fluctuation); then, it leaves the horizon to become a superhorizon mode, and, finally, it re-enters the comoving horizon in the post-inflationary universe. Large-scale modes (with smaller k) exit the horizon earlier than small-scale modes and re-enter the horizon at a later time in the post-inflationary era. For definiteness, in Figure 2, we represent the behavior of the comoving curvature perturbation ($\mathcal{R}$) [129,130] (see [128] for a textbook discussion), which plays an important role in our discussion.

In fact, apart from providing seeds for the observed cosmic microwave background (CMB) anisotropies at large scales, the dynamics of curvature fluctuations ($\mathcal{R}$) may also be at play for PBH formation, provided that cosmological fluctuations exhibit specific ‘initial conditions’ at small scales. For this purpose, denoting $k=k_{\mathrm{pbh}}\gg k_{\mathrm{cmb}}$ as the comoving momentum associated with PBH formation, we assume that the curvature power spectrum at these small scales acquires an amplification well above the level required to match CMB observations: ${\mathcal{P}}_{\mathcal{R}}\left(k_{\mathrm{pbh}}\right)\gg {\mathcal{P}}_{\mathcal{R}}\left(k_{\mathrm{cmb}}\right)\sim {10}^{-9}$[4] (more on this later). Soon after the end of inflation, i.e., after reheating6, the modes associated with the enhancement (e.g., modes with a comoving size of $k_{\mathrm{pbh}}^{-1}$) become the seeds of density perturbations in the RDU:

$${\mathcal{P}}_{\mathcal{R}}{\left(k_{\mathrm{pbh}}\right)}^{1/2}\sim \frac{\delta \rho }{\rho }.$$

(3)

Since the comoving Hubble scale grows with respect to comoving scales in the RDU, the characteristic scale of perturbations eventually becomes comparable to the comoving horizon at $a=a_{\mathrm{f}}$, where $k_{\mathrm{pbh}}={\left(aH\right)}_{\mathrm{f}}$ (the subscript $\mathrm{f}$ indicates PBH formation time). At this point, gravitational interactions can trigger the collapse of over-dense regions if the latter have a sufficient over-density above a certain collapse threshold (see the discussion below).

Notice that at subhorizon scales, the radiation pressure can overcome the gravitational collapse. Therefore, the production of PBHs effectively occurs around horizon re-entry. This implies that the concept of horizon re-entry is crucial for our understanding of PBH formation as a causal process. In fact, only when the physical wavelength of a perturbation becomes comparable to the causal distance ($1/H$) is gravity able to communicate the presence of an over-density and to initiate gravitational collapse. A schematic diagram that summarizes the discussion presented above is illustrated in Figure 2.

In what comes next, we briefly discuss relevant quantities such as the threshold for collapse and the mass and collapse fractions of PBHs, which are important for the computation of the current PBH abundance.

2.2. The Relevant Quantities for PBH Abundance

• The threshold for collapse

The first criterion of the collapse threshold—defined as ${\delta }_{\mathrm{c}}$ hereafter—for PBH formation was formulated by B. Carr in 1975 using a Jeans-type instability argument within Newtonian gravity [13]. In Carr’s estimate, an over-density in RDU would collapse upon horizon re-entry if the fractional over-density of the perturbation is larger than the sound speed square of density perturbations ($c_s^2$), which is a measure of how fast a pressure wave caused by the over-density can travel from the center to the edge of a local fluctuation. In the RDU, the speed of sound of perturbations satisfies $c_s^2=w=1/3$ so that its square is directly related to EoS. This result implies that a perturbation can collapse to form PBHs if its over-density is larger than the pressure exerted by the radiation pressure7. Implementing general relativistic effects, the second analytical estimate was derived in [132], obtaining ${\delta }_{\mathrm{c}}\simeq 0.4$ during the RDU. These results notwithstanding, further studies in this direction have found that the threshold depends on the initial profile of curvature perturbation [133,134,135,136]8.

At this point, we emphasize that a more precise criterion for the threshold of collapse is provided in terms of the so-called compaction function (which is a measure of excess mass within a given spatial volume; see, e.g., [135,139]) the maximum value of which results in ${\delta }_c$. In this prescription, we note that contrary to the simplistic argument presented in Appendix B, there is no one-to-one correspondence between the threshold value and the peak over-density ((3) $\delta \equiv \delta \rho /\rho$), although they are related; see, e.g., [140]. Keeping this in mind, in the rest of this work, we will continue to refer to the threshold using the ${\delta }_{\mathrm{c}}$ notation as commonly practiced in the literature.

An accurate characterization of the threshold requires a dedicated analysis of the evolution of perturbations in the non-linear regime after horizon re-entry, which can be done with the help of numerical simulations (see, e.g., [141] for a review). Utilizing the peak of the criterion of compaction function as the definition of the threshold (${\delta }_{\mathrm{c}}$), simulations performed in the radiation fluid for different primordial perturbation profiles lead to a range of values (${\delta }_{\mathrm{c}}=0.4-2/3$) depending on the shape of the density peak [135,142]. Although in general, threshold depends on the shape of the primordial perturbation, refs. [142,143] showed the existence of an approximate universal value for the threshold that depends only on the type of ambient fluid and geometry. For a detailed account of the threshold and its history, we refer the reader to a recent review [141].

• The mass of PBHs

The characteristic mass of PBHs can be related to the mass contained within the Hubble horizon at the time of formation ($a=a_{\mathrm{f}}$) through an efficiency factor ($\gamma =0.2$), as suggested by the analytical model developed in [13]:

$$M_{\mathrm{pbh}}^{\left(\mathrm{f}\right)}=\gamma M_H{|}_{a=a_{\mathrm{f}}}=\gamma M_H^{\left(\mathrm{eq}\right)}\left(\frac{M_H^{\left(\mathrm{f}\right)}}{M_H^{\left(\mathrm{eq}\right)}}\right)={\left(\frac{a_{\mathrm{f}}}{a_{\mathrm{eq}}}\right)}^2\gamma M_H^{\left(\mathrm{eq}\right)}.$$

(4)

In this formula, $M_H\left(t\right)\equiv 4\pi \rho \left(t\right)/\left(3H{\left(t\right)}^3\right)$ is the time-dependent horizon mass, where the sub/superscripts “f” and “eq” denote quantities evaluated at the time of PBH formation and matter–radiation equality, respectively; we use the standard relations $H^2\propto \rho \propto a^{-4}$ during the RDU. Noting that the horizon mass at the time of equality is given by $M_H^{\left(\mathrm{eq}\right)}\simeq 2.8\times {10}^{17}{M}_{\odot }$ [90], Equation (4) informs us that PBHs, contrary to astrophysical black holes, can, in principle, span a wide range of masses, depending on their formation time with respect to matter–radiation equality.

Making use of the time-dependent horizon mass as above, we can relate the PBH mass at formation to the characteristic size of the perturbations that leave the horizon during inflation and are responsible for PBH formation. For this purpose, we first rewrite the PBH mass at formation as

$$M_{\mathrm{pbh}}^{\left(\mathrm{f}\right)}={\left(\frac{{\rho }_{{}_{\mathrm{f}}}}{{\rho }_{{}_{\mathrm{eq}}}}\right)}^{1/2}{\left(\frac{H_{\mathrm{eq}}}{H_{\mathrm{f}}}\right)}^2\gamma M_H^{\left(\mathrm{eq}\right)}.$$

(5)

Using the property of entropy conservation ($g_s\left(T\right)T^3{a}^3=\mathrm{constant}$) and the scaling property of the energy density with respect to the temperature of the plasma during the RDU ($\rho \propto g_{*}\left(T\right)T^4$), Equation (5) can be re-expressed as the following relation

$$\begin{array}{cc}\hfill M_{\mathrm{pbh}}^{\left(\mathrm{f}\right)}\left(k_{\mathrm{pbh}}\right)& ={\left(\frac{g_{*}\left(T_{\mathrm{f}}\right)}{g_{*}\left(T_{\mathrm{eq}}\right)}\right)}^{1/2}{\left(\frac{g_s\left(T_{\mathrm{eq}}\right)}{g_s\left(T_{\mathrm{f}}\right)}\right)}^{2/3}{\left(\frac{k_{\mathrm{eq}}}{k_{\mathrm{pbh}}}\right)}^2\gamma M_H^{\left(\mathrm{eq}\right)},\hfill \\ \hfill & \simeq \left(\frac{\gamma }{0.2}\right){\left(\frac{g_{*}\left(T_{\mathrm{f}}\right)}{106.75}\right)}^{-1/6}{\left(\frac{k_{\mathrm{pbh}}}{3.2\times {10}^5{\mathrm{Mpc}}^{-1}}\right)}^{-2}30M_{\odot }.\hfill \end{array}$$

(6)

In the second line, we assume that the effective number of relativistic degrees of freedom in energy density and entropy are equal, i.e., we set $g_{*}\left(T\right)=g_s\left(T\right)$ and take $g_{*}\left(T_{\mathrm{eq}}\right)\simeq 3.38$9 with $k_{\mathrm{eq}}\simeq 0.0104{\mathrm{Mpc}}^{-1}$, accordingly with the latest Planck results [7]. Equation (6) indicates that for masses of PBHs that can be associated with recent LIGO observations ($M_{\mathrm{pbh}}\simeq 30M_{\odot }$), the peak scale of perturbations responsible for PBH formation is much smaller compared to CMB scales ($k_{\mathrm{pbh}}\gg k_{\mathrm{cmb}}$). For subsolar-mass PBHs, the corresponding peak scale for PBH formation becomes progressively smaller. For example, considering the currently allowed sublunar range ($M_{\mathrm{moon}}\simeq 3.7\times {10}^{-8}{M}_{\odot }$) of PBH masses (${10}^{-17}\lesssim M_{\mathrm{pbh}}\left[M_{\odot }\right]\lesssim {10}^{-12}$), which are objects that can account for the totality of dark matter, the range of scales associated with PBH formation is quite small: ${10}^{12}\lesssim k_{\mathrm{pbh}}\left[{\mathrm{Mpc}}^{-1}\right]\lesssim {10}^{15}$. See

Table 1. Range of PBH mass(es) vs. the corresponding wave number(s) ($k_{\mathrm{pbh}}$) (see Equation (6)) of the primordial modes, together with the approximate horizon crossing time measured with respect to e-folding number of the pivot scale $k_{\mathrm{cmb}}=0.002{\mathrm{Mpc}}^{-1}$ leaves the horizon during inflation, $\Delta N\equiv N_{\mathrm{pbh}}-N_{\mathrm{cmb}}>0$ (see Equation (7)). The first row refers to the corresponding quantities for a typical supermassive black hole (SMBH), such as $\mathrm{Sagittarius}{\mathrm{A}}^{*}$ in the center of our galaxy [144]. The third row refers to asteroid-mass PBHs that can still account for a significant fraction (or all) of DM density today [26]. The corresponding mass of the PBHs in terms of the Earth’s mass ($M_{\odot }\simeq 3.33\times {10}^5{M}_{\oplus }$) and the mass of Mount Everest ($M_{\mathrm{Everest}}=8.1\times {10}^{14}\mathrm{kg}\simeq 4.1\times {10}^{-16}{M}_{\odot }$) are given in the last two columns on the right.

M${}_{\mathbf{pbh}}\left[{\mathit{M}}_{\odot }\right]$	$\mathbf{\Delta }\mathit{N}$	${\mathit{k}}_{\mathbf{pbh}}\left[{\mathbf{Mpc}}^{-1}\right]$	${\mathit{M}}_{\mathbf{pbh}}\left[{\mathit{M}}_{\oplus }\right]$	${\mathit{M}}_{\mathbf{pbh}}\left[{\mathit{M}}_{\mathbf{Everest}}\right]$
${10}^6$	14	${10}^3$	${10}^{11}$	${10}^{21}$
${10}^0–{10}^2$	$18–21$	${10}^5–{10}^6$	${10}^5–{10}^7$	${10}^{15}–{10}^{17}$
${10}^{-17}–{10}^{-12}$	$34–40$	${10}^{12}–{10}^{15}$	${10}^{-12}–{10}^{-7}$	${10}^{-2}–{10}^3$

for an easier-to-visualize summary of these considerations.

Elaborating on Equation (6), we can also derive a rough relation between the PBH mass at formation with the the number of e-folds ($N_{\mathrm{pbh}}$) at which the PBH-forming modes leave the horizon during inflation. For this purpose, we first notice that $k_{\mathrm{pbh}}/k_{\mathrm{cmb}}={\left(aH\right)}_{\mathrm{pbh}}/{\left(aH\right)}_{\mathrm{cmb}}$, where the values of the Hubble rate and scale factor should be evaluated at the scales of horizon exit during inflation (see Figure 2). Assuming a roughly constant slow-roll parameter ($ϵ\equiv -\dot{H}/H^2\ll 1$) between the horizon exit time of modes associated with CMB and PBH formation, respectively, we can relate the Hubble and the scale factor as $H_{\mathrm{pbh}}=H_{\mathrm{cmb}}{e}^{-ϵ(N_{\mathrm{pbh}}-N_{\mathrm{cmb}})}$ and $a_{\mathrm{pbh}}=a_{\mathrm{cmb}}{e}^{N_{\mathrm{pbh}}-N_{\mathrm{cmb}}}$, where $N_{\mathrm{pbh}}>N_{\mathrm{cmb}}$ so that we count e-folds forward in time with respect to the horizon exit of the CMB mode10. Using the last two relations, we find $k_{\mathrm{pbh}}/k_{\mathrm{cmb}}\simeq e^{(N_{\mathrm{pbh}}-N_{\mathrm{cmb}})(1-ϵ)}$; once plugged into Equation (6), assuming $k_{\mathrm{cmb}}=0.002{\mathrm{Mpc}}^{-1}$, we find

$$M_{\mathrm{pbh}}^{\left(\mathrm{f}\right)}\left(N_{\mathrm{pbh}}\right)\approx 7.7\times {10}^{17}{M}_{\odot }{e}^{-2(N_{\mathrm{pbh}}-N_{\mathrm{cmb}})(1-ϵ)}\left(\frac{\gamma }{0.2}\right){\left(\frac{g_{*}\left(T_{\mathrm{f}}\right)}{106.75}\right)}^{-1/6}.$$

(7)

Modes that leave the horizon much later compared to CMB scales have $N_{\mathrm{pbh}}-N_{\mathrm{cmb}}\gg 1$; therefore, the exponential in Equation (7) can considerably reduce the overall large normalization, leading to small PBH masses.

• PBH abundance

After discussing possible masses for PBHs and how they depend on the dynamics of inflation, we analyze the notion of abundance of PBHs relative to the energy density of other species. We can compute this quantity during two epochs: today and at PBH formation.

When considering the present-day fraction of PBH density, it is a common practice to relate the PBH abundance to present-day dark matter density, introducing the quantity

$$f_{\mathrm{pbh}}\equiv \frac{{\Omega }_{\mathrm{pbh}}}{{\Omega }_{\mathrm{dm}}},$$

(8)

where for each species (i), we define ${\Omega }_i\equiv {\rho }_{{}_{i,0}}/{\rho }_{{}_{c,0}}$, with subscript “0” denoting quantities evaluated today and ${\rho }_{{}_{c,0}}=3H_0^2{M}_{\mathrm{pl}}^2$ representing the critical density. Planck measurements provide the following value for the dark matter abundance [7],

$${\Omega }_{\mathrm{dm}}{h}^2=0.120\pm 0.001,$$

(9)

in terms of $h=0.6736\pm 0.0054$, which measures the Hubble rate ($H_0$) in units of $100\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$.

We can then relate $f_{\mathrm{pbh}}$ today to the density fraction of PBHs at the epoch of their formation, denoting this quantity with $\beta$. In fact, since we assume that PBH formation takes place during the RDU and sinc after formation, the PBH density scales as ${\rho }_{\mathrm{pbh}}\propto a^{-3}$, we can write

$$\begin{array}{cc}\hfill \beta \equiv \frac{{\rho }_{{}_{\mathrm{pbh}}}}{\rho }{|}_{a=a_{\mathrm{f}}}=\frac{{\rho }_{{}_{\mathrm{pbh},\mathrm{f}}}}{{\rho }_{{}_{\mathrm{pbh},0}}}\frac{{\rho }_{{}_{c,0}}}{{\rho }_{{}_{\mathrm{f}}}}{\Omega }_{\mathrm{dm}}{f}_{\mathrm{pbh}}& =\frac{{\rho }_{{}_{\mathrm{pbh},\mathrm{f}}}}{{\rho }_{{}_{\mathrm{pbh},0}}}\frac{{\rho }_{{}_{c,0}}}{{\rho }_{{}_{\mathrm{eq}}}}\frac{{\rho }_{{}_{\mathrm{eq}}}}{{\rho }_{{}_{\mathrm{f}}}}{\Omega }_{\mathrm{dm}}{f}_{\mathrm{pbh}},\hfill \\ \hfill & \simeq \frac{1}{2}\frac{a_{\mathrm{f}}}{a_{\mathrm{eq}}}{\Omega }_m^{-1}{\Omega }_{\mathrm{dm}}{f}_{\mathrm{pbh}}\hfill \end{array}$$

(10)

where ${\Omega }_m$ is the current matter density in the universe. In (10), we normalize the scale factor today as $a_0=1$ and use the fact that the total energy density evolves as $\rho \propto a^{-4}$ for $a_{\mathrm{f}}<a<a_{\mathrm{eq}}$ and $\rho \left(a_{\mathrm{eq}}\right)=2{\rho }_{{}_{m,0}}{a}_{\mathrm{eq}}^{-3}$. Using the conservation of total entropy ($g_s\left(T\right)T^3{a}^3=\mathrm{constant}$), we can re-express the factor $a_{\mathrm{f}}/a_{\mathrm{eq}}$ appearing in Equation (10) as follows:

$$\begin{array}{cc}\hfill \frac{a_{\mathrm{f}}}{a_{\mathrm{eq}}}& =\frac{T_{\mathrm{eq}}}{T_{\mathrm{f}}}{\left(\frac{g_s\left(T_{\mathrm{eq}}\right)}{g_s\left(T_{\mathrm{f}}\right)}\right)}^{1/3},\hfill \\ \hfill & \simeq 3.17\times {10}^{-9}{\left(\frac{\gamma }{0.2}\right)}^{-1/2}{\left(\frac{g_{*}\left(T_{\mathrm{f}}\right)}{106.75}\right)}^{-1/12}{\left(\frac{M_{\mathrm{pbh}}^{\left(\mathrm{f}\right)}}{M_{\odot }}\right)}^{1/2},\hfill \end{array}$$

(11)

where we make use of Equation (5) to relate $T_{\mathrm{eq}}/T_{\mathrm{f}}$ to the mass of PBHs at formation, and, as before, we assume $g_{*}\left(T\right)=g_s\left(T\right)$. Finally, plugging (11) into (10) and implementing Planck measurements on ${\Omega }_{\mathrm{dm}}$ (9) and ${\Omega }_m$ (see Equation (A18)), we directly relate the abundance of PBHs at formation ($\beta$) to their present-day fraction ($f_{\mathrm{pbh}}$) in terms of the PBH mass at formation:

$$\beta \simeq 1.33\times {10}^{-9}{\left(\frac{\gamma }{0.2}\right)}^{-1/2}{\left(\frac{g_{*}\left(T_{\mathrm{f}}\right)}{106.75}\right)}^{-1/12}{\left(\frac{M_{\mathrm{pbh}}^{\left(\mathrm{f}\right)}}{M_{\odot }}\right)}^{1/2}{f}_{\mathrm{pbh}}.$$

(12)

Therefore, we learn that in the case in which PBH abundance accounts for the total DM density today ($f_{\mathrm{pbh}}\to 1$), the fraction of the total density in the form of PBHs ($\beta$) at the time of their formation takes extremely small values when considering an interesting range of masses ($M_{\mathrm{pbh}}^{\left(\mathrm{f}\right)}$). This situation reflects the fact that PBH formation in the early universe is a very rare event. There is also another way to parameterize $\beta$ in terms of the (relative) number of collapsing regions to form PBHs. This approach is especially useful to relate the PBH abundance to the statistical properties of primordial fluctuations, as we discuss below.

• Collapse fraction of PBHs at formation

The PBH abundance at formation can also be interpreted as the fraction ($\beta$) of local regions in the universe that has a density larger than a certain threshold. The standard treatment for estimating $\beta$ is based on the so-called Press–Schechter model of gravitational collapse, which is widely used in the literature on large-scale structure formation [145]11,

$$\beta \equiv \frac{{\rho }_{{}_{\mathrm{pbh}}}}{\rho }{|}_{a=a_{\mathrm{f}}}\equiv {\int }_{{\Delta }_{\mathrm{c}}}^{\infty }P\left(\delta \right)\mathrm{d}\delta ,$$

(13)

where $P\left(\delta \right)$ is the probability distribution function (PDF), which describes the likelihood that a given fluctuation has an over-density $\delta$; we assume that a perturbation will collapse to form a BH if its amplitude is larger than a critical value (${\Delta }_c$). Notice that ${\Delta }_c$ is the critical density contrast and, in general, cannot be identified with the threshold ${\Delta }_{\mathrm{c}}\ne {\delta }_{\mathrm{c}}$, where the latter can be defined as the peak value of the compaction function and can be related to an averaged density contrast [141]. Note also that an alternative, more accurate approach to compute the PBH abundance requires the use of the PDF of the compaction function (see, e.g., [149]). A discussion of this approach and the actual relation between ${\Delta }_{\mathrm{c}}$ and ${\delta }_{\mathrm{c}}$ would bring us outside the pedagogical purposes of this review. For a detailed account of the use the PDF of the compaction function and the influence of resulting non-linearities on the PBH abundance, we refer the reader to [149] and references therein. For more details on the relation between ${\Delta }_{\mathrm{c}}$ and ${\delta }_{\mathrm{c}}$, see also the review [141]. Keeping an eye on the recent progress in the literature on these issues, we will continue to identify ${\Delta }_{\mathrm{c}}={\delta }_{\mathrm{c}}$ for estimates of $\beta$ in this and the next section and continue to work with a PDF on over-density $\delta$.

Let us assume that the latter follows a Gaussian distribution,

$$P_G\left(\delta \right)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-{(\delta -\mu )}^2/2{\sigma }^2},$$

(14)

where $\mu$ is the mean and ${\sigma }^2$ is the variance of the distribution. In Figure 3, we represent two Gaussian PDFs that have the same mean value and two different variances satisfying ${\sigma }_2^2>{\sigma }_1^2$. As the second distribution is more “spread” with a larger variance compared to the first one, the probability of an over-density larger than the critical threshold (${\delta }_{\mathrm{c}}$) is larger, as is PBH abundance $\beta$ at formation, since the integral in Equation (13) has more support within the integration limits $\delta \in [{\delta }_{\mathrm{c}},\infty ]$. Hence, we can expect that the quantity ($\sigma$) plays an important role in estimating the PBH abundance.

Using (12) and (13), we can estimate the required variance (${\sigma }^2$) of (14) that can give rise to a large population of PBH today, as controlled by the quantity ($\beta$). Focusing on a distribution with zero mean $\mu =0$ in (14) and integrating (13), we have

$$\beta ={\int }_{{\delta }_{\mathrm{c}}}^{\infty }\frac{\mathrm{d}\delta }{\sqrt{2\pi }\sigma }\exp \left(-\frac{{\delta }^2}{2{\sigma }^2}\right)=\frac{1}{2}\mathrm{Erfc}\left(\frac{{\delta }_{\mathrm{c}}}{\sqrt{2}\sigma }\right)\simeq \frac{\sigma }{\sqrt{2\pi }{\delta }_{\mathrm{c}}}\exp \left(-\frac{{\delta }_{\mathrm{c}}^2}{2{\sigma }^2}\right),$$

(15)

where $\mathrm{Erfc}\left(x\right)=1-\mathrm{Erf}\left(x\right)$ is the complementary error function, and in the last equality, we take ${\delta }_{\mathrm{c}}\gg \sigma$. As we will learn shortly, this is a good approximation for all practical purposes. As concrete examples, substituting Equation (15) into Equation (12), we infer that a solar mass PBH population with $f_{\mathrm{pbh}}={10}^{-3}$ requires ${\delta }_{\mathrm{c}}/\sigma \simeq 7$, whereas for a population with $M_{\mathrm{pbh}}^{\left(\mathrm{f}\right)}={10}^{-12}{M}_{\odot }$ and $f_{\mathrm{pbh}}=1$, we need ${\delta }_{\mathrm{c}}/\sigma \simeq 7.9$. Assuming a threshold of ${\delta }_{\mathrm{c}}=0.4$, these results translate into $\sigma \simeq 0.06$ and $\sigma \simeq 0.05$, respectively. Notice also from (15) that the PBH abundance at the time of formation is exponentially sensitive to the variance of the distribution. We will dwell more on this dependence below. However, first, we discuss the implications of these findings in terms of the amplitude of scalar fluctuations generated during the phase of cosmic inflation.

2.3. Relating PBH Properties to Primordial Scalar Fluctuations

We now examine how to relate the notion of PBH abundance to the properties of the comoving curvature fluctuation ($\mathcal{R}$) [129,130], as produced in the early universe by cosmic inflation. $\mathcal{R}$ is conserved on super Hubble scales as the modes evolve from the inflationary phase to the RDU. We start by connecting the amplitude of $\mathcal{R}$ with the fractional over-density ($\delta$), which triggers PBH formation, as we learned in our previous discussion. Working in Fourier space with Taylor expansion at the leading order in a gradient expansion (controlled by the small parameter $k/aH$) and at a linear order in $\mathcal{R}$, one finds [135]:

$$\delta (\overrightarrow{x},t)\simeq \frac{2(1+w)}{(5+3w)}\frac{{\nabla }^2\mathcal{R}\left(\overrightarrow{x}\right)}{{\left(aH\right)}^2}+\cdots ⟹{\delta }_k\simeq -\frac{4}{9}{\left(\frac{k}{aH}\right)}^2{\mathcal{R}}_k,$$

(16)

where $w=1/3$ is used during the RDU, and … denotes terms of higher order ($\mathcal{O}\left({\mathcal{R}}^2\right)$) in the curvature perturbation12. Defining the power spectrum of a Fourier variable ($X_k$) as

$$\langle X_k{X}_{k^{\prime }}\rangle =\frac{2{\pi }^2}{k^3}{\mathcal{P}}_X\left(k\right){\delta }^{\left(3\right)}(\overrightarrow{k}+{\overrightarrow{k}}^{\prime }),$$

(17)

the relation between the power spectrum of over-density and curvature perturbation is then given by

$${\mathcal{P}}_{\delta }\left(k\right)\simeq \frac{16}{81}{\left(\frac{k}{aH}\right)}^4{\mathcal{P}}_{\mathcal{R}}\left(k\right).$$

(18)

In the computation of the density contrast, one should typically use a window function ($\mathcal{W}$) to smooth $\delta (\overrightarrow{x},t)$ on a scale ($R\approx k^{-1}\approx {\left(aH\right)}^{-1}$) (e.g., on scales of size $k_{\mathrm{pbh}}^{-1}$ at horizon re-entry, as shown in Figure 2) relevant to PBH formation. Therefore, the variance of density contrast can be related to the primordial power spectrum as [30,152]13

$${\sigma }^2\left(R\right)\equiv {\langle {\delta }^2\rangle }_R={\int }_0^{\infty }\mathrm{d}\ln q{\mathcal{W}}^2(q,R){\mathcal{P}}_{\delta }\left(q\right)\simeq \frac{16}{81}{\int }_0^{\infty }\mathrm{d}\ln q{\mathcal{W}}^2(q,R){\left(qR\right)}^4{\mathcal{P}}_{\mathcal{R}}\left(q\right)$$

(19)

where $\mathcal{W}$ is the Fourier transform of a real space-window function. Popular choices for $\mathcal{W}$ include a volume-normalized Gaussian or a top-hat window function, the Fourier transforms of are respectively given by

$$\mathcal{W}(k,R)=\exp \left(-\frac{k^2{R}^2}{2}\right),\mathcal{W}(k,R)=\frac{3\sin \left(kR\right)-3kR\cos \left(kR\right)}{{\left(kR\right)}^3}.$$

(20)

When selecting a curvature power spectrum (${\mathcal{P}}_{\delta }$) characterized by a narrow peak around the wave number $k_{\mathrm{pbh}}$, the integral in (19) can be approximated as ${\sigma }^2\sim {\mathcal{P}}_{\delta }\left(k_{\mathrm{pbh}}\right)$. Then, utilizing (18) at horizon entry ($k\simeq aH$) (i.e., at the time of PBH formation), since $81/16\sim 5$, we can roughly relate the variance ($\sigma$) to the primordial curvature power spectrum as

$${\mathcal{P}}_{\mathcal{R}}\left(k_{\mathrm{pbh}}\right)\sim 5{\sigma }^2.$$

(21)

Finally, recall our considerations after Equation (15): a Gaussian PDF of $\delta$ requires $\sigma \simeq 0.06$ ($\sigma \simeq 0.05$) for $M_{\mathrm{pbh}}^{\left(\mathrm{f}\right)}=M_{\odot }$ ($M_{\mathrm{pbh}}^{\left(\mathrm{f}\right)}={10}^{-12}{M}_{\odot }$) to generate a population of $f_{\mathrm{pbh}}={10}^{-3}$ ($f_{\mathrm{pbh}}=1$) today. Hence, we can estimate the amplitude of the scalar power spectrum needed at scales relevant for PBH formation:

$${\mathcal{P}}_{\mathcal{R}}\left(k_{\mathrm{pbh}}\right)\sim 5{\sigma }^2\sim {10}^{-2}\mathrm{for}\mathrm{Gaussian}\mathrm{fluctuations}.$$

(22)

This estimate holds for a wide range of subsolar PBH masses. This implies that we need a very large amplification of the curvature spectrum between large CMB and small PBH formation scales:

$$\Delta {\mathcal{P}}_{\mathcal{R}}\equiv \frac{{\mathcal{P}}_{\mathcal{R}}\left(k_{\mathrm{pbh}}\right)}{{\mathcal{P}}_{\mathcal{R}}\left(k_{\mathrm{cmb}}\right)}\sim {10}^7,$$

(23)

and the task is to produce such amplification in a controllable way by an appropriate inflationary mechanism.

It is worth pointing out that this estimate does not change much for even smaller-mass PBHs with $M_{\mathrm{pbh}}^{\left(\mathrm{f}\right)}<{10}^{-12}{M}_{\odot }$ because the power spectrum has a logarithmic sensitivity to the PBH fraction ($\beta$). In order to see this, we can invert the expression (15) and use (22) to relate the primordial power spectrum of curvature perturbations to $\beta$ as

$${\mathcal{P}}_{\mathcal{R}}\left(k_{\mathrm{pbh}}\right)\sim 5{\sigma }^2\sim \frac{5{\delta }_{\mathrm{c}}^2}{2\ln (1/\beta )}.$$

(24)

Now, as an extreme case, we can consider the smallest PBHs ($M_{\mathrm{pbh}}^{\left(\mathrm{f}\right)}\simeq {10}^{-18}{M}_{\odot }$) that can survive until today (not yet eliminated by Hawking radiation), which have the tightest available observational constraints, restricting their current abundance to $f_{\mathrm{pbh}}\lesssim {10}^{-9}$ [32]. Plugging these values into (12), PBH fraction at formation gives $\beta \lesssim {10}^{-28}$, which, in turn, leads to the constraint ${\mathcal{P}}_{\mathcal{R}}\left(k_{\mathrm{pbh}}\right)\lesssim 6\times {10}^{-3}$ in (24) for a threshold of ${\delta }_{\mathrm{c}}=0.4$. Therefore, we conclude that for Gaussian perturbations and for any PBH mass of interest, the amplitude of the scalar power spectrum relevant for PBH formation requires ${\mathcal{P}}_{\mathcal{R}}\sim {10}^{-2}$ for any potentially observable PBH fraction ($f_{\mathrm{PBH}}$) today. The discussion presented above informs us that a small change in the amplitude of the power spectrum leads to a difference of many orders of magnitude in the fraction of regions collapsing into PBHs, as clearly implied by the exponential dependence of $\beta$ on ${\mathcal{P}}_{\mathcal{R}}$ in (24). Similarly, a small change in the choice of threshold (${\delta }_{\mathrm{c}}$) could lead to very different estimates in terms of $\beta$. For example, focusing on a fixed value of variance (${\sigma }^2\simeq 0.05$) as relevant for PBH formation, $\beta$ (15) can change by various orders of magnitude if we reduce the threshold (${\delta }_{\mathrm{c}}$) by about $20\%$. In fact,

$$\frac{\beta ({\delta }_{\mathrm{c}}=0.4)}{\beta ({\delta }_{\mathrm{c}}=1/3)}\simeq {10}^{-5},$$

(25)

demonstrating how tuned the conditions are for producing a cosmologically interesting population of PBHs.

• Collapse fraction vs. curvature perturbation

While it is customary to use the smoothed density contrast at horizon crossing to estimate the number of collapsing regions, it is also possible to work directly with the comoving curvature perturbation to approximately compute the PBH fraction ($\beta$) at the time of formation. In this case, there is no need of to rely on the smoothing procedure of subhorizon fluctuations provided by the window functions [152,153]. Interestingly, as we will see later, this approach also provides a way to assess the effects of large primordial non-Gaussianity that might be present in some of the PBH-forming inflationary scenarios.

To understand the role of the primordial curvature fluctuation ($\mathcal{R}$), we approximate its variance with the power spectrum ${\sigma }_{\mathcal{R}}^2\approx {\mathcal{P}}_{\mathcal{R}}$. Using the Press–Schechter approach with a Gaussian PDF for the curvature fluctuation spectrum, the fraction of collapsing regions at formation can be estimated as

$${\beta }_G={\int }_{{\mathcal{R}}_{\mathrm{c}}}^{\infty }\mathrm{d}\mathcal{R}\frac{e^{-{\mathcal{R}}^2/\left(2{\sigma }_{\mathcal{R}}^2\right)}}{\sqrt{2\pi }{\sigma }_{\mathcal{R}}}\simeq \frac{1}{2}\mathrm{Erfc}\left(\frac{{\mathcal{R}}_{\mathrm{c}}}{\sqrt{2{\mathcal{P}}_{\mathcal{R}}}}\right),$$

(26)

where ${\mathcal{R}}_{\mathrm{c}}$ is the threshold. To roughly estimate ${\mathcal{R}}_{\mathrm{c}}$, we can assume an almost scale-invariant power spectrum of $\mathcal{R}$ for a logarithmic range of wave numbers relevant to PBH formation. Making use of a Gaussian window function in (19) gives, in this case, ${\sigma }_{\mathcal{R}}^2=\langle {\mathcal{R}}^2\rangle \simeq 8{\mathcal{P}}_{\mathcal{R}}/81$. Finally, plugging the latter into (15), and comparing the resulting expression with (26), we obtain [152]:

$${\mathcal{R}}_{\mathrm{c}}\approx \frac{9}{2\sqrt{2}}{\delta }_{\mathrm{c}}.$$

(27)

For a density threshold of ${\delta }_{\mathrm{c}}=0.4$, the relation presented above gives ${\mathcal{R}}_{\mathrm{c}}\simeq 1.3$, which we set as the fiducial value for the estimates below. Following these considerations and using the formulas derived so far—in particular Equations (12) and (26)—we can then repeat the previous estimates and determine the approximate amplitude of the power spectrum required for PBH formation. The result is that for Gaussian primordial fluctuations, a sensible PBH population today requires ${\mathcal{P}}_{\mathcal{R}}\sim {10}^{-2}$. These findings confirm our earlier results of Equation (22).

• The case of non-Gaussian curvature fluctuations

So far, we have assumed that primordial fluctuations obey Gaussian statistics in order to estimate the amplitude of the power spectrum required for PBH formation. Since PBHs are expected to form through extremely rare large fluctuations (see Figure 3), any small deviations in the shape of the tail of the fluctuation distribution—which essentially depend on the amount of non-Gaussianity (i.e., skewness of the PDF)—can have a significant impact on the PBH abundance [153,154,155,156,157,158,159,160,161]. For the sake of obtaining a lower limit on the amplitude of the PBH-forming curvature power spectrum, we now consider scenarios with large primordial non-Gaussianity. A particularly interesting case of this type occurs if the main source of the curvature perturbation results from a higher-order interaction, where the distribution of $\mathcal{R}$ can be modeled as a ${\chi }^2$ distribution (see, e.g., [79,153,162]):

$$\mathcal{R}=g^2-\langle g^2\rangle ,$$

(28)

where g is a Gaussian random variable ($\langle g\rangle =0$) with variance ${\sigma }_g^2\equiv \langle g^2\rangle$. The PDF of $\mathcal{R}$ in this case can be determined by making a change of variable $P_{\mathrm{NG}}\left(\mathcal{R}\right)=P_G\left(g\right)|\mathrm{d}g/\mathrm{d}\mathcal{R}|$, which takes the following form

$$P_{\mathrm{NG}}\left(\mathcal{R}\right)=\frac{e^{-(\mathcal{R}+{\sigma }_g^2)/2{\sigma }_g^2}}{2\sqrt{2\pi (\mathcal{R}+{\sigma }_g^2)}{\sigma }_g}.$$

(29)

Making a change of variable to a quantity t through the definition ${\sigma }_g^{2}t=\mathcal{R}+{\sigma }_g^2⟶\mathrm{d}\mathcal{R}={\sigma }_g^2\mathrm{d}t$, the fraction of regions in the universe that can collapse to form PBHs can be estimated as

$${\beta }_{\mathrm{NG}}={\int }_{{\mathcal{R}}_{\mathrm{c}}}^{\infty }\mathrm{d}\mathcal{R}{P}_{\mathrm{NG}}\left(\mathcal{R}\right)\simeq \frac{1}{2}\mathrm{Erfc}\left(\sqrt{\frac{1}{2}+\frac{{\mathcal{R}}_{\mathrm{c}}}{\sqrt{2{\mathcal{P}}_{\mathcal{R}}}}}\right),$$

(30)

where in the last step, we approximate the variance as $\langle {\mathcal{R}}^2\rangle =2{\langle g^2\rangle }^2=2{\sigma }_g^4\approx {\mathcal{P}}_{\mathcal{R}}$, in order to express ${\beta }_{\mathrm{NG}}$ in terms of the curvature power spectrum14. We can now compute the amplitude of the curvature power spectrum required for PBH formation when the statistics of fluctuations are strongly non-Gaussian. Using (12), together with (30), a population of solar mass PBHs with $f_{\mathrm{pbh}}={10}^{-3}$ requires ${\mathcal{P}}_{\mathcal{R}}\simeq 1.5\times {10}^{-3}$, whereas for a population of PBHs with $M_{\mathrm{pbh}}^{\left(\mathrm{f}\right)}={10}^{-12}{M}_{\odot }$ and $f_{\mathrm{pbh}}=1$, we find ${\mathcal{P}}_{\mathcal{R}}\simeq 9\times {10}^{-4}\simeq {10}^{-3}$. Hence, we conclude that the required amplitude of the power spectrum is roughly given by

$${\mathcal{P}}_{\mathcal{R}}\left(k_{\mathrm{pbh}}\right)\sim {10}^{-3}\mathrm{for}\mathrm{non}\text{-}\mathrm{Gaussian}\mathrm{fluctuations}.$$

(31)

Compared to the case of Gaussian distributed curvature perturbation (see Equation (22)), we learn that the required amplitude of the power spectrum is reduced by about one order of magnitude. Therefore, a curvature perturbation with a smaller amplitude can produce the same amount of PBH abundance if non-Gaussianity is present (see Equation (28)). In particular, for ${\mathcal{R}}_{\mathrm{c}}^2\gg {\mathcal{P}}_{\mathcal{R}}$, which is typically satisfied by a satisfactory approximation, we can disregard the $1/2$ factor in (30). Compared with Equation (26), the power spectrum required to generate the same collapse fraction of PBHs in both cases can be related as

$${\mathcal{P}}_{{\mathcal{R}}_{\mathrm{NG}}}\simeq \frac{2}{{\mathcal{R}}_{\mathrm{c}}^2}{\mathcal{P}}_{{\mathcal{R}}_{\mathrm{G}}}^2.$$

(32)

To further illustrate these points, in Figure 4, we show the quantity ($\beta$) both for Gaussian and non-Gaussian cases, represented as a function of the curvature power spectrum. We learn from the figure that in the non-Gaussian case, within a phenomenologically relevant interval of $\beta$ (see, e.g., (12)), a given value of the power spectrum leads to a much larger value of $\beta$. We emphasize that we focus on a specific type of non-Gaussian distribution (namely ${\chi }^2$) to estimate the amplitude of the power spectrum required for PBH formation; therefore, the derived results could change for milder cases depending on the amplitude and sign of the non-Gaussianity (i.e., depending on whether the PDF in (26) is positively or negatively skewed) [153].

2.4. Brief Summary and the Path Ahead

Let us summarize the arguments that we have reviewed so far. We computed the required amplitude of the small-scale primordial power spectrum (${\mathcal{P}}_{\mathcal{R}}\left(k_{\mathrm{pbh}}\right)$) required to generate a sizable population of PBHs that can account for all or a fraction of DM density. The typical small scale of PBH formation ($k_{\mathrm{pbh}}$) is related to the BH mass through Equation (6) (see Table 1 for examples). Compared with the power spectrum at large CMB scales, we need a

$$\Delta {\mathcal{P}}_{\mathcal{R}}\equiv \frac{{\mathcal{P}}_{\mathcal{R}}\left(k_{\mathrm{pbh}}\right)}{{\mathcal{P}}_{\mathcal{R}}\left(k_{\mathrm{cmb}}\right)}\sim {10}^6-{10}^7$$

(33)

enhancement in the spectrum amplitude between small and large scales, depending on the statistics obeyed by the primordial curvature perturbation (see Equations (22) and (31)).

We also learned that the PBH abundance is extremely sensitive to the amplitude of the primordial curvature spectrum. Notice that the results we reviewed are derived for the case of PBHs produced during the RDU; if early phase transitions or early phases of non-standard cosmology occur, the corresponding modified equations of state can also considerably influence the properties of the PBH population [143,163]. An interesting example is the QCD phase transition, which can lead to a high peak in the distribution of solar-mass PBHs several orders of magnitude larger than the corresponding values in the RDU [164,165].

There are various opportunities to improve and elaborate on these results. In our considerations, we assumed that PBHs form at a particular mass (Equations (4) and (6)), for example, by means of a sharply peaked primordial power spectrum; moreover, we ignored the effects of clustering [22,166,167,168,169] and mergers [170,171,172] on the PBH abundance and evolution. As shown in [140,173,174], assumptions about the shape of the primordial spectrum may alter the PBH abundance and the corresponding clustering properties. Another topic of debate concerns the use of over-density ($\delta$) versus the curvature perturbation ($\mathcal{R}$) when computing the PBH abundance (see, e.g., [175] for a discussion on these issues). In light of the discussion presented above, we emphasize that the calculations we carried out in this section should be regarded as rough order-of-magnitude estimates, requiring more precise numerical analysis. Furthermore, in discussing the effects of non-Gaussianities on PBH formation, we stress that we computed the corresponding power spectrum only for an extreme example of ${\chi }^2$ non-Gaussian statistics. For a detailed analysis of the impact of primordial non-Gaussianity on PBH formation and abundance, we refer the reader to the general discussion in [176]. An additional phenomenological consequence of PBH-forming inflationary scenarios is the inevitable production of a stochastic gravitational wave (GW) background [87,88]. In fact, an enhanced spectrum of curvature fluctuations, as needed to produce PBH, acts as a source of GW [89]. The characterization of this GW background (see, e.g., [93,94,96]), along with the properties of its scalar sources (see, e.g., [92,97,177]) and the corresponding forecasts for its detection, is an important avenue for the experimental probing of PBH-forming inflationary models. We refer the reader to [105] for a detailed recent review of scalar-induced GW backgrounds.

All the topics mentioned above are being actively developed by the PBH community. The arguments and results we reviewed in this section are sufficient for introducing our specific purpose, which is to review the theoretical foundation of inflationary scenarios leading to PBHs. In the following sections, we discuss different conceptual ideas and concrete inflationary mechanisms for obtaining the enhancement (33) of the curvature power spectrum, as needed to generate PBHs. We focus on the inflationary theory aspects only, without computing the resulting PBH abundance, as well as other phenomenological properties that are already covered in various recent complementary reviews [28,29,30,31,32,33,34].

3. Enhancement of Scalar Perturbations during Single-Field Inflation

We now focus our attention on inflationary scenarios able to lead to PBH formation. As we learned in the previous section, such scenarios are characterized by a significant enhancement in the curvature power spectrum at a scale ($k_{\mathrm{pbh}}$) (which depends on the PBH mass) much smaller with respect to CMB scales ($k_{\mathrm{cmb}}$). The condition to satisfy is Equation (23), which we rewrite here as:

$$\frac{{\mathcal{P}}_{\mathcal{R}}\left(k_{\mathrm{pbh}}\right)}{{\mathcal{P}}_{\mathcal{R}}\left(k_{\mathrm{cmb}}\right)}\sim {10}^7.$$

(34)

We classify inflationary models as single-field (this section) or multifield type (next section) depending on whether the mechanism responsible for the enhancement in the scalar fluctuations relies on a single- or multifield scenario. In general, existing inflationary mechanisms amplify the spectrum of curvature fluctuations by means of significant gradients in the background evolution of fields responsible for inflation. In this section, we phrase our discussion as model-independently as possible, mostly focusing on conceptual aspects of the problem. We aim to discuss the dynamics and the general properties of curvature fluctuations in inflationary models leading to PBHs and refer to specific representative scenarios when necessary.

• Main References: Our discussion in this section is based on [41,42,52,64].

3.1. The Dynamics of Curvature Perturbation

In order to analyze the behavior of the scalar power spectrum in single-field scenarios, we consider the second-order action of scalar perturbations around an inflationary phase of evolution. The background metric corresponds to a (quasi) de Sitter background with a nearly constant Hubble parameter (H). Cosmological inflation is controlled by a slow-roll parameter ($ϵ\equiv \dot{H}/H^2$) satisfying $ϵ\ll 1$, with $ϵ=1$ corresponding to the condition required to conclude the inflationary process. We work with conformal time ($\tau \le 0$) during inflation (see, e.g., [178] for a classic survey of inflationary models).

The dynamics of scalar fluctuations can be formulated in terms of the comoving curvature perturbation ($\mathcal{R}$) [129,130], the second-order action (at the lowest order in derivatives)15 of which takes the following form (see, e.g., [64])

$$S_{\mathcal{R}}^{\left(2\right)}=\frac{1}{2}\int \mathrm{d}\tau {\mathrm{d}}^{3}x\frac{2a^2\left(\tau \right)M^2\left(\tau \right)ϵ\left(\tau \right)}{c_s^2\left(\tau \right)}\left[{{\mathcal{R}}^{\prime }}^2-c_s^2\left(\tau \right){\left(\overrightarrow{\nabla }\mathcal{R}\right)}^2\right].$$

(35)

In this formula, $c_s$ is the sound speed of the curvature perturbation, M is an effective time-dependent Planck mass and $ϵ$ is the aforementioned slow-roll parameter.

First, we present a few initial words to contextualize single-field models aimed to produce PBHs, leading to a set of dynamics of curvature perturbation controlled by action (35). The simplest option to consider is PBH-forming models with unit sound speed and constant Planck mass characterized only by the shape of the potential ($V\left(\varphi \right)$). As mentioned in the Introduction, models in this class require a potential characterized by a flat, plateau-like region (see, e.g., [46,47,48,49,50,51,52,53,55] for a choice of works studying this possibility; we will discuss its implications for the dynamics of curvature perturbations in the next subsection). PBH-forming potentials with the required characteristics can find explicit realizations, for example, in models of Higgs inflation [47,185,186,187],alpha attractors [188,189] and string inflation [51,52,190]. Considering more complex possibilities, PBH-generating models that exploit a time dependence for the sound speed are based on non-canonical kinetic terms for the inflation scalar, such as K-inflation [191,192] (see, e.g., [64,65,66,67,193,194,195,196,197] for concrete examples and Section 3.3 for some of their implications). Finally, scenarios with a time-dependent effective Planck mass can be generated by non-minimal couplings of the inflation scalar with gravity, as in the Horndeski action [198] and its cosmological applications to G-inflation scenarios [199]. Realizations of PBH-forming models in setups with non-minimal couplings belonging to the Horndeski sector include [200,201,202,203]. To the best of our knowledge, early universe models based on the more recent covariant DHOST actions [204,205,206] have not yet been explored in the context of PBH model building.

Interestingly, despite the many distinct concrete realizations, all single-field scenarios rely on a few common mechanisms to enhance the spectrum of curvature fluctuations, which exploit the behavior of background quantities. We are now going to discuss these mechanisms in a model-independent way. We treat M, $c_s$ and $ϵ$ as appearing in action (35) as time-dependent quantities controlled by the single scalar background profile that drives inflation. To start with, it is convenient to redefine the time variable in action (35) to adsorb the time-dependent $c_s$ into a rescaled conformal time and impose an equal-scaling condition of time and space coordinates:

$$\mathrm{d}\overline{\tau }=c_s\mathrm{d}\tau ⟹S_{\mathcal{R}}^{\left(2\right)}=\frac{1}{2}\int \mathrm{d}\overline{\tau }{\mathrm{d}}^{3}x{z}^2\left(\overline{\tau }\right)\left[{\mathcal{R}}^{\prime 2}-{\left(\overrightarrow{\nabla }\mathcal{R}\right)}^2\right],$$

(36)

with a prime indicating a derivative with respect to $\overline{\tau }$, the rescaled conformal time. Importantly, we introduce a so-called time-dependent ‘pump field’ ($z\left(\overline{\tau }\right)$) as

$$z^2\left(\overline{\tau }\right)=\frac{2a^2\left(\overline{\tau }\right)M^2\left(\overline{\tau }\right)ϵ\left(\overline{\tau }\right)}{c_s\left(\overline{\tau }\right)}.$$

(37)

The dynamics of $\mathcal{R}$ are strongly tied to the time dependence of the pump field ($z\left(\overline{\tau }\right)$) and, more generally, to the behavior of the background quantities ($M,ϵ,c_s$) that constitute it.

To analyze the evolution mode by mode, we work in Fourier space and write the Euler–Lagrange mode equation for curvature perturbation derived from the action (36):

$$\frac{1}{z^2\left(\overline{\tau }\right)}{\left[z^2\left(\overline{\tau }\right){\mathcal{R}}_k^{\prime }\left(\overline{\tau }\right)\right]}^{\prime }=-k^2{\mathcal{R}}_k\left(\overline{\tau }\right),$$

(38)

where $k\equiv |\overrightarrow{k}|$ is the magnitude of the wave number that labels a given mode. This is a differential equation involving derivatives along the time direction acting on the function ${\mathcal{R}}_k\left(\overline{\tau }\right)$ depending on both time and momentum (k).

To express its solution, we implement a gradient expansion approach (see, e.g., [41,42,52]), starting from the solution in the limit of small $k/\left(aH\right)$ and including its momentum-dependent corrections, which solve (38) order by order in a $k/\left(aH\right)$ expansion. This approach is particularly suitable for our purpose of describing scenarios in which the size of small-scale curvature fluctuations ($k/\left(aH\right)$, large) differs considerably from that of large-scale curvature fluctuations ($k/\left(aH\right)$, small) (see condition (23)). Indeed, a gradient expansion allows us to better understand the physical origin of possible mechanisms that raise the curvature spectrum at small scales.

The most general solution of Equation (38), up to second order in powers of $k/\left(aH\right)$, is formally given by the following integral equation16

$${\mathcal{R}}_k\left(\overline{\tau }\right)={\mathcal{R}}_k^{\left(0\right)}\left[1+\frac{{\mathcal{R}}_k^{\left(0\right){}^{\prime }}}{{\mathcal{R}}_k^{\left(0\right)}}{\int }_{{\overline{\tau }}_0}^{\overline{\tau }}\frac{\mathrm{d}{\overline{\tau }}^{\prime }}{{\tilde{z}}^2\left({\overline{\tau }}^{\prime }\right)}-k_2{\int }_{{\overline{\tau }}_0}^{\overline{\tau }}\frac{\mathrm{d}{\overline{\tau }}^{\prime }}{\tilde{z}{\left({\overline{\tau }}^{\prime }\right)}^2}{\int }_{{\overline{\tau }}_0}^{\overline{\tau }{}^{\prime }}\mathrm{d}{\overline{\tau }}^{{}^{″}}{\tilde{z}}^2\left({\overline{\tau }}^{{}^{″}}\right)\frac{{\mathcal{R}}_k\left({\overline{\tau }}^{{}^{″}}\right)}{{\mathcal{R}}_k^{\left(0\right)}}\right],$$

(39)

where the sub- and superscripts 0 denote a reference time, and a tilde over a time-dependent quantity indicates that it is normalized with respect to its value at $\overline{\tau }={\overline{\tau }}_0$.

Typically, we are interested in relating the late time curvature perturbation at $\overline{\tau }$ to the same quantity computed at some earlier time (${\tau }_0$). For this purpose, it is convenient to identify ${\overline{\tau }}_0$ as the time coordinate evaluated soon after horizon crossing and ${\mathcal{R}}_k^{\left(0\right)}$ as the mode function computed at ${\overline{\tau }}_0$. In order to enhance the spectrum of curvature fluctuations at small scales (recall the PBH-forming condition of Equation (23)), we can envisage two possibilities. One option is to exploit the structure of Equation (39), making sure that its contributions within the square parentheses become increasingly important as time proceeds after modes leave the horizon. In this way, we generate a sizable scale dependence for ${\mathcal{R}}_k\left(\overline{\tau }\right)$ after horizon crossing, with the possibility of amplifying the small-scale curvature spectrum. Alternatively, we can design methods that lead to significant scale dependence already at horizon crossing, i.e., for the quantity ${\mathcal{R}}_k^{\left(0\right)}$, which then maintains its value at superhorizon scales. We explore both these options in Section 3.2 and Section 3.3, respectively.

To develop a quantitative discussion, it is convenient to introduce the so-called slow-roll parameters as

$$\eta \equiv \frac{\mathrm{d}\ln ϵ}{\mathrm{d}N},s\equiv \frac{\mathrm{d}\ln c_s}{\mathrm{d}N},\mu \equiv \frac{\mathrm{d}\ln M^2}{\mathrm{d}N}$$

(40)

where, in our definition, we make use of the relation between e-foldings and the time coordinate ($\overline{\tau }$: $\mathrm{d}N=H\mathrm{d}t=(aH/c_s)\mathrm{d}\overline{\tau }$).

In standard models of inflation based on inflationary attractor dynamics, one imposes the so-called slow-roll conditions throughout the entire inflationary period, corresponding to the requirements

$$\eta ,s,\mu ,\frac{\mathrm{d}\eta }{\mathrm{d}N},\frac{\mathrm{d}s}{\mathrm{d}N},\frac{\mathrm{d}\mu }{\mathrm{d}N}\ll 1,$$

(41)

which imply that the pump field always grows in time $z^2\propto {(-\overline{\tau })}^{-2}$ as $\overline{\tau }\to 0$ (see Equation (37)). As a consequence, the second and third terms in the general solution (39) decay as ${(-\tau )}^3$ and ${(-\tau )}^2$, respectively in the late time limit ($(-\tau )\to 0$). Hence, they can be identified as decaying modes17 that rapidly cease to play any role in the dynamics of curvature perturbations. This is a regime of a slow-roll attractor, wherein soon after horizon crossing, the curvature perturbation settles into a nearly-constant configuration (${\mathcal{R}}_k^{\left(0\right)}$), the spectrum of which is almost scale-invariant. In this case, the momentum-dependent terms in Equation (39) do not have the opportunity to raise the curvature spectrum at small scales.

Hence, in order to produce PBHs, we need to go beyond the slow-roll conditions of Equation (41), as first emphasized in [207]. Before discussing concrete ideas to do so, in view of numerical implementations, as well to improve our physical understanding, it is convenient to express the curvature perturbation (Equation (38)) in a way that makes more manifest the role of slow-roll parameters in controlling the mode evolution. We introduce a canonical variable ($v_k$) defined as

$$v_k\left(\overline{\tau }\right)=z\left(\overline{\tau }\right){\mathcal{R}}_k\left(\overline{\tau }\right).$$

(42)

Plugging this definition into Equation (38), we obtain the so-called Mukhanov–Sasaki equation, which reads

$$v_k^{{}^{″}}\left(\overline{\tau }\right)+\left(k^2-\frac{z^{{}^{″}}}{z}\right)v_k\left(\overline{\tau }\right)=0,$$

(43)

where

$$\frac{z^{{}^{″}}}{z}=2{\left(\frac{aH}{c_s}\right)}^2\left[1+\Theta \right].$$

(44)

Expanding the derivatives of z (37) in terms of the slow-roll parameters of Equation (40), we define

$$\begin{array}{cc}\hfill \Theta & \equiv -\frac{1}{2}(ϵ+s)+\frac{3}{4}(\eta -s+\mu )+\frac{1}{8}{(\eta -s+\mu )}^2-\frac{1}{4}(ϵ+s)(\eta -s+\mu )\hfill \\ \hfill & +\frac{1}{4}\left(\frac{\mathrm{d}\eta }{\mathrm{d}N}-\frac{\mathrm{d}s}{\mathrm{d}N}+\frac{\mathrm{d}\mu }{\mathrm{d}N}\right).\hfill \end{array}$$

(45)

Standard slow-roll attractor scenarios correspond to situations in which $\Theta$ is negligibly small; the quantity in Equation (44) then reads $2{\overline{\tau }}^{-2}$, leading to a scale-invariant curvature power spectrum. To break the scale invariance of curvature perturbation, we need to consider a sizable time-dependent $\Theta$. We note that the expression (45) is exact and does not assume any slow-roll hierarchy as Equation (41). Hence, it can be used to study the system beyond slow roll, as we will do in the following section.

3.2. Enhancement through the Resurrection of the Decaying Mode

• The idea

An interesting mechanism to enhance the curvature perturbation at superhorizon scales is suggested by the structure of the integrals within the square parentheses of Equation (39). Suppose that for a brief time interval, a given mode (k) experiences a background evolution during which the pump field (z) rapidly decreases after the horizon exit epoch (${\overline{\tau }}_0$). Then, the ‘decaying’ mode can increase, and the integrals in the parentheses of Equation (39) can contaminate the nearly constant solution (${\mathcal{R}}_k^{\left(0\right)}$), eventually leading to a late-time value (${\mathcal{R}}_k\left(\overline{\tau }\right)\gg {\mathcal{R}}_k^{\left(0\right)}$) on superhorizon scales. This situation signals a significant departure from the attractor slow-roll regime discussed after Equation (41). In fact, in this case, the criterion for the enhancement of the curvature perturbation can be explicitly phrased in terms of the derivative of the pump-field transiently changing sign during a short time interval during inflation:

$$\frac{z^{\prime }}{z}=\frac{aH}{c_s}\left[1+\frac{\eta -s+\mu }{2}\right]<0.$$

(46)

This condition implies that the combination of the slow-roll parameters ($\eta -s+\mu$) should be of order $\mathcal{O}\left(1\right)$ and negative during some e-folds during inflation, violating the slow-roll conditions (41). In particular, we require

$$\eta -s+\mu <-2.$$

(47)

If Equation (46) is satisfied, the slow-roll conditions (41) are not satisfied, and the contributions within parentheses in Equation (39) can increase. Strong time gradients of homogeneous background quantities, which lead to condition (47), can then be converted into a small-scale amplification of the curvature power spectrum. As discussed in [64], the expression (46), along with the considerations mentioned above, generalizes to a time-dependent sound speed and Planck mass according to the arguments first developed in [41,42].

• Model building parameterization of the non-attractor phase

To illustrate a viable model that can generate a seven-order-of-magnitude enhancement required for PBH formation—see Equation (23)—we focus on canonical single-field models ($c_s\to 1,M\to M_{\mathrm{pl}}$ and $s\to 0,\mu \to 0$) in order to simplify our analysis. The background evolution for the single scalar field driving inflation is

$$\ddot{\varphi }+3H\dot{\varphi }+V^{\prime }\left(\varphi \right)=0,$$

(48)

where $V\left(\varphi \right)$ is the scalar potential, and the time derivatives are carried on in coordinate time (t). The non-slow-roll dynamics are controlled by the properties of the potential (V), as we are going to discuss, and by its consequences for the behavior of the inflation velocity ($\dot{\varphi }$).

Since in these scenarios, the pump field can be parameterized purely in terms of the slow-roll parameter ($ϵ$) as $z=a\sqrt{2ϵ}{M}_{\mathrm{pl}}$ (see, e.g., Equation (37)), the linear dynamics of ${\mathcal{R}}_k$ (Equation (39)) are dictated by the first slow-roll parameter, the evolution of which is, in turn, determined by the sign and amplitude of the slow-roll parameter $\eta$. Hence, the criterion required to realize the desired growth in the spectrum can be simply parameterized as a condition of the second slow-roll parameter as $\eta <-2$ in (47).

From a concrete model-building perspective, scalar potentials ($V\left(\varphi \right)$) that can induce this type of dynamics include a characteristic ‘plateau’ within a non-vanishing field range ($\Delta \varphi \ne 0$) [38,39,208]. This property gives rise to phases of transient non-attractor dynamics of ultra slow-roll (USR) [40,43,44] or constant-roll (CR) [45,209,210] evolution, depending on the shape profile of the potential around the aforementioned feature. In particular, for USR evolution, the potential typically has a very flat plateau with $V^{\prime }\simeq 0$, whereas for constant-roll evolution, $V^{\prime }<0$ so that the field climbs a hill by overshooting a local minimum18. As the scalar field, during its evolution, traverses such a flat region with negligible potential gradient, the acceleration term ($\ddot{\varphi }$) is balanced by the Hubble damping term in the Klein–Gordon Equation (48), and the inflation speed is no longer controlled by the scalar potential. This phenomenon significantly changes the values of the inflation velocity ($\dot{\varphi }$) during the transient non-attractor phase and inevitably leads to the violation of one of the slow-roll conditions:

$$\ddot{\varphi }+3H\dot{\varphi }+V^{\prime }\left(\varphi \right)=0⟹\eta =-6-\frac{2V^{\prime }}{\dot{\varphi }H}+2ϵ;$$

(49)

hence, $\eta \simeq -6$ ($\eta <-6$) for transient USR ($V^{\prime }=0$) and CR ($V^{\prime }<0$) phases 19. We emphasize that since the non-slow-roll inflationary era is characterized by a large negative $\eta$ for a brief interval of e-folds, the pump field, as well as the first slow-roll parameter ($ϵ$), quickly decays during this stage as required for the activation of the decaying modes. In fact,

$$\frac{\mathrm{d}\ln ϵ}{\mathrm{d}N}\equiv \eta ⟹z^2\propto ϵ\propto e^{-\left|\eta \right|\Delta N},$$

(50)

where, for simplicity, we assume a constant $\eta$ during the non-attractor phase. For explicit inflationary scenarios that can realize such transient phases in the context of PBH formation, see, e.g., [49,51,52,211]. Nevertheless, it is worth pointing out that, although possible, explicit constructions of suitable inflationary potentials involve a high degree of tuning to render the potentials extremely flat for a small region in the field range and ensure an appropriate transition for the scalar velocity among different epochs (see, e.g., the discussion in [50], as well as the comments at the end of this section).

After this general discussion of model building, in the analysis that follows, we do not need to work with an explicit form of potential ($V\left(\varphi \right)$) to analyze the enhancement through the non-attractor dynamics. Instead, we exploit the general idea we are discussing in a model-independent way and model PBH-forming inflationary scenarios as a succession of distinct phases that connect smoothly with one another, each parameterized by a constant $\eta$ (related approaches are developed in [212,213,214,215,216]). Our perspective encompasses the important features of scenarios based on the idea of transient resurrection of the decaying mode at superhorizon scales, satisfying Equations (46) and (47). In order to capture the transitions among phases, we multiply each phase by the smoothing function [217]:

$$\sigma (N,\Delta )=\frac{1}{2}\left[\tanh \left(\frac{N-N_i}{\Delta }\right)-\tanh \left(\frac{N-N_f}{\Delta }\right)\right],$$

(51)

where N denotes e-folds; $N_i$ and $N_f$ are the e-folding numbers at the beginning and end of the constant $\eta$ phase, respectively; and $\Delta$ signifies the duration of the smoothing procedure. Keeping this smoothing prescription in mind, the inflationary evolution can be divided into three phases:

• Phase I. The initial phase of inflationary evolution is characterized by a standard slow-roll (SR) regime, where $ϵ,\eta \ll 1$ and $ϵ<\eta$ at the pivot scale ($k_{\mathrm{cmb}}=0.05{\mathrm{Mpc}}^{-1}$) (assuming that modes at the pivot scale exit the horizon at the beginning of evolution ($N_{\mathrm{b}}=0$)) in order to match Planck observations [4];
• Phase II. As the scalar field starts to traverse the flat plateau-like region in its potential, its dynamics eventually enter the non-attractor era, lasting a given number of e-folds of evolution. This phase is characterized by a large negative $\eta <-6$, during which the first slow-roll parameter ($ϵ$) decays exponentially:

  $$ϵ\left(N\right)\equiv \exp \left[{\int }_{60}^N\eta \left(N^{\prime }\right)\mathrm{d}{N}^{\prime }\right];$$

  (52)
• Phase III. The final phase of evolution ensures a graceful exit from the non-attractor phase into a final slow-roll epoch, leading to the end of inflation. Since $ϵ$ decays quickly in the non-attractor era, this final phase is characterized by a hierarchy between the slow-roll parameters:

  $$\eta \gg ϵ$$

  (53)

  where $\eta >0$. We typically require a large positive $\eta$ to bring back $ϵ$ from its tiny values at the end of the non-attractor era, towards the value ($ϵ=1$) needed to conclude inflation. To capture this behavior accurately, we split the final phase of evolution into two parts, parameterizing $\eta$ as

  $$\eta \left(N\right)={\eta }_{\mathrm{III}}^{\left(1\right)}{\sigma }_1(N,\Delta )+{\eta }_{\mathrm{III}}^{\left(2\right)}{\sigma }_2(N,\Delta ).$$

  (54)
  The relevant parameter choices to model the dynamics can be found in the third column of

  Table 2. Parameter choices that characterize the background evolution of $\eta$ smoothed by the function (51) in each phase. Note that the final phase of evolution is divided into two in orders to accommodate the end of inflation with $ϵ=1$ at $N_{\mathrm{end}}=60$.

  	Phase I	Phase II	Phase III
  $\eta$	$0.02$	$-6.30$	0.30	3.00
  $N_i$	$0.00$	$33.2$	35.7	55.0
  $N_j$	$33.2$	$35.7$	55.0	65.0
  $\nu$	$1.00$	$0.50$	1.00	2.00

  .
  We note that our choice of $\eta$ in the initial stage of the Phase III and in Phase II is not a coincidence; in most of the single-field modes, there exists a correspondence that relates the $\eta$ values in Phase II and Phase III, i.e., ${\eta }_{\mathrm{III}}=-6-{\eta }_{\mathrm{II}}$, which is a consequence of Wands’ duality [218]. We will elaborate the consequence of this correspondence in the context of the power spectrum below, in particular for modes that exit the horizon as the background evolves from Phase II to Phase III.

In light of the discussion presented above, we can characterize the full background evolution using the Hubble hierarchy in (52) and $H\left(N\right)=H_{\mathrm{end}}\exp [-{\int }_{60}^Nϵ\left(N^{\prime }\right)\mathrm{d}{N}^{\prime }]$, where $H_{\mathrm{end}}$ denotes the Hubble rate at the end of inflation, where $N_{\mathrm{end}}=60$.

For a representative set of parameter choices (see Table 2), in Figure 5, we show an example of background evolution in which we plot $ϵ,\eta$ and $z^{\prime }/z$. The right panel of the figure makes manifest that the background evolution leads to $z^{\prime }/z<0$ for a short interval of e-folds ($N=33–35.7$), as highlighted by the red region in the plot. In accordance with our discussion so far, this behavior is appropriate for triggering a significant enhancement in the power spectrum of curvature perturbation through the resurrection of the decaying mode.

• Numerical analysis

Having obtained the background evolution, we are ready to describe mode evolution to obtain the power spectrum of curvature perturbation towards the end of inflation20:

$${\mathcal{P}}_{\mathcal{R}}(k,N_{\mathrm{end}})=\frac{k^3}{2{\pi }^2}|\frac{v_k\left(N_{\mathrm{end}}\right)}{z\left(N_{\mathrm{end}}\right)}{|}^2,$$

(55)

where to study the evolution of curvature perturbations, we make use of the canonical variable ($v_k$) and consider the Mukhanov–Sasaki system of Equations (43)–(45) after setting $s=\mu =0$. In general, it is not possible to find full analytic solutions for this system of equations, and a numerical analysis is needed. However, as we will explain soon, interesting properties of the resulting curvature spectrum can be derived and understood analytically. We implement the numerical procedure explained in detail in the technical Appendix C, which solves the Mukhanov–Sasaki equation with Bunch–Davies initial conditions, and we provide a Python code that reproduces our numerical findings21. The resulting power spectrum is represented in Figure 6; it manifestly grows in amplitude towards small scales, exhibiting a peak at around $k_{\mathrm{peak}}\simeq {10}^{12}{\mathrm{Mpc}}^{-1}\gg k_{\mathrm{cmb}}$. Notice that the spectrum grows as $k^4$ towards its peak and is characterized by a dip preceding the phase of steady growth [212]. We will have more to say soon about these features.

Interestingly, for the system under consideration, the bulk of the enhancement can be attributed to the active dynamics of the would-be ‘decaying modes’, i.e., the second and third terms of Equation (39). To show this explicitly, we study superhorizon solution of the curvature perturbation in Appendix C by applying the Formula (A38), which is a special case of Equation (39) applied to the canonical single-field scenario we discuss here. For a grid of wave numbers that exit the horizon during the initial slow-roll era, the amplitude of the power spectrum obtained in this way is shown by blue dots in Figure 6. The accuracy of these locations with respect to the full numerical result (black solid line) confirms our expectation that decaying modes in (39) play a crucial role for the enhancement of the curvature perturbation for this scenario. In the right panel of Figure 6, we zoom in on the growth and the subsequent decay of power spectrum following the peak.

• The features of the spectrum: analytic considerations

Besides the numerical findings presented above, we can derive general analytic results for the spectrum of curvature fluctuations in scenarios activating the would-be decaying modes through a brief non-attractor era.

We start by noticing that for modes that leave the horizon during the initial slow-roll stage (leftmost light blue region in Figure 6), the spectrum shows characteristic features such as the presence of a dip, followed by an enhancement parameterized by a spectral index of $n_s-1=4$ during the bulk of the growth [95,212,224,225,226]. The dip is physically caused by a disruptive interference between the ‘constant’ mode of curvature fluctuation at superhorizon scales and the ‘decaying’ mode that is becoming active and ready to contribute to the enhancement of the spectrum. The position and depth of the dip are analytically calculable in terms of other features of the spectrum, at least in a limit of short duration of the non-slow-roll epoch. It is found that the position of the dip in momentum space is proportional to the inverse fourth root of the enhancement of the spectrum, and the depth of the dip is proportional to the inverse square root of the enhancement of the spectrum [226]. These relations are valid for any single-field models that enhance the spectrum through a brief deviation from the standard attractor era, including cases with a time-varying sound speed and Planck mass. They are accompanied by consistency conditions on the squeezed limit of non-Gaussian higher-order point functions around the dip [227,228,229], as expected in single-field scenarios22.

While in the considerations of the previous paragraph, we considered modes leaving the horizon during the first stage of slow-roll evolution, we can also derive analytic results for what happens during the non-attractor epoch. In fact, for modes that exit the horizon deep in the non-attractor era (light green region in the middle of Figure 6) and the following final slow-roll era, the spectrum behaves as expected in a standard slow-roll phase, with spectral index

$$n_s-1=-2ϵ-\eta \simeq -{\eta }_{\mathrm{III}}=6+{\eta }_{\mathrm{II}}.$$

(56)

(Recall that the latin numbers $\mathrm{II}$ and $\mathrm{III}$ relate to the phases of evolution; see Equations (52) and (54).) This behavior is a manifestation of the duality invariance of perturbation spectra within distinct inflationary backgrounds, called Wands’ duality (see, e.g., [218,231]). Wands’ duality can be understood by noticing that the structure of the Mukhanov–Sasaki Equation (43) is unchanged by a redefinition of the pump field that leaves the combination $z^{{}^{″}}/z$ invariant:

$$\begin{array}{ccc}\hfill z\left(\overline{\tau }\right)& \to & \tilde{z}\left(\overline{\tau }\right)\equiv z\left(\overline{\tau }\right)\left(c_1+c_2{\int }^{\overline{\tau }}\frac{d\tilde{\tau }}{z\left(\tilde{\tau }\right)}\right)\Rightarrow \frac{{\tilde{z}}^{{}^{″}}}{\tilde{z}}=\frac{z^{{}^{″}}}{z}\hfill \end{array}$$

(57)

where $c_{1,2}$ are arbitrary constants. If $z\left(\overline{\tau }\right)$ controls a phase of the slow-roll attractor, $\tilde{z}\propto 1/\overline{\tau }$, a dual phase whose pump field ($z\propto {\overline{\tau }}^2$) as given by Equation (57) describes a non-attractor era. Although the statistics of the canonical variable ($v_k$) are identical in the two regimes, the amplitude of the curvature perturbation spectrum (${\mathcal{R}}_k$) increases in the non-attractor epoch. In scenarios in which the parameter $\eta$ is considerably larger than the other slow-roll parameters, Wands’ duality (57) analytically prescribes the relation (56), in agreement with the numerical findings plotted in Figure 6. Subtleties can arise in joining attractor and non-attractor phases, since consistency conditions can be violated [232] due to the effects of boundary conditions at the transitions between different epochs. All these considerations are relevant for our topic, given the sensitivity of PBH formation and properties on the shape of the spectrum near the peak.

For further detailed accounts of the characterization of the interesting features in the power spectrum of PBH-forming single-field scenarios, we refer the reader to [95,212,214,215,217,224,225,226,233].

• Stochastic inflation and quantum diffusion

While, so far, we have focused on the predictions of the second-order action (35), non-linearities and non-Gaussian effects can play an important role in the production of PBHs, as we learned in Section 2.3. For the case of ultra-slow-roll (USR) models based on non-attractor phases of inflation, there are sources of non-Gaussianity associated with stochastic effects during inflation.

The stochastic approach to inflation pioneered by Starobinsky [234] constitutes a powerful formalism for describing the evolution of coarse-grained superhorizon fluctuations during inflation. It is based on a classical (but stochastic) Langevin equation, which reads in canonical single-field inflation (N is the number of e-folds, and we assume constant sound speed and Planck mass):

$$\frac{d\varphi }{dN}=-\frac{V^{\prime }}{3H^2}+\frac{H}{2\pi }\xi \left(N\right).$$

(58)

Here, $\varphi$ represents a coarse-grained version of superhorizon scalar fluctuations; $V^{\prime }$ is the derivative of the inflationary potential, which leads to a deterministic drift for the coarse-grained superhorizon mode; and $\xi$ is a source of stochastic noise acting on long-wavelength fluctuations caused by the continuous kicks of modes that cross the cosmological horizon and pass from sub- to superhorizon scales during inflation.

Besides the physical insights that it offers, the inflationary stochastic formalism [234,235,236,237,238,239,240,241,242] offers the opportunity to obtain accurate results for the probability distribution function controlling coarse-grained superhorizon modes beyond any Gaussian approximation. As a classic example, by solving the Fokker–Planck equation associated with (58), the seminal work [238] analytically obtained the full non-Gaussian distribution functions for certain representative inflationary potentials, going beyond the reach of a perturbative treatment of the problem.

Returning to the discussion of a USR inflationary evolution for PBH scenarios, we can expect that stochastic effects can be very relevant in this context (see, e.g., [241,243,244,245,246,247,248,249,250,251]). In fact, since the amplitude of scalar fluctuations are amplified, the stochastic noise can increase relative to that under slow-roll inflation. Moreover, during USR, the derivative of the potential ($V^{\prime }=0$), the classical drift is absent, and the stochastic evolution is driven by stochastic effects only. Various researchers have studied this topic by solving the stochastic evolution equations [252,253,254,255] and found that non-Gaussian effects can change the predictions of PBH formation, depending on the duration of the USR phase. In fact, the stochastic noise modifies the tails for the curvature probability distribution function, which decays with an exponential (instead of a Gaussian) profile23 and, consequently, tends to overproduce PBHs. Refs. [252,253,254,255] set constraints on the duration of the USR phase, which (depending on the scenario) can last, at most, a few e-folds before overproducing PBHs. There has been increased activity surrounding these subjects, and we refer the reader to the aforementioned literature for details on the state of the art with respect to this important topic.

3.3. Growth in the Power Spectrum When the Decaying Modes Are Slacking

• Slow-roll violation without triggering decaying modes

We learned in the previous subsections that a possible way to enhance the spectrum of fluctuations at small scales with respect to its large-scale counterpart is to amplify the k-dependent corrections to the constant-mode solution (${\mathcal{R}}_k^{\left(0\right)}$) within the parentheses of Equation (39).

However, as we anticipated in the paragraph following Equation (39), we can also design scenarios in which an enhanced time dependence of the slow-roll parameters leads to a scale-dependent curvature power spectrum at horizon crossing, even without exciting the decaying mode at superhorizon scales. The idea is to still make sure that the pump field ($z\left(\tau \right)$) increases with time—hence, conditions (46) and (47) are not satisfied, the decaying mode remains inactive and the terms within parentheses in Equation (39) can be neglected. However, at the same time, each individual slow-roll parameter changes considerably during a short time interval during inflation. The derivatives of slow-roll parameters can be large; they can contribute significantly to the quantity ($\Theta$) controlling the Mukhanov–Sasaki equation, and they can influence the scale dependence of the curvature spectrum at horizon crossing (see Equations (43) and (45)).

We start this section by setting the stage and derive formulas to describe this possibility. We then present an explicit realization of this scenario. It is convenient to work with the canonical variable ($v_k$) defined through Equation (42), and solve the Mukhanov–Sasaki system in the form of the set of Equations (43)–(59). We assume that the pump field (z) is monotonic and always increasing with time, and we identify the sound horizon of fluctuations as $aH/c_s\simeq \sqrt{z^{{}^{″}}/\left(2z\right)}$24. We can identify two asymptotic regimes for each mode k: (i) an early-time regime, when each mode is deep inside the horizon; and (ii) a late-time regime, when the modes are stretched to become a superhorizon. On the one hand, in the former regime, the modes satisfy $k^2\gg z^{{}^{″}}/z$ and behave as the standard vacuum fluctuations in Minskowski space time

$$v_k\left(\overline{\tau }\right)=\frac{e^{-ik\overline{\tau }}}{\sqrt{2k}},\mathrm{for}{k}^2\gg \frac{z^{{}^{″}}}{z}.$$

(59)

On the other hand, later during inflation, the fluctuations are stretched outside the horizon, entering the second regime and eventually satisfying $k^2\ll z^{{}^{″}}/z$, with a solution given by

$$v_k\left(\overline{\tau }\right)\simeq {\mathcal{C}}_{1,k}z+{\mathcal{D}}_{2,k}z\int \frac{\mathrm{d}{\overline{\tau }}^{\prime }}{z^2\left({\overline{\tau }}^{\prime }\right)}+\mathcal{O}\left(k^2\right),$$

(60)

where the finite $k^2$ corrections to this solution can be derived in a similar fashion as in (39). Recall that we are now interested in attractor background configurations; hence, we can neglect the last two terms in Equation (60) that rapidly decay. Shortly after horizon crossing, the canonical variable settles into the solution $v_k=z{\mathcal{C}}_{1,k}$. Using the field redefinition (42), we can identify the constant mode as the curvature perturbation at late times (${\mathcal{C}}_{1,k}={\mathcal{R}}_k={\mathcal{R}}_k^{\left(0\right)}$). In order to determine its expression, we match the solutions sometime around horizon crossing ($\overline{\tau }={\overline{\tau }}_0$), and we obtain

$$|{\mathcal{C}}_{1,k}{|}^2={\left|{\mathcal{R}}_k^{\left(0\right)}\right|}^2=\frac{1}{2k}\frac{1}{z{\left({\overline{\tau }}_0\right)}^2}=\frac{1}{2k}\frac{c_s}{a^2{M}^{2}2ϵ}{|}_{\overline{\tau }={\overline{\tau }}_0}.$$

(61)

The horizon-crossing time can be conveniently expressed as a leading contribution in a WKB approximation [64]:

$$k^2=\frac{z^{{}^{″}}}{2z}{|}_{{\overline{\tau }}_0}={\left(\frac{aH}{c_s}\right)}^2(1+\Theta ){|}_{{\overline{\tau }}_0},$$

(62)

with $\Theta$ given in Equation (45). Collecting there results, we can write the late-time power spectrum for curvature fluctuations as

$${\mathcal{P}}_{\mathcal{R}}\left(k\right)\equiv \frac{k^3}{2{\pi }^2}{\left|\frac{v_k\left(\overline{\tau }\right)}{z}\right|}^2=\frac{k^3}{2{\pi }^2}{\left|{\mathcal{R}}_k^{\left(0\right)}\right|}^2=\frac{H^2}{8{\pi }^2ϵc_s{M}^2}\left(1+\Theta \right){|}_{{\overline{\tau }}_0}.$$

(63)

From (63), we observe that rapid changes in the background quantities ($ϵ,c_s,M$) and slow-roll parameters constituting the quantity ($\Theta$) of Equation (45) as a function of ${\overline{\tau }}_0$ can then translate into a scale-dependent amplification of the power spectrum. As we will see, this situation leads to a scale-dependent enhancement in the power spectrum realized through the ‘constructive interference’ of the time-dependent background parameters ($ϵ,c_s,M$).

• An explicit realization

We now review a possible realization of this scenario, closely following the discussion presented in [64]. We focus on the generalized single-field framework discussed in Section 3.1, and we consider background dynamics that include simultaneous pronounced dips in the time-dependent profiles for the parameters $ϵ,c_s,M^2$. We then study the power spectrum by numerically solving numerically the Mukhanov–Sasaki equation for curvature perturbations (see Appendix C), and we compare the result with the analytical expressions discussed in Section 3.1 (see, e.g., (63)).

To analyze a representative scenario in this category, we parameterize the three time-dependent quantities ($X=\{ϵ,c_s,M^2\}$) as [64]

$${\log }_{10}X\left(N\right)=\left(n_{\mathrm{a}}-n_{*}\right){\tanh }^2\left(\frac{N-N_{*}}{\sigma }\right)+n_{*}.$$

(64)

Each of these quantities tends to ${10}^{n_{\mathrm{a}}}$ asymptotically away from $N_{*}$ in both directions, $|N-N_{*}|\gg \sigma$ and becomes equal to ${10}^{n_{*}}$ at $N_{*}$ while staying around the neighborhood of this value for a number of e-folds determined by $\sigma$. Notice that we are interested in features of the inflationary dynamics that affect modes at scales much smaller than the CMB pivot scale. Hence, we can assume $N_{*}\gg 1$, where modes associated with CMB scales leave the horizon at $N\simeq 0$, while inflation ends at $N_{\mathrm{end}}=60$. A representative set of parameter choices that leads to a localized decrease in the slow-roll parameters is presented in

Table 3. Parameter choices that characterize the background evolution of the time-dependent parameters ($ϵ$, $c_s$, $\tilde{M}=M/M_{\mathrm{pl}}$).

X	${\mathit{n}}_{\mathbf{a}}$	${\mathit{n}}_{*}$	${\mathit{N}}_{*}$	$\mathit{\sigma }$
$ϵ$	$-2.0$	$-4.0$	$36.0$	$2.0$
$c_s$	$0.0$	$-2.0$	$35.5$	$2.5$
${\tilde{M}}^2$	$0.0$	$-3.0$	$34.8$	$4.0$

. The resulting background evolution and the behavior of $z^{\prime }/z$ are shown in the left and right panels of Figure 7, respectively. We observe from the right panel of the figure that $z^{\prime }/z>0$ is always satisfied during inflation, suggesting (as expected) that the decaying modes do not grow over time. Using the parameterization (64), we then numerically solve the Mukhanov–Sasaki equation, following Appendix C.

Our results for ${\mathcal{P}}_{\mathcal{R}}(k,N_{\mathrm{end}})$ are presented in Figure 8. In the left panel, we notice that the expression (63) accurately describes the behavior of the power spectrum towards its peak, confirming our expectation that the constant growing mode is responsible for the enhancement. Notice the absence of a dip proceeding the growth (which characterizes the scenario shown in Figure 6). This is in agreement with our interpretation presented in the previous section; the dip is due to disruptive interference between ‘growing’ and ‘decaying’ modes, while in this context, the decaying mode is not active. On the other hand, as shown in the right panel of the figure, the power spectrum exhibits oscillations for scales following the peak. As we discussed in Note 24, this is due to multiple horizon crossing of modes within certain momentum scales, leading to excited states and oscillations in the spectrum. We illustrate this phenomenon in Figure 9—for two neighboring modes labeled by red and purple dots in Figure 8—where we plot k versus $\sqrt{|z^{{}^{″}}/z|}$ as a function of e-folds. As shown in the right panel of Figure 9, although these modes are neighboring, they exhibit non-trivial behavior before their final horizon exit ($N\simeq 40$) so that their asymptotic values ($N\to N_{\mathrm{end}}$) differ considerably, giving rise to the sizable modulations in the late-time power spectrum (see Figure 8), as discussed in [64].

Besides [64], the ideas discussed in this subsection for enhancing the spectrum at small scales are further realized in [65]. Conceptually similar frameworks include a sound speed resonance scenario proposed in [66]—which can be explicitly realized within a DBI inflation model [67]—and a parametric resonance model discussed within the single-field canonical inflation framework [257].

3.4. Brief Summary

We find it remarkable that many distinct single-field models of inflation built with the aim of producing PBHs share common features in the properties of the resulting power spectrum because the enhancement of the curvature spectrum at small scales is the result of few mechanisms common to several scenarios. We identified the idea of resurrecting the decaying mode of curvature fluctuations at super Hubble scales by increasing the absolute value of the slow-roll parameter ($\eta$)—see Section 3.2—and the idea of having very rapid changes in the values of slow-roll parameters (keeping them relatively small)—see Section 3.3.

The resulting properties of the curvature spectrum might be essential for testing single-field models of inflation at small scales, thanks to their predictions for the PBH population properties, as well as for the spectrum of gravitational waves induced at the second order in perturbations [87,88]. In fact, the latter is a very interesting topic that is relevant and well-studied for single-field models; we refer the reader to [105] for a detailed review.

Nevertheless, the existing explicit scenarios of single-field inflation, which are able to generate an appreciable PBH abundance, typically suffer from severe fine-tunings on the choices of the parameters characterizing their Lagrangians, for example, in producing sufficiently flat plateau regions of their potential and in ensuring regular transitions between attractor and non-attractor eras during the inflation process (see e.g., the discussion in [50] and [258] for a scenario that can partially ameliorate the tuning involved). Moreover, as we learned in Section 2.2, the resulting PBH population can be very sensitive to the details of the spectrum profile and to the presence of primordial non-Gaussianities. The latter should, to a certain extent, be generated in many concrete single-field realizations (see Section 3.2) and render particularly delicate the task of making precise theoretical predictions. In fact, large non-Gaussianities might not only change the predictions for PBH production but can also impose restrictions on single-field model building due to large one-loop corrections [259,260,261,262] (see [263]).

It is then interesting to try to consider PBH inflationary scenarios in contexts in which more than one field acquires dynamics during inflation. The hope is to find qualitatively new ideas for the production of PBHs or, alternatively, more natural realizations of known mechanisms in order to amplify the primordial curvature spectrum at small scales. It is possible that multiple fields affect the predictions on the statistics of curvature fluctuations with respect to single-field models, with important implications for PBH formation. We review this topic in the following section.

4. Enhanced Primordial Power Spectrum in Multifield Models

Despite the remarkable agreement of single-field, slow-roll inflation with the CMB and the LSS data [4], the fundamental nature of inflation continues to elude us. In particular, it is very important to know whether additional fields took part in the dynamics of inflation, besides the scalar driving cosmic expansion. In fact, while current cosmological observations do not provide hints of isocurvature fluctuations associated with extra inflationary degrees of freedom, it might also be the case that additional fields play a role during epochs of inflation that are not well probed by current large-scale surveys. The formation of PBHs occurring at small scales can be sensitive to isocurvature fluctuations and represent a valuable probe of inflationary multifield dynamics. Multifield inflation can offer new possibilities for the production of PBHs, exploiting novel, distinctive ways of converting large gradients of background quantities into the spectrum of curvature fluctuations. Moreover, when used in tandem with the techniques discussed in the previous section, multifield models might alleviate some of the fine-tuning issues the occur in single-field scenarios (see the discussion in Section 3.4). For these reasons, in this section, we review a selected set of theoretical frameworks aimed at enhancing the curvature spectrum at small scales25 by using tools that specifically involve more than one field during inflation, in particular the tachyonic behavior of field dynamics in some extra sectors during inflation.

First, in Section 4.1, we review PBH scenarios based on axions interacting with gauge fields during inflation, exploiting a mechanism based on particle production during inflation. This possibility is well-motivated by particle physics constructions, and it has been thoroughly explored in the literature; we make efforts to carefully review the state of the art, discussing opportunities and challenges for PBH formation in this context. Then, in Section 4.2, we review multifield inflationary scenarios in which the rise in the spectrum is a result of large isocurvature perturbations induced by a curved inflationary trajectory traversing a rich multifield moduli space.

• Main References: In the discussion presented in Section 4.1, we benefit from the works presented in [79,80,81,82,85,86], while the material presented in Section 4.2 is based on [72,73].

4.1. Enhanced Scalar Perturbations from Axion-Gauge Field Dynamics

Axions [269,270] are pseudoscalar particles theoretically introduced in studies of particle physics theories beyond the Standard Model. First motivated by scenarios addressing the strong CP problem in terms of a Peccei–Quinn mechanism [271,272], they find natural realizations in string theory (see, e.g., [273]), as well as many useful applications in cosmology—see, e.g., [274] for a comprehensive review. The possibility of using the physics of axions to produce PBHs is very interesting, given their strong theoretical and experimental motivations that enable connections between particle physics, quantum gravity and cosmology. In this section, we review several representative ideas and their realizations in this context.

Axion fields, denoted here as $\chi$, correspond to Goldstone bosons of Peccei–Quinn (PQ) global symmetry spontaneously broken at a scale dubbed f. Their Goldstone nature equips axions with a shift symmetry $\chi \to \chi +$ constant, which is valid at all orders in perturbations. The shift symmetry makes them interesting inflation candidates, as their potential is naturally very flat [275]. The PQ symmetry suffers from a chiral anomaly, which breaks the shift symmetry and assigns axions a potential ($V\left(\chi \right)$) through non-perturbative contributions in field theory, or monodromy effects in string theory (see, e.g., [68] for a review). The chiral anomaly implies that axions interact with gauge vectors through higher-dimensional operators. This property is essential for our arguments. The typical structure of an axion Lagrangian considered for cosmological purposes reads

$$\mathcal{L}=-\frac{1}{2}{\partial }_{\mu }\chi {\partial }^{\mu }\chi -V\left(\chi \right)-\frac{1}{4}{F}_{\mu \nu }{F}^{\mu \nu }-\frac{g_{\mathrm{cs}}}{4f}\chi F_{\mu \nu }{\tilde{F}}^{\mu \nu }.$$

(65)

The axion $\chi$ is equipped by a kinetic term and a potential term ($V\left(\chi \right)$). The potential is often a polynomial function of the axion field, at least for small values of $\chi$26. As anticipated, the Lagrangian (65) also includes a $U\left(1\right)$ gauge field ($A_{\mu }$) with field strength ($F_{\mu \nu }$), coupling to the axion through a dimension-5 gauge-preserving operator

$${\mathcal{L}}_{\mathrm{int}}=-\frac{g_{\mathrm{cs}}}{4f}\chi F_{\mu \nu }{\tilde{F}}^{\mu \nu },$$

(66)

where f is the PQ symmetry-breaking scale, and $g_{\mathrm{cs}}$ is a dimensionless coupling constant. The dual field strength appearing in Equation (66) is ${\tilde{F}}^{\mu \nu }\equiv {\eta }^{\mu \nu \rho \sigma }{F}_{\rho \sigma }/\left(2\sqrt{-g}\right)$, and the alternating symbol $\eta$ is 1 for even permutations of its indices and $-1$ otherwise (starting with ${\eta }^{0123}=1$). The dimension-5 operator (66) is especially important for our arguments. As we review below, this axion–vector coupling leads to a tachyonic instability in the gauge sector, which exponentially enhances one of the helicities of the gauge vector modes and produces a large amount of gauge quanta [280]. The energy adsorbed in the production of vector modes is carried away from the axion kinetic energy, which slows down its evolution; this is good news for axion inflationary models, since the axion can slowly roll, even along steep potentials [281]. Moreover, the inverse decay of enhanced gauge quanta into axion fluctuations (schematically, $A+A\to \chi$) can amplify curvature perturbations [282,283] and lead to a large curvature spectrum at small scales at a level able to produce PBHs [79]. We refer the reader to [284] for a review on axion inflation written one decade ago, which includes many of the topics mentioned above. We now focus on reviewing the consequences of these ideas for scenarios amplifying curvature fluctuations, in addition to covering opportunities and challenges put forward in the most recent literature on this subject.

• Amplification of gauge field fluctuations

When considering a dynamical axion field with a time-dependent background profile ($\overline{\chi }\left(t\right)$), the dimension-5 axion-vector interaction of Equation (66) leads to a copious production of vector fluctuations. This property can be understood by considering the equation of motion (EoM) for the gauge field mode functions (See Appendix D) [281]

$${\partial }_x^2{A}_{\pm }+\left(1\pm \frac{2\xi }{x}\right)A_{\pm }=0,\xi \equiv -\frac{g_{\mathrm{cs}}\dot{\overline{\chi }}}{2Hf}$$

(67)

where $A_{\pm }$ corresponds to the gauge vector polarizations, and we define a dimensionless time variable ($x\equiv -k\tau =k/\left(aH\right)$), as well as the effective dimensionless coupling ($\xi$) between the spectator axion and the gauge field. Without loss of generality, we work within the following conditions

$$\xi >0\mathrm{and}\dot{\overline{\chi }}<0.$$

(68)

Hence, the axion rolls along its potential from large positive to small values of $\overline{\chi }\ge 0$.

Notice from (67) that the dimension-5 operator of Equation (66) introduces a time-dependent mass term in the dispersion relation of the $\mathrm{U}\left(1\right)$ field, which changes sign depending on the helicity ($\lambda =\pm$). This property reflects the parity-violating nature of the dimension-5 interaction. When the gauge modes are deep inside the horizon ($(x=k/(aH)\gg 1)$), the time-dependent correction term is negligible, and both vector polarizations obey a standard dispersion relation as in Minkowski space. However, as the modes stretch outside the horizon, the correction becomes dominant for $x=k/\left(aH\right)\lesssim 2\xi$, leading to an instability for one of the circular polarizations of the gauge fields. In our conventions (68), only the $A_{-}$ state experiences a tachyonic instability, while $A_{+}$ is unaffected. The dynamics of the axion-like field controlled by the axion velocity ($\dot{\overline{\chi }}\ne 0$), therefore induces a significant production of helical vector fields.

The nature of gauge field production and its consequences as a source for scalar perturbations are sensitive to the scalar potential ($V\left(\chi \right)$), since this quantity determines the profile of the axion background velocity ($\dot{\overline{\chi }}$). The dynamics of the curvature perturbations, as generated through the axion-gauge field dynamics, depend on whether we identify the axion field as the inflation or whether it belongs to some extra spectator sector during inflation. In what comes next, we arrange our discussion so as to clearly distinguish between these possibilities. We consider the following scenarios:

1. 

Section 4.1.1 As a first possibility, we identify the axion ($\chi$) with $\varphi$ that drives inflation: $\chi =\varphi$. Then, the order parameter ($\dot{\overline{\varphi }}$) controlling the gauge-field production increases over time, generating scalar perturbations through an inverse vector decay: $\delta A+\delta A\to \delta \varphi$. We refer to this scenario as “Smooth Axion Inflation”; models that can be considered in this category are studied in [79,80,81,82]. Such scenarios suffer from dynamical instabilities associated with large back-reaction effects from the gauge fields on the axion evolution.

2. 

Section 4.1.2 In certain axion-inflation models, the axion potential has special features located far away from the field range corresponding to scales affecting CMB physics. A sudden increase in the axion inflation velocity occurs at their location, with enhanced scalar perturbations amplified at small PBH scales through inverse vector decay. This possibility was first discussed in [81], while further developments were studied in [83,86]. We refer to scenarios producing PBHs by exploiting localized particle production as “Bumpy Axion Inflation”, and we analyze how they address the instabilities mentioned above.

3. 

Section 4.1.3 A final possibility corresponds to gauge field production in a hidden sector through the dynamics of an axion spectator field. This case also allows one to address the aforementioned dynamical instabilities, given that the back-reaction effects from the vector sector can be placed under control. Depending on the shape of the spectator axion potential, such a scenario can lead to localized peaks in the scalar curvature power spectrum [81,85]. We refer to this scenario as “Spectator axion-gauge field dynamics”.

We build our discussion mainly in terms of representative, concrete examples, presenting the main results and relegating technical details to appendices.

4.1.1. Smooth Axion Inflation

In the first scenario, we identify the axion ($\chi$) as the inflaton $\varphi$: $\chi \to \varphi$ to study the behavior of scalar perturbations in a setup described by Lagrangian (65). We assume that the profile for the scalar potential ($V\left(\varphi \right)$) is sufficiently flat to support inflation. In this case, the effective coupling ($\xi$), as defined in Equation (67), adiabatically increases during the inflationary process:

$$\xi \equiv g_{\mathrm{cs}}\sqrt{\frac{ϵ}{2}}\frac{M_{\mathrm{pl}}}{f},$$

(69)

where $ϵ$ is the standard slow-roll parameter. The amplified gauge-field mode function can be analytically expressed as [281] (recall that the axion speed has a negative sign; see Equation (68)):

$$A_{-}(\tau ,k)\simeq \frac{1}{\sqrt{2k}}{\left(\frac{-k\tau }{2\xi }\right)}^{1/4}\exp (\pi \xi -2\sqrt{-2\xi k\tau }),\xi \equiv -\frac{g_{\mathrm{cs}}\dot{\overline{\varphi }}}{2Hf}.$$

(70)

The slowly changing time-dependent parameters ($\xi$ and H) in Equation (70) are evaluated at the epoch of horizon crossing. The amplification of gauge field modes is maximized when the size of the mode is comparable to the horizon, $-k\tau \sim \mathcal{O}\left(1\right)$.

In fact, the analytic expression in (70) is valid within the interval ${\left(8\xi \right)}^{-1}\ll -k\tau \ll 2\xi$. For the values $(\xi \sim \mathcal{O}(1)$) we are be interested in, this range corresponds to a phase during which the gauge modes grow and then remain frozen to their maximal value before being diluted away by universe expansion. We proceed by discussing the time evolution of the relevant quantities in this setup. For concreteness, we focus on a representative example. For a more detailed account of the gauge field amplification by the slowly rolling scalar, we refer the reader to Section 2.1 of [285] or Appendix B of [286].

• Background evolution in a concrete example

As the effective coupling ($\xi$) of Equation (70) increases during inflation, the enhanced vector modes eventually lead to a sizable back-reaction on the background evolution of the axionic inflation field. This fact can significantly affect the homogeneous dynamics of the system.

The back-reaction is mainly controlled by the vector-dependent friction term in the equation of motion for the homogeneous inflation field [281]. This phenomenon implies that the gauge field amplification occurs at the expense of the inflation velocity ($\dot{\overline{\varphi }}$). We consider a situation characterized by a transition between standard slow-roll dynamics in the early stages of inflation and a new attractor regime at late stages, when the gauge field enhancement dominates the Hubble friction [287] (see the cautionary remark towards the end of this section).

The aforementioned friction effect induced by the gauge-mode production can be analyzed considering modified Klein–Gordon and Friedmann equations:

$$\begin{array}{cc}\hfill \ddot{\overline{\varphi }}+3H\dot{\overline{\varphi }}& +V^{\prime }\left(\overline{\varphi }\right)=\frac{g_{\mathrm{cs}}}{f}\langle \overrightarrow{E}·\overrightarrow{B}\rangle ,\hfill \\ \hfill 3H^2{M}_{\mathrm{pl}}^2& =\frac{1}{2}{\dot{\overline{\varphi }}}^2+V\left(\overline{\varphi }\right)+\frac{1}{2}〈{\overrightarrow{E}}^2+{\overrightarrow{B}}^2〉,\hfill \end{array}$$

(71)

where we introduced ‘electric’ and ‘magnetic’ field contributions27, as discussed in detail in Appendix D (see, in particular, Equation (A63)). Using the solutions (70) for gauge-field modes, the expectation values for the electric and magnetic fields can be computed analytically [281,286], finding (see Appendix D)

$$\langle \overrightarrow{E}·\overrightarrow{B}\rangle \simeq 2.1\times {10}^{-4}\frac{H^4}{{\xi }^4}{e}^{2\pi \xi },{\rho }_A\equiv \frac{1}{2}\langle {\overrightarrow{E}}^2+{\overrightarrow{B}}^2\rangle =1.4\times {10}^{-4}\frac{H^4}{{\xi }^3}{e}^{2\pi \xi }.$$

(72)

The quantity (${\rho }_A$) corresponds the total energy density contained in the gauge field fluctuations. Using the evolution Equation (71), we can identify the conditions corresponding to a small back-reaction of the vector modes into the background evolution. They are $3H|\dot{\overline{\varphi }}|\gg g_{\mathrm{cs}}\langle \overrightarrow{E}·\overrightarrow{B}\rangle /f$ and ${\rho }_A\ll 3H^2{M}_{\mathrm{pl}}^2$, which can be expressed as

$$\begin{array}{cc}\hfill {\xi }^{-3/2}{e}^{\pi \xi }& \ll 79\frac{\dot{\overline{\varphi }}}{H^2},⟶\mathrm{negligible} \mathrm{back}\text{-}\mathrm{reaction} \mathrm{on} \varphi \mathrm{equation},\hfill \\ \hfill {\xi }^{-3/2}{e}^{\pi \xi }& \ll 146\frac{M_{\mathrm{pl}}}{H},⟶\mathrm{negligible} \mathrm{back}\text{-}\mathrm{reaction} \mathrm{on} \mathrm{the} \mathrm{Friedmann} \mathrm{equation}.\hfill \end{array}$$

(73)

Assuming that inflation starts in a standard slow-roll regime, the relation $\dot{\varphi }=\sqrt{2ϵ}H{M}_{\mathrm{pl}}\ll H{M}_{\mathrm{pl}}$ implies that the first condition in (73) is more demanding than the second. When these conditions are not met, we enter into an inflationary phase characterized by a strong back-reaction of gauge modes, which goes beyond the standard slow-roll inflationary attractor.

To concretely illustrate the back-reaction of the gauge field production on the homogeneous background evolution of inflation, we focus on a modified version of the linear monodromy-type potential [276] that interpolates between $V\propto \varphi$ and $V\propto {\varphi }^2$ from large to small field values around the global minimum:

$$V\left(\varphi \right)=\lambda M_{\mathrm{pl}}^4\left[\sqrt{1+\frac{{\varphi }^2}{M_{\mathrm{pl}}^2}}-1\right],$$

(74)

where $\lambda$ is a dimensionless parameter that fixes the overall amplitude of the potential. Assuming negligible back-reaction in the beginning of inflation, $\lambda$ can determined by enforcing the standard normalization of the power spectrum in the slow-roll regime for a given choice of ${\varphi }_{\mathrm{in}}$ using

$${\mathcal{P}}_{\mathcal{R}}\left(k_{\mathrm{cmb}}\right)=\frac{H^2}{8{\pi }^2ϵM_{\mathrm{pl}}^2}\simeq \frac{V\left({\varphi }_{\mathrm{in}}\right)}{24{\pi }^2{ϵ}_V\left({\varphi }_{\mathrm{in}}\right)M_{\mathrm{pl}}^4}\equiv 2.1\times {10}^{-9}.$$

(75)

We then numerically solve Equation (71) utilizing (72) by setting the coupling between the inflation and the gauge fields to its maximally allowed value28 by Planck: $g_{\mathrm{cs}}{M}_{\mathrm{pl}}/f=48$ (see e.g., [289]). Note that for a given initial value of the scalar field, this number can be translated to an initial value of the effective coupling, which we label as ${\xi }_{\mathrm{cmb}}$:

$${\xi }_{\mathrm{cmb}}\simeq \frac{g_{\mathrm{cs}}{M}_{\mathrm{pl}}}{f}\frac{1}{2{\varphi }_{\mathrm{in}}}.$$

(76)

For ${\varphi }_{\mathrm{in}}=10.5M_{\mathrm{pl}}$ (${\xi }_{\mathrm{cmb}}=2.5$), the results of the numerical computation are shown in Figure 10. In the left panel, we show the inflation profile as a function of e-folds until the end of inflation, where $\dot{H}=-H^2$. We can clearly observe that in the absence of coupling to the gauge fields ($g_{\mathrm{cs}}=0$) (black dashed line with $\xi =0$), inflation lasts for about ≈48 e-folds. Turning on the interactions with gauge fields at the beginning of inflation (solid orange line with ${\xi }_{\mathrm{cmb}}=2.5$), the back-reaction induced by the vector fluctuations become noticeable around ≈36 e-folds, which, in turn, has the overall effect of extending the duration of inflation for about ≈12 e-folds compared to the standard slow-roll case ($\xi =0$).

These findings are also supported by the right panel of Figure 10, where we compare the standard Hubble damping term ($3H\dot{\varphi }$) with the friction term ($g_{\mathrm{cs}}\langle \overrightarrow{E}·\overrightarrow{B}\rangle /f$) induced by the gauge field production as a function of e-folds until the end of inflation. We observe that the standard Hubble friction dominates the dynamics at early times, but as the effective coupling ($\xi$) increases, the dynamics gradually evolve into a regime dominated by the back-reaction of the produced gauge quanta. For the parameter choices we adopt, the onset of this regime starts around $N_{\mathrm{br}}\simeq 36$ e-folds, and at around $N\simeq \mathcal{O}\left(50\right)$ the system enters into a strong back-reaction regime wherein the dynamics become completely controlled by the gauge modes [281].

Therefore, in any viable model of axion inflation, back-reaction effects induced by gauge field production become important only towards the latest stages of inflation. For example, for a rough estimate of the orders of magnitude involved, we can adopt a linear potential ($V\propto \varphi$). Assuming the dynamics are initially in the slow-roll regime, we have $\xi \simeq g_{\mathrm{cs}}\sqrt{2{ϵ}_V}{M}_{\mathrm{pl}}/f\propto 1/\varphi$ and ${\varphi }_{\mathrm{cmb}}/{\varphi }_{\mathrm{br}}={\xi }_{\mathrm{br}}/{\xi }_{\mathrm{cmb}}\simeq 2.2–1.9$. For typical initial field values that can sustain enough inflationary evolution, we have ${\varphi }_{\mathrm{cmb}}\sim \mathcal{O}\left(10\right)M_{\mathrm{pl}}$, which implies that ${\varphi }_{\mathrm{br}}\sim \mathcal{O}\left(5\right)M_{\mathrm{pl}}$. Hence, back-reaction becomes dominant in the latest stages but before the end of inflation corresponding to ${\varphi }_{\mathrm{end}}\sim \mathcal{O}\left(1\right)$.

• Scalar fluctuations sourced by gauge fields

As inflation progresses, particle production becomes increasingly efficient. Hence, the vector modes start to act as an additional source of inflation fluctuations ($\delta \varphi$) due to the presence of trilinear coupling ($\delta {\mathcal{L}}_{\mathrm{int}}\propto \delta \varphi F\tilde{F}$) (see Appendix D.2 and references therein). In particular, a tachyonic amplification of the vector fields leads—at the second order in perturbations—to an enhancement in the scalar sector. Since these contributions to the scalar fluctuations are statistically independent of their vacuum counterpart, the resulting curvature power spectrum acquires two separate contributions (see Appendix E) that can be parameterized as [79,281]

$$\begin{array}{cc}\hfill {\mathcal{P}}_{\mathcal{R}}& ={\mathcal{P}}_{\mathcal{R}}^{\left(\mathrm{v}\right)}+{\mathcal{P}}_{\mathcal{R}}^{\left(\mathrm{s}\right)},\hfill \\ & \simeq \frac{H^2}{8{\pi }^2ϵM_{\mathrm{pl}}^2}+{\left(\frac{g_{\mathrm{cs}}}{f}\frac{\langle \overrightarrow{E}·\overrightarrow{B}\rangle }{3H\beta \dot{\overline{\varphi }}}\right)}^2{\mathcal{F}}^2.\hfill \end{array}$$

(77)

The first term in the second line denotes the standard vacuum contribution to the curvature power spectrum. The second accounts for the vector part, with the time-dependent quantities ($\beta$ and $\mathcal{F}$) defined in Equations (A81) and (A99). In the weak back-reaction regime, the power spectrum of the curvature perturbation is dominated by the vacuum contribution in the first term of Equation (77). As inflation progresses, the effective coupling between the inflation and the vector fields increases; the source contribution to the power spectrum starts to kick in, and the power spectrum grows as ${\mathcal{P}}_{\mathcal{R}}\propto {\langle \overrightarrow{E}·\overrightarrow{B}\rangle }^2\sim {\mathrm{e}}^{4\pi \xi }$ (the system is still in the weak back-reaction regime, where $\beta \simeq 1$). Then, the scalar fluctuations grow; they start influencing the gauge fields sources and lead to an increase in the damping factor ($\beta \sim \langle \overrightarrow{E}·\overrightarrow{B}\rangle /3H\dot{\overline{\varphi }}$). Therefore, the sourced contribution to the power spectrum eventually saturates towards late times during inflation. In fact, the time-dependent factor ($\mathcal{F}$) is introduced to capture the modified definition of the curvature perturbation in the strong back-reaction regime. As discussed in detail in [81], this factor is responsible for an order-one correction to the power spectrum amplitude towards the end of inflation.

We summarize these arguments with a plot (Figure 11) of the total power spectrum with respect to the e-folding number. We notice that as the effective coupling ($\xi )$ between the scalar and the gauge fields increases along the inflationary trajectory, the power spectrum acquires a large contribution from the vector source. As desired, it grows towards smaller scales, eventually saturating towards a constant value. Since the source of the peak in the power spectrum originates from a trilinear coupling between the inflation and the gauge fields, the statistics of the curvature perturbation at those scales is non-Gaussian (see also [290]). The values of the curvature spectrum amplitude at the peak appearing in Figure 11 are sufficiently large to generate a sizable population of PBHs (recall the discussion in Section 2.3).

• Cautionary remarks: an instability in the inflation dynamics

We warn the reader that the findings reviewed in this section have been recently questioned. In particular, the authors of [290,291,292,293,294] studied the axion-gauge field dynamics for different choices of potentials and parameters, focusing on the strong back-reaction regime. By implementing sophisticated numerical techniques, these works go beyond the constant-velocity approximation assumed in our discussion above. Their findings indicate that after entering the strong back-reaction regime, the inflation velocity is characterized by oscillations of increasing amplitude, in sharp contrast to analytical studies approximating $\dot{\overline{\varphi }}$ (and $\xi$) as constant. Finally, the recent analytic study reported in [295] fully supports the numerical results. The physical source of this instability seems to be a delayed response of the vector source to changes in the axion velocity ($\dot{\overline{\varphi }}$) when the system enters a strong back-reaction regime. This phenomenon causes the aforementioned oscillations with increasing amplitude in the inflation velocity.

It is not clear at the moment whether these secular effects can be tamed by adding ingredients to the setup discussed in this section. A potential way out is to modify the inflationary dynamics in such a way as to push the strong vector back-reaction regime towards an epoch that does not affect the predictions of PBH formation in the scalar sector. This possibility can be achieved in scenarios in which the gauge field production is localized, leading to a pronounced peak in the scalar power spectrum before entering into a strong back-reaction regime (see Section 4.1.2). Such a way out does not avoid the instability but at least moves it towards epochs that do not affect our predictions for the production and properties of curvature perturbations. In Section 4.1.3, we instead discuss the possibility of avoiding a strong vector back-reaction regime altogether in scenarios in which the axion is a spectator field and does not drive inflation. Appropriate choices of the axion potential nevertheless lead to an enhancement of the scalar sector at a level sufficient for producing PBHs. All the opportunities we review next are possible in scenarios of inflation including axion–vector couplings, which have very rich ramifications and are well studied, given their particle physics motivations.

4.1.2. Bumpy Axion Inflation

In this section, motivated by the instability arguments discussed above, we consider scenarios with localized production of gauge fields able to enhance the curvature spectrum and produce PBHs before entering any strong back-reaction regime. As we will see, the mechanism is based on the ideas discussed in the previous Section 4.1.1 but also uses tools for the activation of the decaying mode that we introduced in Section 3.2.

For definiteness, we focus on a representative example, discussing its properties in the main text and referring to appendices for technical details. We assume that the shift symmetry of the axion is broken both by non-perturbative instanton corrections, as well as by a non-periodic monomial term in the potential29 (see, e.g., [276,278,279,297,298,299,300,301]),

$$V\left(\varphi \right)=\frac{1}{2}{m}^2{\varphi }^2+{\Lambda }^4\frac{\varphi }{f}\sin \left(\frac{\varphi }{f}\right),$$

(78)

where $m,\Lambda ,f$ are parameters of unit mass dimension. In fact, potentials such as those mentioned above characterized by modulations on top of a smooth profile are well-motivated and studied in theoretical constructions of axion inflation models in a variety of situations. For these reasons, we find it interesting to review their useful applications for our aim of producing PBHs.

We are interested in obtaining noticeable modulations of the functional form of the potential within the regime ${\Lambda }^4\lesssim m^2{f}^2$ (which we refer to as the “bumpy” regime in what follows). In particular, we exploit the aforementioned non-perturbative corrections to generate a roller-coaster profile to the otherwise smooth potential, where plateau-like regions are connected to each other by steep cliffs (see Figure 12). While plateau-like regions are suitable to sustain long enough inflation—and to generate nearly scale-invariant fluctuations at CMB scales—when reaching the cliff-like regions, the scalar velocity speeds up intermittently.

During such brief stages, a localized production of gauge fields occurs. To illustrate this phenomenon explicitly, we rely on a numerical procedure along the same lines as that discussed with respect to single-field models of inflation. We numerically solve the background Equation (71) in the bumpy regime of the potential (78), neglecting possible back-reaction effects induced by gauge fields30. The resulting background evolution during inflation is represented in Figure 13 in terms of the first two slow-roll parameters ($ϵ\left(N\right),\eta \left(N\right)$) for the parameter choices used in Figure 12.

Adopting the initial condition ${\varphi }_{\mathrm{in}}\simeq 4.8M_{\mathrm{pl}}$, the scalar field completes a total of 60 e-foldings of its evolution when traversing a single bump-like region, around which the slow-roll parameter ($ϵ$) reaches its maximal value of ${ϵ}_{*}=0.455$ at $N_{*}=35.8$. Following the stage of maximal velocity for the scalar field, the background dynamics enter into a short non-attractor phase with $\eta \lesssim -6$, during which $\varphi$ slows down before the dynamics re-enter into a final inflationary slow-roll stage ($\eta \ll 1$). Since the inflationary background proceeds through an epoch of slow-roll violation, in our calculations, we compute the vacuum contribution to the scalar power spectrum using the numerical methods discussed in Section 3.2 Appendix C. In fact, as anticipated above, in developing this example, we make use both of ideas based on the conversion into curvature fluctuations of gauge-field modes (as in the previous section) and on the activation of the would-be decaying mode of the inflation scalar sector (as discussed in Section 3.2).

Since while traveling the plateau-like regions, the scalar velocity is small, the only contributions to the scalar power spectrum arise from the cliff-like region of the potential (see Figure 12). Around the cliff, the scalar field velocity and the effective coupling ($\xi$) acquire the following profile [86]:

$$\xi \left(N\right)=\frac{g_{\mathrm{cs}}\delta }{1+{\delta }^2{\left(N-N_{*}\right)}^2},⟶{\xi }_{*}\equiv \xi \left(N_{*}\right)=g_{\mathrm{cs}}\delta ,$$

(79)

where ${\xi }_{*}$ denotes the maximal value of the effective coupling. The dimensionless parameter ($\delta \simeq \sqrt{2/3}(M_{\mathrm{pl}}/f)(m/H)$) determines the number of e-folds during which the parameter $\xi$ maintains its maximal value, i.e., $\Delta N\sim {\delta }^{-1}$. For $\delta \propto m$, the duration of such an epoch is inversely proportional to the restoring force provided by the potential; hence, it inversely depends on the mass parameter (m). A time-dependent profile for the effective coupling gives rise to a scale-dependent amplification of the gauge fields whose growing solutions correspond to the real part of the mode functions (see Equation (A61)), using Equation (A62) (see Appendix D.1 for more details).

As expected, the resulting sourced contribution to the power spectrum induced by the $\delta A+\delta A\to \delta \varphi$ inherits the scale dependence of the gauge fields, and the curvature power spectrum can be parameterized as [86],

$${\mathcal{P}}_{\mathcal{R}}\left(k\right)={\mathcal{P}}_{\mathcal{R}}^{\left(\mathrm{v}\right)}\left(k\right)\left[1+\frac{H^2}{64{\pi }^2{M}_{\mathrm{pl}}^2}{f}_{2,\mathcal{R}}\left({\xi }_{*},\frac{k}{k_{*}},\delta \right)\right].$$

(80)

The vacuum contribution (${\mathcal{P}}_{\mathcal{R}}^{\left(\mathrm{v}\right)}$) can be numerically calculated following some of the methods explained in Section 3. Due to the presence of a brief non-attractor phase, it is amplified by the activation of the would-be scalar decaying mode (see Section 3.2). The dimensionless function ($f_{2,\mathcal{R}}$) characterizes the contributions from the gauge fields and has a log-normal shape [86],

$$f_{2,\mathcal{R}}\left({\xi }_{*},\frac{k}{k_{*}},\delta \right)\simeq f_{2,\mathcal{R}}^c\left[{\xi }_{*},\delta \right]\exp \left[-\frac{1}{2{\sigma }_{2,\mathcal{R}}^2\left[{\xi }_{*},\delta \right]}{\ln }^2\left(\frac{k}{k_{*}{x}_{2,\mathcal{R}}^c\left[{\xi }_{*},\delta \right]}\right)\right].$$

(81)

As suggested by the expression (81), the information about the location, width and height of the sourced signals in (80) depends on the motion of $\varphi$ through the step-like features of its wiggly potential, particularly through ${\xi }_{*}$ and $\delta$ dependence of the functions ($x_{2,\mathcal{R}}^c,{\sigma }_{2,\mathcal{R}},f_{2,\mathcal{R}}^c$). We learn from (81) that the sourced signal is maximal at $k=k_{*}{x}_{2,\mathcal{R}}^c$, where it is elevated to $f_{2,\mathcal{R}}^c$, whereas ${\sigma }_{2,\mathcal{R}}$ controls the width of the signal. For the background evolution presented in Figure 13, the fitting functions describing the height, width and location of the sourced signal is calculated in [85] using the semianalytical techniques developed in [302].

Collecting all the results above, we present in Figure 14 the full power spectrum in Equation (80) in terms of the e-folding number by converting the scale dependence to e-fold dependence using the horizon crossing condition ($k=a\left(N\right)H\left(N\right)$). Due to non-trivial background evolution during which the slow-roll conditions are briefly violated, the vacuum contribution to the power spectrum (dotted lines) acquires a non-trivial shape, which shares the characteristic features of the single-field scenarios we discussed in Section 3.2. Notice in the presence of the dip due to the interference between growing and decaying modes of the curvature perturbation that immediately precedes a phase of steady growth of the spectrum.

However, in the present example, the duration of the non-attractor era is not sufficient to provide a prominent peak in the power spectrum solely given the presence of vacuum fluctuations of the curvature perturbation. At this point, the additional source provided by the gauge fields come to the rescue, giving rise to extra-steep growth, leading to the required power to generate a sizable peak at wave numbers $k_{\mathrm{peak}}=k_{*}{x}_{2,\mathcal{R}}^c\simeq 1.5\times {10}^{13}{\mathrm{Mpc}}^{-1}$ corresponding to peak PBH mass at the time of formation $M_{\mathrm{PBH}}\simeq 2.2\times {10}^{-13}{M}_{\odot }$ (see Equation (6)). Hence, in the presence of couplings to the Abelian gauge fields, the roller-coaster architecture of the potential provides an assisted amplification mechanism of the power spectrum. The curvature perturbation sourced by the vector fields and its vacuum counterpart help each another to generate a sufficiently pronounced peak to produce PBHs. Interestingly, the profile of the curvature spectrum as a function of momenta, in particular its growth rate and the details of the peak, are quite different from the corresponding predictions of single-field inflation, as investigated in Section 3. These differences make the two scenarios distinguishable.

In conclusion, the localized production scenario we presented in this section might be considered a possible way to address the instability issues associated with the background of the smooth-axion inflation discussed in Section 4.1.1. In fact, for the model we considered in this section, back-reaction effects become prominent for ${\xi }_{*}>15.6$—see [86] for details—a much larger value than the phenomenological value (${\xi }_{*}\simeq 10.5$) needed for the amplification of the power spectrum. Hence, we make predictions in a safe region of cosmological evolution.

Although this is promising in terms of building workable models of PBH production in axion inflation, the issue regarding the stability of the background is likely to reappear towards the end of bumpy axion inflation—in particular towards its final stages, when $ϵ\to 1$, where $\xi \equiv g_{\mathrm{cs}}\sqrt{ϵ/2}(M_{\mathrm{pl}}/f)\simeq \mathcal{O}\left(10\right)$. In the next section, we discuss an extension of the model presented here, in which a spectator axion sector generates a localized scalar enhancement through its coupling to an Abelian gauge sector. Since in this model, the dynamical evolution of the axion field stops long before inflation ends, we can avoid a strong back-reaction regime that leads to instabilities of the background configuration.

4.1.3. Spectator Axion-Gauge Field Dynamics

We continue to discuss possible methods to locally enhance the spectrum of scalar fluctuations in axion-gauge field models, with an eye to solving the dynamical instability problems mentioned above. We focus on scenarios in which the axion is a spectator field, i.e., does not directly drive inflation. We call $\sigma$ the spectator axion field and couple it to an Abelian gauge sector. The Lagrangian is

$${\mathcal{L}}_{\mathrm{m}}={\mathcal{L}}_{\inf }-\frac{1}{2}{\partial }_{\mu }\sigma {\partial }^{\mu }\sigma -U\left(\sigma \right)-\frac{1}{4}{F}_{\mu \nu }{F}^{\mu \nu }-\frac{g_{\mathrm{cs}}}{4f}\sigma F_{\mu \nu }{\tilde{F}}^{\mu \nu },$$

(82)

where ${\mathcal{L}}_{\inf }=-{(\partial \varphi )}^2/2-V\left(\varphi \right)$ is the inflation Lagrangian, and $V\left(\varphi \right)$ is its potential. We are interested in an inflationary setup whereby the spectator sector provides a subleading contribution to the total energy density during inflation. This implies that the energy densities of the scalar fields in the model (82) satisfy ${\overline{\rho }}_{\sigma }\ll {\overline{\rho }}_{\varphi }$. Assuming small back-reaction from the gauge field fluctuations (${\rho }_A\ll {\overline{\rho }}_{\sigma }$) (more on this later), the Friedmann equation simplifies to

$$3H^2{M}_{\mathrm{pl}}^2\simeq {\overline{\rho }}_{\varphi }+{\overline{\rho }}_{\sigma }⟶3H^2{M}_{\mathrm{pl}}^2\simeq V\left(\overline{\varphi }\right),$$

(83)

so that the inflationary background is completely dictated by the inflation potential. We assume that $V\left(\varphi \right)$ is flat enough to support sufficient inflation, but otherwise, we leave it unspecified, as the fine details of the inflation dynamics are not essential for what follows. Since we are working in a weak back-reaction system, we can avoid instabilities in the inflation sector such as that reviewed in Section 4.1.1.

From now on, we consider representative examples of axion potentials and their corresponding dynamics. The choice of our examples is motivated by their interesting ramifications for PBH production and for their motivations from high-energy physics. If the spectator axion ($\sigma$) is displaced from its global minimum, it can be dynamically active during inflation. Thanks to the perturbative shift symmetry of the axion sector, we expect the axion potential to be nearly flat; hence, the axion dynamics should occur in a slow-roll regime ($\ddot{\overline{\sigma }}\ll 3H\dot{\overline{\sigma }}$). In such a case, we can neglect back-reaction effects induced by gauge field fluctuations31, and the background evolution of the spectator axion can be studied through the following equation,

$$\ddot{\overline{\sigma }}+3H\dot{\overline{\sigma }}+U^{\prime }\left(\overline{\sigma }\right)=0,⟶3H\dot{\overline{\sigma }}+U^{\prime }\left(\overline{\sigma }\right)\simeq 0.$$

(84)

To realize localized gauge field production through the last term in (82), we consider two class of transiently rolling spectator axion models characterized by the following potentials [85,302]

$$U\left(\sigma \right)=\{\begin{array}{cc}{\Lambda }^4\left[1-\cos \left(\frac{\sigma }{f}\right)\right],\hfill & \hfill \left(\mathrm{M}1\right),\\ {\mu }^3\sigma +{\Lambda }^4\left[1-\cos \left(\frac{\sigma }{f}\right)\right]\mathrm{and}{\Lambda }^4\lesssim {\mu }^{3}f\hfill & \hfill \left(\mathrm{M}2\right),\end{array}$$

(85)

where $\mu$ and $\Lambda$ are parameters of mass dimension one.

The first model (M1) features a standard (discrete) shift symmetric potential akin to natural inflation [275], where the size of the axion modulations is set by $\Lambda$. In this model, the axion dynamics span an interval bounded by the maximum ($\sigma =\pi f$) and the minimum ($\sigma =0$) of the potential. Therefore, for large and small field values (early and late times, respectively), the axion rolls with small velocities. However, $\dot{\overline{\sigma }}$ obtains relatively large values at an intermediate time when $\sigma$ passes through an inflection point (${\sigma }_{*}=\sigma \left(t_{*}\right)$) with $U^{{}^{″}}\left({\sigma }_{*}\right)=0$, where the slope of the potential ($U^{\prime }\left(\sigma \right)$) becomes maximal.

In the second model (M2), the axion field range is extended via a linear term [276,278] proportional to a soft symmetry breaking mass parameter ($\mu$). We assume that the axion ($\sigma$) can probe the corresponding bumpy potential32 in the ${\Lambda }^4\lesssim {\mu }^{3}f$ regime, similar to the axion inflation model we discussed in Section 4.1.3. In the plateau-like region(s) and towards the global minimum ($\sigma =0$)33, the spectator axion acquires very small velocities, but obtains a relatively large transient peak when the slope of the potential ($U^{\prime }\left(\sigma \right)$) becomes maximal in the cliff-like region(s), in particular at the inflection point(s) denoted by * in Figure 15.

The background evolution ($\sigma$) can be analytically derived in the slow-roll regime (84), starting from the expressions for the scalar potentials of Equation (85). For typical field ranges dictated by these potentials34, the spectator axion velocity ($\dot{\overline{\sigma }}$) and the effective coupling strength ($\xi =-g_{\mathrm{cs}}\dot{\overline{\sigma }}/\left(2Hf\right)$) obtain a peaked time-dependent profile [85,302]:

$$\xi \left(N\right)=\{\begin{array}{cc}\frac{2{\xi }_{*}}{{\mathrm{e}}^{\delta (N-N_{*})}+{\mathrm{e}}^{\delta (N_{*}-N)}},\delta \equiv \frac{{\Lambda }^4}{3H^2{f}^2}\mathrm{and}{\xi }_{*}\equiv \frac{g_{\mathrm{cs}}\delta }{2}\hfill & \hfill \left(\mathrm{M}1\right)\\ \frac{{\xi }_{*}}{1+{\delta }^2{(N-N_{*})}^2},\delta \equiv \frac{{\mu }^3}{3H^{2}f}\mathrm{and}{\xi }_{*}\equiv g_{\mathrm{cs}}\delta \hfill & \hfill \left(\mathrm{M}2\right)\end{array}$$

(86)

where $N_{*}$ denotes the e-fold number when the axion passes through the inflection point (see Figure 15), and ${\xi }_{*}$ is the maximal value of the effective coupling parameter at $N_{*}$. As in the bumpy axion inflation model, the width of the time-dependent peak in $\xi$ depends on the dimensionless ratio ($\delta$), which characterizes the mass of the spectator axion in its global minimum ($\delta \approx m_{\sigma }^2/H^2$). For a heavier axion (corresponding to larger $\delta$), the restoring force towards the global minimum becomes very relevant; the axion $\sigma$ traverses the inflection point faster, the result being a sharper peak in $\xi$. In other words, $\delta$ controls the acceleration ($\dot{\xi }/\left(\xi H\right)=\ddot{\overline{\sigma }}/\left(\dot{\overline{\sigma }}H\right)\sim \delta$) of $\sigma$ as it rolls down on its potential. Notice that given our slow-roll approximation, we require $\delta \ll 3$.

The peaked structure of the $\xi$ profile controls a critical scale ($k_{*}=a_{*}{H}_{*}$) characterizing the equation of motion (A56) corresponding to the scale of momenta leaving the horizon at an epoch when the axion velocity is maximal (i.e, when $\xi ={\xi }_{*}$). Since the mass of the $\mathrm{U}\left(1\right)$ field in (A56) is as tachyonic as possible around this scale, it gives rise to a scale-dependent growth of the gauge field fluctuations, where only modes with a size comparable to the horizon size at $N=N_{*}$, i.e $k\sim \mathcal{O}\left(1\right)a_{*}{H}_{*}$, are efficiently amplified. In Appendix D.1, we briefly review the parametric amplification of the gauge field modes corresponding to both models (M1,2) we are presenting in this section; for more details, we refer the reader to [85,285,302].

In the spectator axion models we are focusing on, the impact of the vector field production on the visible scalar fluctuations is encoded only indirectly by the presence of gravitational interactions [303]. In fact, although we consider a matter Lagrangian (82) where the spectator axion-gauge field sector is decoupled from the visible inflation sector, when integrating the non-dynamical lapse and shift metric fluctuations, we find a mass mixing between inflation ($\delta \varphi$) and spectator axion ($\delta \sigma$) fluctuations. This phenomenon introduces the possibility that the gauge fields influence the curvature perturbation ($\mathcal{R}\simeq H\delta \varphi /\dot{\overline{\varphi }}$) through a succession of inverse decays: $\delta A+\delta A\to \delta \sigma \to \delta \varphi \propto \mathcal{R}$ (see Appendix E).

Therefore, the dynamics of this sourced contribution can be understood by first studying the influence of particle production on the spectator axion fluctuations ($\delta \sigma$) and then by computing their conversion to curvature perturbation ($\delta \varphi \propto \mathcal{R}$). We sketch some of the steps in Appendix D.3; more detailed computations can be found in [85,302].

Taking into account the vacuum fluctuations within the inflation sector, the total power spectrum of curvature perturbation in the spectator axion-gauge field model can be expressed as [85,302],

$${\mathcal{P}}_{\mathcal{R}}\left(k\right)={\mathcal{P}}_{\mathcal{R}}^{\left(\mathrm{v}\right)}\left(k\right)+{\left[{ϵ}_{\varphi }{\mathcal{P}}_{\mathcal{R}}^{\left(\mathrm{v}\right)}\left(k\right)\right]}^2\sum _i{f}_{2,\mathcal{R}}^{\left(i\right)}\left({\xi }_{*},\frac{k}{k_{*}},\delta \right),$$

(87)

where ${\mathcal{P}}_{\mathcal{R}}=H^2/\left(8{\pi }^2{ϵ}_{\varphi }{M}_{\mathrm{pl}}^2\right)$ with ${ϵ}_{\varphi }\equiv {\dot{\overline{\varphi }}}^2/\left(2H^2{M}_{\mathrm{pl}}^2\right)$ and the sum over i denotes contributions from each particle production location, e.g., for model M1, the sum is over only one of such regions, while for M2, there are two of them (due to our assumptions about the initial conditions of $\sigma$ in the bumpy regime of the potential in (85)), as the spectator rolls along the potential towards its global minimum. Equation (87) teaches us that the part of the spectrum sourced by the gauge fields has an extra slow-roll suppression (${ϵ}_{\varphi }$), in particular with respect to the direct coupling model we discussed in the previous section (see Equation (80)). However, the part that parameterizes the scale dependence of the sourced contribution, i.e., $f_{2,\mathcal{R}}$ in (87), exhibits the same log-normal scale dependence of the bumpy axion model in (81), the peak height, width and location ($f_{2,\mathcal{R}}^c,{\sigma }_{2,\mathcal{R}},x_{2,\mathcal{R}}^c$) of which are calculated in [85,302] in terms of ${\xi }_{*}$ for phenomenologically interesting values of $\delta$ that characterize the background evolution of the spectator axion.

To understand the full-scale dependence of the power spectrum (and, in particular, its vacuum contribution (${\mathcal{P}}_{\mathcal{R}}^{\left(\mathrm{v}\right)}$)), we need to specify the scalar potential ($V\left(\varphi \right)$) in the inflationary sector. Instead of specifying $V\left(\varphi \right)$, we take a phenomenological approach to determine the important set of parameters summarizing the inflationary dynamics. For this purpose, first notice that assuming the effects introduced by the rolling of $\sigma$ is negligible at CMB scales (early times during inflation), we have $n_s-1\simeq 2{\eta }_{\varphi }-6{ϵ}_{\varphi }$ and $r\simeq 16{ϵ}_{\varphi }$. Considering the latest Planck results, we adopt $r\lesssim {10}^{-2}$ at CMB scales [4] to obtain ${ϵ}_{\varphi }\simeq 6.25\times {10}^{-4}$. However, the observed value of the spectral tilt gives $n_s-1\simeq -0.035$ [4]. Furthermore, neglecting higher-order slow-roll parameters, we assume that ${ϵ}_{\varphi }$ remains constant throughout the inflation. These simplifying approximations enable us the describe the total power spectrum in the multifield scenarios we described above. Denoting the e-fold dependence of the Hubble parameter during slow-roll inflation as $H\left(N\right)=H_{\mathrm{cmb}}{e}^{-{ϵ}_{\varphi }N}$ (where the subscript $\mathrm{cmb}$ denotes the time at which the CMB pivot scale ($k_{\mathrm{cmb}}=0.05{\mathrm{Mpc}}^{-1}$) exits the horizon), we can turn the the scale dependence of the vacuum scalar power spectrum into e-fold dependence as

$${\mathcal{P}}_{\mathcal{R}}^{\left(\mathrm{v}\right)}\left(k\right)\simeq {\mathcal{A}}_s{\mathrm{e}}^{-\left(1-{ϵ}_{\varphi }\right)\left(1-n_s\right)N},$$

(88)

where ${\mathcal{A}}_s\equiv {\mathcal{P}}_{\mathcal{R}}^{\left(\mathrm{v}\right)}\left(k_{\mathrm{cmb}}\right)=2.1\times {10}^{-9}$, where ${ϵ}_{\varphi }\simeq 6.25\times {10}^{-4}$ and $1-n_s\simeq 0.035$, as discussed above. Collecting all the information we presented above, we plot in Figure 16 the total scalar power spectrum (Equation (87)) for a representative set of parameter choices describing both of the spectator axion-gauge field models (M1 and M2).

We learn from Figure 16 that since the spectator dynamics do not significantly influence the inflationary background, the vacuum contribution (dotted lines) has a smooth red tilted power law form, which should be contrasted with the bumpy axion inflation model we discussed in the previous section. However, the transient roll of the spectator axion ($\sigma$) before settling to its global minimum triggers gauge field production as it probes the inflection point(s) on its potential (see Figure 15). This phenomenon generates an additional Gaussian bump in the scalar power spectrum. Notice that for larger $\delta$ values, the width of the corresponding bump decreases because $\dot{\overline{\sigma }}\propto \xi$ is maximal for a shorter time interval, during which it can affect fewer gauge field modes. This explains why the sourced signals in the second model (M2) has a narrower width. As in the direct coupling cases with the gauge fields we considered earlier, the peak in the power spectrum originates from a cubic term, i.e., the last term in (82); therefore, the statistics of the curvature perturbation are expected to be non-Gaussian. This explains why the exemplifying scenarios we present in Figure 16 with a peak power of ${\mathcal{P}}_{\mathcal{R}}\sim {10}^{-3}$ are sufficient to generate large populations of PBHs at masses corresponding to the peak scales in Equation (7).

To close this section on axion-vector inflationary systems with the aim of producing PBHs, we can conclude that this possibility requires a fair amount of complex model building. Nevertheless, the physics of axion reached so far have a high degree of theoretical sophistication, thanks to the input and motivations from particle physics, quantum gravity and cosmology. The distinctive predictions we explored with respect to the properties of the curvature power spectrum, in particular its profile as a function of momenta and its shape around the peak, make these scenarios distinguishable from single-field inflation, with relevant ramifications for PBH phenomenology. It is certainly worthwhile to explore and apply existing efforts to explore their consequences with respect to the physics of PBH formation.

4.2. Strong Turns in the Multiscalar Field Space

An alternative approach to enhance scalar fluctuations at small scales is to exploit the dynamics of multifield inflation. Multifield inflation is a topic that has been thoroughly studied in the inflationary literature with a variety of motivations (see e.g., [304] for a review). In reviewing the applications to PBH production, we focus on a general class of non-linear sigma models that are well-motivated when embedding inflation high-energy physics. Although we go through representative concrete examples, the essence of the mechanism is again associated with a transient tachyonic instability in the scalar sector—along a direction orthogonal to the inflationary one—which is converted through suitable couplings into an amplification of curvature perturbations.

The generic two-derivative Lagrangian describing the dynamics of such a system of scalar fields (${\varphi }^I$) minimally coupled to gravity is given by ($S\equiv \int {\mathrm{d}}^{4}x\sqrt{-g}{\mathcal{L}}_m$),

$${\mathcal{L}}_m=-\frac{1}{2}{G}_{IJ}\left(\varphi \right){\partial }_{\mu }{\varphi }^I{\partial }^{\mu }{\varphi }^J-V\left(\varphi \right),$$

(89)

where the fields may interact through the potential ($V\left(\varphi \right)$) and the field field space metric ($G_{IJ}\left(\varphi \right)$). For an FLRW background characterized by the scale factor ($a\left(t\right)$) and Hubble parameter ($H\left(t\right)$), the equation of motion of the homogeneous fields is described by

$${\mathcal{D}}_t{\dot{\overline{\varphi }}}^I+3H{\dot{\overline{\varphi }}}^I+G^{IJ}{V}_{,J}=0,$$

(90)

where the time field space covariant derivative of any field space vector ($A^I$) is defined as ${\mathcal{D}}_t{A}^I={\dot{A}}^I+{\Gamma }_{JK}^I{\dot{\overline{\varphi }}}^J{A}^K$. Considering a two-field case for simplicity, the background trajectory can be split into adiabatic ($e_{\sigma }^I={\dot{\overline{\varphi }}}^I/\dot{\sigma }$) and entropic ($e_s^I$) field bases that are orthogonal to each other. Here, $\sigma \equiv {\left(G_{IJ}{\dot{\overline{\varphi }}}^I{\dot{\overline{\varphi }}}^J\right)}^{1/2}$, which is related to the Hubble slow-roll parameter as $ϵ\equiv -\dot{H}/H^2={\dot{\sigma }}^2/\left(2H^2{M}_{\mathrm{pl}}^2\right)$.

We find it worth mentioning that the vanilla multifield slow-roll trajectory corresponds to fields following standard gradient flow of the potential ($3H{\dot{\overline{\varphi }}}^I\simeq -V^{,I}$) corresponding to a field space trajectory that is approximately geodesic with ${\mathcal{D}}_t{\dot{\overline{\varphi }}}^I\simeq 0$ (see, e.g., (90)). However, a sufficiently long phase of inflation only requires a small acceleration along the unit vector tangent to the inflationary trajectory ($e_{\sigma I}{\mathcal{D}}_t{\dot{\overline{\varphi }}}^I\simeq 0$), and the acceleration pointing in the perpendicular direction—which parameterizes the “bending” of the inflationary trajectory—is generically not constrained and can therefore be large. The level of “bending” can be described by the dimensionless parameter (${\eta }_{\perp }$) so that the orthogonal field bases evolve as

$${\mathcal{D}}_t{e}_{\sigma }^I=H{\eta }_{\perp }{e}_s^I,{\mathcal{D}}_t{e}_s^I=-H{\eta }_{\perp }{e}_{\sigma }^I,$$

(91)

where ${\eta }_{\perp }=-e_s^I{V}_{,I}/\left(H\dot{\sigma }\right)$ measures the acceleration of the trajectory perpendicular to its velocity [305,306]. In combination with the Hubble rate, a “turning rate” associated with the trajectory can be defined as $\Omega \equiv {\eta }_{\perp }H$. Therefore, an inflationary background with a strong turn (${\eta }_{\perp }\gg 1$) satisfies $\Omega \gg H$. As we will show below, this is the regime of interest for the amplification of the scalar fluctuations in this multifield setup.

For this purpose, we consider the linear fluctuations around the background we described above, which is given by the following action [306,307,308]

$$S^{\left(2\right)}=\int {\mathrm{d}}^{4}x{a}^3\left[M_{\mathrm{pl}}^2ϵ\left({\dot{\mathcal{R}}}^2-\frac{{(\partial \mathcal{R})}^2}{a^2}\right)+2\dot{\sigma }{\eta }_{\perp }\dot{\mathcal{R}}{Q}_s+\frac{1}{2}\left({\dot{Q}}_s^2-\frac{{(\partial Q_s)}^2}{a^2}+m_s^2{Q}_s^2\right)\right],$$

(92)

where $\mathcal{R}$ is the comoving curvature perturbation35, and $Q_s$ is the instantaneous entropy perturbation, the mass of which is given by

$$m_s^2=V_{;ss}-H^2{\eta }_{\perp }^2+ϵH^2{M}_{\mathrm{pl}}^2{R}_{\mathrm{fs}}$$

(93)

were $V_{;ss}=e_s^I{e}_s^J{V}_{;IJ}$ is the projection of the covariant Hessian of the potential along the entropic direction, and $R_{\mathrm{fs}}$ is the field space curvature. As should be clear from the action (92), the dynamics of the fluctuations in this two-field model are completely dictated by three functions of time ($N=\ln \left(a\right)$), namely the Hubble rate ($H\left(N\right)$), the turning rate (${\eta }_{\perp }\left(N\right)$) and the mass ($m_s^2\left(N\right)$) of the entropy mode. To demonstrate the amplification in the power spectrum via strong turns in curved multifield space, we follow the effective approach presented in [73] to consider a field space that undergoes a strong turn (${\eta }_{\perp }\gg 1$) around an e-folding number ($N_f$) during inflation corresponding to a time well after CMB scales exited the horizon and assume a featureless Hubble rate ($H\left(N\right)$). We also assume that the time dependence of entropic mass is mainly controlled by the bending parameter, taking $m_s^2=(b-{\eta }_{\perp }^2)H^2$, where b is chosen as a constant for simplicity.

The behavior of the bending parameter (and that of the entropic mass ($m_s^2$)) for typical strong turn cases is shown in Figure 17. An important quantity that determines the behavior of the fluctuations is the duration of the turn, which leads to either broad (left panel) or rapid turn (right panel) cases, as shown in the figure. Although the dynamics of the scalar perturbations are qualitatively different for broad and sudden (strong) turn cases, their behaviors share some common characteristics that can be attributed to the transient nature of the strong turn in field space. First of all, the modes that are still deep in the horizon at the end of the turn do not feel the presence of the feature; hence, their dynamics are of standard single-field type, where the power spectrum is given by

$${\mathcal{P}}_{\mathcal{R}}^{\left(0\right)}\left(k\right)=\frac{H^2}{8{\pi }^2ϵM_{\mathrm{pl}}^2}{|}_{k=aH},$$

(94)

with a slight red-tilted scale dependence. For this purpose, we utilize the phenomenological model we discussed in the previous section to parameterize the scale dependence (time dependence) of the power spectrum as in Equation (88), with a constant $ϵ\simeq 6.25\times {10}^{-4}$ and $1-n_s\simeq 0.035$, along with $H_{\mathrm{cmb}}\simeq {10}^{-5}$. On the other hand, large scales that leave the horizon well before the turn experience typical multifield dynamics, resulting in a transfer of power from entropic to curvature perturbations and enhancement in the power spectrum [304]. The corresponding amplification is sensitive to the value of the entropic mass ($m_s$) from the time of horizon exit until the turn. In the following, we assume a sizable $m_s$ (i.e., b) before the turn so that the entropic fluctuations have decayed sufficiently before the onset of the turn so that the power spectrum is given by the standard result (94) for these scales. For scales that exit the horizon around the time of the turn, the entropic fluctuations exhibit transient (exponential) growth due to their tachyonic mass ($m_s^2<0$) around the turn, which is eventually transferred to the curvature perturbation through the kinetic coupling in (92) proportional to the bending (${\eta }_{\perp }$) parameter. The exponential amplification of ${\mathcal{P}}_{\mathcal{R}}$ associated with these scales is a typical observational feature of strong turns, which we discuss briefly, focusing on broad and sharp turn cases separately.

• Broad turns

For a strong turn in field space that lasts a sufficiently long time so that modes of interest (i.e., modes that are enhanced) have already exited the Hubble radius at the end of the turn, the dynamics of the curvature perturbation can be described by the single-field effective theory with an imaginary sound speed parameterizing the growth of the fluctuations [309,310]. From the perspective of two-field description, imaginary sound speed is just a manifestation of the transient instability induced by the strong non-geodesic motion in field space and the associated transient negative entropic mass ($m_s^2$) [311]. Therefore, the dynamics of mode functions are mainly characterized by two time scales: (i) $\tilde{N}$, corresponding to entropic mass crossing and the onset of instability; and (ii) $\overline{N}$, denoting the sound horizon crossing of fluctuations after which the curvature perturbation ($\mathcal{R}$) becomes frozen:

$$\frac{k}{a\left(\tilde{N}\right)}=\left|m_s\left(\tilde{N}\right)\right|,\&\frac{k\left|c_s\right|}{a\left(\overline{N}\right)}=H\left(\overline{N}\right),$$

(95)

where $|c_s|$ denotes the absolute value of the imaginary speed of sound describing the instability/growth in the fluctuations ($\mathcal{R}$) in the single EFT language. Since the background dynamics change mildly between the time of entropic mass crossing ($\tilde{N}$) and freeze-out ($\overline{N}$) in the broad turn case, the exponential enhancement in the power spectrum can be parameterized as [73,311]

$${\mathcal{P}}_{\mathcal{R}}\left(k\right)={\mathcal{P}}_{\mathcal{R}}^{\left(0\right)}\left(k\right)\exp \left[\pi {\eta }_{\perp }(2-\sqrt{3+b/{\eta }_{\perp }^2})\right]{|}_{{\tilde{N}}_k},$$

(96)

for all k modes that satisfies $k=a\left(\tilde{N}\right)\left|m_s\right|$, while $m_s^2<0$. Notice that exponential enhancement in the power spectrum with respect to the base spectrum is proportional to the bending parameter (${\eta }_{\perp }\gg 1$).

In Figure 18, we show the power spectrum profile as a function of e-folds (left panel) using the analytic Formula (96) for Gaussian bending profiles (${\eta }_{\perp }={\eta }_{\perp }^{\max }{e}^{{(N-N_f)}^2/\left(2{\Delta }^2\right)}$) with $({\eta }_{\perp }^{\max },{\Delta }^2)=(20,10)$ (orange curve, Model 1) and $({\eta }_{\perp }^{\max },{\Delta }^2)=(20,2)$ (black curve, Model 2) and $N_f=30$. For both models, the range of scales and times affected by the turn where $m_s^2<0$ are highlighted by red and correspond to $2.4\times {10}^8<k\left[{\mathrm{Mpc}}^{-1}\right]<1.2{\times }^{15}$ (Model 1) $1.7\times {10}^{10}<k\left[{\mathrm{Mpc}}^{-1}\right]<1.6\times {10}^{13}$ (Model 2). As should be clear from the expression (96), the power spectrum reaches its maximal value for the scale for which entropic mass crossing overlaps with the peak of the bending parameter, i.e., when $N_f={\tilde{N}}_{k_p}$ corresponding to

$$k_p\simeq k_f{|{\left({\eta }_{\perp }^{\max }\right)}^2-b|}^{1/2},$$

(97)

where $k_f=a\left(N_f\right)H\left(N_f\right)\simeq 5.24\times {10}^{11}{\mathrm{Mpc}}^{-1}$ is the scale that exits the horizon at $N_f$. To determine the slope of the power spectrum on the rise, one can compute the spectral index with respect to the base power spectrum to obtain

$$\left(n_s-1\right)-{\left(n_s-1\right)}_0\simeq \pi (2-\sqrt{3})\frac{N_{\mathrm{f}}-\tilde{N}}{{\Delta }^2+N_{\mathrm{f}}-\tilde{N}}{\eta }_{\perp }$$

(98)

which shows that stronger and/or less broad turns in field space lead to a steeper power spectrum. In particular, the expression (98) tells us that for broad and strong turn in field space, the slope of the power spectrum towards its peak can be much larger than the single-field models, where $n_s-1\lesssim 4$ [95,212,217,224].

We would like to note that for scales affected by the turn, the analytic profile presented in (96) is symmetric around the peak; however, numerical computations carried out in [73] show that the power spectrum typically exhibits a much quicker fall of behavior for scales following the peak compared to those preceding it. This can be understood by first recalling that the modes that are enhanced are those that go through entropic mass crossing (95) while the entropic mass is negative ($m_s^2<0$). The right panel of Figure 18 can then guide us for an intuitive understanding for the expected asymmetry in the power spectrum following its peak because a larger range of modes enjoys entropic mass crossing before the peak of $|m_s|$36 than afterwards. We therefore conclude that the expression (96), along with the power spectrum shown in the left panel of Figure 18, only characterizes the behavior of the fluctuations qualitatively, and for more accurate results, numerical methods are required, as shown in [73].

• Sudden turns

When the duration of the turn is shorter (a statement that we will make more precise below), the scales that are maximally enhanced are still those that experience entropic mass crossing during the turn; therefore, they are of order $k\sim k_f{\eta }_{\perp }^{\max }$ (see, e.g., Equation (97)). However, in contrast to the broad turn case, these modes get caught in the (rapid turn) feature while they are still deep inside the horizon, i.e., they satisfy $k>aH$ at the end of the turn. This situation effectively generates an initial excited state for these modes that can be studied analytically in the regime of sudden and strong (constant) turns (${\eta }_{\perp }$) [72,312]. Along with a localized exponential amplification, the resulting power spectrum of curvature perturbation exhibits order-one oscillations in k, which is a common characteristic of sharp features [313,314], with a frequency set by the time of the turn. We will review these features below by focusing on the analytic model presented in [312]37. For this purpose, we consider an inflationary trajectory with a top-hat38sudden turn profile (see, e.g., the right panel in Figure 17) characterized by three parameters: $N_f$, the central time of its location (measured with respect to the horizon exit time of the CMB pivot scale); its duration ($\delta$) and typical (constant) value (${\eta }_{\perp }^{\max }$) around $N_f$:

$${\eta }_{\perp }\left(N\right)={\eta }_{\perp }^{\max }\left[\theta \left(N-\left(N_f-\frac{\delta }{2}\right)\right)-\theta \left(N-\left(N_f+\frac{\delta }{2}\right)\right)\right],$$

(99)

along with an entropic mass profile during the turn that is given by

$$\frac{m_s^2}{H^2}=(\xi -1){\eta }_{\perp }^2\left(N\right)$$

(100)

where $\xi <1$ is a constant parameter. Note from the parameterization of the bending parameter (99) that the sudden turn case we consider corresponds to $\delta >\ln \left({\eta }_{\perp }^{\max }\right)$.

In this setup, one can determine the power spectrum generated by sharp turns by studying the scattering problem of the two-field system using WKB methods [72]. In particular, this procedure includes matching operators (and their derivatives) describing the perturbations from the IN region, where the Bunch–Davies vacuum is assumed, passing through the turn and to the OUT region, resulting in an excited “initial” state as sketched in the right panel of Figure 17. The behavior of the entropic and curvature perturbations in both regions (IN and OUT) are standard, where the fields are dynamically decoupled from each other, while the turn sandwiched between the IN and OUT regions results in their exponential growth. Following this procedure, an analytic expression for the power spectrum was derived in [312] for scales satisfying $k/k_f<\sqrt{1-\xi }{\eta }_{\perp }$ as:

$$\frac{{\mathcal{P}}_{\mathcal{R}}\left(k\right)}{{\mathcal{P}}_{\mathcal{R}}^{\left(0\right)}}=\underset{\mathrm{envelope}}{\underbrace{\frac{e^{2{\eta }_{\perp }\delta S}}{2S^2(1+X+\sqrt{X(1+X)})}}}\times {\sin }^2\left(e^{-\delta /2}\kappa {\eta }_{\perp }+\arctan (\kappa /S)\right)$$

(101)

where

$$S\left(k\right)=\sqrt{\sqrt{4{\kappa }^2+\frac{{(3+\xi )}^2}{4}}-\left({\kappa }^2+\frac{(3+\xi )}{2}\right)},\mathrm{and}X\left(k\right)=\frac{{(3+\xi )}^2}{16{\kappa }^2}$$

(102)

where $\kappa \equiv k/\left(k_f{\eta }_{\perp }\right)$ and ${\eta }_{\perp }$ denote the typical constant value of the bending parameter corresponding to ${\eta }_{\perp }^{\max }$ in the top-hat profile (99).

For the choice of $\xi =-3$—which corresponds to an effectively massless entropy mode on superhorizon scales (see Section 2 of [312] for a detailed discussion on this)—the profile of the curvature power spectrum39 for two representative parameter choices characterizing the turn is shown in Figure 19. As described by the analytic expression (101), the power spectrum features characteristic oscillations (${\sin }^2$ term) modulated by an envelope that is described as highlighted in the equation. The peak value of the envelope (see, e.g., the blue dashed lines in Figure 19) controls the maximal enhancement of the power spectrum with respect to the baseline power spectrum (94) (we assume it to be scale invariant for simplicity) which will occur roughly at the maximum of the function $S\left(k\right)$ defined in (102). Arguably, the oscillatory pattern of the power spectrum around the peak is much more interesting. In particular, the ${\eta }_{\perp }\delta \ll e^{-\delta /2}{\eta }_{\perp }$${\sin }^2$ term changes much faster than the envelope. On the other hand, in a sharp turn regime $e^{-\delta /2}{\eta }_{\perp }\gg 1$ so that the first term in the argument of the ${\sin }^2$ term also changes much faster than the second. These two observations imply that oscillations in the power spectrum are periodic to a considerable degree. The period of the maxima in the oscillations occurs roughly in $\Delta \kappa \approx \pi e^{\delta /2}/{\eta }_{\perp }$, corresponding to a linear frequency of

$$\Delta k\approx \pi e^{\delta /2}{k}_f⟶\omega \equiv \frac{2\pi }{\Delta k}\approx \frac{2e^{-\delta /2}}{k_f}.$$

(103)

Therefore, the scale ($k_f$) that exits the horizon at the center of the feature roughly set the frequency of oscillations in the power spectrum.

In summary, strong sudden turns in field space lead to a power spectrum that exhibits order-one oscillations, the amplitude of which is modulated by an exponentially enhanced envelope (101). The oscillations are fast, making their peaks almost periodic with a frequency given by (103). We note that these oscillations do not have an impact on the PBH mass spectrum ($\beta \left(M\right)$) (see, e.g., (15)) because the calculation of the variance (${\sigma }^2\left(M\right)$) (of the density contrast) involves a smoothing procedure over scales comparable to the width of peaks in the power spectrum (${\mathcal{P}}_{\mathcal{R}}$) [73]. Another issue of interest from the perspective of model building is the influence of non-Gaussianity. The periods of strong turns during inflation are known to generate non-Gaussianity of the flattened type [310,316,317]. The implications of this on the PBH distribution and model building (i.e., the required peak amplitude of ${\mathcal{P}}_{\mathcal{R}}$) is not known and subject to future research. We therefore would like to emphasize that in the strong turn examples we discussed in this section (see Figure 18 and Figure 19), the peak amplitudes we provided are not chosen based on a particular bias on the basis of non-Gaussianity present in these models.

5. Outlook

Primordial black holes (PBHs), if they exist, can shed light on long-standing questions with respect to the nature of dark matter and on the mechanisms driving cosmic inflation. They can provide distinctive sources of gravitational waves that are potentially detectable with current or forthcoming gravitational wave experiments. For these reasons, the physics of PBHs offer promising opportunities for collaboration between cosmologists, astronomers and gravitational wave scientists.

In this review, we focused on theoretical aspects of PBH inflationary model building. We learned that generating PBHs from inflation is difficult but possible. We reviewed conceptual ideas for amplifying the primordial curvature spectrum at a level sufficient to trigger black hole formation. These mechanisms find realizations within single-field inflation through a conversion of pronounced gradients of homogeneous quantities into curvature perturbations or within multifield inflation, where curvature fluctuations are instead amplified through appropriate couplings with additional sectors, which are characterized by tachyonic instabilities. The required tuning of model parameters, or the degree of model sophistication required to realize these ideas can be demanding. However, it is certainly a worthwhile effort, given that experimental probes of PBHs are sensitive to the details of the curvature power spectrum, for example, through the properties of the resulting PBH population or through an induced stochastic gravitational wave background sourced by curvature fluctuations. In fact, different categories of models lead to distinct, potentially distinguishable predictions for the statistics of primordial fluctuations. Hence, a detailed analysis relating theoretical scenarios to cosmological and gravitational wave probes offers new precise tests of inflationary mechanisms complementary to traditional ones associated with the physics of the cosmic microwave background and of the large-scale structures of our universe.

Although much theoretical work has been done so far by the community, much more is needed to further investigate and clarify different aspects of the physics of inflation leading to PBH. At the level of model building, it will be important to clarify and address challenges associated with a severe tuning of model parameters or with the dynamical stability of PBH models requiring non-attractor, non-slow-roll phases of inflationary evolution. It will also be crucial to continue to characterize the rich and subtle properties of primordial fluctuations in PBH models of inflation with large enhancements of the primordial power spectrum at small scales. The to-do list includes a deeper analytic understanding of the properties of the curvature spectrum profile, as well as non-linearities and non-Gaussianities and their consequences for observable quantities. Such theoretical analysis will have relevant ramifications for designing appropriate cosmological probes of the physics of PBHs. Hopefully, a concerted effort of theory and experiments motivated by a deeper understanding of PBH physics will allow us to set new bounds or possibly make new discoveries with respect to the mechanisms driving inflation and the nature of dark matter.