Primordial black holes as a dark matter candidate

By considering the Jeans' length and using Newtonian gravity, reference [25] found that a PBH will form if the density contrast δ ≡ δρ/ρ in the comoving slice ⁴ , evaluated when a given scale enters the horizon exceeds a critical, or threshold, value, δ_c, given by ${\delta }_{\text{c}}={c}_{\text{s}}^{2}$. ⁵ Here, c_s is the sound speed, which is equal to $1/\sqrt{3}$ during radiation domination. The mass of the resulting PBH is of order the horizon mass at that time, M_PBH ∼ M_H, where

$$\begin{equation}{M}_{\text{H}}\sim \frac{t}{G}\sim 1{0}^{15}\enspace \text{g}\left(\frac{t}{1{0}^{-23}\enspace \text{s}}\right).\end{equation} \tag{ 1 }$$

See equation (8) in section 2.1.4 below for a more precise expression for M_H. A PBH formed at around the QCD phase transition (t ∼ 10⁻⁶ s) would have mass of order a Solar mass (M_⊙ = 2 × 10³⁰ g). These initial analytic PBH formation calculations were supported by numerical simulations a few years later [30].

Reference [25] also calculated the initial (i.e. at the time of formation) abundance of PBHs

$$\begin{equation}\beta \left({M}_{\text{H}}\right)\equiv \frac{{\rho }_{\text{PBH}}}{{\rho }_{\text{tot}}}={\int }_{{\delta }_{\text{c}}}^{\infty }P\left(\delta \right)\mathrm{d}\delta \sim \sigma \left({M}_{\text{H}}\right)\mathrm{exp}\left(-\frac{{\delta }_{\text{c}}^{2}}{2{\sigma }^{2}\left({M}_{\text{H}}\right)}\right),\end{equation} \tag{ 2 }$$

where the final step assumes that the probability distribution of primordial density perturbations, P(δ), is Gaussian with variance ${\sigma }^{2}\left({M}_{\text{H}}\right)\ll {\delta }_{\text{c}}^{2}$. See section 2.1.4 below for an expansion on this calculation, including a more precise definition of the mass variance, σ(M_H), in equation (6). Since the abundance of PBHs, β, depends exponentially on the typical size of the fluctuations and the threshold for collapse, uncertainties in these quantities lead to large (potentially orders of magnitude) uncertainty in β.

There has been extensive work on the criterion for PBH formation in the past decade. For a more detailed review see the introduction of reference [31] and section 2.1 of reference [32].

Reference [28] calculated the density threshold analytically finding, in the comoving slice, δ_c = 0.41 for radiation domination. The threshold for collapse depends significantly on the shape of the density perturbation [33]. For a Gaussian density field the shape of rare peaks depends on the power spectrum [34]. Consequently the threshold for PBH formation, and hence their abundance, depends on the form of the primordial power spectrum [35]. Recently reference [31] has studied a wide range of shapes and found that the threshold is lowest (δ_c ≈ 0.41) for broad shapes where pressure gradients play a negligible role and highest (δ_c ≈ 0.67) for peaked shapes where pressure gradients are large. The lower limit agrees well with the analytic estimate of the threshold in reference [28]. Subsequently reference [36] has shown that the criterion for collapse is, to an excellent approximation, universal when expressed in terms of the average of the compaction function (which quantifies the gravitational potential) inside the radius at which it is maximised. If the abundance of PBHs is calculated using peaks theory (see section 2.1.4) rather than using equation (2) (or its more accurate form equation (5) below) then the peak amplitude, rather than the average value, of the perturbation is required and this is also calculated in reference [31].

Deviations from spherical symmetry could in principle affect the threshold for collapse and hence the abundance of PBHs. However the ellipticity of large peaks in a Gaussian random field is small [34, 37], and numerical simulations have recently shown that the effect on the threshold for collapse is negligibly small [38]. While non-Gaussianity of the primordial perturbations can have a significant effect on the abundance of PBHs (see section 2.1.4), its effect on the threshold for collapse alone is also relatively small, of order a few percent [39].

In principle it is also possible to calculate the abundance of PBHs using the primordial curvature perturbation (which is introduced in section 2.3.1) rather than the density contrast. However, as emphasised in reference [40] (see also reference [32]), perturbations on scales larger than the cosmological horizon can not (because of causality) affect whether or not a PBH forms. Consequently care must be taken if using the curvature perturbation to ensure that super-horizon modes do not lead to unphysical results.

It was realised in the late 1990s that, due to near critical gravitational collapse [41], the mass of a PBH depends on the amplitude of the fluctuation from which it forms [42, 43]:

$$\begin{equation}{M}_{\text{PBH}}=\kappa {M}_{\text{H}}{\left(\delta -{\delta }_{\text{c}}\right)}^{\gamma },\end{equation} \tag{ 3 }$$

where the constants γ and κ depend on the shape of the perturbation and the background equation of state [43, 44]. Numerical simulations have verified that this power law scaling of the PBH mass holds down to ${\left(\delta -{\delta }_{\text{c}}\right)}^{\gamma }\sim 1{0}^{-10}$ [45, 46]. For PBHs formed from Mexican hat-shaped perturbations during radiation domination γ = 0.357 and κ = 4.02 [45].

The fraction of the energy density of the Universe contained in regions overdense enough to form PBHs is usually calculated, as in Press–Schechter theory [47], as ⁶

$$\begin{equation}\beta \left({M}_{\text{H}}\right)=2{\int }_{{\delta }_{\text{c}}}^{\infty }\frac{{M}_{\text{PBH}}}{{M}_{\text{H}}}P\left(\delta \left(R\right)\right)\;\mathrm{d}\delta \left(R\right).\end{equation} \tag{ 4 }$$

Assuming that the probability distribution of the smoothed density contrast at horizon crossing, δ(R), is Gaussian with mass variance σ(R) and that all PBHs form at the same time (i.e. at the same value of M_H), with the same mass M_PBH = αM_H, then the PBH MF is monochromatic and

$$\begin{equation}\beta \left({M}_{\text{H}}\right)=\sqrt{\frac{2}{\pi }}\frac{\alpha }{\sigma \left(R\right)}{\int }_{{\delta }_{\text{c}}}^{\infty }\mathrm{exp}\left(-\frac{{\delta }^{2}\left(R\right)}{2{\sigma }^{2}\left(R\right)}\right)\mathrm{d}\delta \left(R\right)=\alpha \enspace \mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(\frac{{\delta }_{\text{c}}}{\sqrt{2}\sigma \left(R\right)}\right).\end{equation} \tag{ 5 }$$

The mass variance is given by [48, 49]

$$\begin{equation}{\sigma }^{2}\left(R\right)=\frac{16}{81}{\int }_{0}^{\infty }{\left(kR\right)}^{4}{W}^{2}\left(kR\right){\mathcal{P}}_{\mathcal{R}}\left(k\right){T}^{2}\left(kR/\sqrt{3}\right)\frac{\mathrm{d}k}{k},\end{equation} \tag{ 6 }$$

where W(kR) is the Fourier transform of the window function used to smooth the density contrast on a comoving scale R, ${\mathcal{P}}_{\mathcal{R}}\left(k\right)$ is the power spectrum of the primordial comoving curvature perturbation (see e.g. reference [50]) and T(y) is the transfer function which describes the evolution of the density perturbations on subhorizon scales:

$$\begin{equation}T\left(y\right)=3\frac{\left(\mathrm{sin}\enspace y-y\enspace \mathrm{cos}\enspace y\right)}{{y}^{3}}.\end{equation} \tag{ 7 }$$

The appropriate window function to use for PBH formation is not known and the relationship between the amplitude of the power spectrum and σ(R) (and hence the abundance of PBHs formed) depends significantly on the choice of window function [51]. For a locally scale-invariant power spectrum with amplitude ${\mathcal{P}}_{\mathcal{R}}\left(k\right)={A}_{\text{PBH}}$, one finds σ²(R) = bA_PBH with b = 1.1, 0.09 and 0.05 for real-space top-hat, Gaussian and k-space top-hat window functions respectively [51]. The horizon mass, M_H, within a comoving radius, R, during radiation domination is given by (cf appendix A of reference [52])

$$\begin{equation}{M}_{\text{H}}=5.6{\times}1{0}^{15}{M}_{\odot }{\left(\frac{{g}_{\star ,\text{i}}}{106.75}\right)}^{-1/6}{\left(R{k}_{0}\right)}^{2}.\end{equation} \tag{ 8 }$$

This expression has been normalised to a fiducial comoving wavenumber, k₀ = 0.05  Mpc⁻¹, corresponding to the cosmic microwave background (CMB) pivot scale and assumes that the initial effective degrees of freedom for entropy and energy density are equal and denoted by g_⋆,i. Excursion set (or peaks) theory [34], which uses the heights of peaks in the density field rather than their averaged value, can also be used to calculate the PBH abundance [40, 53, 54].

In the case of critical collapse, where the mass of a PBH depends on the size of the fluctuation from which it forms (see section 2.1.3), even if all PBHs form at the same time they have an extended MF [42]. While the MF is peaked close to the horizon mass, it has a significant low mass tail.

If the mass of each PBH remains constant and mergers are negligible (see section 3.4 for discussion of the latter) then the PBH density evolves with the scale factor, a, as ρ_PBH ∝ a⁻³. During matter domination the fraction of the total density in the form of PBHs remains constant, while during radiation it grows proportional to a. For a monochromatic MF the present day (at t = t₀) PBH density parameter is given by [55]

$$\begin{equation}{{\Omega}}_{\text{PBH}}\equiv \frac{{\rho }_{\text{PBH}}\left({t}_{0}\right)}{{\rho }_{\text{c}}\left({t}_{0}\right)}=\left(\frac{\beta \left({M}_{\text{H}}\right)}{1.1{\times}1{0}^{-8}}\right){\left(\frac{h}{0.7}\right)}^{-2}{\left(\frac{{g}_{\star ,\text{i}}}{106.75}\right)}^{-1/4}{\left(\frac{{M}_{\text{H}}}{{M}_{\odot }}\right)}^{-1/2},\end{equation} \tag{ 9 }$$

where ρ_c is the critical density for which the geometry of the Universe is flat, h is the dimensionless Hubble constant, H₀ = 100  h km s⁻¹  Mpc⁻¹, and again it is assumed that the initial effective degrees of freedom for entropy and energy density are equal. Accretion of gas onto multi-Solar mass PBHs at late times can increase their mass significantly, and this needs to be taken into account when translating constraints on the present day PBH abundance into constraints on the initial PBH abundance [56]. Constraints arising from this accretion are discussed in section 3.6.

If the primordial power spectrum has a finite width peak, then PBHs can be formed at a range of times and in this case the spread in formation times (or equivalently horizon masses at the time of formation) also needs to be taken into account. Often (e.g. references [17, 52, 57–59]) the PBH MF is calculating by binning the primordial power spectrum by horizon mass, calculating the MF for each bin, and then summing these MFs. The resulting MFs can often be well fit by a lognormal distribution [60, 61]. The accurate calculation of the MF resulting from a broad power spectrum is, however, an outstanding problem. In principle a region which is over-dense when smoothed on a scale R₁ could also be over-dense when smoothed on a scale R₂ > R₁, and hence the original PBH with mass M₁ is then subsumed within a PBH of mass M₂ > M₁. This general situation in structure formation is known as the 'cloud in cloud' problem. Reference [54] argued that for a broad power spectrum the probability that a PBH with mass M₁ is subsumed within a PBH with mass M₂ ≫ M₁ is small. For work towards an accurate calculation of the PBH MF, see e.g. references [62, 63].

During phase transitions the pressure is reduced and consequently the threshold for PBH formation, δ_c, is reduced (e.g. references [25, 28]) and PBHs form more abundantly. In particular (if the primordial power spectrum is close to scale-invariant on small scales) the QCD phase transition leads to enhanced formation of Solar mass PBHs [64–66] and other phase transitions may lead to enhanced formation of PBHs with other masses [67].

It has long been realised that since PBHs form from the high amplitude tail of the density perturbation probability distribution, non-Gaussianity can have a significant effect on their abundance [68, 69]. Reference [70] presents a path-integral formulation for calculating the PBH abundance (in principle) exactly in the presence of non-Gaussianity. In practice this is non-trivial as the resulting expression depends on all of the n-point correlation functions. In many PBH-producing inflation models quantum diffusion is important (see section 2.3) and in this case the high amplitude tail of the probability distribution is exponential rather than Gaussian [71, 72]. It should be noted that even if the underlying curvature perturbations are Gaussian, the non-linear relationship between density and curvature perturbations inevitably renders the distribution of large density perturbations non-Gaussian [73–75].

The rare high peaks in the density field from which PBHs form are close to spherically symmetric. Therefore the torques on the collapsing perturbation, and the resulting angular momentum, are small. Consequently PBHs are formed with dimensionless spin parameters, $a=\vert \mathbf{\text{S}}\vert /\left(G{M}_{\text{PBH}}^{2}\right)$ where S is the spin, of order 0.01 or smaller [76, 77]. Note however that accretion of gas at late times may increase the spin of massive (M_PBH  ≳  30M_⊙) PBHs [78].

As PBHs are discrete objects there are Poissonian fluctuations in their distribution [79]. The initial clustering of PBHs was first studied in references [79, 80], which found that PBHs would be formed in clusters. More recently it has been shown that if the primordial curvature perturbations are Gaussian and have a narrow peak then the PBHs are not initially clustered, beyond Poisson [81–83]. Reference [54] has argued that the initial clustering is also small for broad spectra. However local non-Gaussianity ⁷ can lead to enhanced initial clustering [85–87].

In most models of inflation the accelerated expansion is driven by a scalar field known as the inflaton. The Friedmann equation for the expansion of a Universe dominated by a scalar field ϕ with potential V(ϕ) is

$$\begin{equation}{H}^{2}=\frac{8\pi }{3{m}_{\text{pl}}^{2}}\left[\frac{1}{2}{\dot {\phi }}^{2}+V\left(\phi \right)\right],\end{equation} \tag{ 10 }$$

and the evolution of the scalar field is governed by the Klein–Gordon equation

$$\begin{equation}\ddot {\phi }+3H\dot {\phi }+\frac{\mathrm{d}V}{\mathrm{d}\phi }=0.\end{equation} \tag{ 11 }$$

The dynamics of inflation are often studied using the (potential) slow-roll parameters _V and η_V ⁸

$$\begin{equation}{{\epsilon}}_{V}\equiv \frac{{m}_{\text{pl}}^{2}}{16\pi }{\left(\frac{{V}^{\prime }}{V}\right)}^{2},\qquad {\eta }_{V}\equiv \frac{{m}_{\text{pl}}^{2}}{8\pi }\left(\frac{{V}^{\prime \prime }}{V}\right),\end{equation} \tag{ 12 }$$

where ' denotes derivatives with respect to ϕ. Accelerated expansion occurs when _V ≲ 1. In the slow-roll approximation (SRA), _V and η_V are both much less than one, and the $\ddot {\phi }$ term in the Friedmann equation (equation (10)) and the $\ddot {\phi }$ term in the Klein–Gordon equation (equation (11)) are negligible. In this regime the power spectrum of the primordial curvature perturbation, ${\mathcal{P}}_{\mathcal{R}}\left(k\right)\equiv \left({k}^{3}/2{\pi }^{2}\right)\langle \vert {\mathcal{R}}_{k}\vert \rangle$, is given by

$$\begin{equation}{\mathcal{P}}_{\mathcal{R}}\left(k\right)\approx \frac{8}{3{m}_{\text{pl}}^{4}}\frac{V}{{{\epsilon}}_{V}},\end{equation} \tag{ 13 }$$

where V and _V are to be evaluated when the scale of interest exits the horizon during inflation, k = aH. It is common to use a Taylor expansion of the spectral index to parameterise the primordial power spectrum on cosmological scales (k_cos ≈ (10⁻³–1)  Mpc⁻¹): ⁹

$$\begin{equation}{\mathcal{P}}_{\mathcal{R}}\left(k\right)={A}_{\text{s}}{\left(\frac{k}{{k}_{0}}\right)}^{\left({n}_{\text{s}}\left(k\right)-1\right)},\end{equation} \tag{ 14 }$$

where k₀ is the pivot scale about which the expansion is carried out and

$$\begin{equation}{n}_{\text{s}}\left(k\right)={n}_{\text{s}}{\vert }_{{k}_{0}}+\frac{1}{2}{\left.\frac{\mathrm{d}{n}_{\text{s}}}{\mathrm{d}\mathrm{ln}\enspace k}\right\vert }_{{k}_{0}}\mathrm{ln}\left(\frac{k}{{k}_{0}}\right)+\cdots \enspace .\end{equation} \tag{ 15 }$$

In the SRA, ${n}_{\text{s}}{\vert }_{{k}_{0}}-1=2{\eta }_{V}-6{{\epsilon}}_{V}$. Primordial tensor modes, which manifest as GWs, are also generated by inflation and in the SRA the ratio of the amplitudes of the tensor and scalar power spectra on cosmological scales (known as the 'tensor to scalar ratio') is given by $r\equiv {\mathcal{P}}_{\text{t}}\left({k}_{0}\right)/{\mathcal{P}}_{\mathcal{R}}\left({k}_{0}\right)=16{{\epsilon}}_{V}$.

The amplitude and scale dependence of the primordial perturbations on cosmological scales are now accurately measured. In particular a scale-invariant power spectrum (n_s = 1) is excluded at high confidence. From Planck 2018, combined with other CMB and large scale structure data sets [105]:

$$\begin{equation}\enspace \mathrm{ln}\left(1{0}^{10}\;{A}_{\text{s}}\right)=3.044{\pm}0.0014,\end{equation} \tag{ 16 }$$

$$\begin{equation}{n}_{\text{s}}{\vert }_{0.05\enspace {\text{Mpc}}^{-1}}=0.9668{\pm}0.0037,\end{equation} \tag{ 17 }$$

$$\begin{equation}r{< }0.063,\end{equation} \tag{ 18 }$$

where 1σ errors on measured parameters and a 95% upper confidence limit on r are stated. The upper limit on r leads to a relatively tight constraint on the slope of the potential in the region that corresponds to cosmological scales, _V < 0.0039, and various inflation models are tightly constrained, or excluded (e.g. V(ϕ) ∝ ϕ²) [105]. There are also constraints on smaller scales from spectral distortions of the CMB and induced GWs (see section 3.8 for a more detailed discussion). However the current constraints on the amplitude of the power spectrum on small scales are fairly weak. The COBE/FIRAS limits on spectral distortions require $\mathcal{P}\left(k\right)\lesssim 1{0}^{-4}$ [106, 107] for k ≈ (10–10⁴)  Mpc⁻¹ and pulsar timing array (PTA) limits on GWs require $\mathcal{P}\left(k\right)\lesssim 1{0}^{-2}$ for k ≈ (10⁶–10⁷)  Mpc⁻¹ [59, 108, 109]. However a future PIXIE-like experiment could tighten the spectral distortion constraint to $\mathcal{P}\left(k\right)\lesssim 1{0}^{-8}$ [110] and SKA and LISA will improve the current induced GW constraints over a wide range of smaller scales [59, 108].

On cosmological scales the amplitude of the primordial curvature perturbation power spectrum has been measured to be A_s = 2.1 × 10⁻⁹ [105]. If the power spectrum were completely scale invariant ¹⁰ then using equation (5), and assuming δ_c = 0.5 and σ² = A_s, which is sufficient for a rough calculation, the initial mass fraction of PBHs formed during radiation domination would be β ≈ erfc(7000) ≈ exp[(−7000)²]/7000, i.e. completely negligible. Conversely if all of the DM is in PBHs with M_PBH ∼ M_⊙ then, from equation (9), the initial PBH mass fraction must be β ∼ 10⁻⁸. From equation (5) this requires ${\delta }_{\mathrm{c}}/\left(\sqrt{2}\sigma \right)\sim 4$, and hence the mass variance on the corresponding scale must be σ ∼ 0.1. This requires the amplitude of the primordial power spectrum on this scale to be A_PBH ∼ 0.01, 7 orders of magnitude larger than its measured value on cosmological scales. We note here that since β depends exponentially on the amplitude of the perturbations, fine-tuning is required to achieve an interesting (i.e. neither negligible nor unphysically large, Ω_PBH ≫ 1) abundance of PBHs. Furthermore to produce PBHs of a particular mass the peak in the power spectrum must occur on a specific scale, given by equation (8).

Figure 1 compares the amplitude of the power spectrum required on small scales to form PBH DM (${\mathcal{P}}_{\mathcal{R}}\left(k\right)\sim 1{0}^{-2}$) ¹¹ with the current measurements on cosmological scales from the CMB temperature angular power spectrum [105] and the Lyman-α forest [111], and current and future constraints on smaller scales [59, 108, 110], which were mentioned above ¹² . The steepest growth of the power spectrum which can be achieved in single field inflation [59, 112] (see section 2.3.2 for discussion) is also shown. As discussed in more detail in section 3.8, the current constraints already indirectly exclude M_PBH/M_⊙ ≳ 10³ and 10⁻² ≲ M_PBH/M_⊙ ≲ 1 for PBH DM formed from the collapse of large inflationary density perturbations during radiation domination. Significant improvements in these indirect probes in the future will cover most of the mass range for PBHs formed via this mechanism.

Zoom In Zoom Out Reset image size

As we saw in equation (13), in the SRA the power spectrum is inversely proportional to _V and hence the slope of the potential, so that reducing the slope increases the amplitude of the perturbations. The potential then has to steepen again so that _V becomes greater than 1 and inflation ends. This can be achieved by inserting a plateau or inflection point feature in the potential (see section 2.3.2). Alternatively large perturbations can be generated in multi-field models. In this case typically a different field is responsible for the perturbations on small scales than on cosmological scales. This effectively decouples the constraints on cosmological scales from the requirements for generating large perturbations.

In this subsection we discuss various ways of generating large PBH forming perturbations in single-field inflation models, namely features in the potential, running of the inflaton mass, hilltop models, and the reheating period at the end of inflation.

The possibility of producing PBHs by inserting a plateau in the potential was explored in the 1990s by Ivanov et al [98]. More recently reference [113] showed that a PBH-forming peak in the power spectrum could be produced with a potential possessing an inflection point. Such a potential can be generated from a ratio of polynomials [113] or for a non-minimally coupled scalar field with a quartic potential [114]. Reference [115] showed that if a quintic potential is fine-tuned so that a local minimum is followed by a small maximum (so that the field is slowed down but can still exit this region) a sufficiently large peak in the power spectrum can be produced.

As shown in reference [116], for the power spectrum of canonical single-field models to grow by the required seven orders of magnitude the SRA has to be violated. In the limit where the potential is flat a phase of non-attractor ultra-slow roll (USR) inflation occurs [117, 118], where the evolution of the inflaton (equation (11)) is governed by the expansion rate, rather than the slope of the potential. In this case the standard calculation of the power spectrum is not valid. In particular the curvature perturbations grow rapidly on superhorizon scales, rather than remaining constant, or 'frozen out'. As emphasised by e.g. references [114, 115, 119], in models with an inflection point, or shallow local minimum, a numerical calculation using the Sasaki–Mukhanov equations [120, 121] is required to accurately calculate the position and height of the resulting peak in the power spectrum. Also in this regime quantum diffusion (where quantum kicks of the inflaton field are larger than the classical evolution) occurs and has a significant effect on the probability distribution of the primordial perturbations, and hence the number of PBHs formed [69, 71, 72, 122].

Reference [59] studied the fastest possible growth of the primordial power spectrum that can be achieved in principle in single-field inflation with a period of USR, finding $\mathcal{P}\left(k\right)\propto {k}^{4}$. Reference [112] subsequently showed that $\mathcal{P}\left(k\right)\propto {k}^{5}{\left(\mathrm{log}\enspace k\right)}^{2}$ can be achieved for a specific form of the pre-USR inflationary expansion. These constraints on the growth of the power spectrum are important because to form an interesting abundance of PBHs the amplitude has to grow significantly from its value on cosmological scales, while evading the constraints from spectral distortions on intermediate scales [59], see figure 1.

The running mass inflation model has a potential V(ϕ) = V₀ + (1/2)m²(ϕ)ϕ², where the inflaton mass, m, depends on its value [123, 124]. The resulting power spectrum can grow sufficiently for PBHs to form while satisfying constraints on cosmological scales [125, 126]. However this model is not complete; it relies on a Taylor expansion of the potential around a maximum and does not contain a mechanism for ending inflation (see discussion in e.g. section IV of reference [116]).

Inflation can alternatively be studied using the hierarchy of (Hubble) slow-roll parameters rather than the potential. It is possible to 'design' a functional form for (N), where N is the number of e-foldings of inflation, which satisfies all of the observational constraints and produces PBHs [126]. The corresponding potential has a 'hill-top' form [126, 127], with inflation occurring as the field evolves away from a local maximum, towards a minimum with V(ϕ) ≠ 0. However, as for running mass inflation, an auxiliary mechanism is required to terminate inflation, so this is not a complete single-field inflation model.

The reheating era at the end of inflation, where the inflaton oscillates around a minimum of its potential and decays, may offer a mechanism for generating PBHs [128, 129]. Reference [130] has shown that during oscillations about a parabolic minimum perturbations are enhanced sufficiently by a resonant instability for PBHs to be produced.

In this subsection we discuss some multi-field scenarios which can generate large PBH-producing fluctuations: hybrid inflation, double inflation and a curvaton field. Quantum diffusion is also often important in multi-field models, e.g. references [58, 99, 100, 131].

The most commonly studied two-field model in the context of PBH production is hybrid inflation (e.g. references [58, 99, 100, 132]). In hybrid inflation one of the fields, ϕ, initially slow-rolls while the accelerated expansion is driven by the false-vacuum energy of a second scalar field ψ. At a critical value of ϕ there is a phase transition, with ψ undergoing a waterfall transition to a global minimum and inflation ending. Around the phase transition quantum fluctuations are large, and a spike in the power spectrum on small scales is generated, leading to a large abundance of light PBHs [99, 100]. For some parameter values, however, the waterfall transition can be 'mild' so that there is a 2nd phase of inflation as the ψ field evolves to the minimum of its potential [58]. In this case during the initial stage of the waterfall transition when both fields are important, isocurvature perturbations are generated, leading to a broad peak in the curvature perturbation power spectrum. Perturbations on cosmological scales are generated during the initial phase of inflation and can be consistent with CMB observations.

In double inflation [58, 101, 102, 131] there are two separate periods of inflation, with perturbations on cosmological scales being generated during the first period, and those on small scales during the second ¹³ . Hybrid inflation models with a mild waterfall transition, as discussed above, fall into this class.

A curvaton is a field which is dynamically unimportant during inflation and acquires isocurvature perturbations, with adiabatic perturbations being generated when it later decays [133]. If the inflaton is responsible for the perturbations on cosmological scales, while the curvaton generates small-scale perturbations, it is easier to produce large PBH-producing perturbations than in standard single field models (where the inflaton is responsible for perturbations on all scales) [103, 134]. A specific example is the 'axion-like curvaton' [134].

Stellar microlensing occurs when a compact object with mass in the range 5 × 10⁻¹⁰  M_⊙ ≲ M_CO ≲ 10 M  _⊙ crosses the line of sight to a star, leading to a temporary, achromatic amplification of its flux [184]. The duration of the microlensing event is proportional to ${M}_{\text{CO}}^{1/2}$, therefore the range of masses constrained depends on the cadence of the microlensing survey. The EROS-2 survey of the magellanic clouds (MC) found f_CO ≲ 0.1 for masses in the range 10⁻⁶ ≲ M_CO/M_⊙ ≲ 1. The constraint weakens above M_CO ≈ 1M_⊙, reaching f_CO ≲ 1 for M_CO ≈ 30M_⊙ [12]. A MACHO collaboration search for long duration (>150 days) events places a similar constraint on f_CO, for 1 ≲ M_CO/M_⊙ ≲ 30 [161]. Uncertainties in the density and velocity distribution of the DM have a non-negligible effect on the MC microlensing constraints [185–187], and they would also be changed significantly if the compact objects are clustered [187].

Tighter constraints (f_CO ≲ 10⁻² for M_CO ∼ 10⁻³ M_⊙ weakening to f_CO ≲ 0.1 for M_CO ∼ 10⁻⁵ M_⊙ and 10⁻² M_⊙) have been obtained [162] using the OGLE microlensing survey of the Galactic bulge [188]. The OGLE data also contain 6 ultra-short (0.1–0.3 days) events which could be due to free-floating planets, or PBHs with M_CO ∼ (10⁻⁴–10⁻⁶)M_⊙ and f_CO ∼ 0.01–0.1 [162].

A high cadence optical observation of M31 by Subaru HSC [189] constrains f_CO ≲ 10⁻² for 5 × 10⁻¹⁰ M_⊙ ≲ M_CO ≲ 10⁻⁸ M_⊙ weakening to f_CO ≲ 1 at M_CO ∼ 5 × 10⁻¹² M_⊙ and 5 × 10⁻⁶ M_⊙ [163]. The constraints are weaker than initially found, due to finite source and wave optics effects. For M_PBH ≲ 10⁻¹⁰ M_⊙, the Schwarzschild radius of the PBH is comparable to, or less than, the wavelength of the light and wave optics effects reduce the amplification [190]. Furthermore the stars in M31 that are bright enough for microlensing to be detected are typically larger than assumed in reference [189], further weakening the constraint [182, 191]. There are also significantly weaker constraints for M_CO ≈ (10⁻⁷–10⁻⁹)M_⊙ from a search for low amplification microlensing events in Kepler data [192].

When a background star crosses a caustic in a Galaxy cluster it is magnified by orders of magnitude [193]. Microlensing by stars or other compact objects in the cluster can lead to short periods of further enhanced magnification. On the other hand if a significant fraction of the DM within the cluster is composed of compact objects then the overall magnification is reduced (see e.g. reference [194] and references therein). Icarus/MACS J1149LS1 is the first such microlensing event discovered (serendipitously) and is consistent with microlensing by an intracluster star of a source star at a redshift of 1.5 [195, 196]. This leads to a constraint f_CO < 0.08 for 10⁻⁵ < M_CO/M_⊙ < 10 from the compact object population not reducing the magnification [164]. For more massive compact objects a constraint ${f}_{\text{CO}}{< }0.08{\left({M}_{\text{CO}}/10{M}_{\odot }\right)}^{1/3}$ arises from assuming the microlensing event was caused by a dark compact object rather than a star [164]. Both of these constraints have an uncertainty of order a factor of 2, from the uncertainty in the lens-source transverse velocity.

Microlensing by compact objects in the lens Galaxy leads to variations variation in the brightness of multiple-image quasars [197]. Using optical data from 24 lensed quasars, reference [198] finds that (20 ± 5)% of the mass of the lens Galaxies is in compact objects (including stars) with mass in the range 0.05 < M_CO/M_⊙ < 0.45. This is consistent with the expected stellar component, with only a small contribution allowed from dark compact objects, however no constraint on f_CO is stated ¹⁵ .

The effects of gravitational lensing on the magnification distribution of type 1a supernovae (SNe1a) depend on whether or not the DM is smoothly distributed [200]. If the DM is in compact objects with M_CO ≳ 10⁻² M_⊙, then most SNe1a would be dimmer than if the DM were smoothly distributed, while a few will be significantly magnified [165]. Using the JLA and Union 2.1 SNe samples, reference [165] find a constraint f_CO ≲ 0.4, for all M_CO ≳ 10⁻² M_⊙. They argue that this result is robust to the finite size of SNe and peculiar SNe. The former is contrary to reference [201] who argue that the constraint does not apply for M ≲ 3M_⊙.

Strong gravitational lensing of fast radio bursts (FRBs) by COs with M_CO  ≳  (10–100)M_⊙ would lead to two images, separated by a measurable (ms) time delay [202]. No such signal has been seen in the ∼100 FRBs observed to date, which leads to a constraint f_CO ≲ 0.7 for M_CO ≳ 10³ M_⊙, weakening to f_CO ≲ 1 for M_CO ∼ 10² M_⊙ [203].

Reference [204] proposed that asteroid mass compact objects could be probed by femtolensing of gamma-ray bursts (GRBs), specifically via interference fringes in the frequency spectrum due to the different phases of the 2 (unresolved) images during propagation. Using data from the Fermi GRB monitor, reference [205] placed constraints on compact objects in the mass range 5 × 10¹⁷ g ≲ M_CO < 10²⁰ g. However reference [206] demonstrated that most GRBs are too large to be modelled as point-sources, and furthermore that wave optics effects [207] need to be taken into account. Consequently there are in fact no current constraints on f_CO from femtolensing.

Two-body interactions lead to the kinetic energies of different mass populations within a system becoming more equal. In a system made of stars and more massive compact objects, the stars gain energy and their distribution will expand. Ultra-faint dwarf Galaxies (UFDGs) are particularly sensitive to this effect, due to their high mass to luminosity ratios. Reference [169] found that the observed sizes of UFDGs place a constraint f_CO ≲ 0.002–0.004 for M_CO ∼ 10⁴ M_⊙, weakening with decreasing mass to f_CO ≲ 1 for M_CO ∼ 10M_⊙. The uncertainty comes from uncertainties in the velocity dispersion of the compact objects and the assumed initial radius of the stellar distribution. Tighter constraints may be obtained from the survival of the star cluster near the centre of the dwarf Galaxy Eridanus II, depending on its age [169].

Reference [230] used the projected stellar surface density profile of Segue 1, and in particular the absence of a ring feature, to show that f_CO < 0.06 (0.20) for M_CO = 30 (10)M_⊙ at 99.9% confidence. Subsequent Fokker–Planck simulations of dwarfs composed of stars and compact object DM have not found such a feature however [231], and have also found a slightly slower growth of the stellar component than in reference [169]. Consequently their constraints [231] are slightly weaker than references [169, 230]: f_CO < 1 for M_CO > 14M_⊙. Reference [232] has subsequently shown that collectively UFDGs exclude f_CO = 1 for the mass range (1–100)M_⊙ for both delta-function and lognormal MFs. While some low mass UFDGs have stellar populations that are individually consistent with f_CO ∼ 1 for M_CO ∼ 1M_⊙ (see also reference [231]), most would be 'puffed up' too much.

The energy of wide binary stars is increased by multiple encounters with compact objects, potentially leading to disruption [233]. The separation distribution of wide binaries can therefore be used to constrain the abundance of compact objects [234]. Radial velocity measurements are required to confirm that wide binaries are genuine, as otherwise erroneously tight constraints can be obtained from spurious binaries [235].

Using the 25 wide binaries in their catalogue [236] that spend the least time in the Galactic disk (and are hence least affected by encounters with the stars therein) reference [170] finds f_CO ≲ 0.1 for M_CO ≳ 70M_⊙ with the limit weakening with decreasing mass to f_CO ≲ 1 for M_CO = 3M_⊙. Tighter constraints may be possible using data from Gaia, however radial velocity follow-up will be needed to confirm that candidate binaries are genuine (cf reference [237]).

Radiation emitted due to gas accretion onto PBHs can modify the recombination history of the Universe, and hence affect the anisotropies and spectrum of the CMB [238, 239]. There are significant theoretical uncertainties in the accretion rate and also the ionizing effects of the radiation [240, 241].

Reference [242] argues that, contrary to the spherical accretion assumed in previous work, an accretion disk should form, in which case the resulting constraints are significantly tightened. Formation of (non-PBH) DM halos around PBHs tightens the constraints for M_PBH ≳ 10M_⊙ [171]. For spherical (disk) accretion f_PBH ≲ 1 for M_PBH ∼ 10 (1)M_⊙, tightening with increasing PBH mass to f_PBH < 3 × 10⁻⁹ at M_PBH ∼ 10⁴ M_⊙ [171]. There are also model-dependent constraints on PBHs with M_PBH ≳ 10M_⊙ from their effects on the 21 cm spectrum as measured by EDGES [172].

Accretion of interstellar gas onto M_PBH > M_⊙ PBHs in the MW would lead to observable x-ray and radio emission [243]. Comparing predictions from numerical simulations of gas accretion onto isolated moving compact objects with known x-ray and radio sources in the Chandra and VLA Galactic centre surveys leads to a constraint f_PBH ≲ 10⁻³ for M_PBH ∼ (30–100)M_⊙ [173]. Reference [244] uses the observed number density of compact x-ray objects to place a similar constraint on f_PBH, valid up to M_PBH ∼ 10⁷ M_⊙. Reference [174] places a constraint f_PBH ≲ 10⁻⁴ for M_PBH ∼ 10³ M_⊙, weakening to f_PBH ≲ 1 for M_PBH ∼ M_⊙ and 10⁷ M_⊙, from gas heating in dwarf Galaxies.