Primordial Black Holes as a dark matter candidate

The detection of gravitational waves from mergers of tens of Solar mass black hole binaries has led to a surge in interest in Primordial Black Holes (PBHs) as a dark matter candidate. We aim to provide a (relatively) concise overview of the status of PBHs as a dark matter candidate, circa Summer 2020. First we review the formation of PBHs in the early Universe, focusing mainly on PBHs formed via the collapse of large density perturbations generated by inflation. Then we review the various current and future constraints on the present day abundance of PBHs. We conclude with a discussion of the key open questions in this field.


I. INTRODUCTION
There is convincing evidence from astronomical and cosmological observations that ≈ 85% of the matter in the Universe is in the form of cold, nonbaryonic dark matter (DM), see Ref. [1] for a historical review. The study of Primordial Black Holes (PBHs), black holes formed via the collapse of large overdensities in the early Universe, dates back to the 1960s and 70s [2,3]. It was realised early on that PBHs are a potential DM candidate [3,4]. As they form before matter-radiation equality, PBHs are non-baryonic. While PBHs are thought to evaporate via Hawking radiation [5,6], those with initial mass M PBH 5 × 10 14 g have a lifetime longer than the age of the Universe [7,8]. On cosmological scales PBH DM would behave like particle DM, however on galactic and smaller scales its granularity can have observable consequences.
The MACHO collaboration's 2-year Large Magellanic Cloud microlensing results [9] in the mid-late 1990s generated a wave of interest in PBH DM. They observed significantly more events than expected from known stellar populations. This excess was consistent with roughly half of the Milky Way (MW) halo being in 0.5 M compact objects, with astrophysical compact objects excluded by baryon budget arguments [10]. With subsequent data sets the allowed halo fraction decreased somewhat [11,12] (see Sec. III C 1). However, many of the ideas and models for producing PBH DM date back to this time.
In 2016 LIGO-Virgo announced the discovery of gravitational waves from mergers of tens of Solar mass black holes [13]. The possibility that these BHs could be primordial rather astrophysical [14][15][16], has led to a 2nd, larger, wave of interest in PBH DM. At that time Ref. [17] carried out a comprehensive review of PBH DM, highlighting several mass windows where PBHs could make up all of the DM. Subsequently there have been significant refinements of observational constraints on the abundance of PBHs. New constraints have been proposed, while some existing constraints have been weakened, or even removed completely. There have also been significant developments in the theoretical calculations of PBH formation.
Here we aim to provide a relatively concise overview of the current (Summer 2020) status of PBHs as a dark matter candidate, 1 aimed at readers outside the field. For a comprehensive recent review of constraints on the abundances of PBHs of all masses, with an extensive reference list, see Ref. [20]. Reference [21] provides a recent overview of various potential observational consequences of PBHs, including dark matter. For a detailed review (circa 2018) of observational constraints on non-evaporated PBHs, PBH formation from large perturbations produced by inflation and PBH binaries as a source of gravitational waves see Ref. [22]. Reference [23] covers formation mechanisms while Ref. [24] focuses on future electromagnetic probes of PBHs as dark matter.
In Sec. II, we review the formation of PBHs, focusing mainly on the collapse of large density perturbations produced by inflation. In Sec. III, we overview the various constraints on the present day abundance of PBHs, including potential future constraints. Finally we conclude with a summary of the current status and open questions in Sec. IV. Throughout our intention is not to mention all (or even most) papers published on a topic, but to briefly describe the origin of a calculation, model or constraint and to summarise the current status.

II. PBH FORMATION
The most commonly considered PBH formation mechanism is the collapse, during radiation domination, of large (adiabatic) density perturbations generated by a period of inflation in the very early Universe. We first overview the formation of PBHs via the collapse of density perturbations during radiation domination (Sec. II A) and during matter domination (Sec. II B). We then discuss how large perturbations can be generated by inflation (Sec. II C). Finally in Sec. II D we briefly review other PBH formation mechanisms.

A. Collapse of density perturbations during radiation domination
In Sec. II A 1 we will first review the pioneering calculations by Carr [25] of the threshold for PBH formation and the resulting PBH mass and abundance. We then look at refinements to the calculation of the threshold for PBH formation (Sec. II A 2), the mass of an individual PBH (Sec. II A 3), the PBH abundance, including their mass function (Sec. II A 4) and finally spin and clustering (Sec. II A 5). For a more extensive review of PBH formation from large inflationary perturbations, see Ref. [22]. For a detailed recent study of the relationship between the PBH abundance and the primordial power spectrum, see Ref. [26].

Original calculation
By considering the Jeans' length and using Newtonian gravity, Ref. [25] found that a PBH will form if the density contrast δ ≡ δρ/ρ when a given scale enters the horizon exceeds a critical value, δ c , given by δ c = c 2 s . 2 Here, c s is the sound speed, which is equal to 1/ √ 3 during radiation domination. The mass of the resulting PBH is of order the horizon mass at that time, See Eq. (8) in Sec. II A 4 below for a more precise expression for M H . A PBH formed at around the QCD phase transition (t ∼ 10 −6 s) would have mass of order a Solar mass (M = 2 × 10 30 g). These initial analytic PBH formation calculations were supported by numerical simulations a few years later [28].
Reference [25] also calculated the initial (i.e. at the time of formation) abundance of PBHs where the final step assumes that the probability distribution of primordial density perturbations, P (δ), is Gaussian with typical size σ(M H ) δ c . See Sec. II A 4 below for an expansion on this calculation, including a more precise definition of σ(M H ) in Eq. (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 β.

Threshold for PBH formation
There has been extensive work on the problem of the PBH formation threshold in the past decade. For a more detailed review see the introduction of Ref. [29] and Sec. 2.1. of Ref. [30].
Reference [31] calculated the threshold analytically finding, in the comoving slice, 3 δ c = 0.41 for radiation domination. The threshold for collapse depends significantly on the shape of the density perturbation [32]. For a Gaussian density field the shape of rare peaks depends on the power spectrum [33]. Consequently the threshold for PBH formation, and hence their abundance, depends on the form of the primordial power spectrum [34]. Recently Ref. [29] 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 Ref. [31]. Subsequently Ref. [35] has shown that the threshold 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 Sec. II A 4) rather than using Eq. (2) (or its more accurate form Eq. (5) below) then the peak amplitude, rather than the average value, of the perturbation is required and this is also calculated in Ref. [29].
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 [33,36], and numerical simulations have recently shown that the effect on the threshold for collapse is negligibly small [37]. While non-Gaussianity of the primordial perturbations can have a significant effect on the abundance of PBHs (see Sec. II A 4), its effect on the threshold for collapse alone is also relatively small, of order a few percent [38].
In principle it is also possible to calculate the abundance of PBHs using the primordial curvature perturbation (which is introduced in Sec. II C 1) rather than the density contrast. However, as emphasised in Ref. [39] (see also Ref. [30]), 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 don't lead to unphysical results.
The mass variance is given by [47,48] where W (kR) is the Fourier transform of the window function used to smooth the density contrast on a comoving scale R, P R (k) is the power spectrum of the primordial comoving curvature perturbation (see e.g. Ref. [49]) and T (y) is the transfer function which describes the evolution of the density perturbations on subhorizon scales: 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 [50]. For a locally scale-invariant power spectrum with amplitude P R (k) = A PBH , one finds σ 2 (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 [50]. The horizon mass, M H , within a comoving radius, R, during radiation domination is given by (c.f. Appendix A of Ref. [51]) This expression has been normalised to a fiducial comoving wavenumber, k 0 = 0.05 Mpc −1 , 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 [33], which uses the heights of peaks in the density field rather than their averaged value, can also be used to calculate the PBH abundance [39,52,53].
In the case of critical collapse, where the mass of a PBH depends on the size of the fluctuation from which it forms (see Sec. II A 3), even if all PBHs form at the same time they have an extended mass function (MF) [41]. 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 Sec. III D for discussion of the latter) then the PBH density evolves with the scale factor, a, as ρ PBH ∝ a −3 . 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 mass function the present day (at t = t 0 ) PBH density parameter is given by [54] where ρ c is the critical density for which the geometry of the Universe is flat, h is the dimensionless Hubble constant, H 0 = 100 h km s −1 Mpc −1 , 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 [55]. Constraints arising from this accretion are discussed in Sec. III F. 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. Refs. [17,51,[56][57][58]) the PBH mass function is calculating by binning the primordial power spectrum by horizon mass, calculating the mass function for each bin, and then summing these mass functions. The resulting mass functions can often be well fit by a lognormal distribution [59,60]. 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 1 could also be over-dense when smoothed on a scale R 2 > R 1 , and hence the original PBH with mass M 1 is then subsumed within a PBH of mass M 2 > M 1 . This general situation in structure formation is known as the 'cloud in cloud' problem. Reference [53] argued that for a broad power spectrum the probability that a PBH with mass M 1 is subsumed within a PBH with mass M 2 M 1 is small. For work towards an accurate calculation of the PBH MF, see e.g. Refs. [61,62].
During phase transitions the pressure is reduced and consequently the threshold for PBH formation, δ c , is reduced (e.g. Refs. [25,31]) 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 [63][64][65] and other phase transitions may lead to enhanced formation of PBHs with other masses [66].
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 affect on their abundance [67,68]. Reference [69] presents a pathintegral 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 Sec. II C) and in this case the high amplitude tail of the probability distribution is exponential rather than Gaussian [70,71]. 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 [72][73][74].

Spin and clustering
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 = |S|/(GM 2 PBH ) where S is the spin, of order 0.01 or smaller [75,76]. Note however that accretion of gas at late times may increase the spin of massive (M PBH 30M ) PBHs [77].
As PBHs are discrete objects there are Poissonian fluctuations in their distribution [78]. The initial clustering of PBHs was first studied in Ref. [79], 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 [80][81][82]. Reference [53] has argued that the initial clustering is also small for broad spectra. However local non-Gaussianity 5 can lead to enhanced initial clustering [84][85][86].

B. Collapse of density perturbations during matter domination
It is usually assumed that the evolution of the Universe is radiation dominated from the end of inflation up until matter-radiation equality at t eq = 1.7×10 12 s. However it is possible that there could be a period of matter domination (with equation of state p ≈ 0) prior to Big Bang Nucleosynthesis due to, for instance, long-lived particles dominating the Universe and then decaying (see e.g. Ref. [87][88][89][90] for discussion in the context of PBH formation).
The criteria for PBH formation during matter domination are significantly different than during radiation domination. During matter domination the density contrast grows as δ ∝ a, so that in principle small perturbations can grow sufficiently to form a PBH. However a perturbation has to be (close) to spherical and homogeneous for a PBH to form [87]. The modified expansion history of the Universe also modifies the relationship between the initial PBH mass fraction, β, and the present day PBH density parameter, Ω PBH , Eq. (9).
The fraction of horizon-sized regions which collapse to form a PBH can be written as the product of the fraction of regions which separately satisfy the inhomogeneity and anisotropy criteria: β = β inhom × β aniso [91]. If a perturbation is not sufficiently spherically symmetric it will collapse to form a pancake or cigar. Khlopov and Polnarev originally found β aniso ≈ 0.02σ 5 , where σ is the mass variance as defined in Eq. (6). Reference [92] revisited this calculation numerically. Their result, which is valid for all σ, can be approximated by β aniso ≈ 0.056σ 5 for σ 0.01. Polnarev and Khlopov argued that for a PBH to form, a fluctuation must collapse to within its Schwarzschild radius before a caustic can form at its centre, and found that the fraction of regions which satisfy this criteria is given by β inhom ≈ σ 3/2 [93]. Taking into account the finite propagation speed of information β inhom ≈ 3.7σ 3/2 (for σ 1) [91]. The final result (for σ 1) is β ≈ 0.21σ 13/2 [91], and if σ 0.05, PBHs form more abundantly during matter domination than during radiation domination.
As angular momentum plays a significant role in PBH formation during matter domination, they are formed with large spins: a 0.5, with the exact value depending on the duration of the period of matter domination and also σ [94]. Since PBHs formed during matter domination don't form from the high peaks of the density field, local primordial non-Gaussianity would lead to smaller initial clustering than for formation during radiation domination, however it can still be much larger than the Poisson shot noise [95].

C. Generation of large primordial perturbations by inflation
Inflation is a period of accelerated expansion (ä > 0, where˙denotes a derivative with respect to time) in the early Universe, originally proposed to solve various problems with the standard Big Bang. It also provides a mechanism for generating primordial perturbations, via quantum fluctuations of scalar fields. For a comprehensive and comprehensible introduction to inflation see Baumann's lecture notes [96]. In Sec. II C 1 we review the aspects of inflation that are relevant to PBH formation and briefly discuss the requirements for generating large, PBH-producing perturbations. We then look at PBH producing single (Sec. II C 2) and multi (Sec. II C 3) field inflation models in more detail. A significant number of PBH-producing inflation models have been proposed. We will not attempt a detailed study of all possible models, but instead focus on examples of the types of models that can produce large primordial perturbations from a phenomenological point of view. In many cases the initial ideas (a plateau in the potential of a single field [97], hybrid inflation [98,99], double inflation [100,101] and a spectator field [102]) were proposed in the 1990s, motivated by the excess of LMC microlensing events observed by the MACHO collaboration [9]. In recent years, motivated by the LIGO-Virgo discovery of massive BH binaries, these models have been revisited and refined, taking into account theoretical and observational developments in the intervening decades.

Introduction
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 and the evolution of the scalar field is governed by the Klein-Gordon equation The dynamics of inflation are often studied using the (potential) slow-roll parameters V and η 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φ term in the Friedmann equation (Eq. (10)) and thë φ term in the Klein-Gordon equation (Eq. (11)) are negligible. In this regime the power spectrum of the primordial curvature perturbation, P R (k) ≡ (k 3 /2π 2 ) |R k | , is given by 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 where k 0 is the pivot scale about which the expansion is carried out and In the SRA, n s | k0 −1 = 2η V −6 V . Primordial tensor modes, which manifest as gravitational waves, 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 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 [104]: 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 (φ) ∝ φ 2 ) [104]. There are also constraints on smaller scales from spectral distortions of the CMB and induced gravitational waves (see Sec. III H 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 P(k) 10 −4 [105,106] for k ≈ (10 − 10 4 ) Mpc −1 and Pulsar Timing Array (PTA) limits on gravitational waves require P(k) 10 −2 for k ≈ (10 6 −10 7 ) Mpc −1 [58,107]. However a future PIXIE-like experiment could tighten the spectral distortion constraint to P(k) 10 −8 [108] and SKA and LISA will improve the current induced gravitational wave constraints over a wide range of smaller scales [58].
On cosmological scales the amplitude of the primordial curvature perturbation power spectrum has been measured to be A s = 2.1 × 10 −9 [104]. If the power spectrum were completely scale invariant 8 then using Eq. (5), and assuming δ c = 0.5 and σ 2 = 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) 2 ]/7000, i.e. completely negligible. Conversely if all of the DM is in PBHs with M PBH ∼ M then, from Eq. (9), the initial PBH mass fraction must be β ∼ 10 −8 . From Eq. (5) this requires δ c /( √ 2σ) ∼ 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 FIG. 1. Constraints on the primordial power spectrum PR(k) from the CMB temperature angular power spectrum [104] (red line), from the Lyman-α forest [109] (blue), CMB spectral distortions [105,106] (green) and pulsar timing array limits on gravitational waves [58] (magenta). Potential constraints from a future PIXIE-like spectral distortion experiment [108] (blue) and limits on gravitational waves from SKA and LISA (magenta) [108] are shown as dotted lines. In each case the excluded regions are shaded. The spectral distortion and induced gravitational wave constraints depend on the shape of the primordial power spectrum, a k 4 growth followed by a sharp cut-off has been assumed here [58]. The approximate amplitude, PR(k) ∼ 10 −2 , required to form an interesting number of PBHs is shown as a black line (see text for details). The dotted black line shows the steepest growth possible in a single field inflation model [58,110] (see Sec. II C 2). Adapted from Refs. [58] and [108]. 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 Eq. (8). Figure 1 compares the amplitude of the power spectrum required on small scales to form PBH DM (P R (k) ∼ 10 −2 ) 9 with the current measurements on cosmological scales from the CMB temperature angular power spectrum [104] and the Lyman-α forest [109], and current and future constraints on smaller scales [58,108], which were mentioned above. 10 The steepest growth of the power spectrum which can be achieved in single field inflation [58,110] (see Sec. II C 2 for discussion) is also shown. As discussed in more detail in Sec. III H, the current constraints already indirectly exclude M PBH /M 10 3 and 10 −2 M PBH /M 1 for PBH DM formed from the collapse of large 9 The required amplitude is in fact scale dependent. Lighter PBHs form earlier and therefore their density relative to the total density grows for longer. Consequently the initial PBH density corresponding to a given present day density is smaller, and hence the amplitude of the power spectrum required for PBHs to make up all of the DM is smaller. However this scale dependence is smaller than the uncertainties discussed in Sec. II A, and therefore we do not include it here. 10 To translate from k to M H we use Eq. (8), from Ref. [51], which agrees with the conversion given in Ref. [26]. However we note a discrepancy with Ref. [58].
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. As we saw in Eq. (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 Sec. II C 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.

Single-field models
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. [97]. More recently Ref. [111] 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 [111] or for a non-minimally coupled scalar field with a quartic potential [112]. Reference [113] 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 Ref. [114], 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 [115,116], where the evolution of the inflaton (Eq. (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. Refs. [112,113,117], in models with an inflection point, or shallow local minimum, a numerical calculation using the Sasaki-Mukhanov equations [118,119] 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 [68,70,71,120].
Reference [58] 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 P(k) ∝ k 4 . Reference [110] subsequently showed that P(k) ∝ k 5 (log k) 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 [58], see Fig. 1.
The running mass inflation model has a potential V (φ) = V 0 +(1/2)m 2 (φ)φ 2 , where the inflaton mass, m, depends on its value [121,122]. The resulting power spectrum can grow sufficiently for PBHs to form while satisfying constraints on cosmological scales [123,124]. 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. Sec. IV of Ref. [114]).
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 [124]. The corresponding potential has a 'hill-top' form [124,125], 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 [126,127]. Reference [128] has shown that during oscillations about a parabolic minimum perturbations are enhanced sufficiently by a resonant instability for PBHs to be produced.

Multi-field models
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. Refs. [57,98,99,129].
The most commonly studied two-field model in the context of PBH production is hybrid inflation (e.g. Refs. [57,98,99,130]). 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 [98,99]. 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 [57]. 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 [57,100,101,129] 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. 11 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 [131]. If the inflaton is responsible for the perturbations on cosmological scales, while the curvaton generates small-scale perturbations, it is easier to produce large PBHproducing perturbations than in standard single field models (where the inflaton is responsible for perturbations on all scales) [102,132]. A specific example is the 'axion-like curvaton' [132].

D. Other formation mechanisms
There are a variety of other early Universe processes which can produce large, PBH-forming overdensities. These include bubble collisions, cosmic string loop-or cusp-collapse, domain wall collapse and scalar condensate fragmentation.
First order phase transitions occur through the formation of bubbles. If these bubbles collide, PBHs with mass of order the horizon mass can form [133][134][135]. However a non-negligible abundance of PBHs is only formed if the bubble formation rate is fine-tuned so that bubble collisions occur, but the phase transition doesn't complete instantaneously. Recently Ref. [136] has studied PBH formation from the collapse of bubbles nucleated during inflation.
Cosmic strings are topological defects which may form during phase transitions in the early Universe [137]. A network of cosmic strings is formed which quickly reaches a stable scaling solution, in which loops with size smaller than the Hubble radius are constantly being produced via long string interactions and loop self-intercommutations. The loops oscillate and if a loop contracts under its own tension to become smaller than its Schwarzschild radius a PBH can form [138,139]. The loop collapse probability is independent of time, and the mass of the PBH formed is proportional to the typical loop mass, which is proportional to the horizon mass. Consequently the PBHs formed from loop collapse have an extended mass function of the form dn PBH /dM PBH ∝ M −5/2 PBH [140]. The fraction of loops that collapse, f , is not well known. Numerical simulations [141] have found f = 10 4.9±0.2 (Gµ) 4.1±0.1 for large tensions Gµ ∼ 10 −(2−3) . 12 Critical phenomena also arise in this PBH formation mechanism, with the PBH mass scaling as a power law of the difference between the loop radius and the Schwarzschild radius [143]. It has recently been argued that PBHs would form more abundantly from the collapse of cosmic string cusps [144].
Large closed domain walls, produced during a second order phase transition, can collapse to form PBHs [145,146]. The PBHs have masses that depend on the parameters of the field which undergoes the phase transition, typically with a significant spread, and can be clustered.
A scalar field with a sufficiently flat potential (such as the multiple flat directions found in supersymmetric generalizations of the Standard Model of particle physics) forms a coherent condensate at the end of inflation. This condensate typically fragments into lumps, such as oscillons or Q-balls. These lumps can come to dominate the Universe, and have large density fluctuations which can produce PBHs [147,148]. These PBHs are smaller (compared with the horizon mass) than those formed via the collapse of density perturbations during radiation domination and can have larger spin [148].
Baryogenesis scenarios with spontaneous breaking of charge symmetry during inflation generate high density regions that can collapse to form PBHs after the QCD phase transition [149]. In this case the PBHs formed have a lognormal mass function [149], centered at ∼ 10 M or higher [150].

III. ABUNDANCE CONSTRAINTS
PBH DM has a wide range of potentially observable effects. In this section we review the constraints (ordered roughly by increasing PBH mass) on the abundance of PBHs from evaporation (Sec. III A), interactions with stars (Sec. III B), gravitational lensing (Sec. III C), gravitational waves from mergers of PBH binaries (Sec. III D), dynamical effects (Sec. III E), the consequences of accretion (Sec. III F) and large scale structure (Sec. III G). We then review indirect constraints which apply if PBHs are formed via the collapse of large density perturbations during radiation domination (Sec. III H) and potential future constraints (Sec. III I). For more detailed descriptions of the physics behind the constraints, including key equations, see Ref. [22]. Figure 2 shows all of the current limits discussed in the text, grouped by the type of constraint (evaporation, lensing, gravitational waves, dynamical effects, and accretion), while Fig. 3 provides an overview, showing the envelope of each type of constraint. Where different constraints arise from different assumptions on e.g. modeling and backgrounds, we aim to show the most conservative. Code for plotting the constraints is available online at github.com/bradkav/PBHbounds. We restrict our attention to PBHs with M PBH 10 7 M which could, in principle, constitute the DM halos of small dwarf galaxies. There are various constraints on the abundance of more massive PBHs, for an overview, see Ref. [20]. All limits quoted assume that the PBHs have a delta-function mass function and do not form clusters. We discuss the application of delta-function constraints to extended mass functions in Sec. III J. As discussed in Sec. III D understanding the late time clustering of PBHs is an outstanding challenge. In this section we use 'PBH' to denote limits which apply specifically to PBHs and 'CO' to denote limits which apply to any compact object.

A. Evaporation
PBHs with initial mass M PBH < M ≈ 5 × 10 14 g have completed their evaporation by the present day [7,8]. The emission from slightly more massive PBHs (M < M PBH 10 17 g) is sufficient that limits on their evaporation products can be used to constrain their abundance [172].
From the extragalactic gamma-ray background [17,54], f PBH 2 × 10 −8 (M PBH /M ) (3+ ) , where ∼ 0.1 − 0.4 parameterizes the energy dependence of the observed gamma-ray intensity: I obs ∝ E −(1+ ) γ [54]. This constraint can be tightened by a factor of O(10) by taking into account the contribution of known astrophysical sources, such as blazars [173]. There are also similar constraints from the damping of CMB anisotropies due to energy injection during recombination [151,152] and also from heating of neutral hydrogen, as probed by the EDGES measurements of 21cm absorption [153].
Constraints on the e ± flux from Voyager 1 lead to a similar limit on the contribution of PBHs to the local dark matter density f PBH < 0.001 for M PBH = 10 16 g [154], with the constraint on f PBH varying with M PBH in a similar way to the gamma-ray constraint. Again subtraction of backgrounds, in this case supernova remnants and pulsar wind nebulae, leads to constraints that are tighter by ∼ 2 orders of magnitude [154]. Positrons produced by PBHs will also annihilate and contribute to the flux of the 511 keV line [174]. The SPI/INTEGRAL limits on this line lead to limits on the PBH fraction which are similar to the gamma-ray limit for 10 16 g M PBH 10 17 g [155,156]. There are somewhat tighter constraints (that exclude f PBH = 1 up to M PBH ≈ 2 × 10 17 g) from INTEGRAL measurements of the Galactic diffuse flux in the MeV range [157].

B. Interactions with stars
Asteroid mass PBHs can potentially be constrained by the consequences of their capture by, and transit through, stars [175][176][177][178]. See Ref. [178] for detailed recent calculations and discussion.
As a PBH passes through a star it loses energy by dynamical friction, and may be captured. A captured PBH will sink to the centre of the star and also accrete matter, potentially destroying the star. A large capture probability requires a large DM density and low velocity dispersions. Stellar survival constraints have been applied to globular clusters [175]. However, as emphasised by Ref. [176], (most) globular clusters are not thought to have a high DM density. Moreover, Ref. [178] argues that the survival of stars does not in fact constrain the PBH abundance, but that the disruption of stars may lead to constraints, if the observational signatures are worked out (see Ref. [179] for work in this direction).

Stellar microlensing
Stellar microlensing occurs when a compact object with mass in the range 10 24 g M CO 10 34 g crosses the line of sight to a star, leading to a temporary, achromatic amplification of its flux [180]. The duration of the microlensing event is proportional to M 1/2 CO , 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 10 −6 M CO /M 1, with the constraint weakening with increasing M CO to 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 [158]. Uncertainties in the density and velocity distribution of the dark matter have a non-negligible effect on the MC microlensing constraints [181][182][183], and they would also be changed significantly if the compact objects are clustered [183].
A high cadence optical observation of M31 by Subaru HSC [185] constrains f CO 10 −2 for 10 24 g M CO 5×10 25 g weakening to f CO 1 at M CO ∼ 10 23 g and 10 28 g [160]. The constraints are weaker than initially found, due to finite source and wave optics effects. For M PBH 10 −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 [186]. Furthermore the stars in M31 that are bright enough for microlensing to be detected are typically larger than assumed in Ref. [185], further weakening the constraint [160,178]. There are also significantly weaker constraints for M CO ≈ 10 24−26 g from a search for low amplification microlensing events in Kepler data [187].
When a background star crosses a caustic in a galaxy cluster it is magnified by orders of magnitude [188]. 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. Ref. [189] 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 [190,191]. This leads to a constraint f CO < 0.08 for 10 −5 < M CO /M < 10 from the compact object population not reducing the magnification [161]. For more massive compact objects a constraint f CO < 0.08(M CO /10 M ) 1/3 arises from assuming the microlensing event was caused by a dark compact object rather than a star [161]. Both of these constraints have an uncertainty of order a factor of 2, from the uncertainty in the lens-source transverse velocity.

Quasar microlensing
Microlensing by compact objects in the lens galaxy leads to variations in the brightness of images of multiply lensed quasars [192]. Using optical data from 24 lensed quasars, Ref. [193] 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. 13

Type 1a supernovae
The effects of gravitational lensing on the magnification distribution of type 1a supernovae (SNe1a) depend on whether or not the DM is smoothly distributed [195]. If the DM is in compact objects with M CO 10 −2 M , then most SNe1a would be dimmer than if the DM were smoothly distributed, while a few will be significantly magnified [162]. Using the JLA and Union 2.1 SNe samples, Ref. [162] find a constraint f CO 0.4, for all M CO 10 −2 M . They argue that this result is robust to the finite size of SNe and peculiar SNe. The former is contrary to Ref. [196] who argue that the constraint does not apply for M 3M .

Femtolensing
Reference [197] 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 Gamma-ray Burst Monitor, Ref. [198] placed constraints on compact objects in the mass range 5 × 10 17 g M CO < 10 20 g. However Ref. [199] demonstrated that most GRBs are too large to be modelled as point-sources, and furthermore that wave optics effects [200] need to be taken into account. Consequently there are in fact no current constraints on f CO from femtolensing.

D. Gravitational waves from mergers
In the late 1990s Nakamura et al. [201,202] studied the formation of Solar mass PBH DM binaries in the early Universe, when pairs of PBHs may be close enough to decouple from the Hubble expansion before matter-radiation equality. Three-body interactions would impart a small angular momentum on the PBHs, leading to the formation of highly eccentric binaries. If these binaries survive unaffected to the present day then the gravitational waves (GWs) resulting from their coalescence could be detected by LIGO. 14 In fact the merger rate would be several orders of magnitude larger than measured by LIGO-Virgo [13], which places a tight constraint, f PBH < O(10 −3 ) for 10 M PBH /M 300 [16,163,203,204]. A dedicated LIGO-Virgo search for sub-Solar mass mergers constrains f PBH < O(10 −1 ) down to M PBH ∼ 0.2 M . There are also similar constraints from the stochastic gravitational wave background [165,[205][206][207] of such mergers.
If PBHs don't constitute all of the DM, then during matter domination stellar mass (and more massive) PBHs accrete halos of particle dark matter with a steep density profile [208,209]. 15 These DM mini-halos affect the dynamical evolution of the PBH-binaries [203,204], however this has a relatively small effect on the merger rate and resulting constraints [163].
If PBHs make up a significant fraction of the DM then PBH clusters form not long after matter-radiation equality [79,[212][213][214]. While distant three-body interactions are expected to have little impact on isolated PBH binaries [204,215], three-body interactions in PBH clusters could significantly affect the eccentricity distribution and hence the predicted merger rates [216][217][218]. Merger rates may also be increased (leading to stronger constraints) in scenarios where PBHs are formed with large initial clustering (as discussed in Sec. II A 5) [82,219], though Ref. [220] argue that clustering instead should weaken constraints. The clustering of PBHs today, and in particular the survival of binaries, is an outstanding problem.
E. Dynamical

Dwarf galaxies
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 [166] found that the observed sizes of UFDGs place a constraint f CO 0.002 − 0.004 for M CO ∼ 10 4 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 [166].
Reference [221] 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 dark matter have not found such a feature however [222], and have also found a slightly slower growth of the stellar component than in Ref. [166]. Consequently their constraints [222] are slightly weaker than Refs. [166,221]: f CO < 1 for M CO > 14 M . Reference [223] has subsequently shown that collectively UFDGs exclude f CO = 1 for the mass range (1 − 100)M for both delta-function and lognormal mass functions. While some low mass UFDGs have stellar populations that are individually consistent with f CO ∼ 1 for M CO ∼ 1M (see also Ref. [222]), most would be 'puffed up' too much.

Wide binaries
The energy of wide binary stars is increased by multiple encounters with compact objects, potentially leading to disruption [224]. The separation distribution of wide binaries can therefore be used to constrain the abundance of compact objects [225]. Radial velocity measurements are required to confirm that wide binaries are genuine, as otherwise erroneously tight constraints can be obtained from spurious binaries [226].
Using the 25 wide binaries in their catalogue [227] that spend the least time in the Galactic disk (and are hence least affected by encounters with the stars therein) Ref. [167] 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 (c.f. Ref. [228]).

z > 0
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 [229,230]. There are significant theoretical uncertainties in the accretion rate and also the ionizing effects of the radiation [231,232].
Reference [233] 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) dark matter halos around PBHs tightens the constraints for M PBH 10M [168]. For spherical (disk) accretion f PBH 1 for M PBH ∼ 10 (1) M , tightening with increasing PBH mass to f PBH < 3 × 10 −9 at M PBH ∼ 10 4 M [168]. There are also model-dependent constraints on PBHs with M PBH 10 M from their effects on the 21-cm spectrum as measured by EDGES [169].

Present day
Accretion of interstellar gas onto M PBH > M PBHs in the Milky Way would lead to observable X-ray and radio emission [234]. 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 −3 for M PBH ∼ (30 − 100) M [170]. Reference [235] uses the observed number density of compact X-ray objects to place a similar constraint on f PBH , valid up to M PBH ∼ 10 7 M . Reference [171] places a constraint f PBH 10 −4 for M PBH ∼ 10 3 M , weakening to f PBH 1 for M PBH ∼ M and 10 7 M , from gas heating in dwarf galaxies.

G. Large scale structure
If massive PBHs make up a significant fraction of the DM, then Poisson fluctuations in their number density enhance the matter power spectrum at small scales [78,236], which can be probed by observations of the Lyman-α forest [78]. Using the latest data from MIKE/HIRES and high-resolution hydrodynamical simulations Ref. [237] finds a conservative limit f co (100M /M PBH ).

H. Indirect constraints
In this subsection we look at constraints on the amplitude of large primordial perturbations, which lead to indirect constraints on the abundance of PBHs formed via the collapse of large density perturbations during radiation domination (Sec. II A). These constraints do not apply to PBHs formed via other mechanisms (see Sec. II D). As discussed in Sec. II A, there are large uncertainties in the calculation of the abundance of PBHs formed from a given primordial power spectrum.
First order scalar perturbations generate tensor perturbations at second order [238,239]. If the density perturbations are sufficiently large then the amplitude of the resulting 'scalar induced gravitational waves' (SIGWs) is larger than that of the GWs generated by the primordial tensor perturbations. Constraints on the energy density of stochastic GWs, from e.g. Pulsar Timing Arrays, therefore limit the abundance of PBHs formed via the collapse of large density perturbations [240]. These constraints depend on the shape of the primordial power spectrum, and also the assumed probability distribution of the density perturbations, and are therefore (inflation) model dependent [241][242][243]. Models which produce a broad peak in the primordial power spectrum are most tightly constrained [242,243]. For PBHs forming from large density perturbations during radiation domination, Ref. [58] finds f PBH < 1 for 10 −2 M PBH /M 1. Reference [107] finds, using data from NANOGrav, f PBH < 10 −23 for M PBH = 0.1M and f PBH < 10 −6 for 0.002 < M PBH /M < 0.4. However this calculation makes approximations which have a huge effect on the constraint on f PBH (including setting the PBH formation threshold equal to unity, and σ 2 = A). There are also tight constraints on the abundance of light, M PBH ∼ 10 13−15 g, PBHs from limits on SIGWs from LIGO [244]. Such light PBHs are expected to have evaporated by the present day, however if Hawking evaporation is not realised in nature they would be stable and otherwise viable as DM.
The amplitude of the primordial density perturbations can also be constrained by the CMB spectral distortions caused by the dissipation of large perturbations [245]. Limits on CMB spectral distortions from COBE/FIRAS exclude f PBH = 1 for PBHs with M PBH /M 10 3 formed from large density perturbations with a Gaussian distribution [246]. 16 Masses down to M ∼ 10 3 M can similarly be excluded via the effects of dissipation on Big Bang Nucleosynthesis [247].

I. Future constraints
In this subsection we discuss potential future constraints on PBH DM. We start with direct constraints in (roughly) increasing order of PBH mass probed, followed by indirect constraints.
Planned space observatories, such as e-ASTROGAM and ASTRO-H, will allow a lower and more precise measurement of the flux of the isotropic gamma-ray and X-ray backgrounds. This will allow improved constraints to be placed on PBHs in the mass range M PBH ∼ 10 16−18 g via the products of their evaporation [173].
A small subset of GRBs with fast variability have small sizes and are therefore suitable targets for femtolensing [199]. A sample of 100 such GRBs with well-measured redshifts could be used to probe PBHs with 10 17 g M CO 10 19 g. Proposals to measure the lensing parallax of GRBs -the relative brightness measured by multiple telescopes at large spatial separations -suggest that this approach could be sensitive to the entire unconstrained range M PBH ∼ 10 17−23 g [248,249]. Very high cadence observations of white dwarfs in the LMC, by e.g. LSST, could reduce the minimum mass probed by microlensing observations by a factor of a few [186]. However, as discussed in detail in Ref. [178], diffraction and finite source size effects make it difficult to extend the range of masses probed by optical microlensing below M CO ∼ 10 22 g. Microlensing of X-ray pulsars, which have small source sizes, can avoid these restrictions and long observations by X-ray telescopes with large effective areas (e.g. AstroSat, LOFT) of SMC X-1 and other X-ray binaries could probe the range 10 18 g M CO 10 22 g [250].
Per-cent level constraints on f CO for 10 −4 M CO /M 0.1 could be achieved from future Fast Radio Burst (FRB) detections, via the phase difference they produce between unresolved images (as in femtolensing) [251]. Pulsar timing arrays can detect the gravitational redshift and acceleration induced by passing compact objects [252,253]. Planned sub-Hertz gravitational wave observatories such as LISA and DECIGO will be sensitive to extreme mass ratio binaries, composed of astrophysical supermassive black holes and compact objects in the range 10 −6 M CO /M 1, down to the level of around f CO ∼ 10 −3 [254,255]. Via various types of searches, SKA will be able to probe f PBH ≈ 1 over the entire mass range (10 −12 − 100) M , with potentially sub-percent level constraints in some mass regions [253].
Detailed observations of caustic-crossing events (combined with improved theoretical modelling) could place very tight constraints on compact objects with planetary, stellar and larger masses [189,191]. Constraints on frequencydependent gravitational lensing dispersions by future gravitational wave detectors could probe PBHs with M 0.1M via their effects on the matter power spectrum [256]. Measurements of the 21cm power spectrum by HERA and SKA will potentially improve cosmological accretion constraints on PBHs with M PBH > M by an order of magnitude [257]. Accurate astrometric surveys can probe compact objects via the time-dependent weak lensing of stars ("astrometric microlensing") [258]. By the end of its lifetime Gaia could place sub per-cent level constraints on f CO for M CO > 10 M via the large anomalous angular velocities and accelerations produced by close encounters [259]. Similar constraints could be placed on stellar and planetary mass COs by Gaia via non-repeating proper motion anomalies (dubbed 'blips'), with a future survey such as Theia capable of even tighter constraints [259]. Strong gravitational lensing of FRBs by COs with M CO ≈ (20−100)M would lead to two images, separating by a measurable (ms) time delay [260]. Upcoming experiments, such as CHIME and SKA, should lead to per-cent level constraints on f CO in this mass range. Strong lensing of GRBs by compact objects with 10 M CO /M 1000M leads to superimposed images which could be detected by a future GRB observatory, leading to per-cent level constraints [261]. Gravitational lensing of gravitational waves by compact objects with M CO /M 5M would produce fringes which could lead to sub per-cent level constraints from current gravitational wave detectors [262] and 3rd generation experiments like the Einstein Telescope [263].
Future space-based laser interferometers will be able to indirectly constrain PBHs produced by the collapse of large density perturbations in the mass range 10 20−26 g, via induced gravitational waves [240,264]. A PIXIE-like CMB spectral distortion experiment could similarly constrain M PBH M [108].
J. Application of constraints to extended mass functions The constraints described above are typically calculated assuming a delta-function (or monochromatic) mass function. However, as discussed in Sec. II, in many cases PBHs are expected to be produced with an extended mass function.
Carr et al. [265] devised a method for applying constraints calculated assuming a delta-function mass function to specific extended mass functions, without explicitly recalculating the constraint from scratch. The dependence of a given astrophysical observable, A[ψ(M )], on the mass function, ψ(M ), defined so that f PBH = ψ(M ) dM can be expanded as where A 0 is the background contribution and the functions K j (M ) encode how the underlying physics of the observable depend on the PBH mass function. In many cases observations place a bound on a single observable: A[ψ(M )] < A exp . For instance in the case of stellar microlensing the observable is the number of events. Also, for many constraints, PBHs with different masses contribute independently to the constraint so that K j (M ) = 0 for j ≥ 2. In this case the constraint for a delta-function mass function as a function of mass, f max (M ), can be translated into a constraint on a specified extended mass function using: This procedure has to be implemented separately for each observable, and different constraints combined in quadrature.
As emphasised in Ref. [265] there are some caveats in the application of this method. The PBH MF can evolve, due to mergers and accretion, so that the the initial MF is not the same as the MF at the time the constraint applies. For some constraints, e.g. GWs from mergers, the higher order terms, K j (M ) for j ≥ 2, are not zero and a detailed calculation is required for each mass function [206,266]. In some cases, such as the effects of PBH accretion on the CMB, the observable is not a single quantity, and if this is not taken into account artificially tight constraints will be obtained. This method has also been expounded on by Bellomo et al. [267]. They showed that for any specific mass function, for each observable the effects are equivalent to a delta-function MF with a particular 'equivalent mass'. They also emphasised that when considering MFs with extended tails, for instance a lognormal distribution, care should be taken not to apply constraints beyond the limits of their validity.
When applied to extended mass functions the constraints are effectively 'smeared out' [265]. Consequently, when multiple constraints are considered, extended mass functions are more tightly constrained than the delta-function mass function, and small mass windows between constraints are closed [59,265] IV. SUMMARY The LIGO-Virgo discovery of gravitational waves from multi-Solar mass black holes has led to a resurgence of interest in primordial black holes (PBHs) as a dark matter candidate. Consequently there have been significant improvements in the theoretical calculations of PBH formation and also the observational constraints on their abundance. In this final section we summarise the current status and highlight key open questions.
The most popular PBH formation mechanism is the collapse of large density perturbations, generated by a period of inflation in the early Universe, during radiation domination. There have been significant recent improvements in our understanding of the threshold for collapse, δ c , and its dependence on the shape of the perturbations (see Sec. II A 2). To produce an interesting number of PBHs the primordial power spectrum must grow by ∼ 7 orders of magnitude from its measured value on cosmological scales. This can be achieved in single-field models with an inflection point or a shallow local minimum in the potential (Sec. II C 2), or in some multi-field models (Sec. II C 3), and there has been significant recent work revisiting these models and refining calculations. The standard calculation of the abundance and mass function of PBHs (Sec. II A 4) assumes that the primordial density perturbations have a Gaussian probability distribution. However this assumption is not valid for large, PBH-producing perturbations. An accurate calculation of the non-Gaussian probability distribution, which for many models needs to take into account quantum diffusion, is an open issue. Detailed calculations of the mass function and initial clustering of PBHs produced by broad power spectra are also an outstanding problem.
There are a wide range of different observational constraints. The lightest stable PBHs are constrained by the products of the early stages of their evaporation (Sec. III A), planetary and stellar mass and heavier PBHs are constrained by gravitational lensing observations (Sec. III C), with stellar mass PBHs also being constrained by their dynamical effects (Sec. III E), the consequences of accretion onto them (Sec. III F) and gravitational waves from their mergers (Sec. III D). The abundance of PBHs formed by the collapse of large density perturbations is also constrained indirectly via constraints on the amplitude of the power spectrum (Sec. III H). In the past few years there has been significant activity on PBH abundance constraints, with new constraints being proposed and existing constraints being revisited. In some cases (e.g. interactions with stars, Sec. III B, or GRB femtolensing, Sec. III C 4) 'old' constraints have been shown not to hold. While individual constraints may have significant modelling uncertainties the Solar mass region is now subject to multiple constraints. If the constraints are taken at face value, PBHs with masses in the planetary to multi-Solar mass range can only make up a subdominant fraction of the DM. However the robustness of this conclusion depends on the late-time clustering of the PBH population, which remains unclear. The asteroid mass region (10 17 g M PBH 10 21 g) remains open. This largely reflects the difficulty of detecting such light compact objects.
Primordial Black Holes, in particular asteroid mass ones, remain a viable dark matter candidate. Further improvements in the theoretical calculations of the production and evolution of PBHs are required to reliably predict the abundance and properties of PBHs from a given model. Even so, it seems clear that a cosmologically interesting number of PBHs can only be produced in specific models of the early Universe, and often fine-tuning is required. However, whether or not PBHs are a significant component of the DM is a question that has to be answered observationally. Novel ideas are needed here to either detect or rule out the remaining open parameter space.