A Brief Review on Primordial Black Holes as Dark Matter

Primordial black holes (PBHs) represent a natural candidate for one of the components of the dark matter (DM) in the Universe. In this review, we shall discuss the basics of their formation, abundance and signatures. Some of their characteristic signals are examined, such as the emission of particles due to Hawking evaporation and the accretion of the surrounding matter, effects which could leave an impact in the evolution of the Universe and the formation of structures. The most relevant probes capable of constraining their masses and population are discussed.


INTRODUCTION
The hypothesis of the formation of black holes (BHs) in the early Universe was first suggested in 1967 by Zeldovich and Novikov (Zel'dovich and Novikov, 1967), and independently by Hawking in 1971(Hawking, 1971. Soon after, the possibility that primordial black holes (PBHs) could account for at least part of the DM became obvious (Chapline, 1975;Meszaros, 1975). At that time, the DM question started to be outlined as one of the fundamental problems in cosmology [see, e.g., the reviews in Sanders (2010); Bertone and Hooper (2018)]. DM (partly) composed by PBHs constitutes an exciting possibility, presenting an enormous number of observable signatures which can constrain its parameter space, as shall be detailed in Sec. 6. The variety of phenomenological effects produced by PBHs allows placing stringent bounds on the abundance of PBHs, usually indicated by the energy fraction of DM as PBHs, f PBH Ω PBH /Ω DM , with Ω PBH and Ω DM the ratios of energy densities and critical density. Moreover, since PBHs are usually expected to be formed before nucleosynthesis, BBN constraints on the baryon abundance do not apply to them, as they do not intervene in the nucleosynthesis of elements, and thus, can be regarded as non-baryonic DM (Carr and Kühnel, 2020). Several recent reviews are devoted to discuss PBHs in great detail (see, e.g., Sasaki et al., 2018;Carr and Kühnel, 2020;Green and Kavanagh, 2021).
Shortly after the first detection of gravitational waves from a merger of ∼ 30 M ⊙ BHs by LIGO (Abbott et al., 2016), the question whether these could be of primordial nature was raised (Bird et al., 2016). Analysis of posterior data from the gravitational wave detectors LIGO and Virgo showed that the detected mergers are compatible with the hypothesis of their components being of primordial nature, although there is no strong preference over stellar BHs (Sasaki et al., 2016;Clesse and García-Bellido, 2017;Clesse and García-Bellido, 2018;Garcia-Bellido et al., 2020;De Luca et al., 2020b;De Luca et al., 2021;Wong et al., 2021). PBHs with a lognormal mass function have been claimed to better fit data than BHs from astrophysical origin (Dolgov et al., 2020), although this is in contrast to the results of Hall et al. (2020), and a mixed population seems compatible or even favored (Hütsi et al., 2021).
Unlike stellar BHs, formed from the collapse of a massive star, which can present masses only above ∼ 3 M ⊙ [the Tolman-Oppenheimer-Volkoff limit (Oppenheimer and Volkoff, 1939;Tolman, 1939)], PBHs could be produced with any mass. Thus, a positive measurement of a BH with a mass lower than ∼ 3 M ⊙ would be a confirmation of the existence of primordial, non-stellar BHs (Clesse and García-Bellido, 2018). PBHs could also conform intermediate mass BHs, with masses between ∼ 10 2 M ⊙ and 10 5 M ⊙ , too massive to be originated from a single star. It is the case of the merger event of BHs with masses ∼ 60 M ⊙ and ∼ 80 M ⊙ , producing a remnant BH of ∼ 150 M ⊙ , in a so far mostly unobserved range of masses (Abbott et al., 2020). Furthermore, PBHs could constitute the seeds for super massive black holes (SMBHs), present in the nuclei of most galaxies, with masses ranging from 10 5 M ⊙ to 10 10 M ⊙ , and already existing at redshifts z > 6 ( Carr and Silk, 2018). Such massive objects can hardly be produced from accreting stellar remnant BHs (Volonteri, 2010). However, the existence of massive enough PBHs may act as seeds for the SMBHs, from which they could have grown by accretion.
PBH formation is already present in standard cosmologies, although extremely unlikely. However, their production usually requires some exotic inflationary scenarios or physics beyond the Standard Model (BSM) in order to obtain a large enough abundance. The typically considered formation mechanism of PBHs arises from the direct collapse of primordial fluctuations, whose power is enhanced at small scales as a consequence of some inflationary potential, as we shall further comment on below. There are, however, other scenarios which naturally predict a population of PBHs as a result of phase transitions in the early Universe and by the collapse of topological defects (Hawking et al., 1982;Hawking, 1989;Polnarev and Zembowicz, 1991;Garriga et al., 2016;. Hence, the existence of PBHs would provide valuable hints about the still unknown physics of the very early Universe (see, e.g., Polnarev and Khlopov, 1985;Khlopov, 2010;Belotsky et al., 2014), and may allow to probe high-energy scales and supersymmetric theories (Ketov and Khlopov, 2019).
In this review, the most relevant aspects of PBHs as DM are briefly discussed, such as the mechanism of formation in Section 2, the initial abundance and mass distribution in Section 3, the process of accretion in Section 4 and other PBH features, which could leave imprints on different observables, in Section 5. Current observational constraints on their population are summarized in Section 6, and we conclude in Section 7.

FORMATION AND CONDITIONS OF COLLAPSE
The mass of a BH which collapsed in the early Universe depends on its formation time. A BH can be characterized by an extremely dense amount of matter in a very compact region, i.e., within the wellknown Schwarzschild radius, R S 2GM PBH /c 2 ∼ 3M PBH /M ⊙ km. Thus, the mean density inside that region can be estimated as ρ S M PBH /(4πR 3 S /3) ∼ 10 18 (M/M ⊙ ) − 2 g cm −3 . On the other hand, the mean density of the universe in the radiation era scales with time as ρ c ∼ 10 6 (t/s) − 2 g cm −3 . In order to have PBH formation, densities at least of the order of the mean inside the BH horizon, ρ c ∼ ρ S , are needed. Therefore, the mass of the resulting PBHs should be of the order of the horizon mass at that time, i.e., the mass within a region of the size of the Hubble horizon, M PBH ∼ M H (Carr, 1975), which is defined as PBHs with masses of ∼ M ⊙ x2 × 10 33 g would have been formed at around the QCD phase transition, at t ∼ 10 − 6 s. Since the PBH mass is roughly given by the mass within the horizon, it means that fluctuations entering the horizon can collapse into PBHs. A detailed calculation shows that M PBH cM H , where the proportionality factor c depends on the details of gravitational collapse, and gets values lower than 1. Early estimates showed that it can be approximated as cxc 3/2 s x0.2, with c s 1/3 the sound speed at the radiation epoch (Carr, 1975) More refined results show that the PBH mass is given by a scaling relation with the overdensity δ, M PBH κ(δ − δ c ) α (Niemeyer and Jedamzik, 1998;Niemeyer and Jedamzik, 1999), where κ and α are order unity constants, which depend on the background cosmology and on the shape of the perturbation (Niemeyer and Jedamzik, 1998;Niemeyer and Jedamzik, 1999;Musco et al., 2005), and δ c xc 2 s is the collapse threshold [see, e.g., Escrivà et al., 2021;Musco et al., 2021 for recent accurate computations]. Figure 1 depicts a sketch of the process of PBH formation.
Since PBHs are formed when fluctuations cross the horizon by the time of formation, t f , their mass can be related to the wavelength of perturbations. When the mode of wavenumber k crosses the horizon, the condition a(t f ) H(t f ) k holds. Since the mass of PBHs is proportional to the horizon mass at the moment of formation, M PBH ∝ cH −1 , at the radiation dominated era (Sasaki et al., 2018), Hence, probing a given scale k could constrain a PBH population of its corresponding mass. Furthermore, an enhancement in the power spectrum around that scale would result in a large number of PBHs of such masses. In this review, we shall only consider PBHs formed during the radiation era. Those produced before inflation ends would have been diluted due to their negligible density during the inflationary accelerated expansion. PBHs formed during the matter-dominated era, or in an early matter domination era previous to the radiation era, have also been considered in the literature, and may have different imprints, since the conditions of collapse are less restrictive, and could start from smaller inhomogeneities (see, e.g., Green and Kavanagh, 2021).
variance σ 2 (M) at a mass scale M, the initial abundance, defined as β(M H ) ρ PBH (t f )/ρ tot (t f ), is given by (Sasaki et al., 2018;Green and Kavanagh, 2021) where in the last equality, σ(M H ) ≪ δ c has been assumed. For the standard cosmological scenario with an initial scale invariant power spectrum, at CMB scales, the amplitude of the fluctuations is around σ(M PBH ) ∼ 10 − 5 , leading to β ∼ 10 − 5 exp(−10 10 ), which is completely negligible (Green, 2014). Therefore, in order to have a relevant population of PBHs, larger values of the initial power spectrum are needed. On the other hand, the assumption of a Gaussian distribution may not be consistent with enhanced fluctuations and the presence of PBHs, so deviations are unavoidable (Franciolini et al., 2018;De Luca et al., 2019b) [except for specific inflation models presenting an inflection point Atal and Germani (2019)]. Non-gaussianities could have a great impact on the initial fraction and lead to a larger population, as well as leaving detectable signatures in gravitational waves (Cai et al., 2019). Finally, note that the above formula follows the simple Press-Schechter approach, whereas to account for the non-universal nature of the threshold, the use of peak theory provides more accurate results (Green et al., 2004;Germani and Musco, 2019). Its validity, however, is limited to relatively narrow power spectra, while broader spectra require the use of non-linear statistics (Germani and Sheth, 2020). Nonetheless, although the initial fraction β is a very small quantity, since matter and radiation densities scale differently with redshift [ ∝ (1 + z) 3 and ∝ (1 + z) 4 respectively], the PBH contribution can become relevant at current times. The fraction β(M) can be related to the current density parameters of PBH and radiation, Ω PBH and Ω c , as (Carr et al., 2010), For initial fractions as low as β ∼ 10 − 8 of solar mass BHs, the fraction of energy in PBHs could be, thus, of order unity.
Depending on the specific mechanism of formation, a population of PBHs with different masses could be generated. The specific shape of the enhancement of fluctuations determines the mass distribution function. Sharp peaks in the power spectrum imply approximately monochromatic distributions. For instance, chaotic new inflation may give rise to relatively narrow peaks (Yokoyama, 1998;Saito et al., 2008). However, inflation models with an inflection point in a plateau of the potential (García-Bellido, 2017; García-Bellido and Ruiz Morales, 2017), or hybrid inflation (Clesse and García-Bellido, 2015), predict, instead, extended mass functions, which can span over a large range of PBHs masses. In this review, we focus on monochromatic distributions for simplicity, although it is possible to translate constraints to extended mass functions, which can be very stringent despite the fact of having more parameters to fit (Carr et al., 2017;Bellomo et al., 2018). Nonetheless, even if there are bounds that exclude f PBH ∼ 1 in the monochromatic case, there are choices for the mass function which allow PBHs to constitute most or all of the DM abundance (Lehmann et al., 2018;Ashoorioon et al., 2020;Sureda et al., 2020).

ACCRETION ONTO PBHS
One of the consequences of the existence of PBHs with greater impact on different observables is the process of accretion. Infalling matter onto PBHs would release radiation, injecting energy into the surrounding medium, and strongly impacting its thermal state, leaving significant observable signatures. The physics of accretion is highly complex, but one can attempt a simplified approach considering the spherical non-relativistic limit, following the seminal work by Bondi (1952). In this framework, the BH is treated as a point mass surrounded by matter, embedded in a medium which tends to constant density far enough from the BH. The relevant scale is the so-called Bondi radius, defined by where c s,∞ is the speed of sound far enough from the PBH. At distances around ∼ r B , the accretion process starts to become FIGURE 1 | Sketch of the formation of PBHs from overdensities for three different successive moments. When fluctuations larger than a critical threshold δ c ∼ c 2 s enter the horizon, i.e., their wavelength λ 2π/k (which characterizes the size of the perturbation) is of the order of the Hubble horizon (aH) − 1 , the overdense region collapses and a PBH is produced. As can be seen in Eqs 1, 2, longer modes (large λ, low k) enter the horizon later and lead to more massive PBHs.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org May 2021 | Volume 8 | Article 681084 relevant, and the velocity of the accreted matter reaches v ∼ c s,∞ , when the density is still close to the boundary value, ρ ∞ . Taking into account the velocity of the BH relative to the medium, v rel , one can write the accretion rate as (Bondi, 1952) where the parameter λ is the dimensionless accretion rate and depends on the equation of state, but it is order ∼ 1 for the cases of interest (Ali-Haïmoud and Kamionkowski, 2017). In the early Universe, v rel can be estimated as the baryon-DM relative velocity, computed in linear theory. Its root-mean-square value is approximately constant before recombination, dropping linearly with 1 + z at later times (Ali-Haïmoud and Kamionkowski, 2017).
On the other hand, real BHs spin and form an accreting disk. Thus, the spherical symmetric case may not be applicable. Even though PBH spins are expected to be small, an accretion disk would form if the angular momentum is large enough to keep matter orbiting at Keplerian orbits at distances much larger than the innermost stable orbits, which are roughy given by the Schwarzschild radius (Agol and Kamionkowski, 2002). Applying this criterion, the formation of accretion disks around PBHs has been suggested to occur if the condition (Poulin et al., 2017b). This is satisified for M PBH ≳ M ⊙ and large enough abundances at the epoch of CMB decoupling or at later times. Some of the results outlined above are still valid if the dimensionless accretion rate λ is suppresed by roughly two orders of magnitude, after accounting for viscosity effects and matter outflows through jets (Poulin et al., 2017b).
The matter falling onto BHs is greatly accelerated, which gives rise to radiative emission of high-energy photons by processes such as bremsstrahlung. The luminosity of accreting BHs is proportional to its accretion rate, and can be written as (Xie and Yuan, 2012) where ϵ( _ M PBH ) denotes the radiative efficiency, which is in general a complicated function of the accretion rate _ M PBH , and depends upon the geometry, viscosity and other hydrodynamical considerations.
If the accretion disk is optically thin, most of the energy released through viscous dissipation is radiated away, and the luminosities obtained can be close to the Eddington luminosity, L Edd , explaining the extreme brightness of many far Active Galactic Nuclei (AGN) (Shakura and Sunyaev, 1973). However, PBHs, like nearby astrophysical BHs, may radiate in a much less efficient way, through the Advective-Dominated-Accretion-Flow (ADAF) (Ichimaru, 1977;Rees et al., 1982;Narayan and Yi, 1994;Narayan and Yi, 1995) (see, e.g., Yuan and Narayan, 2014 for a review). In this scenario, the dynamics is ruled by advective currents, forming a hot thick disk. Most of the emitted energy is deposited in the same accretion disk, heating it up. Thus, only a small portion of energy is released to the surrounding medium, the radiative process being inefficient. In the ADAF scenario, the efficiency function can be fitted by a broken power-law formula, with the slopes and amplitudes dependent on the mass range and the specific modeling of viscosity effects (Xie and Yuan, 2012).
Finally, the energy emitted in the accretion processes is deposited through different channels into the medium. The energy deposition rate for each channel reads Poulin et al., 2017b) where the subscript c denotes the channel in which energy is deposited, namely: ionization, heating, or atomic excitations (where the Lyα transitions are the most relevant ones). The energy deposition factors f c (z) quantify the fraction of energy which goes to the different channels, and has been computed numerically (Slatyer, 2016). This energy injection would affect different types of observables, as we briefly outline below.

OTHER GENERIC FEATURES
A noteworthy phenomenon of BHs is that of evaporation. Due to quantum effects in curved spacetimes, BHs may emit particles at their event horizon, as was noticed in Hawking (1974). The emitted radiation would have a nearly thermal black body spectrum, with a temperature given by (Carr et al., 2010) which is known as Hawking temperature. Due to this emission, BHs would slowly lose mass until completely evaporate. The lifetime of a PBH of initial mass M is (Lopresto, 2003) Thus, the lower the PBH mass, the earlier it evaporates. Those with masses of 10 15 g or below would have already evaporated by now, having lifetimes shorter than the age of the Universe (Page, 1976), so they cannot contribute to the current DM abundance. These evaporation products or the effects they produce in different observables can be search for in a variety of experiments, probing different mass ranges. Detailed computations of the emitted spectra can be performed by codes such as BlackHawk Arbey and Auffinger (2019).
Another important feature is that of clustering. If fluctuations are originally Gaussian distributed and around a relatively narrow peak, PBHs are not expected to be originated in clusters, being initially randomly distributed on small scales (Ali-Haïmoud, 2018; Desjacques and Riotto, 2018). However, either primordial non-gaussianities or a broad peak in the power spectrum could lead to a significant initial clustering (Ballesteros et al., 2018;Suyama and Yokoyama, 2019)  On another hand, since PBHs would be formed from the collapse of high density peaks relatively spherically symmetric, their torques and angular momentum are expected to be small (Chiba and Yokoyama, 2017;Mirbabayi et al., 2020). It is usually quantified with the dimensionless spin parameter, S S/(GM 2 PBH ), where S is the spin. Estimations of for PBHs show that it is a small quantity, equal or lower than 0.01 (De Luca et al., 2019a). In contrast, astrophysical BHs are expected to have substantially larger spins, since angular momentum must be conserved during the collapse of their stars of origin, which are often rotating. Hence, the spin can serve as a good proxy to distinguish the nature of a population of BHs. The measurement of low spin parameters could represent a hint for the detection of PBHs. The latest Bayesian analyses of LIGO/Virgo mergers suggest that low values of the spin parameter are strongly preferred by data, regardless of the assumed priors (Garcia-Bellido et al., 2020). Note, however, that the PBH mass and spin depend on the accretion mechanism and their time evolution is correlated (De Luca et al., 2020c).
Furthermore, due to the discrete nature of PBHs, a Poisson shot noise contribution to the matter power spectrum, constant in wavenumber, P sn (k) ∝ f 2 PBH n −1 PBH ∝ f PBH M PBH , would be expected (Afshordi et al., 2003;Gong and Kitajima, 2017). PBHs fluctuations give rise to isocurvature modes (Afshordi et al., 2003;Chisholm, 2006;Inman and Ali-Haïmoud, 2019), and thus, affect only at scales smaller than the horizon at the epoch of matter-radiation equality (Peacock, 1999). Therefore, this leads to an enhancement in the matter power spectrum, increasing the population of low-mass halos, which can be constrained by large scale structure and Lyα forest analyses. This effect is different from the one induced by other non-CDM candidates, such as warm DM or fuzzy DM, which suppress small scale fluctuations, washing out small structures. This contribution becomes relevant for low-mass halos not large enough to cool and collapse to form stars, which are commonly known as minihalos. It has been argued that this enhancement may produce a non-negligible 21 cm signal from the neutral hydrogen in minihalos Kitajima, 2017, Gong andKitajima, 2018). However, consistently accounting for the heating of the IGM due to PBH accretion increases the Jeans mass and suppresses the minihalo 21 cm signal, making it almost negligible (Mena et al., 2019).

OBSERVATIONAL CONSTRAINTS ON PBHS AS DM
PBHs can impact cosmology and astrophysics in a wide range of ways, leaving different observational effects which allow to constrain their properties. In this section, we review the most important bounds on the current fraction of PBHs as DM, f PBH Ω PBH /Ω DM , for a wide range of masses M PBH , for monochromatic mass functions. A collection of limits from the different probes is depicted in Figure 2. For a more comprehensive list of constraints, see Green and Kavanagh (2021); Carr and Kühnel (2020).

Evaporation
Since BHs emit energy due to Hawking radiation, those with a lifetime shorter than the age of the Universe must have disintegrated nowadays, a fact which excludes PBHs with M PBH < M * x4 × 10 14 g to form part of the current DM (Page, 1976). Moreover, PBHs with masses small enough, although still present, should emit a strong c ray and cosmic ray background which could be observed. Absence of its detection strongly constrains the range of masses M PBH ≲ 10 17 g. In particular, the maximum fraction allowed is f PBH ≲ 2 × 10 − 8 (M PBH /M * ) 3+ϵ , with ϵ ∼ 0.1 − 0.4 (Carr et al., 2016). Comparable limits have been found from INTEGRAL and COMPTEL observations of the Galactic Center (DeRocco and Graham, 2019; Laha, 2019; FIGURE 2 | Compilation of constraints on the PBH fraction (with respect to DM) as a function of the PBH mass, assuming a monochromatic mass function. The different probes considered are: impact of PBH evaporation (red) on the extragalactic c-ray background (Carr et al., 2010) and on the CMB spectrum (Clark et al., 2017); non-observation of microlensing events (blue) from the MACHO (Alcock et al., 2001), EROS (Tisserand et al., 2007), Kepler (Griest et al., 2014), Icarus (Oguri et al., 2018), OGLE (Niikura et al., 2019b) and Subaru-HSC (Croon et al., 2020) collaborations; PBH accretion signatures on the CMB (orange), assuming spherical accretion of PBHs within halos (Serpico et al., 2020); dynamical constraints, such as disruption of stellar systems by the presence of PBHs (green), on wide binaries (Monroy-Rodríguez and Allen, 2014) and on ultra-faint dwarf galaxies (Brandt, 2016); power spectrum from the Lyα forest (cyan) (Murgia et al., 2019); merger rates from gravitational waves (purple), either from individual mergers (Kavanagh et al., 2018;Abbott et al., 2019) or from searches of stochastic gravitational wave background (Chen and Huang, 2020). Gravitational waves limits are denoted by dashed lines, since they could be invalidated (Boehm et al., 2021). Dotted brown line corresponds to forecasts from the 21 cm power spectrum with SKA sensitivities (Mena et al., 2019) and from 21 cm forest prospects (Villanueva-Domingo and Ichiki, 2021). Figure created with the publicly available Python code PBHbounds (Kavanagh, 2019 Coogan et al., 2020;Laha et al., 2020). Furthermore, bounds from the isotropic X-ray and soft gamma-ray background have also been recently updated (Ballesteros et al., 2020;Iguaz et al., 2021). Additionally, data on the diffuse supernova neutrino background at Super-Kamiokande are also able to set constraints (Dasgupta et al., 2020).

Microlensing
If a compact object crosses the line of sight of a star, it may produce a so-called microlensing effect, which implies a transient and achromatic amplification of its flux. The range of masses of the objects which can produce it span from 5 × 10 − 10 to ∼ 100 M ⊙ (Paczynski, 1986;Green and Kavanagh, 2021). The non-detection of these events leads to bounds on the maximum abundance of PBHs about f PBH ≲ 0.01 − 0.1 by the MACHO (Alcock et al., 2001) and EROS (Tisserand et al., 2007) surveys in the Large and Small Magellanic Clouds, the Subaru Hyper Suprime Cam (HSC) in M31 (Andromeda) (Niikura et al., 2019a) and the Optical Gravitational Lensing Experiment (OGLE) in the Galactic bulge (Niikura et al., 2019b). Nonetheless, the existence of Earth-mass PBHs (M PBH ∼ 10 − 5 M ⊙ ) with a fraction f PBH ∼ 0.03 could explain the observed excess of six microlensing events found in the OGLE data (Mróz et al., 2017), which is consistent with other constraints in this range of PBHs masses (Niikura et al., 2019b). Although this may constitute a hint of their existence, it cannot be regarded as a detection of PBHs, since these microlensing observations could also be explained by free-floating planets. There are some caveats regarding the results of the MACHO collaboration (Hawkins, 2015), since the limits reported are model dependent and could be biased by the assumption of an over-massive halo. Moreover, the results of the MACHO and EROS projects have been found to be statistically incompatible. Therefore, these bounds are not completely reliable, and PBHs could not be definitely ruled out within this range of masses. Variations of the lensing effect in type Ia Supernovae (Zumalacárregui and Seljak, 2018) and gamma-ray bursts (Barnacka et al., 2012) have also been proposed to constrain the PBH parameter space, although these limits have been later invalidated (Garcia-Bellido et al., 2017;Katz et al., 2018).

Gravitational Waves
The observation of BH mergers by LIGO and Virgo collaborations can be employed to constrain the allowed number of PBHs. Demanding that the predicted merger rates of PBH binaries cannot exceed the ones measured by gravitational waves, tight upper bounds of f PBH ≲ 0.01 have been found for PBHs masses between 1 and 300 M ⊙ (Ali Kavanagh et al., 2018). The nonobservation of a stochastic gravitational wave background of mergers expected from a population of PBHs has also been used to constrain their abundace (Chen and Huang, 2020). However, to derive these limits, BHs have been treated as Schwarzschild BHs, while it would be more appropiate to use cosmological BH solutions embedded in a FLRW metric, such as the Thakurta metric (Thakurta, 1981). This implies a timedependent mass, and that PBH binaries created before galaxy formation would have merged at much earlier times, allowing to obtain merger rates consistent with the LIGO observations, and completely avoiding these constraints (Boehm et al., 2021).

Dynamical Constraints
Due to two-body interactions, kinetic energies of systems of different masses usually become balanced and match. If a stellar system also has an additional MACHO population, its stars would gain kinetic energy and, due to the virial theorem, the system would expand. Therefore, the presence of PBHs would dynamically heat star clusters, making them larger and with higher velocity dispersions, leading to an eventual dissolution into its host galaxy. Populations with high mass to luminosity ratios are more sensitive to this effect, as happens with ultra faint dwarf galaxies (UFDW), which would be disrupted by the presence of PBHs. Making use of these effects, tight bounds have been obtained at f PBH ∼ 10 − 3 for M PBH ∼ 10 4 M ⊙ , weakening at lower masses down to f PBH ≲ 1 for M PBH ∼ 10 M ⊙ (Brandt, 2016). In a similar way, wide binary stellar systems may be perturbed by compact objects, potentially being disrupted after multiple encounters. The separation distribution of wide binaries restricts the PBH fraction from f PBH ≲ 1 for M PBH x3 M ⊙ down to f PBH ≲ 0.1 at M PBH ≳ 70 M ⊙ (Monroy-Rodríguez and Allen, 2014).

CMB
Radiation emitted either by accretion or from Hawking evaporation may affect the CMB spectrum in two ways: producing spectral distortions and modifying temperature anisotropies. The energetic radiation can enhance the ionization rate, delay recombination and shift the peaks of the CMB anisotropy spectrum, as well as induce more diffusion damping. The polarization spectrum could also be modified, since the increase of the free electron fraction would increase the Thomson optical depth and enhance the reionization bump at large angular scales. Although early CMB analyses (Ricotti et al., 2008) found very stringent bounds on the allowed abundance of accreting PBHs, later works revisited these computations and found much milder constraints (Horowitz, 2016;. On another hand, while the former constraints rely on the assumption of spherical accretion, accreting disks have been argued to be more realistic for PBHs, resulting in tighter limits (Poulin et al., 2017b). Taking into account that PBHs could be immersed in DM halos with higher densities than the background, their accretion rates would be increased, also leading to more stringent constraints (Serpico et al., 2020). CMB limits from accretion are currently the most stringent ones for masses ≳ 10 M ⊙ . The main caveat is their dependence on some details of the accretion mechanisms, such as the effective velocity and the accretion rate, which may not be very well understood yet. On the other hand, the energy injection from PBHs evaporation would produce anisotropies and spectral distortions in the CMB spectrum, which would also limit the maximum abundance, leading to similar costraints to those obtained from the extra-galactic c-ray background commented above (Clark et al., 2017;Poulin et al., 2017a;Acharya and Khatri, 2020 spectral distortions can also be produced by other means, such as the diffusion of photons due to Silk damping at small scales. This allows the translation of constraints on spectral distortions from FIRAS to stringent upper bounds on the PBH abundance, for masses M PBH > 10 5 M ⊙ (Nakama et al., 2018).

Lyα Forest
The discrete nature of PBHs would lead to a shot-noize contribution to the matter power spectrum, enchancing small scale fluctuations. Observations of the Lyα forest, which traces matter distribution at small scales, have been employed to extract limits on the maximum allowed fraction of PBHs (Afshordi et al., 2003). The shot-noize contribution to the power spectrum depends on the joint product f PBH M PBH , for which the upper bound f PBH M PBH ≤ 60 M ⊙ has been obtained (Murgia et al., 2019). The drawback of this method is on the priors of the reionization modeling and, as any Lyα forest analysis, is model dependent.

cm Cosmology
The 21 cm line signal from the hyperfine structure of the hydrogen is highly sensitive to the thermal state of the IGM, and thus, energy injection from PBH accretion or evaporation may leave strong observable signatures. The first claimed measurement of a global absorption dip by the EDGES collaboration (Bowman et al., 2018) may lead to competitive bounds on the PBH abundance, either from accretion processes (Hektor et al., 2018) or from energy injection by evaporation (Clark et al., 2018;Halder and Banerjee, 2021;Halder and Pandey, 2021). It must be noted, however, that the EDGES signal has not been confirmed yet by other experiments, and it has been argued that it could be explained by alternative mechanisms (Hills et al., 2018;Bradley et al., 2019;Sims and Pober, 2020). On the other hand, although the cosmological 21 cm power spectrum has not been detected yet, forecasts with future experiments such as HERA and SKA have shown that 21 cm power spectrum data from these two radiotelescopes could potentially improve the bounds up to f PBH < 10 − 2 − 10 − 6 for masses above M ⊙ (Mena et al., 2019). The 21 cm forest observed as absorption troughs in the spectra of radio sources at z ∼ 10 − 15 could also provide similar limits on the PBH abundance, due to the Poisson shot noise and to the accretion heating effect (Villanueva-Domingo and Ichiki, 2021).

CONCLUSION
The extremely rich physics involved in the formation, evolution and distribution of PBHs implies a large number of observable effects which allow probing them. A myriad of constraints are present for a large range of PBHs masses. In recent years some of these limits, as those from microlensing, femtolensing, CMB accretion or BH mergers, have been revisited. Some of these bounds have been significantly weakened or have even disappeared, opening windows in the parameter space where PBHs could still form a substantial part of the DM, if not all. On the other hand, future experiments with better sensitivities may be able to reach yet unexplored regions in the parameter space and tighten up current limits. New probes, such as the 21 cm line pursued in radio interferometers like SKA, will present a promising and powerful way to proof or refute the existence of PBHs formed in the early universe and their potential contribution to the DM in the Universe.

AUTHOR CONTRIBUTIONS
PV-D is the main author of the text and the figures. OM and SP-R have substantially collaborated in the editing process.