$$ \langle m_{\text{PBH}}^2 \rangle = \frac{\int m_{\text{PBH}}^2 \mathrm{d}n_{\text{PBH}}}{n_{\text{PBH}}} = \frac{\int m_{\text{PBH}} f(m_{\text{PBH}}) \mathrm{d} \ln m_{\text{PBH}}}{\int \frac{f(m_{\text{PBH}})}{m_{\text{PBH}}} \mathrm{d} \ln m_{\text{PBH}}} \quad (105) $$

where $n_{\text{PBH}}$ denotes the total PBH number density. The function $C$ encodes the transition between small and large $\bar{N}$ limits. A good approximation is given by (Hütsi et al. 2021c)

$$ C \approx \frac{f_{\text{PBH}}^2 \langle m_{\text{PBH}}^2 \rangle}{\sigma_m^2 \langle m_{\text{PBH}} \rangle^2} \times \left\{ \left[ \frac{\Gamma(29/37)}{\sqrt{\pi}} U\left(\frac{21}{74}, \frac{1}{2}, \frac{5f_{\text{PBH}}^2}{6\sigma_m^2}\right) \right]^{-74/21} - 1 \right\}^{-1}, \quad (106) $$

where $\Gamma$ is the Euler function and $U$ is the confluent hypergeometric function.

The value of $S_1$ as a function of $f_{\text{PBH}}$ has been represented in Fig. 20 for different cases. Assuming $\bar{N} = 2$, as expected for a close-to-monochromatic distribution, the value of $S_1$ ranges from $10^{-2}$ when $f_{\text{PBH}} \simeq 10^{-3}$ to 0.1 when $f_{\text{PBH}} \simeq 1$. In particular, one can also obtain a maximum value $S_1^{\text{max}}$ that is independent of the mass distribution when $\bar{N} \to 0$. This maximal value is around unity for $f_{\text{PBH}} \gtrsim 0.1$ and goes down to $10^{-1}$ for lower values. In this limit, one recovers the results from Chen and Huang (2018), Ali-Haïmoud et al. (2017) that only take into account the effect of matter inhomogeneities.

Finally, we would like to comment on the calculation of $S_1$ using different definitions of the PBH mass function, which in particular modify the way to compute $\langle m \rangle$ and $\langle m^2 \rangle$ depending on whether the mass distribution is defined with respect to the PBH number density $n$ or the PBH density $\rho$, or if the normalisation changes. The possible definitions, the way they are related, and the corresponding calculation rules are summarized in Table 2 from Escrivà et al. (2023). It is important to do these calculations carefully because they can significantly impact the resulting rates, and because some inconsistencies were found in the literature. As pointed out in Escrivà

<figure><img src="image_22.png" /><figcaption>Fig. 20 Suppression factor $f_{\text{sup}}$ and its possible contributions $S_1$ and $S_2$ in different cases, in terms of the PBH abundance. When PBH masses are specified, it assumes a PBH mass function from a nearly-scale invariant primordial power spectrum of curvature fluctuations and including QCD-induced features. Image</figcaption></figure>
et al. (2023), there was a conversion factor $m/\langle m \rangle$ not considered in Franciolini et al. (2022e) which led to an inconsistency with Hütsi et al. (2021c) which uses exactly the same definition and normalization of the PBH mass distribution, but a different merger rate formula. There was also a typo in the Eq. (A4) of Cesse and García-Bellido (2022) that has an incorrect mass dependence even if the correct formula was used for the calculations and paper figures.

Following these rules for an extended mass function, one gets that the abundance of light PBHs can strongly impact $\langle m \rangle$ and the resulting values of $\bar{N}$ that typically increases. In turn, $S_1$ and the resulting value of $f_{sup}$ can be strongly suppressed due to the exponential dependence in $\bar{N}$, and their exact value is dependent on the chosen lower mass cut for the PBH mass distribution. Physically, this suppression comes from the numerous light black holes, reducing the value of the average black hole mass, at proximity of the binary that are considered as being able to disrupt it when one uses the presented conservative approach. Such an issue was pointed out in Gow et al. (2020), Cesse and García-Bellido (2022), Escrivà et al. (2023) and suggests that the rates of early binaries should probably be revised in the case of a very broad mass distribution. In absence of a clear solution, the approach of Cesse and García-Bellido (2022), Escrivà et al. (2023) was to consider the limit obtained for a monochromatic mass function, motivated by the QCD-induced peak, and assume that $\bar{N} \approx 2$. In Franciolini et al. (2022e), which used another averaging procedure, there was a strong dependence on the high-mass PBH distribution and a high-mass cut-off was introduced that may naturally arise from the transition in the primordial power spectrum. The different approaches however lead to different merger rates and different possible interpretations of the LIGO/Virgo observations, as well as different limits on the possible PBH abundance.

<table><caption>**Table 2** Different definitions of the normalized PBH mass distribution proposed in various references with their conversion, their normalisation rule and the corresponding $\langle m \rangle$ and $\langle m^2 \rangle$</caption><thead><tr><th></th><th>Escrivà et al. (2023) Cesse and<br>García-Bellido (2022)* Bagui and<br>Cesse (2022)</th><th>Raidal et al. (2019)<br>Franciolini et al. (2022e)*<br>Hall et al. (2020)</th><th>Hütsi et al.<br>(2021c)</th><th>Kocsis et al.<br>(2018)</th></tr></thead><tbody><tr><td>$f \equiv \frac{1}{\rho_{PBH}} \frac{d\rho_{PBH}}{d\ln m}$</td><td></td><td>$\psi_1 \equiv \frac{1}{\rho_{PBH}} \frac{d\rho_{PBH}}{dm}$</td><td>$\psi_2 \equiv \frac{1}{n_{PBH}} \frac{dn_{PBH}}{d\ln m}$</td><td>$\psi_3 \equiv \frac{1}{n_{PBH}} \frac{dn_{PBH}}{dm}$</td></tr><tr><td>$f = m\psi_1$</td><td></td><td>$\psi_1 = f/m$</td><td>$\psi_2 = f\langle m \rangle/m$</td><td>$\psi_3 = \langle m \rangle f/m^2$</td></tr><tr><td>$= m\psi_2/\langle m \rangle$</td><td></td><td>$= \psi_2/\langle m \rangle$</td><td>$= \langle m \rangle \psi_1$</td><td>$= \langle m \rangle \psi_1/m$</td></tr><tr><td>$= m^2\psi_3/\langle m \rangle$</td><td></td><td>$= m\psi_3/\langle m \rangle$</td><td>$= m\psi_3$</td><td>$= \psi_2/m$</td></tr><tr><td>$\int f d\ln m = 1$</td><td></td><td>$\int \psi_1 dm = 1$</td><td>$\int \psi_2 d\ln m = 1$</td><td>$\int \psi_3 dm = 1$</td></tr><tr><td>$\langle m \rangle = \left(\int \frac{f}{m} d\ln m\right)^{-1}$</td><td></td><td>$\left(\int \frac{\psi_1}{m} dm\right)^{-1}$</td><td>$\int m\psi_2 d\ln m$</td><td>$\int m\psi_3 dm$</td></tr><tr><td>$\langle m^2 \rangle = \langle m \rangle \int m f d\ln m$</td><td></td><td>$\langle m \rangle \int m\psi_1 dm$</td><td>$\int m^2 \psi_2 d\ln m$</td><td>$\int m^2 \psi_3 dm$</td></tr></tbody></table>

The asterisk denotes the references in which an inconsistency has been found (see details in the text). In this review, we considered both $f$ and $\psi_1 \equiv \psi$ with our notations. Table adapted from Escrivà et al. (2023)
### 4.1.3 Suppression from early Poisson-induced clustering

The second factor $S_2(f_{PBH})$ accounts for the effect of binary disruption in PBH clusters that rapidly form after PBH formation if their abundance is large enough (Inman and Ali-Haïmoud 2019), due to their initial Poisson fluctuations. The interaction of binaries with other PBHs in early-Universe clusters typically induces a modification of their semi-major axis and eccentricity, which leads to longer merger timescales and a suppression of the merger rate (Trashoras et al. 2021; Jedamzik 2020, 2021). In order to have a conservative estimate of the merger rate, Vaskonen and Veermäe (2020) (and De Luca et al. 2020a includes an analytical modelling of clustering evolution) computed the fraction of binaries which never enter in dense enough environments. This fraction $S_2$ today can be approximated by (Hütsi et al. 2021c)

$$S_2 \approx \min(1, 9.6 \times 10^{-3} f_{\text{PBH}}^{-0.65} E^{0.03 \ln^2 f_{\text{PBH}}}). \quad (107)$$

The value of $S_2$ as a function of $f_{PBH}$ has been represented and compared to $S_1$ and the resulting total suppression factor $f_{sup}$ in Fig. 20. When it is considered at a time $t$, a good approximation is to replace $f_{PBH}$ by $f_{PBH}[t(z)/t_0]^{0.44}$, valid for $z \le 100$ (Hütsi et al. 2021c).

### 4.1.4 Merging rates of perturbed binaries

Numerical studies find that binaries are likely to be perturbed, i.e. early binaries whose semi major axis and eccentricity are modified but remain bound after the interactions with surrounding perturbers after their formation. This can happen in two ways: (1) The initial configuration contains a third PBH close to the initial PBH pair that will very likely collide with the forming binary; (2) Binaries are absorbed by dense PBH clusters that form early in the history of the Universe. A conservative calculation of the possible merging rates of PBH binaries that have been disrupted in early-forming halos has been proposed in (Vaskonen and Veermäe 2020). Moreover, let us note that estimates indicate that most PBH binaries are not perturbed in typical DM halos larger than 10 PBHs, and that nearly all initial binaries might be disrupted within the age of the Universe, since the fraction of disrupted binaries is relatively high at $z = 1100$ when $f_{PBH} \gtrsim 0.1$.

For a monochromatic mass function, the merger rate of perturbed binaries at time $t$ is given by (Vaskonen and Veermäe 2020)

$$R_p(t) \propto f_{\text{PBH}}^{\frac{144}{259} + \frac{47}{37} \gamma} t^{\frac{5}{7} - 32} m_{\text{PBH}}^{\frac{5}{7} - 32}, \quad (108)$$

where $\gamma \in [1, 2]$, based on the angular momentum distribution for perturbed binaries in the early Universe. The main result of this analysis is that the rate of perturbed binaries can exceed the rate of the unperturbed ones discussed in the previous section, when $f_{PBH} \gtrsim 0.1$. However, when $f_{PBH} \lesssim 0.02$, the rate of perturbed binaries should be subdominant, at least for peaked mass functions. But it is still unclear how this formula and results are generalized to any broad PBH mass distribution.
4.1.5 Link with observations of compact binary coalescences

In Fig. 20, we represented the values of $S_1$, $S_2$ and the resulting suppression factor $f_{sup}$ as a function of $f_{PBH}$ for the choice of the broad mass function of Clesse and García-Bellido (2022), at the mean masses of the exceptional GW events GW190425, GW190814 and GW190521.<fnref n="7" /> We also calculated $S_1^{\text{max}}$ in the limit $\bar{N} \ll \min(C, 1)$ that is independent of the two binary component masses.

For an arbitrary broad mass distribution, one critical difference with respect to a monochromatic or a relatively sharp log-normal mass function comes from the large number density of tiny black holes, implying that $\langle m_{\text{PBH}} \rangle \ll M_{\odot}$ and $\langle m_{\text{PBH}}^2 \rangle / \langle m_{\text{PBH}} \rangle^2 \ll 1$. This implies that generically $\bar{N} \gg 1$. If the analytical prescription is strictly applied, this leads to a huge exponential suppression of the EB merging. However, this rate suppression is very likely overestimated because PBHs that are much lighter than the binary components are probably not able to bring enough energy to the system in order to disrupt it.

Instead, it is probably more realistic to consider only the disruption by nearby PBHs whose mass is similar to the mean of the binary component masses. By doing so, one gets $\langle m_{\text{PBH}} \rangle \sim (m_1 + m_2)/2$ and $\langle m_{\text{PBH}}^2 \rangle / \langle m_{\text{PBH}} \rangle^2 \sim 1$. In such a case, one gets $\bar{N} \approx 2$ (just as in the monochromatic case) and the suppression factor obtained when $f_{\text{PBH}} \simeq 1$ becomes independent of the mass, slightly below $S_2^{\text{max}}$, depending on the exact value of $\bar{N}$. Different values of $\bar{N}$ (but still of the same order) typically lead $f_{sup}$ to be of order $S_2^{\text{max}} \times \mathcal{O}(0.1 - 1)$. In particular, for $\bar{N} = 2$ and $f_{\text{PBH}} = 1$, one gets $f_{sup} \simeq 0.002$. We have chosen this value to compute the EB merging rate for two typical mass functions impacted by the thermal history. They are represented on Fig. 21. One can see that these merging rates can in principle accommodate some of the LIGO/Virgo GW events and are significantly enhanced when the binary includes a PBH from the QCD peak.

4.1.6 Limitations

Even if the merging rates given by Eq. (101) are consistent with N-body simulations for monochromatic and log-normal mass functions in the stellar-mass range, there are still a series of uncertainties that could limit the applicability of this formalism. Hereafter we list some of these possible limitations.

1. So far, no N-body simulations have been performed on PBH models that include thermal features, whereas it is proven that they have an important impact for stellar-mass PBHs. In general, broad-mass functions could see also deviations from those estimations, in particular for binaries with very low mass ratios.

2. Eq. (101) including the effect described in Eq. (108) should be improved by more appropriate calculation of the distribution of the binary parameters as

<footnote n="7">⁷ We note that the parameter estimation of GW events characterised by a small signal-to-noise ratio at the LVKC detectors may still be sensitive to the priors assumed in the analysis (see e.g. Vitale et al. 2017). Therefore, the interpretation of some events may change if PBH motivated priors are assumed (Bhagwat et al. 2021), in particular when focusing on the binary mass ratio and effective spin.</footnote>
Model 1: Carr 2019, $n_s = 0.97$

<figure><img src="image_23.png" /></figure>

Model 2: De Luca 2020, $n_s = 1$

<figure><img src="image_24.png" /><figcaption>**Fig. 21** Expected differential merging rates $dR/d\ln m_1 d\ln m_2$ of PBH of masses $m_1$ and $m_2$, for the two mass models represented on Fig. 13 (top panels: Model 1, bottom panels: Model 2), for the two considered binary formation channel: primordial binaries (see Eq. (101)) on the left panels, and tidal capture in halos (see Eq. (109)) on the right panels. Image reproduced with permission from Auclair et al. (2023), copyright by the author(s). Image reproduced with permission</figcaption></figure>

interactions have occurred, encoded in the parameter $\gamma$ (Vaskonen and Veermäe 2020). Analytical prescriptions for non-equal mass binaries should also be further developed and checked.

3. The fraction of PBH binaries that are disrupted as a function of the mass of an infalling PBH is uncertain. This question is particularly relevant for binaries with strong differences in mass components, and for broad mass functions with lots of small PBHs that would dominate the calculation of $\bar{N}$, without being massive enough to disrupt a massive binary. The same question arises in the case of sub-dominant intermediate-mass binaries that would hardly be disrupted by a dominant population of stellar-mass or asteroid-mass PBHs.

4. The effect of heavy PBHs in the tail of a broad mass distribution would be to seed PBH clusters, an effect that would be superimposed to the Poisson clustering, but that has been neglected so far in merging rate computations.

5. Slightly different results and another possible dependence in $f_{\text{PBH}}$ have been obtained in Kocsis et al. (2018) using analytical methods.

6. The suppression due to PBH clusters depends on the exact clustering history, on the cluster relaxation time, on the evaporation as well as on internal dynamics (see e.g. De Luca et al. 2020a).
7. PBH clustering at formation may either boost or suppress the merger rates of early binaries, but this effect is very model dependent and still debated (see e.g. Raidal et al. 2017; Young and Byrnes 2020; De Luca et al. 2021d for different models and conclusions).

8. Dark matter spikes around PBHs should affect the merging time of PBH binaries and enhance their merger rates, which typically leads to more stringent limits on their abundance inferred from GW observations (Kavanagh et al. 2018; Pilipenko et al. 2022; Jangra et al. 2023). Nevertheless the whole process might be more complicated with various possible regimes and dynamics, boosting or suppressing the merger rates and leading to model-dependent results, as noticed in Jangra et al. (2023).

9. It has recently been claimed in Boehm et al. (2021) that subtle general relativistic effects may highly suppress this PBH binary formation channel, but this result has been disputed in De Luca et al. (2020a), Hütsi et al. (2021a, 2021b). This problem is related to the question of which metric is physically relevant to describe the PBH environment (Takhurta metric, generalized Mc-Vittie, other).

In general, one should not forget that early binaries can be impacted by their environment during the whole cosmic history, and this environment has a complex evolution, influenced by the clustering after matter-radiation equality, accretion, dynamical heating, etc. Therefore one should be cautious and strong claims relying on these merging rates are probably still premature. Nevertheless, Eq. (101) likely gives good estimations for some models, at least in some regimes. To that end, they are relevant for estimating not only the merging rates but also the resulting SGWB background, based on our current (but rapidly evolving) knowledge.

## 4.2 Late binaries

### 4.2.1 General formula

The second binary formation channel is through dynamical capture in dense environments. We will give a particular emphasis to the case in which PBHs comprise the totality of the dark matter in the universe and investigate the dynamical formation in dense PBH halos. As any other DM candidate, PBHs are expected to form halos during the cosmic history. The clustering properties typically determine the overall merging rate. For a generic PBH mass function $f(m_{\text{PBH}})$, an effective formula has been proposed (Carr et al. 2021b),

$$R_{\text{LB}} \approx R_{\text{clust}} f_{\text{PBH}}^2 f(m_1) f(m_2) \frac{(m_1 + m_2)^{10/7}}{(m_1 m_2)^{5/7}} \text{yr}^{-1} \text{Gpc}^{-3}, \quad (109)$$

where $R_{\text{clust}}$ is a scaling factor that depends on the PBH clustering properties and velocity distribution. This formula comes from the two-body capture in a cluster and assumes that the time it takes for the binary to merge is much shorter than cosmological times.
Hereafter, we present the estimation for $R_{\text{clust}}$ in the case of Poisson clustering. We then discuss the implications in the context of the present LIGO/Virgo observations. Finally, we present the possible limitations of these approaches and how different models could lead to different merging rates.

Poisson clustering introduces a natural clustering scale, leading to $R_{\text{clust}}$ in the range between $10^2$ and $10^3$ for solar-mass PBHs and $f_{\text{PBH}} \simeq 1$ (De Luca et al. 2020a). Other clustering scenarios (e.g. clusters at formation) could induce even larger values of $R_{\text{clust}}$, fixed by the model as well as the evolution of these clusters (dynamical heating, mergers, etc.) through the cosmic history.

In general, one can consider a simplified picture where the PBH halo density is described by the local density contrast $\delta_{\text{loc}}$. The number density of PBHs can be given as $n(m_{\text{PBH}}) \equiv \delta_{\text{loc}}\bar{\rho}_{\text{DM}}/m_{\text{PBH}}$, with the mean dark matter cosmological energy density $\bar{\rho}_{\text{DM}} = \Omega_{\text{DM}}\rho_{\text{c}}$, where $\Omega_{\text{DM}} \simeq 0.25$ is the density parameter for dark matter. In this case, $R_{\text{clust}}$ was estimated as (Clesse and García-Bellido 2017a; García-Bellido et al. 2022)

$$R_{\text{clust}} \approx 3.6h^4 \left(\frac{\Omega_{\text{DM}}}{0.25}\right)^2 \left(\frac{\delta_{\text{loc}}}{10^8}\right) \left(\frac{v_0}{10 \text{ km/s}}\right)^{-11/7}, \quad (110)$$

where $v_0$ is the velocity of PBHs and $h$ is the Hubble parameter today in units of $H_0 = 100 \text{ h km/s/Mpc}$. If one assumes $v \sim O(1-10)\text{km/s}$ and local densities $\delta_{\text{loc}} \sim 10^8$ typical to the ones of ultra-faint-dwarf galaxies, as one may also expect for Poisson clusters of size around $\sim 10-100\text{pc}$, one typically get $R_{\text{clust}} \sim 10^2-10^3$. This provides an effective formula but in realistic scenarios, several effects could impact the exact value of $R_{\text{clust}}$ such as the dynamical heating, merger and disruption history of PBH clusters, how extended is the PBH mass function leading to different efficiencies of the Poisson clustering effect, the radial distribution of PBHs inside a cluster and their possible mass segregation for extended mass functions.

### 4.2.2 Merging rates of binaries induced by three-body interactions

PBH binary formation is also expected to take place in the Poisson induced PBH small scale structure at high redshift. After the matter-radiation equality, the seeding effect of Poisson fluctuations in the PBH distribution give rise to the formation of small scale clusters (see Inman and Ali-Haïmoud 2019 for a N-body simulation and De Luca et al. 2020a for an analytical treatment). In this clusters, three-body interactions may efficiently produce binaries which can then merge in a time-scale comparable to the age of the Universe and be visible with present day GW experiments (Franciolini et al. 2022c) (at odds with mergers induced by GW capture which merge promptly). In this sense, this formation mechanism is similar to the one taking place in the early Universe out of binary directly decoupling from the Hubble flow.

The rate density for three bodies to interact and form PBH binaries can be expressed as (Rodriguez et al. 2022)
$$
\begin{align}
\gamma_{3b}(\eta \ge \eta_{\min}) &= \frac{3^{9/2} \pi^{13/2}}{2^{25/2}} \eta_{\min}^{-11/2} (1+2\eta_{\min})(1+3\eta_{\min}) \times \frac{n_{\text{PBH}}^3 (Gm_{\text{PBH}})^5}{\sigma_v^9} \tag{111} \\
&\approx 3.8 \times 10^{-2} \text{Gyr}^{-1} \text{pc}^{-3} \left(\frac{n_{\text{PBH}}}{\text{pc}^{-3}}\right)^3 \left(\frac{m_{\text{PBH}}}{30M_{\odot}}\right)^5 \left(\frac{\sigma_v}{\text{km/s}}\right)^{-9}, \nonumber
\end{align}
$$

where $n_{\text{PBH}}$ is the number density of PBHs. We also defined the hardness ratio $\eta = Gm_{\text{PBH}}/a\sigma_v^2$ as the binding energy of a binary with size $a$ normalized to the average cluster kinetic energy, where $\sigma_v$ is the PBH velocity dispersion within the clusters. Also, in order to require the efficiency of binary formation to be close to unity and PBH binaries resulting from this channel to be in the region of parameter space where binaries are hard, one should set $\eta_{\min} = 5$ as used in the literature (see e. g. Morscher et al. 2015; Rodriguez et al. 2022).

One can compute the merger rate density of binaries produced by 3-body interactions by integrating the binary formation rate over the PBH cluster distribution throughout the evolution of Universe and by multiplying by the fraction of binaries merging within the remaining time window (see Franciolini and Pani 2022 for more details). Depending on the assumed cluster core size and eccentricity distribution $e \equiv \sqrt{1-j^2}$ for the 3-body induced binaries (Ambartsumian 1937; Heggie 1975; Stone and Leigh 2019; Raidal et al. 2019; Vaskonen and Veermäe 2020) (i.e. $f(j) \simeq j^\gamma$ with $\gamma = 1$ for thermal or $\gamma = 0$ for super-thermal), this channel was found to give at least a contribution comparable to the one from dynamical capture discussed in the following section (Franciolini and Pani 2022).

One can also estimate how this contribution scales with the PBH abundance. Following Franciolini and Pani (2022), this was estimated to be $R^{3b} \propto f_{\text{PBH}}^{2+(1+\gamma)/2}$, showing this channel may only be relevant, compared to the EB contribution, in case of large values of the PBH abundance $f_{\text{PBH}}$.

### 4.2.3 Link with observations of compact binary coalescences

If PBHs significantly contribute to the dark matter, late binaries merging rates are enhanced around $30-100M_\odot$, where most compact binary mergers have been observed. Assuming Poisson clustering and $R_{\text{clust}} = 450$, the corresponding merging rates for the PBH mass function of Fig. 13 that includes thermal features have been represented on Fig. 21. Typically one can obtain the LIGO/Virgo rates in this range with values of $f_{\text{PBH}}^2 \times f(m_{\text{PBH}})^2 \sim 10^{-3}$. At the solar-mass scale, the rates for early and late binaries compete at a similar level and it is possible to get the LIGO/Virgo rates inferred for the latest events observed. These are canonically interpreted as neutron star binary mergers, even though the absence of the observation of tidal effects may allow for a primordial interpretation, with the exception of GW170817 for which the electromagnetic counterpart has also been observed.
### 4.2.4 Limitations

For realistic, extended mass functions, a series of effects can either boost or suppress the merging rate of late binaries and make it a rather complex and model-dependent process, subject to large uncertainties (see e.g. Chisholm 2006, 2011; Ali-Haïmoud 2018; De Luca et al. 2020a; Inman and Ali-Haïmoud 2019; Clesse and García-Bellido 2022, 2017b; Moradinezhad Dizgah et al. 2019; Belotsky et al. 2014; Young and Byrnes 2020; De Luca et al. 2021f; Padilla et al. 2021; Suyama and Yokoyama 2019; Ballesteros et al. 2018 for studies on PBH clustering). Some of these limitations or effects are listed below:

1. The Poisson clustering is modified by the mass function and could be dominated by other mass scales than the dominant PBH masses. Therefore, one in principle needs to consider the whole mass spectrum in order to derive the importance of Poisson clustering.

2. The mass function also impacts the process of dynamical heating and dilution of sub-halos inside larger halos, which is essential to derive the natural clustering scale and the $R_{\text{clust}}$ parameter.

3. Poisson clustering should induce a broad range of cluster masses, evolving and interacting with each others in possibly complex ways through dynamical heating, dilution, collisions, tidal disruptions, etc, resulting in possibly complex time dependence and various possible PBH local environments.

4. The existence of intermediate-mass and supermassive PBHs leads to a seed effect that should also impact the clustering properties.

5. Clustering at formation, obtained e.g. for non-Gaussian models (Atal et al. 2020), also have an uncertain impact.

6. For extended mass function, the role of mass segregation in a cluster should be considered. In average, more massive PBHs falling at the center should have a higher local density than lighter PBHs, leading to an additional mass dependence in the merger rates that was not considered so far in most analysis.

### 4.3 Hyperbolic encounters

In the scenario of *clustered* PBH of Clesse and García-Bellido (2015), it is expected that a large fraction of BH encounters will not end up producing bound systems, which would then inspiral, but rather produce a single scattering event, via an hyperbolic encounter. This could happen, e.g. if the relative velocity or relative distance of the two PBHs is high enough that capture is not possible. The emission of GWs in close encounters of compact bodies has been studied in Kocsis et al. (2006), O’Leary et al. (2009), Capozziello et al. (2008), De Vittori et al. (2012), García-Bellido and Nesseris (2017, 2018), Gröbner et al. (2020), Mukherjee et al. (2021). These events generate *bursts* of gravitational waves, which can be sufficiently bright to be detected at distances up to several Gpc. For clustered PBHs, the waveform and characteristic parameters of the GW emission in hyperbolic encounters are different to those of the inspiralling binaries, and both provide complementary information
that can be used to determine the evolved mass distribution of PBHs as a function of redshift, as well as their spatial distribution.

Hyperbolic encounters are single scattering events where the majority of the energy is released near the point of closest approach, and have a characteristic peak frequency which is a function of only three variables (see Fig 22): the impact parameter *b*, the eccentricity *e* and the total mass of the system $M_{tot} = m_1 + m_2$. Furthermore, the duration of such events is of the order of a few milliseconds to several hours, depending on those parameters. The case of inspiralling and merging PBH has been studied extensively, see e.g. Clesse and García-Bellido (2017a, 2017b), and estimated the production of a few tens of events/year/Gpc³ in the range of $m_{PBH} \sim O(10-100) M_{\odot}$. In García-Bellido and Nesseris (2017, 2018) it was shown that, within the parameter space of the clustered PBH scenario (Clesse and García-Bellido 2015, 2017b), we can expect a similar but somewhat lower rate of GW burst events in the millisecond range (Fig. 22).

The waveforms of the GW emission in hyperbolic encounters are very different from those of the inspiralling binaries, since the majority of the energy is released near the point of closest approach, generating a burst of GWs with a characteristic “tear-drop” shape of the emission in the time-frequency domain (García-Bellido and Nesseris 2017, 2018). The burst has a characteristic peak frequency

$$f_{\text{peak}} = 0.32 \, \text{mHz} \times \frac{\beta(e+1) \, \text{AU}}{(e-1) \, b}, \qquad (112)$$

which corresponds to the maximum GW emission and depends only on the impact parameter *b*, the total mass of the system $M_{tot}$ and the eccentricity of the hyperbolic orbit $e = \sqrt{1 + b^2 v_0^4 / G^2 M_{tot}^2}$, where $v_0$ is the asymptotic relative velocity of the encounter and $\beta \equiv v_0/c$. The maximum strain amplitude and power of the GW burst is given by

$$h_{\text{c}}^{\text{max}} = 3.24 \times 10^{-23} \frac{R_{\text{s}}(\text{km})}{d_{\text{L}}(\text{Gpc})} \frac{q\beta^2 g_{\text{max}}}{(1+q)^2}, \qquad (113)$$

$$P_{\text{max}}(e) = 5.96 \times 10^{26} L_{\odot} \frac{q^2 \beta^{10} (e+1)}{(1+q)^4 (e-1)^5}, \qquad (114)$$

where $L_{\odot}$ is the solar luminosity, $R_{\text{s}}$ is the Schwarzschild radius associated to the

<figure><img src="image_25.png" /><figcaption>**Fig. 22** The scattering of one BH of mass $m_2$ on another of mass $m_1$ induces the emission of gravitational waves which is maximal at the point of closest approach, $r_p$. Image reproduced with permission</figcaption></figure>
total mass, $q \equiv m_1/m_2$ is the binary mass ratio, $m_1 = qm_2 \ge m_2$ and $g_{\max} = 2\sqrt{18(e+1)+5e^2/(e-1)}$.

GW signals in the LISA range could be generated by close hyperbolic encounters (CHE) of an IMBH and a SMBH as expected in the centers of galaxies, as well as from encounters of two SMBHs that could occur during galaxy collisions at low redshift. In the first case, an IMBH of mass $m_2 = 10^3 M_\odot$ and a SMBH of mass $m_1 = 10^6 M_\odot$, with an impact parameter $b = 1$ AU and velocity $v_0 = 0.05 c$ gives an eccentricity parameter of $e = 1.031$ with a duration of the event of 440 s. The maximum stress amplitude would be $h_c^{\max} = 1.02 \times 10^{-19}$ at a distance $d_L = 1$, with the peak at frequency $f_{\text{peak}} = 1.05$ mHz, well within the sensitivity band of LISA. In a hyperbolic encounter between two SMBHs of equal masses $m_1 = m_2 = 10^6 M_\odot$ with impact parameter $b = 10$ AU and relative velocity $v_0 = 0.015c$, the eccentricity is low, $e = 1.01$, and the stress amplitude is huge, $h_c^{\max} = 2.22 \times 10^{-17}$ at $f_{\text{peak}} = 1.51 \times 10^{-4}$Hz, which falls again in the LISA sensitivity range.

The event rate can be estimated by the cross-section of a CHE event, which is given by $\sigma = \pi b^2 = \pi (GM_{\text{tot}}/v_0^2)^2 (e^2 - 1)$. In the picture of the local density contrast $\delta_{\text{loc}}$, the total event rate for a generic PBH mass function is given by (García-Bellido et al. 2022)

$$ \frac{d\tau^{\text{CHE}}}{dm_1 dm_2} \approx 6.1 \times 10^{-8} \text{ yr}^{-1} \text{Gpc}^{-3} h^4 \left(\frac{\Omega_{\text{DM}}}{0.25}\right)^2 \left(\frac{\delta_{\text{loc}}}{10^8}\right) \frac{f(m_1)f(m_2)}{m_1 m_2} \frac{M^2}{m_1 m_2} \frac{e^2 - 1}{(v_0/c)^3}. \quad (115) $$

These CHE events are very common in dense clusters formed soon after recombination, while the PBH scatter off each other and puff-up the cluster. Some of these events lose so much energy in their GW emission that they end up in bounded systems. Many also are responsible for the disruption of previously formed binaries. All in all, there is a significant production of gravitational waves from CHE that can contribute to a SGWB early on. In the next Sections we will explore such a stochastic background from CHE events across the history of the universe, and will see that it can be significant for a wide PBH mass function, even in the LISA range of frequencies.

## 4.4 How to distinguish PBH from astrophysical binaries

In this Section we elaborate on which observable may be used to assess or rule out the primordial origin of an individual merger event. We follow the comprehensive description given in Franciolini et al. (2022b). A systematic strategy to use these discriminators is summarized in the flowchart of Fig. 23, based on the predictions of the standard PBH scenario where binaries are assembled in the early Universe. Notice that most of predictions would still apply to a much larger class of models, see Franciolini et al. (2022b) for more details.

One can identify two “smoking-gun” signatures of the primordial scenario, which are high redshift mergers and sub-solar masses. A general prediction of the PBH model is a merger rate density which grows monotonically with redshift (Ali-
<figure><img src="image_26.png" /><figcaption>**Fig. 23** Schematic flowchart representing how to systematically rule out or potentially assess the primordial origin of a binary merger. These criteria are based on measurements of the redshift *z*, eccentricity *e*, tidal deformability Λ, component masses *m*, and dimensionless spin χ. Each arrow indicates if the condition in the box is met (green) or violated (red), while the marks indicate: ✓ likely to be a PBH binary; ✗ cannot to be a PBH binary; ? may be a PBH binary. Image reproduced with permission from Franciolini et al. (2022b), copyright by APS Image reproduced with permission</figcaption></figure>

Haïmoud et al. 2017; Raidal et al. 2019; De Luca et al. 2020d). Focusing on early binaries, one expects $R^{\text{EB}}(z) \propto (t(z)/t(z=0))^{-34/37}$, extending up to redshifts $z \gtrsim \mathcal{O}(10^3)$. On the contrary, astrophysical-origin mergers are generically not expected to occur at $z \gtrsim 30$, even though the redshift corresponding to the epoch of first star formation is still poorly known. To give a conservative estimate, theoretical calculations and cosmological simulations suggest star formation does not precede $z \sim 40$ (Schneider et al. 2000, 2002, 2003; Bromm 2006; Tornatore et al. 2007; Trenti and Stiavelli 2009; de Souza et al. 2011; Koushiappas and Loeb 2017b; Mocz et al. 2020; Liu and Bromm 2020a). Additionally, the time delay between Pop III star formation and BBH mergers was estimated to be around $\mathcal{O}(10)$ Myr (Kinugawa et al. 2014, 2016; Hartwig et al. 2016; Belczynski et al. 2017; Inayoshi et al. 2017; Liu and Bromm 2020b, a; Kinugawa et al. 2020; Tanikawa et al. 2021; Singh et al. 2022a), implying merger redshifts $z \gtrsim 30$ to represent smoking guns for primordial binaries (Koushiappas and Loeb 2017b; De Luca et al. 2021c; Ng et al. 2022a). Observations of such distant events may be characterised by large measurement uncertainties on the inferred luminosity distance due to the low SNR at 3 G detectors (Ng et al. 2022a; Franciolini et al. 2022b; Martinelli et al. 2022; Ng et al. 2023), including LISA. However, it was recently shown that focusing on constraining the merger rate evolution at redshift larger than $\mathcal{O}(10)$ may allow to constrain PBH populations up to abundances as low as $f_{\text{PBH}} \approx 10^{-5}$ (Ng et al. 2022b) in the solar mass range, even accounting for the contamination of Pop III binaries.
Sub-solar BHs are not expected in standard stellar evolution.<fnref n="8" /> Other compact objects like white dwarfs, brown dwarfs, or exotic compact objects (Cardoso and Pani 2019) (e.g. boson stars; Guo et al. 2019b) may be misinterpreted as sub-solar PBHs. Distinguishing PBHs from other compact objects requires taking into account tidal disruption and tidal deformability measurements. Less compact objects like brown and white dwarfs are expected to be tidally disrupted well before the contact frequency. This leads to a different GW signature, see Franciolini et al. (2022b) and references therein for more details on the observability of these effects.

Another key prediction of the primordial model involves the eccentricity *e* of PBH binaries. While formed with large eccentricity at high redshift, PBH binaries then have enough time to circularize before the GW signal can enter the observation band of current and future detectors. Therefore, observing a non-zero eccentricity would rule out the interpretation as a primordial binary formed in the early Universe, while it may still be compatible with a PBH binary formed in the late-time universe (Cholis et al. 2016; Wang and Nitz 2021). Finally, the spin distribution of PBHs inherits characteristic mass-spin correlations induced by accretion effects, as discussed in Sect. 3.10. Using this criterion for determining the possible primordial nature of individual GW events would require reducing uncertainties on the accretion model. However, searching for signatures of this correlation can also be performed on a population level. This was recently done in Franciolini and Pani (2022) focusing on the GWTC-3 catalog, finding a similar correlation (indistinguishable at this stage from the one present in dynamically formed astrophysical binaries) may be present in the data.

In Sect. 6, dedicated to correlating GWs with LSS, we will mention further possibilities of distinguishing astrophysical and primordial merger populations through measurements of the event bias.

## 4.5 Summary

The calculation of PBH merger rates or encounter rates has been a very active line of research during the last few years, with some important twists. Two main binary formation channels have been envisaged: early binaries that formed quickly in the radiation era when two PBHs formed sufficiently close to each other and late binaries that can form inside PBH clusters, when two PBHs pass sufficiently close to each other to form a bound system. Today, it is still unclear which binary formation channel provides the most important merger rates on the mass scales probed by LVK, especially since these rates depend on the exact shape of the PBH mass function, on the dark matter fraction made of PBHs, on the clustering properties, on accretion, possible hierarchical mergers, on binary disruptions, etc. It is also debated if the limits on those rates coming from the latest GW observations favor or disfavor PBHs as the dark matter.

We tried to list all these limitations and we now attempt to summarize the current status: at first, the merger rates of early binaries were found to importantly surpass the

<footnote n="8">⁸ See, however, Shandera et al. (2018), Kouvaris et al. (2018) for models in which subsolar BHs are born out of dark sector interactions or capture of asteroidal mass PBHs (Oncins et al. 2022).</footnote>
ones of late binaries in clusters, at least in the standard scenario of cosmological
structure formation. However, on the one hand, N-body simulations have revealed
several mechanisms that importantly suppress those rates if PBHs constitute a
significant fraction of the dark matter, due to binary disruption by other PBHs, by
matter fluctuations, and above all by the PBH clusters seeded by the inevitable Pois-
son fluctuations in the initial PBH distribution. On the other hand, the same effects
should be responsible for a boost of the merger rates of late binaries in these PBH
clusters, at a level comparable to that of early binaries on solar-mass scales, if PBHs
significantly contribute to the dark matter. An important difference between early and
late binary merger rates is the mass dependence, roughly going like $1/m_{\text{PBH}}^{32/37}$ for
early binaries and $1/m_{\text{PBH}}^{11/21}$, by taking into account the halo mass distribution, for late
binaries (Bird et al. 2016). Another difference is the redshift evolution of those rates.
In both cases, however, it seems very disfavored that all the dark matter is made of
$30M_{\odot}$ black holes. The scenario where the PBH mass function peaks at the solar
mass scale, where the QCD transition should have favored PBH formation, may
represent an alternative possibility, provided the overall abundance is below the
current observational constraints in that mass range, see in particular Petač et al.
(2022), Gorton and Green (2022), Juan et al. (2022), De Luca et al. (2022b),
Franciolini et al. (2022e).

# 5 Stochastic backgrounds

## 5.1 From second-order curvature fluctuations

If PBHs are generated by the collapse of large density perturbations, they are
unavoidably associated to the emission of induced GWs at second order by the same
scalar perturbations, due to the intrinsic nonlinear nature of gravity (Tomita 1975;
Matarrese et al. 1994; Acquaviva et al. 2003; Mollerach et al. 2004; Ananda et al.
2007; Baumann et al. 2007a).<fnref n="9" /> The origin of these multimessenger signals is
therefore not coming from PBHs themselves but from the underlying cosmological
curvature perturbations (including ones that do not lead to PBHs. As a consequence,
the conclusions will be robust to changes in the astrophysical evolution of PBHs
(accretion and merger history) and clustering effects (Ünal et al. 2021). The
phenomenological implications have been investigated in various contexts also
associated to PBHs (Saito and Yokoyama 2009; Bugaev and Klimai 2010; Saito and
Yokoyama 2010; García-Bellido et al. 2017; Ando et al. 2018; Bartolo et al.
2018, 2019c, b; Clesse et al. 2018; Ünal 2019; Wang et al. 2019; Domènech 2020;
Domènech et al. 2020, 2022; Pi and Sasaki 2020; Ragavendra et al. 2021; Fumagalli

<footnote n="9">⁹ We stress that we do not refer here to the GWs produced only in the overdense regions that contract at horizon crossing due to the gravitational force (De Luca et al. 2019b), but from everywhere in the Universe, due to the general increase of the power of the density perturbations at small scales. See also here Lin et al. (2020), Chen et al. (2021a), Kawai and Kim (2021), Lin et al. (2023), Papanikolaou et al. (2022), Ivanov et al. (2022), Zhang (2022), Yi (2023), Cheong et al. (2022), Feng et al. (2023), Zhang et al. (2022), Arya et al. (2024), Tzerefos et al. (2023), Papanikolaou et al. (2024c) for studies within modified gravity setups.</footnote>
et al. 2021a, 2022a; Yuan et al. 2019) (see recent reviews by Yuan and Huang 2021a; Domènech 2021). This will be also a conclusive test for the intriguing proposal that the seeds of the SMBHs are formed by PBHs (Duechting 2004; Belotsky et al. 2014; Clesse and García-Bellido 2015; Nakama et al. 2016; García-Bellido et al. 2016).

Such a source of GWs is present at all scales, but due to the necessary enhancement of the scalar power spectrum responsible for the generation of PBHs around the characteristic scale $k_*$, it will peak in the frequency range around $f_* = k_*/2\pi$, and can thus reach the detection sensitivities of many GW experiments like LISA. Since the emission mostly occurs when the corresponding scales cross the horizon, one can relate the GW frequency to the PBHs mass $m_{PBH}$ as (see for example Saito and Yokoyama 2010; García-Bellido et al. 2017)

$$f \approx 6 \, \text{mHz} \sqrt{\gamma} \left( \frac{m_{\text{PBH}}}{10^{-12} M_{\odot}} \right)^{-1/2}, \qquad (116)$$

where the parameter $\gamma$ describes the efficiency of the gravitational collapse into PBHs.

There are several current and future experiments searching for a SGWB in various frequency ranges. In the ultra-low frequency range (around nHz), the observations at Pulsar Timing Array (PTA) experiments like NANOGrav (Agazie et al. 2023b, c), EPTA (in combination with InPTA) (Antoniadis et al. 2024, 2023a, b), PPTA (Rardon et al. 2023a; Zic et al. 2023; Reardon et al. 2023b) and CPTA (Xu et al. 2023), give rise to the most stringent constraints on the GWs abundance. Future experiments like SKA (Dewdney et al. 2009) (see also Moore et al. 2015) will greatly improve the sensitivity.

In the LVK frequency range, an additional constraint has been set by the non-observation of a stochastic background after O1-O2 runs (Abbott et al. 2019a). These searches can be translated into a constraint on the amplitude of the comoving curvature perturbation at the corresponding scales (Inomata and Nakama 2019; Ünal et al. 2021). Those bounds are also affecting the maximum allowed PBH fraction of DM with the hypothesis that they originate from the collapse of density perturbations. Detailed studies with the LVK data affecting the mass range $[10^{-20}, 10^{-18}] M_{\odot}$ are reported in Kapadia et al. (2020). Also, in Cai et al. (2019b), the dependence of the result on non-Gaussianities is also investigated, finding that local non-Gaussianity can for example alleviate the bounds (see Sect. 5.3 for details). This is possible for two main reasons: (i) Non-Gaussianity can allow PBHs to be produced more efficiently for the same power spectra, as the tail of the non-Gaussian probability distribution function of the perturbations has more area/probability for the perturbations greater than the threshold. (ii) Due to several contractions of the curvature field within non-Gaussian correlators, there is a higher symmetry factor compared to the Gaussian case. Hence, non-Gaussian perturbations can produce a large fraction of PBHs with a smaller amplitude power spectrum, and therefore with a smaller amount of sourced GWs.

SKA (Moore et al. 2015; Zhao et al. 2013) is the next-generation PTA experiment that will probe primordial perturbations very sensitively (Byrnes et al. 2019; Inomata and Nakama 2019; Kalaja et al. 2019; Gow et al. 2020). It has been recently shown
in Ünal et al. (2021) that PTA-SKA combined with CMB distortions will robustly test the PBHs from inflationary fluctuations, namely they will detect the stochastic GW background or distortion signatures and possibly make an extraordinary discovery or constrain the PBHs heavier than a solar mass for 13 orders of magnitude (1 − 10$^{13}$M$_{⊙}$) to a completely negligible amount $^{p_{PBH}}/p_{DM}$ < 10 − $^{10}$.

Finally, the LISA experiment will be able to provide insights for the intermediate frequencies around the mHz, such that the peak of the signal falls within the LISA sensitivity band for PBH masses around m$_{PBH}$ ∼ O (10 − 15 − 10 − $^{8}$)M$_{⊙}$ (Saito and Yokoyama 2009; García-Bellido et al. 2017; Bartolo et al. 2019c; Cai et al. 2019a; Ünal 2019; Bhaumik and Jain 2021).

In the following, we review the procedure to compute the induced SGWB spectrum from primordial scalar perturbations. First of all, let us define the scalar and tensor perturbations in the Newtonian gauge as

$$
ds^2 = a^2 \left\{ -(1 + 2\Psi) d\eta^2 + \left[ (1 - 2\Psi) \delta_{ij} + \frac{1}{2} h_{ij} \right] dx^i dx^j \right\}, \quad (117)
$$

where we neglected the anisotropic stress. The GWs emission is captured by the
equation of motion for the GWs<fnref n="10" /> as

$$
h''_{ij} + 2Hh'_{ij} - \nabla^2 h_{ij} = -16T^{lm}_{ij} \left[ \Psi \partial_l \partial_m \Psi + 2 \partial_l \Psi \partial_m \Psi - \partial_l \left( \frac{\Psi'}{H} + \Psi \right) \partial_m \left( \frac{\Psi'}{H} + \Psi \right) \right], \quad (118)
$$

where derivatives are taken with respect to conformal time η, ℋ ≡ a'/a is the
conformal Hubble parameter, and the source is evaluated assuming a radiation-
dominated (RD) epoch,<fnref n="11" /> The solution of the first order equation of motion in a
radiation-dominated Universe relates the Bardeen’s potential Ψ to the gauge
invariant comoving curvature perturbation through (Lyth and Riotto 1999)

$$
\Psi(\eta, \vec{k}) = \frac{2}{3} T(k\eta)\zeta(\vec{k}), \quad \text{where} \quad T(z) = \frac{9}{z^2} \left[ \frac{\sin(z/\sqrt{3})}{z/\sqrt{3}} - \cos(z/\sqrt{3}) \right]. \tag{119}
$$

One can decompose the tensor field in Fourier space as

$$
h_{ij}(\eta, \vec{x}) = \int \frac{d^3 k}{(2\pi)^3} \sum_{\lambda=R,L} h_{\lambda}(\eta, \vec{k}) e_{ij,\lambda}(\hat{k}) E^{\lambda\cdot\hat{k}}, \quad (120)
$$

where $h_\lambda$ are the two helicity modes. Introducing the dimensionless variables $x = p/k$ and $y = |\vec{k} - \vec{p}|/k$, the GWs emitted take the form, see for example (Espinosa et al. 2018),

<footnote n="10">¹⁰ We do not consider the free-streaming effect of neutrinos on the GW amplitude (Weinberg 2004).</footnote>

<footnote n="11">¹¹ The projector $\mathcal{T}_{ij}^{lm}$ to the transverse and traceless component is defined in Fourier space in terms of the polarisation tensors $e_{ij}^l(\vec{k})$ the chiral basis (Lr) as $\tilde{\mathcal{T}}_{ij}^{lm}(\vec{k}) = e_{ij}^l(\vec{k}) \otimes e^{lm}(\vec{k}) + e_{ij}^R(\vec{k}) \otimes e^{Rlm}(\vec{k})$. The normalisation adopted is $e_{ij,\lambda}(\vec{k})e_{ij,\lambda'}^*(\vec{k}) = \delta_{\lambda\lambda'}$.</footnote>
$$h_{\lambda}(\eta, \vec{k}) = \frac{4}{9k^3\eta} \int \frac{d^3p}{(2\pi)^3} E_{\lambda}^*(\vec{k}, \vec{p})\zeta(\vec{p})\zeta(\vec{k} - \vec{p})[\mathcal{I}_c(x, y) \cos(k\eta) + \mathcal{I}_s(x, y) \sin(k\eta)], \quad (121)$$

where $E_{\lambda}(\vec{k}, \vec{p}) \equiv E_{ij,\lambda}(\vec{k})\vec{p}_i\vec{p}_j$, and two functions $\mathcal{I}_{c,s}$ are (Espinosa et al. 2018; Kohri and Terada 2018)

$$\mathcal{I}_c(x, y) = -36\pi \frac{(s^2 + d^2 - 2)^2}{(s^2 - d^2)^3} \theta(s - 1) \quad (122)$$

$$\mathcal{I}_s(x, y) = -36 \frac{(s^2 + d^2 - 2)}{(s^2 - d^2)^2} \left[ \frac{(s^2 + d^2 - 2)}{(s^2 - d^2)} \log \frac{(1 - d^2)}{|s^2 - 1|} + 2 \right], \quad (123)$$

with

$$d \equiv \frac{1}{\sqrt{3}}|x - y|, \quad s \equiv \frac{1}{\sqrt{3}}(x + y), \quad (d, s) \in [0, 1/\sqrt{3}] \times [1/\sqrt{3}, +\infty). \quad (124)$$

The definition of the energy density associated to GWs is (Misner et al. 1973; Maggiore 2000; Flanagan and Hughes 2005)

$$\rho_{\text{GW}} = \frac{M_{\text{Pl}}^2}{4} \langle \dot{h}_{ab}(t, \vec{x}) \dot{h}_{ab}(t, \vec{x}) \rangle_T, \quad (125)$$

where the angular brackets denote a time average performed on a timescale $T$ much greater than the GW phase oscillations ($Tk_i \gg 1$) but much smaller than the cosmological time ($TH \ll 1$).

Assuming that the scalar perturbations $\zeta$ are Gaussian (see the next sections for a possible relaxation of this assumption), one finds

$$\begin{aligned} \langle \rho_{\text{GW}}(\eta, \vec{x}) \rangle &\equiv \rho_c(\eta) \int d\ln k \, \Omega_{\text{GW}}(\eta, k) \\ &= \frac{2\pi^4 M_{\text{P}}^2}{81\eta^2 a^2} \int \frac{d^3k_1 d^3p_1}{(2\pi)^6} \frac{1}{k_1^4} \frac{[p_1^2 - (\vec{k}_1 \cdot \vec{p}_1)^2 / k_1^2]^2}{p_1^3 |\vec{k}_1 - \vec{p}_1|^3} \mathcal{P}_{\zeta}(p_1) \mathcal{P}_{\zeta}(|\vec{k}_1 - \vec{p}_1|) \quad (126) \\ &\quad [\mathcal{I}_c^2(\vec{k}_1, \vec{p}_1) + \mathcal{I}_s^2(\vec{k}_1, \vec{p}_1)], \end{aligned}$$

as a function of the curvature perturbation power spectrum $\mathcal{P}_{\zeta}$. In the first line, we have defined the fractional GW energy density for log interval $\Omega_{\text{GW}}$, in terms of the critical energy density of a spatially flat Universe, $\rho_c = 3H_0^2 M_{\text{P}}^2$.

The current GW abundance can then be obtained as
$$
\Omega_{\text{GW}}(\eta_0, k) = \frac{a_f^4 \rho_{\text{GW}}(\eta_f, k)}{\rho_r(\eta_0)} \Omega_{r,0} = \frac{g_*(\eta_f)}{g_*(\eta_0)} \left( \frac{g_{*s}(\eta_0)}{g_{*s}(\eta_f)} \right)^{4/3} \Omega_{r,0} \Omega_{\text{GW}}(\eta_f, k), \quad (127)
$$

in terms of the present radiation energy density fraction $\Omega_{r,0}$ if the neutrinos were massless and the effective degrees of freedom for entropy density $g_{*s}$. For the frequencies of interest, using the relation with the temperature of the universe at the time of horizon crossing $f \simeq 8$ mHZ($g_*/10)^{1/4}$ ($T/10^6$ GeV), one can show that the emission of GWs takes place at $\eta_f$ (with corresponding scale factor $a_f$) well before the time at which top quarks start annihilating, above which we can assume a radiation-dominated Universe with constant effective degrees of freedom.

As one can see, the GW's abundance depends on the curvature perturbation power spectrum. One possibility is represented by a monochromatic power spectrum with support at the momentum scale $k_★$

$$
\mathcal{P}_{\zeta}(k) = A_s k_{\star} \delta(k - k_{\star}), \quad (128)
$$

enhanced with respect to the power spectrum on large CMB scales, for which the current GWs abundance can be computed analytically as (see also Saito and Yokoyama 2010; Bugaev and Klimai 2010)

$$
\begin{equation}
\begin{split}
\Omega_{\text{GW}}(\eta_0, k) ={}& \frac{\Omega_{r,0} A_s^2 g_*(\eta_0)}{15552 g_*(\eta_0)} \left( \frac{g_{*S}(\eta_0)}{g_{*S}(\eta_\mathrm{f})} \right)^{4/3} \left( \frac{4k_\star}{k} - \frac{k}{k_\star} \right)^2 \\
& \qquad \times \theta(2k_\star - k) \left[ \mathcal{I}_c^2 \left( \frac{k_\star}{k}, \frac{k_\star}{k} \right) + \mathcal{I}_s^2 \left( \frac{k_\star}{k}, \frac{k_\star}{k} \right) \right],
\end{split}
\tag{129}
\end{equation}
$$

where $\theta$ is the Heaviside step function.

A more realistic spectrum is provided by a log-normal shape with width $\sigma$

$$
\mathcal{P}_{\zeta}(k) = A_{\zeta} \exp\left(-\frac{\ln^2(k/k_{\star})}{2\sigma^2}\right), \quad (130)
$$

whose abundance is shown by the black line of Fig. 24, compared to the monochromatic shape in red. For illustrative purposes, the power spectra peak scale has been chosen to be the one at which LISA has its maximum sensitivity, i.e. $f_{\star} = f_{\text{LISA}} \sim 3.4$ mHZ where $f_{\star} \equiv k_{\star}/2\pi$, corresponding to PBH with masses around $10^{-12} M_{\odot}$, while the parameters were chosen in order to have PBHs which account for the totality of the DM. One can see that the curves fall well within the sensitivity curves of LISA. It is thus clear that if PBHs of such masses, produced from the collapse of sizeable primordial fluctuations, represent the totality (or a fraction) of the dark matter, LISA will be able to measure the GWs sourced during the PBH formation time for a wide variety of PBH masses.

The spectral shapes present differences at high frequencies, because of the absence of the sharp cut-off at $2f_{\star}$ characteristic of the Dirac delta case, and at lower frequencies, since physical spectra would typically give a white-noise ($\propto f^3$) behaviour (Espinosa et al. 2018; Cai et al. 2020), while the slow fall-off of the monochromatic case is an unphysical effect of assuming such a power spectrum.
<figure><img src="image_27.png" /><figcaption>**Fig. 24** Induced GWs spectra for both examples considered in Eqs. (128) and (130), for the parameters choice $A_s = 0.033$, $A_\zeta = 0.044$ and $\sigma = 0.5$. For comparison, we also show the estimated sensitivity for LISA (Amaro-Seoane et al. 2017), following the proposed design (4y, 2.5 Gm of length, 6 links). The PBH mass corresponding to the characteristic frequency is depicted on the top horizontal axis, according to Eq. (116). Image reproduced with permission</figcaption></figure>

Furthermore, in the latter a resonant effect at $f \sim 2f_\star/\sqrt{3}$ produces the spike, see for example Ananda et al. (2007).

As one can appreciate from Eq. (121), the non-linear coupling to the curvature perturbation naturally leads to an intrinsically non-Gaussian GW signal imprinted in phase correlations. However, it has been shown in Bartolo et al. (2019b, 2019c) (see also Bartolo et al. 2018; Margalit et al. 2020), that the coherence is washed out by the propagation of the waves in the perturbed Universe mainly due to time delay effects originated from the presence of large scale variations of gravitational potential. This is simply a consequence of the central limit theorem applied to a number $N \sim (k_\star\eta_0)^2 \gg 1$ of independent lines of sight (Bartolo et al. 2019b). Possible small deformations smearing the GW spectrum can also arise from similar effects (Domcke et al. 2020).

## 5.2 Gauge invariance of the SGWB spectrum

The fact that tensor modes are generated at second-order in perturbation theory raises the issue of the possible gauge dependence of the result commonly computed in the Newtonian gauge. Indeed, second order tensor modes, contrary to the first order ones, are not gauge invariant. This issue has been recently highlighted in the literature (Hwang et al. 2017; Domènech and Sasaki 2018; Gong 2022; Tomikawa and Kobayashi 2020; Lu et al. 2020).

A physical observable is however not dependent on gauge choices by definition and it can be identified by understanding how the measurement in GW experiments is performed. The description of the detector's response to the GW signal can be best performed in the so-called transverse-traceless (TT) frame (Maggiore 2008; De Luca et al. 2020b), particularly in space based experiments like LISA. The TT frame is
defined as the one where the coordinates are fixed with the position of the mirrors
and the effect of the passing GWs is captured by the delay of the arrival times
between the arms of the experiment. The peculiarity of the LISA experiment is due to
the relation between the characteristic frequency observable and the arm length L,
which is $f_{\text{LISA}}L \sim O(1)$. The TT choice is optimal because it avoids keeping track of
the otherwise large corrections $O(fL)$ appearing if expanding around a locally
inertial reference frame. Indeed, the projected sensitivity curves for the interferometer
LISA are provided in such a frame.

One should remark that, during the period in which the source of the tensor mode
is active and therefore the tensor field is coupled to the scalar perturbations at second-
order in perturbation theory, it retains a gauge dependence and cannot be identified
with the freely-propagating GW, which can be treated linearly and is therefore gauge
invariant.

In the case of the production of tensor degrees of freedom from scalar perturbation
during a radiation-dominated phase of the Universe, the source is active only when
the scalar perturbations re-enter the horizon. Once produced, the GWs effectively
decouple from the second order source and behave like linear perturbations of the
metric in the late-time limit well within the cosmological horizon. Therefore, the
initial gauge dependence of the result is lost (De Luca et al. 2020b; Inomata and
Terada 2020; Yuan et al. 2020) and does not affect the spectrum shown in Fig. 24.

5.3 Impact of primordial non-Gaussianity on scalar induced SGWB

As discussed above, the existence of amplified primordial fluctuations at small scales
can lead to PBHs and also accompanying scalar induced (secondary) SGWB.
Depending on the amplification/enhancement scale, this background can be detected
at different frequencies with various GW experiments including LISA, PTA-SKA,
DECIGO and BBO. Since the physical scales for PBHs of mass less than million
solar mass correspond to very small sub-galactic distances
(perturbation wavelength <kpc), the nonlinearity of scalar fluctuations in the current
Universe prevents us from extracting primordial information about such small scales.
However, the freely propagating SGWB sourced by primordial curvature fluctuations
could be an excellent probe for these scales. The spectral shape of the scalar induced
(or secondary) SGWB can be influenced by the statistical properties of the scalar
density perturbations. Therefore, the detection and spectral analysis of the induced/
secondary GW background can also shed light on the properties of the inflationary
era.

Since the SGWB is sourced directly by the primordial curvature fluctuations in the
horizon re-entry, it is natural to expect to observe some statistical properties of these
scalar fluctuations in the SGWB spectrum. The SGWB includes a 4-pt function of the
primordial curvature perturbations ($P_h \propto \langle h^2 \rangle \propto \langle \zeta^4 \rangle$ see also Eqs. (121) and (126)),
hence for a Gaussian case, it can be reduced to only powers of a 2-pt function (see
left panel of Fig. 25). In the presence of non-Gaussianity, a 4-pt correlator can not be
written in terms of power spectra exactly.
Curvature perturbations with non-Gaussianity that can be expressed in terms of quadratic order Gaussian curvature perturbations contribute to GWs with 5 more diagrams given in Figs. 25 and 26. Note that the interaction could be of local form, or might include time and position derivatives. Hence, the topologies of these five diagrams are independent of the interaction type and universal for curvature perturbations with the form

$$ \zeta = \zeta_g + O \cdot \zeta_g^2. \qquad (131) $$

Here $O$ is an operator acting on a quadratic order Gaussian field. It can be constant or depend on space and time. For two powers of interaction, we have 2 diagrams in Fig. 25

* Hybrid, since the power spectrum of one leg is tree level and the other one is the interaction;

* Walnut, due to the walnut shape. Note that the walnut topology diagram can be obtained by two different contractions.

For 4 powers of interaction, there are 3 diagrams given in Fig. 26

* Reducible, since this diagram can be reduced to a 1 loop diagram;

* Planar and Non-planar, respectively, as their topology can or cannot be drawn on a plane.

The properties of the spectrum have been studied in a number of references. For instance, Nakama et al. (2017) estimated the order of magnitude of the SGWB in the presence of non-Gaussianity. In García-Bellido et al. (2017), the SGWB spectrum was computed in the presence of large non-Gaussianities, and in such a case the main contribution to the SGWB comes from the 3 diagrams given in Fig. 26. The SGWB spectrum in the presence of small non-Gaussianity has been computed in Ünal (2019) with the diagrams in Figs. 25 and 26. Small non-Gaussianities have also been studied in Atal and Domènech (2021), Cai et al. (2019a).

It is important to note that Walnut diagram topology has two distinct contributions from two distinct contractions. This missing point in previous studies has been identified and computed carefully in Adshead et al. (2021). With these results, the whole diagrams have been completed up to quartic order in interaction ($f_{NL}^4$) and quadratic order in scalar density perturbations. The scalar induced SGWB with scalar non-Gaussianities has been worked out in Yuan and Huang (2021b) up to third order

<figure><img src="image_28.png" /><figcaption>Fig. 25 Nonvanishing diagrams for $O(f_{NL}^0)$ (left) and $O(f_{NL}^2)$: center and right, called Hybrid and Walnut diagrams due to their topology and shape. Images reproduced with permission from Ünal (2019), copyright by APS Image reproduced with permission</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>**Fig. 26** Diagrams of $O(f_{nl}^4)$, called (from left to right) Reducible, Planar and non-Planar due to their topological properties. Images reproduced with permission from Ünal (2019), copyright by APS; see also García-Bellido et al. (2017) Image reproduced with permission</figcaption></figure>

in density perturbations. The contribution of higher order in perturbation theory do not change the peak amplitude, but since higher order correlation includes convolution of higher number of propagators, the maximum frequency of external legs grow linearly with this number. Therefore, instead of sharp decrease in the high frequency regime, there are corrections at high frequencies from higher-order terms which slows down the decrease of the spectrum.<fnref n="12" />

The following results about the imprints of non-Gaussianity on the secondary/
induced GW spectrum have been found:

* Non-Gaussianity can enhance the production efficiency of the SGWB. Namely, for the same amplitude of Gaussian or non-Gaussian scalar fluctuations, the resulting SGWB has a larger amplitude for non-Gaussian than for Gaussian fluctuations.

* Because non-Gaussianity is a result of interactions, this leads to a convolution in the momentum space and to a broader distribution of scalar fluctuations, with a peak at frequencies about $O(1)$ higher, especially if the curvature fluctuations have a narrow spectrum.

* In the limit of large non-Gaussianities, the diagrams in Fig. 26 are important and dominate the signal. For small and mild non-Gaussianities, the signal is dominated by the diagrams in Fig. 25.

Some of the results above can be visualized in Fig. 27, where $\Omega_{GW}$ indicates the total GW spectrum, $\Omega_{GW}^g$ indicates GWs from Gaussian scalar perturbations, and $\Omega_{GW}^{\text{NG}}$ GWs from non-Gaussian scalar perturbations.<fnref n="13" /> Due to the interactions of the scalar perturbations, the high frequency regime of the SGWB is dominated by the NG scalar perturbations. This fact can lead to a second peak whose amplitude is set by the strength of the interactions, ie. the non-Gaussianity parameter $f_{\text{NL}}$. If the NG contribution is considerable, the SGWB leads to two peaks, as shown in the right panel of Fig. 27. The second peak, which results from interactions/non-Gaussian perturbations, is at a frequency about $O(1)$ times higher with respect to the first peak from Gaussian perturbations. It can be seen on the left panel of the same Figure that

<footnote n="12">¹² Explicitly, $\zeta + \bar{\zeta} \rightarrow h$ interaction allows GW to have highest frequency at $2k_{\text{UV}}$ due to momentum conservation, $k_{\text{UV}}$ is the highest momentum in scalar perturbations, on the other hand for $\zeta + \bar{\zeta} + \zeta \rightarrow h$ interaction, GW has the highest frequency at $3k_{\text{UV}}$.</footnote>

<footnote n="13">¹³ Recall that induced GWs are sourced by scalar sources at quadratic order hence they are inherently non-Gaussian. By G and NG symbols we indicate the statistical properties of scalar perturbations sourcing GWs.</footnote>
<figure><img src="image_30.png" /><figcaption>**Fig. 27** Second peak produced by primordial non-Gaussian component of curvature perturbations. Image reproduced with permission from Ünal (2019), copyright by APS Image reproduced with permission</figcaption></figure>

even if the NG contribution is tiny, it still dominates the large frequency regime of the GW spectrum and produces an abrupt second bump at smaller amplitudes which can be detected with sensitive next generation GW experiments, such as LISA, Cosmic Explorer, DECIGO and PTA-SKA, and allow probing non-Gaussianity as sensitive as $f_{\text{NL}} \sim 0.5$ (Ünal 2019). This level of non-Gaussianity could be probed even better than using next generation CMB and LSS experiments.

In Garcia-Saenz et al. (2023), the impact of primordial non-Gaussianity on the scalar induced SGWB has been studied beyond the generalized local ansatz (131), paying attention to whether it can self-consistently be important in concrete early-universe scenarios. As pointed out above, the gravitational waves spectrum can be split into two contributions: the one set by the power spectrum, and the one determined by the connected four point function of the curvature perturbation, namely the primordial trispectrum. First, it is shown that this trispectrum-induced gravitational wave spectrum can always be written as a sum of three “channels” contributions, independently of the precise shape of the trispectrum. Second, the question addressed in Garcia-Saenz et al. (2023) is: can the trispectrum be observationally relevant for the scalar induced SGWB in conventional set-ups where non-Gaussianity provides subleading corrections to the Gaussian scalar signal?, i.e. in inflationary scenarios where interactions generating non-Gaussianity maintain perturbative control. Formulating a precise quantitative criterion ensuring perturbative control in strongly scale-dependent theories is difficult and model-dependent. Hence, the focus of Garcia-Saenz et al. (2023) is on scale-invariant theories. In this context, it is shown that neither regular trispectrum shapes peaking in the so-called equilateral configurations, nor local trispectrum shapes diverging in soft momentum limits, can contribute significantly. Indeed, those contributions are always bound to be smaller than an order-one (or smaller) number multiplying the relative one-loop correction to the scalar power spectrum, necessarily much smaller than unity in order for the theory to be under perturbative control. This result is shown to be also valid in a toy-model for the phenomenologically more relevant situation of a scale-dependent scalar spectrum, calling for more in-depth investigations of this question.
5.4 Induced SGWB anisotropies

It is interesting to note that the induced SGWB associated with PBH production<fnref n="14" /> will present some angular anisotropies, which can be quantified by computing the two-point correlation function of the density field in different directions. The LISA experiment will be able to detect anisotropies with an angular resolution of $ℓ \lesssim 15$ (Contaldi et al. 2020) (see also Baker et al. 2021a), which will require the spatial points to be separated by a non-negligible fraction of the present horizon. Since the characteristic scales of the perturbations generating the SGWB are much smaller than those distances, $k_* |\vec{x} - \vec{y}| \gg 1$, and that the emission takes place near horizon crossing with the assumption of Gaussian scalar perturbations, the Equivalence Principle dictates that the anisotropies decay as $(k_* |\vec{x} - \vec{y}|)^{-2}$. Therefore, the anisotropy coarse grained with the resolution of the experiment will be undetectable (Bartolo et al. 2020a). Another possibility, as far as anisotropies imprinted at formation are concerned, is given by the presence of non-Gaussianity which, by providing a correlation between small and large scales, may allow large scale modulation of power and lead to a large-scale anisotropy (Bartolo et al. 2020a). Finally, the propagation of GWs across disconnected regions of the Universe leads to large-scale anisotropies at detection (Alba and Maldacena 2016; Contaldi 2017; Cusin et al. 2017, 2019a; Jenkins and Sakellariadou 2018; Bartolo et al. 2019a, 2020b; Renzini and Contaldi 2019; Bertacca et al. 2020).

Assuming a local, scale-invariant, shape of non-Gaussianity $\zeta = \zeta_g + \frac{3}{5}f_{NL}\zeta_g^2$ and taking into account propagation effects, one can compute the two-point correlation function of the direction dependent GW energy density contrast $\delta_{GW}(\eta, \vec{x}, k) = \Omega_{GW}(\eta, \vec{x}, k)/\Omega_{GW}(\eta, k) - 1$, defined as the relative difference between the direction dependent GW energy density and its monopole contribution, in a spherical harmonics decomposition as (following the notation used in Bartolo et al. 2020a)

$$ \langle \delta_{GW,ℓm} \delta_{GW,ℓ'm'}^* \rangle = \delta_{ℓℓ'} \delta_{mm'} \hat{C}_ℓ(k), \qquad (132) $$

where (Bartolo et al. 2020a)

$$ \sqrt{\frac{\ell(\ell+1)}{2\pi}} \hat{C}_{\ell}(k) \approx \frac{3}{5} \left| 1 + \tilde{f}_{\text{NL}}(k) \right| \left| 4 - \frac{\partial \ln \Omega_{\text{GW}}(\eta, k)}{\partial \ln k} \right| \mathcal{P}_{\zeta_{\text{L}}}^{1/2}, \qquad (133) $$

and $k$ denotes the momentum on small-scales. Here $\mathcal{P}_{\zeta_{\text{L}}}$ denotes the power spectrum at large scales, while the momentum dependent non-Gaussian parameter is defined as

<footnote n="14">¹⁴ The generation of a SGWB in the early Universe may also happen due to other mechanisms, possibly with characteristic anisotropic signatures, see for example García-Bellido et al. (2016), Cook and Sorbo (2012), Guzzetti et al. (2016), Bartolo et al. (2016), Ricciardone and Tasinato (2018), Caprini and Figueroa (2018), Geller et al. (2018), Dimastrogiovanni et al. (2020).</footnote>
$$ \tilde{f}_{\text{NL}}(k) \equiv \frac{8 f_{\text{NL}}}{4 - \frac{\partial \ln \Omega_{\text{GW}}(\eta, k)}{\partial \ln k}}. \qquad (134) $$

The presence of non-Gaussianity in the curvature perturbation, constrained by the Planck collaboration to be $-11.1 \le f_{\text{NL}} \le 9.3$ at 95% C.L. (Akrami et al. 2020a), would also induce a significant variation on large scales of the PBH abundance through the modulation of the power on small scales induced by the long modes. If PBHs compose a non-negligible fraction of the DM, the presence of isocurvature modes in the DM density fluid associated to the non-Gaussianity (Young and Byrnes 2015) is strongly constrained by the CMB observations (Akrami et al. 2020b).

In Fig. 28, we show the expected GW anisotropy for the choice of a monochromatic and log-normal small-scale power spectrum of the curvature perturbation peaked at the LISA maximum sensitivity frequency. The dot-dashed lines identify the corresponding GW abundance at the present epoch evaluated at the peak frequency, while the colored region identifies the range of parameters allowed by the Planck constraints, and we choose the non-linear parameter $f_{\text{NL}} \ge 1/3$ to avoid the failure of the perturbative approach in the computation of the PBH abundance happening at larger negative values, see for example (Young and Byrnes 2013; Yoo et al. 2019). One finds that if a large fraction of the dark matter is composed of PBH, one expects a highly isotropic and Gaussian SGWB, up to propagation effects. On the contrary, the detection of a large amount of anisotropy in the signal associated to the PBH formation scenario would imply that only a small fraction of the DM can be accounted by PBHs (Bartolo et al. 2020a).

## 5.5 From PBH poisson fluctuations

At small scales, the gravitational interaction between individual PBHs (scattering, merging, etc.) results in the production of GWs. In this section, we are interested in the large-scale counterpart of this signal, i.e. at distances much larger than the mean separation length between PBHs, for which PBHs can be described by a pressureless fluid. This fluid is endowed with density perturbations, which can be treated within

<figure><img src="image_31.png" /><figcaption>Fig. 28 Contour plot showing the amount of GWs anisotropy in the parameter space of $f_{\text{PBH}}$ and $f_{\text{NL}}$ allowed by the Planck constraints for the choice of a monochromatic and log-normal small-scale power spectrum, respectively, with peak frequency around the maximum sensitivity of LISA. The dot-dashed lines identify the corresponding present GWs abundance. Image reproduced with permission from Bartolo et al. (2020a), copyright by IOP and SISSA Image reproduced with permission</figcaption></figure>
the framework of cosmological perturbation theory. In particular, these PBH density perturbations induce the production of GWs through second-order effects studied both in the context of general relativity (Papanikolaou et al. 2021; Domènech et al. 2021) as well in the context of modified gravity (Papanikolaou et al. 2022, 2023b).

To be explicit, we assume that PBHs form in a radiation-dominated, homogeneous Universe (at large scales, primordial curvature perturbations provide a contribution to inhomogeneities but they are negligible compared to the ones generated by PBHs, at least in the range of scales we are interested in). One can therefore view PBHs as a dust fluid, formed from the transition of a fraction of the radiation energy density into PBH dust matter. Assuming that PBHs are randomly distributed at formation time (i. e. they have Poisson statistics), the energy density associated to them is inhomogeneous while the total energy density of the background is homogeneous. Consequently, the energy density perturbation of the PBH matter field can be described by an isocurvature Poisson fluctuation. If the initial abundance of PBHs is large enough, PBHs can potentially dominate the Universe energy density content and in that case the isocurvature PBH energy density perturbation in the radiation era is converted into an adiabatic curvature perturbation deep in the PBH dominated era. These early PBH dominated eras can be naturally driven by ultralight PBHs, which evaporate before BBN (García-Bellido et al. 1996; Hidalgo et al. 2012; Martin et al. 2020b; Zagorac et al. 2019) and addressing a plethora of cosmological problems, most importantly that of Hubble tension (Hooper et al. 2019; Papanikolaou 2023b).

At second-order in cosmological perturbation theory, the gravitational potential of these Poisson distributed PBHs induce a SGWB which may be detected by LISA (Papanikolaou et al. 2021; Domènech et al. 2021). More precisely, as found in Papanikolaou et al. (2021), this gravitational potential gives rise to the following power spectrum for Φ:

$$ \mathcal{P}_{\Phi}(k) = \frac{2}{3\pi} \left(\frac{k}{k_{UV}}\right)^3 \left(5 + \frac{4k^2}{9k_d^2}\right)^{-2}, \qquad (135) $$

where Φ is the Bardeen gravitational potential, $k_d$ is the comoving scale exiting the Hubble radius at PBH domination time and $k_{UV} \equiv a/\bar{r}$, where $\bar{r}$ corresponds to the mean PBH separation distance. Note that $k > k_{UV}$ corresponds to distances within the mean separation distance, where the granularity of the PBH matter field and the associated non-linear effects become important, and where the gas of PBHs cannot be described by a fluid anymore. This is why the above expression should be restricted to $k < k_{UV}$ with $k_{UV}$ acting as a UV cutoff (hence the notation). Since $\bar{r} \propto a$, $k_{UV}$ is a fixed comoving scale. From Eq. (135), one can see that $\mathcal{P}_{\Phi}$ is made of two branches: when $k \ll k_d$, $\mathcal{P}_{\Phi} \propto k^3$, while $\mathcal{P}_{\Phi} \propto 1/k$ when $k \gg k_d$. It reaches a maximum when $k \sim k_d$, where $\mathcal{P}_{\Phi}$ is of order $(k_d/k_{UV})^3$.

The energy density of the associated induced GWs can be computed according to the prescription described in Sect. 5.1, where one should use the linear transfer function during a matter (PBH) dominated era. One can show that it peaks at the comoving scale that exits the Hubble radius at the onset of the PBH domination era (Papanikolaou et al. 2021). This scale corresponds to the frequency displayed in Fig. 29 as a function of the fraction of the Universe made of PBHs at the time they
<figure><img src="image_32.png" /><figcaption>**Fig. 29** Frequency at which the GWs induced by a dominating gas of PBHs peak, as a function of their energy density fraction at the time they form, Ω<sub>PBH,f</sub> (horizontal axis), and their mass m<sub>PBH</sub> (colour coding). The region of parameter space that is displayed corresponds to values of m<sub>PBH</sub> and Ω<sub>PBH,f</sub>, such that black holes form after inflation, dominate the Universe content for a transient period and Hawking evaporate before BBN. We also impose that the induced GWs do not lead to a backreaction problem before they evaporate, see Papanikolaou et al. (2021) for more details. For comparison, the frequency detection bands of ET, LISA and SKA are shown. Image reproduced with permission from Papanikolaou et al. (2021), copyright by IOP and SISSA Image reproduced with permission</figcaption></figure>

form, and their mass (assuming an initial monochromatic distribution of PBHs). One
can see that a substantial fraction of parameter space lies within the LISA detection
band. This would provide a unique access to ultralight PBHs, which otherwise
cannot be detected through the PBH merging channel or the second-order induced
GW channel presented in Sect. 5.1.

Let us point out that the calculation of the precise amplitude of the signal is more
complex as it implies to follow its transition from the PBH-dominated era to the
subsequent radiation-dominated (RD) era. This is why Papanikolaou et al. (2021)
was followed up by Inomata et al. (2020), Domènech et al. (2021) where, by
considering a monochromatic PBH mass function and a sudden transition from the
PBH-dominated era to the radiation era after evaporation, it was found that the GW
signal discussed here is not suppressed due to the PBH → RD transition and it can be
detectable from LISA and other GW experiments. In fact, the fate of the signal
through the transition depends on whether the transition is sharp or slow. In
particular, on the one hand, when the transition is gradual, the gravitational potential
oscillates with a decreasing amplitude (Inomata et al. 2019b), suppressing in this way
the signal. On the other hand, if the transition is sudden, the gravitational potential
oscillates very fast with large amplitudes up to the end of the transition and therefore
the GW signal is expected to be enhanced in this case (Inomata et al. 2019a).
It is also important to mention that another calculation was performed in Kozaczuk et al. (2022), where the contribution from scales for which the density contrast becomes non linear during the PBH-dominated era was removed (note that the gravitational potential always remains linear). Assuming again a monochromatic PBH mass function, a reduction was found in the GW signal which removed it from the reach of LISA, while remaining accessible by other GW experiments like BBO. However, to fully access the signal, one should take into account more realistic extended mass functions which will make PBHs evaporate at different times, leading in this way to a longer period of PBH-domination, a fact which can potentially enhance the detectability of the signal discussed here. However, if the transition is gradual due to different PBH evaporation times, this will give an extra suppression as found in Inomata et al. (2019b). One therefore is met with two competitive effects which should be considered together in the case of a calculation with realistic extended mass distributions. Interestingly, as it was shown recently in Papanikolaou (2022a) if one accounts for a cosmologically motivated power-law primordial curvature power spectrum for the generation of an extended PBH mass function, they can find that despite the gradualness of the transition the GW signal can be detectable by LISA and can serve as a novel probe to constrain cosmological parameters of the early Universe. It is worth to mention as well that as pointed out recently in Papanikolaou et al. (2024b), the portal of GWs induced by PBH energy density fluctuations can serve also a novel probe of primordial non-Gaussianities.

5.6 From PBH mergers

5.6.1 Formulation

Overlapped GWs from PBH mergers form a SGWB. The spectral shape of the SGWB, its constraints by GW observations, and its implication to PBH physics have been discussed in the literature (Mandic et al. 2016; Clesse and García-Bellido 2017a; Wang et al. 2018; Raidal et al. 2017; Chen et al. 2019; Wang et al. 2019; De Luca et al. 2020d; Bavera et al. 2022; Bagui and Clesse 2022). For a binary system with a circular orbit, the single source energy spectrum $\frac{dE_{GW}}{d\ln f_r}$ for the inspiral phase is given by

$$
\frac{dE_{GW}}{d \ln f_r} = \frac{\pi^{2/3}}{3} M_c^{5/3} (G f_r)^{2/3}, \qquad (136)
$$

where $f_r = (1+z)f$ is the GW frequency in the source frame with $f$ being the
frequency at the observer. One may use a more accurate fitting formula, which
includes the contributions of the inspiral, merger, and ringdown parts of the BBH
waveform (see (Ajith et al. 2011) for the details). The amplitude of the SGWB is
given by summing up the energy spectrum of each binary system and by taking into
account the merger rate distribution. At the end, the GW spectral abundance is
written as
$$ \Omega_{\text{GW}}(f) = \frac{1}{\rho_c} \int_{0}^{z_{\text{max}}} dz \frac{N(z)}{(1+z)} \frac{dE_{\text{GW}}}{d \ln f_r}, \qquad (137) $$

where $\rho_c$ is the critical density of the Universe, $\rho_c = \frac{3H_0^2}{8\pi G}$. Here $N(z)$ is the number density of GW events at redshift $z$ and can be related to the merger rate in a comoving volume from the previous section $R_{\text{EB/LB}}$

$$ N(z) = \frac{R_{\text{EB/LB}}}{(1+z)H(z)}. \qquad (138) $$

In the case of a broad mass distribution, we have to integrate Eq. (137) over the BH masses, $m_1$ and $m_2$, and we get (Clesse and García-Bellido 2017a)

$$ \Omega_{\text{GW}}(f) = \frac{8\pi^{5/3} G^{5/3}}{9H_0^2} f^{2/3} \int_0^{z_{\text{max}}} dz \frac{1}{H(z)(1+z)^{4/3}} \times \int d \ln m_1 \int d \ln m_2 R_{\text{EB/LB}} M_c^{5/3}(m_1, m_2). \qquad (139) $$

The maximum frequency at the observer is determined by the innermost stable circular orbit which is given by

$$ f_{\text{ISCO}} \approx 4.4 \, \text{kHz} \frac{M_{\odot}}{m_1 + m_2} \frac{1}{1+z}. \qquad (140) $$

This can be translated to the maximum redshift as

$$ z_{\text{max}} = \frac{f}{f_{\text{ISCO}}} - 1. \qquad (141) $$

As described in Sect. 4, there are two binary formation channels, dubbed early and late binaries. Below, we discuss the SGWB for each case.

## 5.6.2 Early binaries

The SGWB from early binaries is obtained by inserting the corresponding merging rate distribution of Eq. (101) into the above expression, Eq. (139). Since the merging rate depends on the redshift $z$, $R_{\text{EB}}(z, m_1, m_2) = R_{\text{EB}}(m_1, m_2) \times (t(z)/t_0)^{-34/37-0.29}$ for $f_{\text{PBH}} = 1$, with $t(z)$ being the elapsed time between the Big Bang and a certain redshift $z$, this redshift-dependent factor can be included in the integral over $z$ (Bagui and Clesse 2022),

$$ \int_{0}^{z_{\text{max}}} \frac{(t(z)/t_0)^{-34/37-0.29}}{(1+z)^{4/3} H(z)} dz = 5.92 H_0^{-1}, \qquad (142) $$

where $z_{\text{max}} \approx 100$, as discussed in Sect. 4.1.3. The integral on the masses can be performed by direct numerical integration or by using a Monte-Carlo method as in Clesse and García-Bellido (2017a), e.g. for a synthetic population of $10^6$ BHs.
Expressing the merging rate in yr⁻¹Gpc⁻³, the chirp mass in solar masses M⊙ and the frequency in Hz, one obtains

$$ \Omega_{\text{GW}}(f)h^2 \approx 6.46 \times 10^{-14} \int d\ln m_1 \int d\ln m_2 \left( \frac{R_{\text{EB}}(m_1, m_2)}{\text{yr}^{-1}\text{Gpc}^{-3}} \right) \left( \frac{f}{\text{HZ}} \right)^{2/3} \left( \frac{M_{\text{c}}}{M_{\odot}} \right)^{5/3}, \quad (143) $$

where the constant h takes the value h = 0.67. The resulting SGWB coming from early PBH binaries following a log-normal mass function (with mean mass μ = 2.5M⊙ and width σ = 1) as well as a mass distribution with thermal features are displayed in Fig. 30, as a function of the frequency. In the log-normal case, one notices that the spectrum reaches its highest point at 100 Hz and is suppressed at lower frequencies, whereas in the broad mass distribution case including QCD-induced features, it peaks at 10⁻⁶, 10⁻³ and around 100 Hz. In the latter case, for ground-based detectors, one observes that the SGWB is above the projections for the LVK and ET detector designs, but just below the current limits imposed by the third observing run of LVK. The next observing runs of LVK may therefore be decisive for the detection of this SGWB. In addition, the detected signal by NANOGrav located within the interval 10⁻⁸ – 10⁻⁷ Hz with ΩGWh² ~ 10⁻⁹ (gray window, see Sect. 7.6.4 for more details) could be interpreted as a background of GWs in the model for late binaries in clusters (see next Subsection). However, in the early binaries model, the amplitude of the SGWB is always lower than the signal detected by NANOGrav, and therefore this model can not explain the observed signal. The detectability with LISA will be discussed in Sect. 8. The dependence of the SGWB on the PBH masses for

<figure><img src="image_33.png" /><figcaption>Fig. 30 The SGWB spectrum ΩGWh² for early PBH binaries with a log-normal mass function (in red, with central mass μ = 2.5M⊙ and width σ = 1) and a broad mass distribution with scalar spectral index n_s = 0.970 (in blue), and fPBH = 1. The numerical spectrum also shows the sensitivities of the ground-based interferometers: the LVK O3 Run, the final LVK and the Einstein Telescope (ET). The sensitivity of future space-based interferometers is also shown (LISA, BBO/DECIGO). The Pulsar Timing Array (PTA) considered here is the Square Kilometer Array (SKA) (Schmitz 2021). The NANOGrav 12.5 signal is represented by the gray square (see Sect. 7.6.4 for more details) Image reproduced with permission</figcaption></figure>
the case of a broad mass function from the thermal history of the Universe is shown in Fig. 31.

### 5.6.3 Late binaries in clusters

In the case of late binary formation channel discussed in Sect. 4.2, the merger rate is given by Eq. (109). Substituting it into Eq. (139), the energy spectrum for the low-frequency inspiral regime can be estimated as

$$ \Omega_{\text{GW}}(f)h^2 \approx 2.41 \times 10^{-14} R_{\text{clust}} \left(\frac{f}{\text{HZ}}\right)^{2/3} \times \int dm_1 dm_2 \frac{f(m_1)f(m_2)(m_1+m_2)^{23/21}}{(m_1 m_2)^{5/7}}, \quad (144) $$

with the masses in solar masses and we use $h = 0.67$. In Fig. 32, we show the SGWB spectrum generated by late PBH binaries in clusters (we have chosen the virial velocity in those clusters to be $v_{\text{vir}} = 5$ km/s, compatible with velocity dispersions in faint dwarf galaxies) for the log-normal mass function (with the central mass $\mu = 2.5M_{\odot}$ and the variance $\sigma = 1$) and the thermal history mass distribution (with $n_s = 0.97$ and no running). Here, the normalization is taken by requiring that the integral of the merger rate at $z = 0$ over the mass range [5, 100]$M_{\odot}$ gives a total rate of $38 \text{ yr}^{-1}\text{Gpc}^{-3}$. Note that the SGWB spectrum at lower frequency is highly uncertain, as it can be suppressed by many factors such as the tilt and running of the primordial curvature perturbation spectrum, and by the eventual presence of a cutoff in the merger rate due to the fact that massive PBHs become isolated (Braglia et al.

<figure><img src="image_34.png" /><figcaption>Fig. 31 Contribution to different PBH masses to the SGWB from early PBH binaries (with $f_{\text{PBH}} = 1$). The color bar indicates the values of the quantity $\log_{10}(\Omega_{\text{GW}}h^2)$, which is represented as a function of the logarithm of the masses $m_1$ and $m_2$ of the PBHs following a mass distribution from the thermal history of the Universe (with $n_s = 0.97$) for LISA frequencies ($10^{-3}$ Hz). Image reproduced with permission from Bagui and Clesse (2022), copyright by Elsevier Image reproduced with permission</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>**Fig. 32** The SGWB spectrum $\Omega_{\text{GW}}h^2$ for late PBH binaries in clusters with $f_{\text{PBH}} = 1$ for a log-normal mass function (in red, with the central mass $\mu = 2.5M_{\odot}$ and the width $\sigma = 1$) and a broad mass distribution (in blue, with $n_s = 0.97$ and no running). Note that here we use Ajith et al. (2011) for the single source GW energy spectrum instead of Eq. (139) Image reproduced with permission</figcaption></figure>

2021b). We can see that, in the case of the log-normal mass function, the spectrum peaks around theISCO frequency of PBH binaries (for $m_{\text{PBH}} = 2.5M_{\odot}$, we get $f_{\text{ISCO}} \sim 900$ Hz), while in the case of the thermal history mass function, we see that

<figure><img src="image_36.png" /><figcaption>**Fig. 33** The mass contribution to the SGWB from late PBH binaries in clusters (with $f_{\text{PBH}} = 1$). The color bar indicates the values of the quantity $\log_{10}(\Omega_{\text{GW}}h^2)$, which is represented as a function of the logarithm of the masses $m_1$ and $m_2$ of the PBHs following a mass distribution from the thermal history of the</figcaption></figure>
the SGWB is produced over a wide range of frequencies because of the wide mass distribution. In Fig. 33, we show which PBH masses contribute to the GW amplitude $\Omega_{GW}$ at the LISA frequency $10^{-3}$ Hz. We can see that $2M_{\odot}$ mass PBHs paring with more massive ones are giving large contribution to the SGWB.

## 5.7 SGWB duty cycle

The SGWB of binary origin may have a characteristic property which could be used to infer the origin of BHs. If the interval between BBH events is larger than the typical duration of the signal, the waveforms do not overlap, and the SGWB becomes strongly non-Gaussian (sometimes referred to as intermittent or popcorn signal).

One of the quantities to characterize such a non-Gaussian characteristics of the SGWB is the so-called astrophysical duty cycle. The duty cycle gives the average number of events present in the frequency bin and can be defined using the event rate $dR/dz$ and the duration of the signal staying in the frequency bin $d\bar{\tau}/df$ as

$$ \frac{dD}{df} = \int dz \frac{dR}{dz} \frac{d\bar{\tau}}{df}. \qquad (145) $$

Here, $d\bar{\tau}/df$ is determined by the chirp mass $M_c^z = (1+z)M_c$ as

$$ \frac{d\bar{\tau}}{df} = \frac{5}{96\pi^{8/3}} (G M_c^z)^{-5/3} f^{-11/3}, \qquad (146) $$

where $M_c^z = (1+z)M_c$. When the duty cycle is larger than unity, the GW events overlap, and the background is in the so-called continuous regime. On the other hand, if the interval between events is comparable or larger than the typical duration of the signal, the duty cycle becomes smaller than unity, and the statistical properties are strongly non-Gaussian. The value of the duty cycle can differ depending on the mass function and the redshift distribution. Thus, the non-Gaussian property could be used to distinguish primordial and astrophysical scenarios (Mukherjee and Silk 2020; Braglia et al. 2023a) together with the spectral shape (Mukherjee and Silk 2021).

## 5.8 From close encounters

Overlapped GWs from PBH close encounters can also form a SGWB (García-Bellido et al. 2022). The formulation is the same as the PBH mergers, Eq. (139). In the case of close hyperbolic encounters (CHE), the energy emitted per logarithmic frequency bin is given by (De Vittori et al. 2012; García-Bellido and Nesseris 2017)

$$ \frac{dE_{GW}^{CHE}}{d \ln f_r} = v \frac{dE_{GW}}{dv} = \frac{4\pi G^3 m_1^2 m_2^2}{45 a^2 c^5 v_0} v^5 F_E(v), \qquad (147) $$

where we have defined $v \equiv 2\pi v_0 f_r$ and $v_0^2 \equiv a^3/GM$, and the semi-major axis $a$ is related to the initial velocity as $a = GM/v_0^2$. The function $F_E(v)$ describes the dependence on eccentricity $e$ and it is given by (García-Bellido and Nesseris 2018)
$$
\begin{align}
v^5 F_E(v) &\approx \frac{12F(v)}{\pi y (1+y^2)^2} e^{-2v \xi(y)}, \\
F(v) &= v^2 (1-y^2 - 3vy^3 + 4y^4 + 9vy^5 + 6v^2y^6), \tag{148} \\
\xi(y) &= y - \tan^{-1}y, \\
y &= \sqrt{e^2 - 1}. \nonumber
\end{align}
$$

Substituting Eqs. (115) and (148) into Eq. (139), we can estimate the peak frequency

$$
f_{\text{peak}} \approx 4.3 \, \text{HZ} \left(\frac{y}{0.01}\right)^{-3} \left(\frac{M}{2000 M_{\odot}}\right)^{1/2} \left(\frac{a}{1 \, \text{AU}}\right)^{-3/2}, \quad (149)
$$

and the peak amplitude

$$
\begin{equation}
\begin{split}
\Omega_{\text{GW}}(f_{\text{peak}}) \approx {}& 0.9 \times 10^{-13} h \left(\frac{\Omega_{\text{m}}}{0.3}\right)^{-1/2} \left(\frac{\Omega_{\text{DM}}}{0.25}\right)^2 \left(\frac{\delta_{\text{loc}}}{10^8}\right) \left(\frac{a}{0.1\,\text{AU}}\right)^{-2} \\
& \times \left(\frac{y}{0.01}\right)^{-5} \frac{m_1}{100M_{\odot}} \frac{m_2}{100M_{\odot}} \frac{m_1+m_2}{200M_{\odot}}.
\end{split}
\tag{150}
\end{equation}
$$

In Fig. 34, we show the SGWB spectrum originating from close hyperbolic encounters, compared with the one from binary PBHs. We find that, for the modest

<figure><img src="image_37.png" /><figcaption>**Fig. 34** Comparison of the SGWB spectrum originating from BBHs and CHEs, both for $\beta = 0$ (solid) and $\beta = 1.28$ (dashed), where $\beta$ is a parameter characterizing the redshift dependence of the merger rate as $\tau^{\text{BBH}} \propto (1+z)^\beta$. For the BBH curves, we take $m_1 = m_2 = 100-300 M_\odot$ and $v_0 = 30\text{km/s}$. The CHE curves correspond to the same range of masses with $a_0 = 5\text{AU}$, $y_0 = 2 \times 10^{-3}$ for frequencies around 10 Hz, and $a_0 = 5 \cdot 10^7 \text{AU}$, $y_0 = 10^{-5}$ in the mHz range. For all cases, we assume log-normal distributions for $a_0, y_0$ and the PBH mass, of respective widths $\sigma_a, \sigma_y = 0.1, 0.5$, as well as $f_{\text{PBH}} = 1$. The bands come from the possible range of parameters $a_0$ and $y_0$. For a smaller fraction of PBHs, the GW spectral amplitude simply scales as $\Omega_{\text{GW}} \propto f_{\text{PBH}}^2$. Image and estimated LISA sensitivity from García-Bellido et al. (2022), copyright by the author(s) Image reproduced with permission</figcaption></figure>
choice of the parameter values, it is difficult to be reached by LISA. This difficulty can easily be understood if we rewrite Eq. (150) in terms of the peak frequency,

$$ \Omega_{\text{GW}}^{\text{CHE}}(f_{\text{peak}}) \approx 3.1 \times 10^{-13} h \times \left(\frac{\Omega_{\text{m}}}{0.3}\right)^{-1/2} \left(\frac{\Omega_{\text{DM}}}{0.25}\right)^2 \left(\frac{\delta_{\text{loc}}}{10^8}\right) \left(\frac{f_{\text{peak}}}{50 \text{ Hz}}\right)^{4/3} \times \left(\frac{y}{0.01}\right)^{-1} \frac{m_1}{100 M_{\odot}} \frac{m_2}{100 M_{\odot}} \left(\frac{m_1 + m_2}{200 M_{\odot}}\right)^{1/3}. \quad (151) $$

We can see that the maximum amplitude $\Omega_{\text{GW}}^{\text{CHE}}(f_{\text{peak}})$ grows with $f_{\text{peak}}^{4/3}$. Thus, the peak amplitude of the background decays significantly at low frequencies. However, we have observed that the peak amplitude tends to get enhanced when we make the log-normal distribution of the mass function wider. Thus, by considering a different PBH mass function, we may be able to find a case where the SGWB from CHEs gives an interesting contribution at the LISA frequency.

It is worth mentioning that it could be tested by third-generation ground-based GW detectors such as Einstein Telescope and Cosmic Explorer. Although the CHE curves are below the BBH curves, loud BBH events can be detected individually and subtracted from the data, which may allow us to probe the other SGWB component below the green curve.

## 5.9 Summary

In this section we have explored the different sources of SGWBs related to PBHs:

* The SGWB induced at second order by the large curvature fluctuations leading to PBH formation. Interestingly, it peaks at PTA frequencies for stellar-mass PBHs, which makes it a very complementary probe to test a possible primordial origin of GW events seen in ground-based detectors. This SGWB could be connected to the recent pulsar timing array observations, see Sect. 7.6.4 for more details. Future GW detectors like LISA, Cosmic Explorer and Einstein Telescope will probe small PBHs, in particular in the asteroid-mass range that remains so far unconstrained (see the section on the limits on the abundance of PBHs).

* The SGWB from the Poisson fluctuations in the initial PBH distribution, even if they were so small that they have totally evaporated. This mechanism could allow to probe the possible PBH formation at very high energy scales, up to the reheating phase and the GUT energy.

* The stochastic background from early binaries, in particular if PBHs have an extended mass function. In such a case the SGWB is boosted by the merging of binaries with extreme mass ratios.

* The stochastic background from late binaries in clusters, which for stellar-mass PBHs is smaller than the SGWB from early binaries at frequencies corresponding to ground-based GW detectors. Because the merger rates are not suppressed with
the PBH mass, it would however be the dominant background at low frequencies, i.e. in the LISA and PTA frequency range.

* The stochastic background from close hyperbolic PBH encounters in clusters, probably undetectable with LISA but in the range of future ground-based GW detectors.

The SGWB is therefore expected to be an interesting complementary signal that will help to distinguish between PBH formation models, mass functions, clustering histories and binary origins.

# 6 Gravitational-waves and large scale structure correlations

## 6.1 Introduction

Primordial black holes, both in the case where they comprise a large part of the dark matter or if they exist only in small numbers, are part of the Large-Scale Structure (LSS) of the Universe, and their observational relation with the other matter-energy components can be generally different with respect to black holes of stellar origin. By the time LISA will be online, there will be data from several large-scale structure experiments, both preceding and overlapping in time with LISA. These include observations from the ground and from space, such as DESI (Aghamousa et al. 2016a, b), the Vera Rubin Observatory (Blum et al. 2022), Euclid (Amendola et al. 2018), SPHEREx (Doré et al. 2014), Nancy Roman (Wenzl et al. 2022), ATLAS (Wang et al. 2019) and the SKAO (Weltman et al. 2020), all together providing galaxy catalogs over a very wide range of scales both in width and deepness, at unprecedented accuracy. It will therefore be paramount to take advantage of such data for PBH science. In this Section we discuss a few possible ways to use gravitational wave observations in combination with the LSS to constrain the existence and abundance of PBHs. LSS observations can be cross-correlated with different types of GW maps, both from the resolved mergers of compact objects and the stochastic backgrounds.

Cross-correlations between GWs from resolved compact object mergers and galaxies have been first investigated in Laguna et al. (2010). Later, it was suggested that the cross-correlation of mergers with EM counterparts could be used to constrain dark energy and modified gravity models (Camera and Nishizawa 2013). After the first detections of gravitational waves, Oguri (2016) investigated the possibility to constrain the distance-redshift relation using such correlations.

It was then shown that the cross-correlation of resolved mergers with LSS catalogs can be used to test primordial black hole scenarios (Raccanelli et al. 2016a; Raccanelli 2017; Scelfo et al. 2018, 2022) and different astrophysical models (Scelfo et al. 2020). More recently, an additional technique for constraining PBH abundances, formation channels and redshift distributions using only GW data has been proposed and studied (Libanore et al. 2022).

The GWxLSS probe became very popular and several different approaches have been suggested and investigated (see e.g., Raccanelli 2017; Scelfo et al. 2018; Cañas
Herrera et al. 2020; Calore et al. 2020), looking at both correlations with resolved sources and the SGWB (Cusin et al. 2018a; Jenkins et al. 2018, 2019a; Bertacca et al. 2020; Yang et al. 2020; Mukherjee and Silk 2020; Alonso et al. 2020). Here we present the main ideas and formalism developed so far; specific predictions for LISA and investigations involving similar analyses for extreme mass-ratio inspirals (EMRIs) and intermediate mass black holes (IMBHs) are currently being developed.

## 6.2 GW-galaxy cross-correlations for resolved events

The cross-correlation of LSS catalogs with maps of gravitational waves coming from the merger of compact objects will provide information on several astrophysical and cosmological parameters and models, and has been investigated in a variety of works from exploring ways to indirectly detect GWs to proper data correlations.

The idea of cross-correlations with the LSS in order to constrain the existence of PBHs and their abundances was first explored in Raccanelli et al. (2016a) and then further investigated in e.g., Raccanelli (2017), Scelfo et al. (2018, 2020). The possibility of using the GWxLSS correlation to test the existence of PBHs comes from the fact that primordial and astrophysical binary black holes trace the LSS in different ways. This reflects into the fact that the bias of the hosts of the binary system will be different when the BHs have different progenitors and formation channels. This bias quantifies the relation between underlying matter fluctuations and in the number of observed sources.

In Bird et al. (2016), it was shown that primordial black holes preferentially merge in lower biased objects and thus have a lower cross-correlation with luminous galaxies. This happens because binary PBHs would preferentially merge in halos with low velocity dispersion, which are low-mass halos and hence with very little or absent star formation; these objects have a bias $b<1$. On the other hand, mergers of compact objects that are the endpoint of stellar evolution naturally happen for the vast majority within star forming-rich halos, which have larger galaxy bias values, $b > 1$. In the linear biasing scheme, $b$ quantifies the relation between the underlying matter distribution and observed sources, $\delta_0 = b\delta_m$. This reasoning is valid for PBH binary systems that form by capture (late binary formation scenario), as these systems form in the local Universe.

However, there is also another way to have mergers of PBH binaries, which we have referred to as *early binary* mergers. In this case, within the standard formation scenario (i.e., the collapse of primordial curvature perturbations at horizon re-entry), PBHs should form in the peaks of the initial matter distribution, and therefore they will trace the dark matter field. As a consequence, such mergers should be unbiased tracers of the underlying density field.

The overall catalog of PBH mergers will consequently be composed of a mixture of early and late binaries, with the total bias being the weighted average of the two.

Therefore, by measuring the amplitude of the angular cross-correlation of galaxy maps with catalogues of compact-object mergers, which directly depends on the bias of the mergers’ hosts, one can statistically probe the abundance of PBH mergers. This information can also be used to discriminate between different astrophysical models (see e.g., Scelfo et al. 2020).
In order to measure the mergers’ hosts bias, one can use the 3D angular power spectrum:

$$
C_{\ell}^{\text{GW}\times\text{LSS}}(z_1, z_2) = \frac{2}{\pi} \int dk \, k^2 P(k) \Psi_{\ell}^{\text{GW}}(k, z_1) \Psi_{\ell}^{\text{LSS}}(k, z_2), \quad (152)
$$

where $P(k)$ is the matter power spectrum and the kernels $\Psi_{\ell}^{X}(z)$ include the obser-
vational window functions and the relevant physical effects:

$$
\Psi_{\ell}^{X}(k, z_i) = \int_{z_j - \Delta z}^{z_i + \Delta z} N^{X}(z) b^{X}(z) D(z) W(z) \Delta_{\ell}^{X}(k, z) \, dz; \quad (153)
$$

here X will be either GW or LSS, and W are observational window functions related to the experiment specifications (i.e., they include information on the redshift distribution and survey geometry and sensitivity). The Δℓ terms are the (gauge-independent) observed overdensities and include effects from intrinsic clustering, peculiar velocity and Doppler contributions, lensing and gravitational potentials containing in this way information information on the underlying cosmological and astrophysical models. They are usually expressed as:

$$
\Delta_{\ell}(k) = \Delta_{\ell}^{\text{den}}(k, z) + \Delta_{\ell}^{\text{rsd}}(k, z) + \Delta_{\ell}^{\text{dop}}(k, z) + \Delta_{\ell}^{\text{len}}(k, z) + \Delta_{\ell}^{\text{gr}}(k, z). \quad (154)
$$

For details on these terms, their derivation and physical meaning, see e.g., Bonvin
and Durrer (2011), Challinor and Lewis (2011), Raccanelli et al. (2016b), Scelfo
et al. (2018).

6.3 GW-galaxy cross-correlations for the stochastic background

In addition to the signal from mergers detected by interferometers, there is another
signal coming from the superposition of many unresolved sources. These will form a
stochastic background of gravitational waves, as discussed in Sect. 5. Obviously,
primordial black holes, if they exist, would contribute to this signal, and their
associated GW backgrounds have been discussed in Sect. 5. This will be particularly
relevant for LISA, as the SGWB is expected to be a product of LISA’s observations.

In the same fashion as for the catalogs of resolved events described above, we can
cross-correlate the stochastic background with the LSS; such correlation will contain
information on the redshift distribution, mass function and clustering behavior of the
sources. Moreover, the observed signal will be the combination of the stochastic GW
background from astrophysical objects, denoted AGWB, the one from primordial
black holes, and the one generated from inflation; disentangling the signal and
distinguishing between the different sources will be of paramount importance to
obtain information about black hole and early Universe physics, and one of the most
promising ways to do so is thanks to the SGWB–LSS correlation.

Given that the GW energy density depends not only on astrophysical properties
but also on cosmological perturbations, it will correlate with other cosmological
probes; moreover, GWs will experience projection effects that will need to be
accounted for in order to observe the signal in the appropriate frame (Bertacca et al.
2020; Bellomo et al. 2022). Some preliminary works on this have been done in the last few years, including forecasts for the cross-correlation signals between GW observatories and future galaxy surveys, as e.g., Euclid and SKA (see, e.g., Contaldi 2017; Cusin et al. 2017, 2018b, a; Jenkins et al. 2019a, 2018; Cusin et al. 2019b, 2020; Jenkins and Sakellariadou 2019; Jenkins et al. 2019b; Bertacca et al. 2020; Pitrou et al. 2020; Mukherjee and Silk 2020; Alonso et al. 2020).

In a similar way to the case of resolved mergers, the most natural observable is the angular cross power spectrum that correlates the energy density of the AGWB with galaxy number counts:

$$C_{\ell}^{\text{AGWB} \times \text{LSS}} = \frac{2}{\pi} \int dk k^2 P(k) \Psi_{\ell}^{\text{AGWB}}(k) \Psi_{\ell}^{\text{LSS}}(k), \qquad (155)$$

where most quantities are the same of Eq. (152), and the dependence on astrophysical models and parameters can be included in the expression for the AGWB energy density. This is usually written as the total GW energy density (denoted by the superscript TOT) per logarithmic frequency $f_0$ and solid angle $\Omega_0$ along the line-of-sight $\hat{\mathbf{n}}$ of a SGWB, defined in e.g., Cusin et al. (2017), Bertacca et al. (2020) as

$$\Omega_{\text{GW}}^{\text{TOT}}(f_0, \hat{\mathbf{n}}) = \frac{f_0}{\rho_c} \frac{d\rho_{\text{GW}}^{\text{TOT}}}{d\Omega_0}; \qquad (156)$$

this contains both a background, which is homogeneous and isotropic, and a direction-dependent contribution; the total relative fluctuation can thus be defined as (similarly to Eq. (132)):

$$\Delta_{\text{AGWB}}^{\text{TOT}}(f_0, \hat{\mathbf{n}}) = \frac{\Omega_{\text{AGWB}}^{\text{TOT}} - \Omega_{\text{AGWB}}^{\text{TOT}}}{\Omega_{\text{AGWB}}^{\text{TOT}}}. \qquad (157)$$

Studies of the cross-correlation between the AGWB with galaxy number counts are presented in Cusin et al. (2018a), Mukherjee and Silk (2020), Cañas Herrera et al. (2020), Alonso et al. (2020), Cusin et al. (2019b); Yang et al. (2020), Bellomo et al. (2022). This allows a tomographic reconstruction of the redshift distribution of sources, and can be also useful for the shot noise characterization (in the case of ground-based detectors) (Alonso et al. 2020). The cross-correlation shows that the combination of galaxy surveys with the AGWB can be a powerful probe for GW physics and can be a robust observational probe for multi-messenger cosmology.

## 6.4 GW × LSS forecasts

Constraints on the presence of PBHs from measurements of the GWxLSS have been forecasted for the case of resolved GW signals, correlated with both galaxy surveys (Raccanelli et al. 2016a; Raccanelli 2017; Scelfo et al. 2018; Bosi et al. 2023) and intensity mapping experiments (Scelfo et al. 2022). Constraints coming from measurements of the stochastic background signal are currently being worked out. In this Section we introduce the formalism for obtaining such forecasts and report some of the results present in literature.
The correlation of two tracers (X, Y) can be estimated from tomographic maps of their number counts and evaluating their N-point statistics; in the case we consider here, this is the angular power spectrum $C_{\ell}^{XY}(z_i, z_j)$ of the tracer X in redshift bin $z_i$ and the tracer Y in $z_j$. Assuming that the noise comes only from the standard shot noise contribution, the total signal extracted from the maps is:

$$ \langle a_{\ell m}^{X}(z_i) a_{\ell' m'}^{Y*}(z_j) \rangle = \delta_{\ell \ell'} \delta_{mm'} \tilde{C}_{\ell}^{XY}(z_i, z_j) = \delta_{\ell \ell'} \delta_{mm'} \left[ C_{\ell}^{XY}(z_i, z_j) + \delta_{XY} \delta_{ij} \mathcal{N}_{\ell}^{X}(z_i) \right], $$

where $\mathcal{N}_{\ell}$ is the noise angular power spectrum.

The observed signal will receive contributions from several physical effects, coming from the cosmological perturbations described above, from density, velocity, lensing and gravity effects. It is important to include all of them to avoid misestimating parameters, best fit values and errors when doing parameter inference, as shown in Bellomo et al. (2022), Bernal et al. (2020).

Forecasts obtained so far have been based on the Multi_CLASS (Bellomo et al. 2020) and GW_CLASS codes (Bellomo et al. 2022). Multi_CLASS is the first public Boltzmann code that allows to compute the angular power spectrum for multiple galaxy (and other tracers) populations. The code lets the user specify specific properties of different tracers, such as its number density redshift distribution and bias, magnification bias and evolution bias parameters, while GW_CLASS calculates the stochastic signal from astrophysical mergers and accounts for all the astrophysical and projection effects contributions.

The first suggestion that this observable can be used to constrain PBHs was presented in Raccanelli et al. (2016a), where after introducing the concept, it is shown that the cross-correlation between galaxies mapped by future radio surveys and GW interferometers could distinguish between primordial and stellar origin of merging black holes. Such possibility depends on the merger rate of primordial black holes and the angular resolution of the GW interferometer. In that work, to obtain the results, a minimum-variance estimator was introduced for the effective correlation amplitude, $A_c \equiv r \times b_{GW}$, where r is the cross-correlation coefficient of Eq. (152). This cross-correlation coefficient parametrises the extent to which two biased tracers of the matter field are correlated (Tegmark and Peebles 1998).

The minimum-variance estimator for the effective correlation amplitude is given by (see e.g. Jeong and Kamionkowski 2012):

$$ \widehat{A_c} = \frac{\Sigma_{\ell} \tilde{C}_{\ell} F_{\ell} / \mathrm{Var}[\tilde{C}_{\ell}]}{\Sigma_{\ell} F_{\ell}^2 / \mathrm{Var}[\tilde{C}_{\ell}]}, \qquad (158) $$

where $\tilde{C}_{\ell}$ is the measured power spectrum and $F_{\ell} \equiv d\tilde{C}_{\ell}/d\widehat{A_c} \propto b_g$. The variance of this estimator is then:

$$ \sigma_{\widehat{A_c}}^2 = \left[ \sum_{\ell} \frac{F_{\ell}^2}{\mathrm{Var}[\tilde{C}_{\ell}]} \right]^{-1}, \qquad (159) $$

which can be used to forecast the measurement error when neglecting that of other parameters.
More generally, the measurement error for specific parameters in a given experiment can be estimated using Fisher analysis. For this case, we write the Fisher matrix as:

$$F_{\alpha\beta} = \sum_{\ell} \frac{\partial C_{\ell}}{\partial \vartheta_{\alpha} \partial \vartheta_{\beta}} \sigma_{C_{\ell}}^{-2}, \qquad (160)$$

where $\vartheta_{\alpha} = \{A_c, b_g\}$; the derivatives of the power spectra $C_{\ell}$ are evaluated at fiducial values $\bar{\vartheta}_{\alpha}$ corresponding to the scenario at hand, and $\sigma_{C_{\ell}}$ are errors in the power spectra.

The results in Raccanelli et al. (2016a) were obtained by computing the $2 \times 2$ Fisher matrix for the parameters $\{A_c, b_g\}$, using a prior on the galaxy bias corresponding to the precision reached by fitting the amplitude of the galaxy auto-correlation function $C_{\ell}^{gg}$, using a value of $\ell_{max} = 200$, which yields a $\sim 10\%$ precision in the measurement of the bias $b_g$. For more details on this analysis, see Raccanelli et al. (2016a).

Subsequently, in Raccanelli (2017), such analysis was updated to include improved merger rate estimates from the LVK collaboration, for the correlation with different possible surveys with the SKA. In Fig. 35 we show forecasts for different experiments, considering different versions of LVK and ET correlated with HI and continuum SKA surveys. The bars show the constraints on $A_c$ for different values of a parameter $R$ that encapsulates uncertainties on the merger rate, which is degenerate with observation times from the gravitational wave interferometers.

Following this, in Scelfo et al. (2018), Libanore et al. (2023) these forecasts are updated by including all projection effects, and inserted in an analysis of a null hypothesis test by comparing two models, one in which the BHs origin is stellar, the other in which it is primordial. One model is assumed as a fiducial model. Then, by computing the signal-to-noise ratio (SNR), it is checked whether the alternative

<figure><img src="image_38.png" /><figcaption>Fig. 35 Forecast errors on the cross-correlation amplitude, $A_c$, for different experiment combinations, varying merger rates and years of observations. Each column corresponds to a GW detector experiment, for merger rates from 1 to 10 Gpc$^{-3}$yr$^{-1}$. Horizontal lines show the expected difference in the cross-correlation between (late binary) PBH and stellar binaries, for different values of $f_{PBH}$. Image reproduced with permission</figcaption></figure>
model can be differentiated from the fiducial one. The null hypothesis is that the model is indistinguishable from the fiducial, which happens for low values of the signal-to-noise ratio.

For instance, the distance of an alternative model from the fiducial can be quantified using a $\Delta\chi^2$ statistics. The $\Delta\chi^2$ is given by the logarithm of a likelihood, and the SNR can be written as:

$$ \left(\frac{S}{N}\right) = f_{\text{sky}} \sum_{2}^{\ell_{\text{max}}} (2\ell + 1) (\mathbf{C}_{\ell}^{\text{Alternative}} - \mathbf{C}_{\ell}^{\text{Fiducial}})^{T} \text{Cov}_{\ell}^{-1} (\mathbf{C}_{\ell}^{\text{Alternative}} - \mathbf{C}_{\ell}^{\text{Fiducial}}), $$

where $\mathbf{C}_{\ell}^{T} = (C_{\ell}^{gg}(z_1, z_1), \dots, C_{\ell}^{gGW}(z_1, z_1), \dots, C_{\ell}^{GWGW}(z_1, z_1), \dots)$, $f_{\text{sky}}$ is the sky fraction covered by both GW and galaxy surveys and $\text{Cov}_{\ell}$ is the covariance matrix, computed from the angular power spectra of the fiducial model.

In Fig. 36, we show the results where the fiducial model is taken to that of stellar BH mergers (see Scelfo et al. 2018). The bars span values between 0.1 and 10 for the factor $r$, (similarly to the case of Fig. 35, this parametrises the uncertainties on the BBH merger rates, which is degenerate with the observation time of the experiment).

After updates on the BBH merger rate and the specifications for the DESI survey (Aghamousa et al. 2016a), forecasts for the SNR for detecting a component of PBHs are shown in Bosi et al. (2023) and a summary plot can be found in Fig. 37. Here the color code is for the SNR, while on the x-axis is the fraction of PBH mergers from the total of observed BBH mergers observed by the Einstein Telescope, and on the y-axis the parameter $r$, as defined above.

<figure><img src="image_39.png" /><figcaption>Fig. 36 Expected signal-to-noise ratio $(S/N)_{\Delta B/B}$ from the Fisher matrix analysis, for cross-correlations between GW observations (from either LVK or Einstein Telescope) and high-redshift LSS surveys (EMU, DESI and SKA), assuming stellar black holes for the fiducial scenario. In each case, the left bar corresponds to late PBH binaries and the right bar to early PBH binaries. Image reproduced with permission from Scelfo et al. (2018), copyright by IOP and SISSA Image reproduced with permission</figcaption></figure>
<figure><img src="image_40.png" /><figcaption>**Fig. 37** SNR for detecting a fraction $\Gamma$ of PBH mergers out of the total BBH mergers observed by the ET, as a function of the parameter $r$, cross-correlated with DESI. The left panel shows the result using just DESI x ET, while the right panel includes Planck priors. Results from Bosi et al. (2023) Image reproduced with permission</figcaption></figure>

These results show that the presence of primordial black holes can be detected by using the cross-correlation with LSS and the fact that these objects cluster differently than astrophysical black hole mergers. Future GW detectors and LSS surveys will therefore allow a detection of a primordial component in the catalog of binary black hole mergers.

The cross-correlation forecasts shown so far are for the correlation of GWs with resolved galaxies; naturally, a similar analysis can be done for observations of intensity mapping. A first investigation for the 21 cm Intensity Mapping (IM) observed with the SKAO is presented in Scelfo et al. (2022), focusing on the IM of the neutral hydrogen (HI) from the proposed 21 cm IM survey with the SKAO and on resolved GW events from the merger of BBHs as detected by the Einstein Telescope (ET) (Sathyaprakash et al. 2012). A cross-correlation signal is expected, as both HI and GWs trace the cosmic density field, and once again they do so in different ways depending on astrophysical and cosmological models. In particular, for constraining the PBH abundance, results were obtained for the predicted IM from the SKAO correlated with GWs from ET.

## 6.5 Summary

Gravitational-wave observations of resolved black hole mergers and of the stochastic background can be cross-correlated with the LSS to constrain the existence and abundance of PBHs. This technique may provide additional information to test and constrain the different PBH scenarios, the PBH abundance within a broad range of masses, as well as binary formation channels. This will be of particular interest for the next generation of LSS surveys, such as DESI, Euclid and SKA, and gravitational-wave detectors like LISA, Einstein Telescope or Cosmic Explorer.
# 7 Current limits

The most stringent limits on $f_{\text{PBH}}$ as a function of the PBH mass for a monochromatic mass function are represented in Fig. 38. Zooms on evaporation-based limits (between $10^{-19}M_{\odot}$ and $10^{-15}M_{\odot}$), on microlensing limits (between $10^{-11}M_{\odot}$ and $1M_{\odot}$) and on the multiple limits above $1M_{\odot}$ are shown in Fig. 39. These astrophysical and cosmological probes of PBHs can be divided in six categories, each related to the physical process at the origin of a PBH signature in observations:

1. Observations impacted by *black hole evaporation*

2. *Microlensing* observations

3. Observational signatures of the black hole *dynamics*

4. Observations impacted by the black hole *accretion*

5. Indirect constraints from *density fluctuations*

<figure><img src="image_41.png" /></figure>

**Fig. 38** We show the most stringent claimed constraints in the mass range of phenomenological interest. They come from the Hawking evaporation producing extra-galactic gamma-ray (EG γ) (Arbey et al. 2020), $e^{\pm}$ observations by Voyager 1 (V $e^{\pm}$) (Boudaud and Cirelli 2019), positron annihilations in the Galactic Center (GC $e^{+}$) (DeRocco and Graham 2019) and gamma-ray observations by INTEGRAL (INT) (Laha et al. 2020) (for other constraints in the ultra-light mass range see also Carr et al. 2010; Ballesteros et al. 2020a; Laha 2019; Poulter et al. 2019; Dasgupta et al. 2020; Laha et al. 2021). We plot microlensing searches by Subaru HSC (Niikura et al. 2019b; Smyth et al. 2020), MACHO/EROS (E) (Alcock et al. 2001; Allsman et al. 2001), Ogle (O) (Niikura et al. 2019a) and Icarus (I) (Oguri et al. 2018). Other constraints come from CMB distortions. In black dashed, we show the ones assuming disk accretion (Planck D in Serpico et al. (2020) and Poulin et al. (2017), from left to right) while in black solid the ones assuming spherical accretion (Planck S in Serpico et al. (2020) and both photo- and collisional ionization in Ali-Haïmoud and Kamionkowski (2017), from left to right). Only Serpico et al. (2020) includes the effect of the secondary dark matter halo in catalysing accretion. Additionally, constraints coming from X-rays (Xr) (Manshanden et al. 2019) and X-Ray binaries (XrB)(Inoue and Kusenko 2017) are shown. Dynamical limits coming from the disruption of wide binaries (WB) (Quinn et al. 2009), survival of star clusters in Eridanus II (EII) (Brandt 2016) and Segue I (SI) (Koushiappas and Loeb 2017a; Stegmann et al. 2020) are also shown. LVK stands for the constraint coming from LVK measurements (Ali-Haïmoud et al. 2017; Raidal et al. 2019; Vaskonen and Veermäe 2020; Wong et al. 2021). Constraints from Lyman-α forest observations (Lα) come from Murgia et al. (2019). We neglect the role of accretion which has been shown to affect constraints on masses larger than $O(10)M_{\odot}$ (De Luca et al. 2020c, d) in a redshift dependent manner. See Carr et al. (2021c) for a comprehensive review on constraints on the PBH abundance. Notice that there are no stringent bounds in the asteroid mass range (Katz et al. 2018; Montero-Camacho et al. 2019) where LISA may constrain PBHs through the search of a second order SGWB Image reproduced with permission
<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>

<figure><img src="image_44.png" /></figure>
► **Fig. 39** Zooms over some claimed limits on the PBH abundance $f_{\text{PBH}}$ for a monochromatic distribution of mass $m_{\text{PBH}}$, in the asteroid-mass range where constraints are dominated by various probes of PBH Hawking evaporation (top panel), in the planetary-mass and low stellar-mass range up to $m_{\text{PBH}} \sim 10M_{\odot}$ coming from microlensing surveys (middle panel), and in the range from stellar-mass up to the supermassive PBHs, from a combination of accretion, dynamical, GW and indirect constraints (bottom panel). The legend indicates the origin of each represented limit. It is worth noticing that all those limits are subject to important uncertainties and can be highly model dependent, moving up and down with different model and astrophysical assumptions. The possible limitations and sources of uncertainties are discussed in the text. Image reproduced with permission

6. Signatures of PBHs in gravitational-wave observations

In this section, we review these current limits on the PBH abundance and comment on possible caveats and model dependences. We also discuss the possible observational evidences for the existence for PBHs, in relation with these limits.

**7.1 PBH evaporation**

We review here the PBH constraints relying on their evaporation. Let us note that the PBH mass is expressed in grams through this Section, and not in solar mass as for all the other probes, for a better consistency with the literature and because these constraints concern PBHs much lighter than a solar mass.

**7.1.1 Big-bang nucleosynthesis (BBN)**

PBHs with masses $10^9 - 10^{13}$ g can affect the abundance of light elements produced during the Big-Bang Nucleosynthesis through Hawking radiation and the related emitted particles, allowing to set constraints on the abundance of low-mass PBHs (Carr et al. 2021c). In particular, for PBHs with masses $m_{\text{PBH}} \approx 10^9 - 10^{10}$ g, Hawking radiated mesons and antinucleons induce extra interconversion processes of protons to neutrons, increasing the neutron-to-proton ratio at the time of freeze-out of the weak interaction (Vainer and Naselskii 1978), and consequently leading to an increase in the final $^4$He abundance (Miyama and Sato 1978). Regarding PBHs with masses $m_{\text{PBH}} \approx 10^{10} - 10^{12}$ g, long-lived high energy hadrons produced out of PBH evaporation, such as pions, kaons and nucleons, remain long enough in the ambient medium to trigger dissociation processes of light elements produced during BBN (Kohri and Yokoyama 2000), reducing the abundance of $^4$He and increasing the one of $d$, $^3$He, $^6$Li and $^7$Li. Finally, for PBHs with $m_{\text{PBH}} \approx 10^{12} - 10^{13}$ g, photons produced from the particle cascade process further dissociate $^4$He, increasing the abundance of light synthesized elements (Keith et al. 2020; Acharya and Khatri 2020). However, it is important to stress that the BBN constraints carry out some uncertainties regarding the baryon-to-photon ratio, the reaction and the decay rates of the elements produced during the BBN processes (Kawasaki et al. 2005a, b).
### 7.1.2 Gamma-rays and neutrinos from Hawking radiation

The lifetime of a PBH due to its evaporation in standard model particles can be approximated as (Carr et al. 2010)

$$ \tau \approx 400 \left( \frac{m_{\text{PBH}}}{10^{10} \text{ g}} \right)^3 \text{s.} \qquad (161) $$

For a PBH with mass $5 \times 10^{14}$ g, its lifetime is close to the age of Universe. Thus, such low-mass PBHs can be probed by looking for bursts of gamma rays. Such searches have been performed by the Water Cherenkov Observatories, Imaging Atmospheric Cherenkov Telescopes (IACTs) and gamma-ray space observatories. So far, no gamma-ray signal that could be due to PBH evaporation has been detected. The current limits are summarized in Table 3. These are expressed in terms of the PBH burst rate densities $\dot{B}_{\text{PBH}}$, describing the number of bursts from PBH evaporation in the local Universe in a unit of volume per unit of time.

The observations of Galactic gamma-ray emission by INTEGRAL (Laha et al. 2020) also put a bound on the mass above which PBHs could constitute the entirety of DM, $m_{\text{PBH}} > 1.2 \times 10^{17}$ g, under some relatively conservative assumptions. New analysis of the diffuse soft gamma-ray emission towards the inner Galaxy measured by INTEGRAL/SPI over 16 years (Berteaud et al. 2022), implementing a spatial template fit of SPI data, an improved instrumental background model, and including the spatial distribution of MeV photons emitted by the evaporating PBHs into the fitting procedure, strengthens this lower bound to $m_{\text{PBH}} \gtrsim 4 \times 10^{17}$ g.

<table><caption>**Table 3** Current limits on the local PBH burst rate density from different gamma-ray observatories</caption><thead><tr><th>Experiment</th><th>B&#x0302;<sub>PBH</sub> Upper limit (pc<sup>-3</sup>yr<sup>-1</sup>)</th><th>Burst duration (s)</th><th>T<sub>obs</sub></th><th>References</th></tr></thead><tbody><tr><td>Milagro</td><td>3.1 &times; 10<sup>5</sup> (99% CL)</td><td>0.001</td><td>1673 d</td><td>Abdo et al. (2015)</td></tr><tr><td>Milagro</td><td>1.2 &times; 10<sup>4</sup> (99% CL)</td><td>0.01</td><td>1673 d</td><td>Abdo et al. (2015)</td></tr><tr><td>Milagro</td><td>5.4 &times; 10<sup>4</sup> (99% CL)</td><td>0.1</td><td>1673 d</td><td>Abdo et al. (2015)</td></tr><tr><td>Milagro</td><td>3.6 &times; 10<sup>4</sup> (99% CL)</td><td>1</td><td>1673 d</td><td>Abdo et al. (2015)</td></tr><tr><td>Milagro</td><td>3.8 &times; 10<sup>4</sup> (99% CL)</td><td>10</td><td>1673 d</td><td>Abdo et al. (2015)</td></tr><tr><td>Milagro</td><td>6.9 &times; 10<sup>4</sup> (99% CL)</td><td>100</td><td>1673 d</td><td>Abdo et al. (2015)</td></tr><tr><td>VERITAS</td><td>2.22 &times; 10<sup>4</sup> (99% CL)</td><td>30</td><td>747 h</td><td>Archambault (2018)</td></tr><tr><td>H.E.S.S</td><td>4.9 &times; 10<sup>4</sup> (95% CL)</td><td>1</td><td>2600 h</td><td>Glicenstein et al. (2013)</td></tr><tr><td>H.E.S.S</td><td>1.4 &times; 10<sup>4</sup> (95% CL)</td><td>30</td><td>2600 h</td><td>Glicenstein et al. (2013)</td></tr><tr><td>H.E.S.S</td><td>5.6 &times; 10<sup>4</sup> (99% CL)</td><td>10</td><td>2860 h</td><td>Tavernier et al. (2020)</td></tr><tr><td>H.E.S.S</td><td>2.0 &times; 10<sup>3</sup> (95% CL)</td><td>120</td><td>4816 h</td><td>Aharonian et al. (2023)</td></tr><tr><td>Fermi-LAT</td><td>7.2<sup>+8.1</sup><sub>-2.4</sub> &times; 10<sup>3</sup> (99% CL)</td><td>1.26 &times; 10<sup>8</sup></td><td>1.26 &times; 10<sup>8</sup> s</td><td>Ackermann et al. (2018)</td></tr><tr><td>HAWC</td><td>3.3<sup>+0.3</sup><sub>-0.1</sub> &times; 10<sup>3</sup> (99% CL)</td><td>0.2</td><td>959 d</td><td>Albert et al. (2020)</td></tr><tr><td>HAWC</td><td>3.5<sup>+0.4</sup><sub>-0.2</sub> &times; 10<sup>3</sup> (99% CL)</td><td>1</td><td>959 d</td><td>Albert et al. (2020)</td></tr><tr><td>HAWC</td><td>3.4<sup>+0.4</sup><sub>-0.1</sub> &times; 10<sup>3</sup> (99% CL)</td><td>10</td><td>959 d</td><td>Albert et al. (2020)</td></tr></tbody></table>
The prospects for detection of gamma-ray signals from PBH evaporation for the future Southern Wide Field of View Gamma-ray Observatory (SWGO) have been estimated in López-Coto et al. (2021): it will be sensitive to local PBH burst rate densities $\dot{B}_{\text{PBH}}$ of the order of $\sim 50 \text{ pc}^{-3}\text{yr}^{-1}$, improving by more than one order of magnitude the limits from current experiments.

### 7.1.3 Positron annihilation in the Galactic Center

PBH limits can be obtained with observations of the 511 keV line from the Galactic Center with INTEGRAL. Positrons produced by evaporating PBHs subsequently propagate throughout the Galaxy. As they annihilate, they contribute to the Galactic 511 keV line. Both monochromatic and log-normal mass distributions have been considered, e.g. giving the corresponding limit of $f_{\text{PBH}} < 1$ for $m_{\text{PBH}} \lesssim 2 \times 10^{17}$ g for a log-normal mass distribution with $\sigma = 0.5$ and $m_{\text{PBH}} < 2 \times 10^{17}$ g with $\sigma = 0.1$ (DeRocco and Graham 2019). Estimates for different background assumptions (Ballesteros et al. 2020a) show that for the low ($\times 0.1$) power law background model, $f_{\text{PBH}}$ can reach 1 up to $m_{\text{PBH}} \sim 5 \times 10^{17}$ g. Estimates for different mass distributions and spins (Dasgupta et al. 2020) give the mass limit $m_{\text{PBH}} \sim 3 - 5 \times 10^{18}$ g for the log-normal mass distribution with $\sigma = 1$, $a^* = 0.5 - 1$ and $f_{\text{PBH}} = 1$. It has been recently shown (Keith et al. 2022) that future MeV-scale gamma-ray observatories e-ASTROGAM or AMEGO would be able not only to detect the Hawking radiation from the Inner Galaxy, but also to precisely measure the abundance and mass distribution of PBHs responsible for the Galactic 511 keV line.

### 7.1.4 $e^{\pm}$ observations by Voyager 1

The measurements of interstellar low-energy (sub-GeV) $e^{\pm}$ flux are used to constrain the PBH fraction of DM in the Galaxy (Boudaud and Cirelli 2019). These electrons and positrons would be due to PBH evaporation. The obtained limits are based on local measurements, independent of the cosmological model. This flux must be shielded for Earth-bounded detectors by the magnetic field of the Sun. However, limits can be obtained from Voyager-1 measurements since it is now beyond the heliopause.

The Hawking radiation flux of $e^{\pm}$ is computed within the fully general diffusion-convection-reacceleration model of propagation with the parameters adjusted to the AMS-02 data. For the log-normal mass distribution of PBHs with a width $\sigma$ larger than 1 and a central value $\mu \lesssim 10^{17}$ g, a PBH fraction larger than 1% of DM is excluded by the Voyager 1 data.

### 7.1.5 Ultra-faint dwarf galaxies

The ultra-faint dwarf galaxies (UFDGs), such as Leo T, also provide some limits on light, evaporating PBHs, through the heating of the interstellar medium gas (Kim 2021). The emitted energetic particles can heat and ionize the gas, in particular the fast electrons can deposit a substantial fraction of their energy and heat the interstellar
gas through the Coulomb interaction. The ultra-faint dwarf galaxy Leo T is of particular interest, not only because it is highly dark matter dominated but also because it has a large amount of neutral hydrogen and a relatively low metal abundance. In Kim (2021), limits are obtained for non-spinning PBHs, with masses between $2 \times 10^{15}$ g and $10^{17}$ g as low as $f_{PBH} \lesssim 10^{-3}$, i.e. more stringent than the limits from galactic and extragalactic emission. The possible spin of PBHs would affect the Hawking radiation, generally increasing the emission of particles, in particular ones with large spins. In comparison, stronger limits are obtained for spinning PBHs (Laha et al. 2021) than for non-spinning ones. In both cases, $f_{PBH} = 1$ is claimed to be excluded for $m_{PBH} \sim (2-2.5) \times 10^{16}$ g.

### 7.1.6 CMB limits on Dark radiation

Evaporating PBHs in the mass range [$3 \times 10^{13}$, $5 \times 10^{16}$] g should have injected an exotic amount of electromagnetic energy in the Universe, which has an impact on the CMB anisotropies. This effect has been computed in Stöcker et al. (2018), using the ExoCLASS code that includes electromagnetic cascade calculations, initially developed in the context of WIMP annihilation. For a monochromatic mass function, the obtained bounds exclude PBHs as the main source of dark matter between $3 \times 10^{13}$ and $2.5 \times 10^{16}$ g, and are up to several orders of magnitude more stringent than the limits from the galactic gamma-ray background between $3 \times 10^{13}$ and $3 \times 10^{14}$ g. A similar analysis was conducted in Acharya and Khatri (2020), leading to limits at slightly lower masses, around $1.1 \times 10^{13}$ g. See also Lucca et al. (2020) for the connection with CMB spectral distortions induced by the energy injection from evaporating PBHs. The non-standard energy injection problem has also been reformulated to account for an extended PBH mass distribution in Poulter et al. (2019). In comparison with a monochromatic mass function, the marginalised and log-normal ones have tighter exclusion bounds for heavier PBH masses, while at low PBH masses these two distributions exhibit weaker bounds than the delta distribution. In the future, experiments dedicated to the detection of the CMB spectral distortions or of the 21 cm signal from the dark ages could improve these limits. Finally, it has been pointed out in Nesseris et al. (2020) that the time-dependence of the matter density, as the PBH density decreases as they lose mass via Hawking radiation, can be formulated as an equivalent early dark energy model with varying equation of state $w(z)$. For instance, a population of ultra-light PBHs, decaying around neutrino decoupling, leads to a dark matter-radiation coupling altering the expansion history of the Universe, and alleviating the $H_0$ tension.

### 7.1.7 Neutrinos

For PBH masses between $10^{15}$ and $10^{17}$ g, black hole evaporation via Hawking radiation leads to sizeable fluxes of MeV neutrinos that could be detected with the next generation of detectors like Hyper-Kamiokande (Abe et al. 2018), JUNO (An et al. 2016), DUNE (Acciarri et al. 2016) and THEIA (Askins et al. 2020). Current experiments like SuperKamiokande already provide some limits on the PBH
abundance (Dasgupta et al. 2020; Bernal et al. 2022), but they are not yet competitive with gamma-ray limits. Those future experiments will have the ability to set complementary limits in the asteroid-mass range. The expected sensitivity of these experiments (Hyper-Kamiokande, JUNO, DUNE and THEIA) with various designs has been calculated in Wang et al. (2021), De Romeri et al. (2021), Bernal et al. (2022) for both monochromatic and log-normal PBH mass functions. Abundances as low as $10^{-6}$ could be probed for PBH with zero spin and even lower if PBHs have significant spin. DUNE and THEIA would also complement each other, because they are respectively sensitive to neutrino and antineutrino fluxes.

Neutrinos resulting from PBHs could also be observed in direct detection experiments, through coherent elastic neutrino-nucleus scattering (see Calabrese et al. 2022c).

High energy extra-galactic neutrinos can be used to put constraints on the local PBH burst rate density $\dot{B}_{\text{PBH}}$ for both non-rotating ($a^* = 0$) and rotating ($a^* = 0.99$) black holes. The 90% C.L. limits from 10-years $\nu_\mu$ IceCube data (Capanema et al. 2021) are presented in Table 4.

### 7.1.8 21 cm signal

The 21 cm line is an electromagnetic radiation spectral line created by an atomic transition between the two hyperfine levels of the 1s ground state of neutral hydrogen. This spectral line is of great interest and importance in cosmology as it seems to be the only known way to probe the dark ages of the Universe from the epoch of recombination to reionization. Observations of 21 cm are known to be very difficult to make, primarily because it is a faint signal and is usually plagued by interferences from various sources. However, a global 21 cm signal was claimed to be detected by the Experiment to Detect the Global Epoch of Reionization Signature (EDGES), which shows an absorption feature with an amplitude of $T_{21} \sim 500$ mK centered at redshift $z \sim 17$, much stronger (by about a factor of 2) than expected in the $\Lambda$CDM model (Bowman et al. 2018). Various explanations, such as interactions between the dark matter particles and baryons, have been put forward to explain this anomalous signal. It must be noticed, however, that EDGES results were later contradicted the SARAS3 radiometer measurements (Singh et al. 2022b) that would

<table><caption>**Table 4** Current limits on the local PBH burst rate density from IceCube (Capanema et al. 2021)</caption><thead><tr><td>Experiment</td><td>B&#x0302;<sub>PBH</sub> Upper limit (pc<sup>-3</sup>yr<sup>-1</sup>)</td><td>Burst duration (s)</td><td>T<sub>obs</sub></td></tr></thead><tbody><tr><td>IceCube</td><td>1.7 &times; 10<sup>7</sup> (90% CL, a<sup>*</sup> = 0)</td><td>10</td><td>10 years</td></tr><tr><td>IceCube</td><td>8.2 &times; 10<sup>6</sup> (90% CL, a<sup>*</sup> = 0)</td><td>10<sup>2</sup></td><td>10 years</td></tr><tr><td>IceCube</td><td>4.1 &times; 10<sup>6</sup> (90% CL, a<sup>*</sup> = 0)</td><td>10<sup>3</sup></td><td>10 years</td></tr><tr><td>IceCube</td><td>1.5 &times; 10<sup>7</sup> (90% CL, a<sup>*</sup> = 0.99)</td><td>10</td><td>10 years</td></tr><tr><td>IceCube</td><td>6.2 &times; 10<sup>6</sup> (90% CL, a<sup>*</sup> = 0.99)</td><td>10<sup>2</sup></td><td>10 years</td></tr><tr><td>IceCube</td><td>3.3 &times; 10<sup>6</sup> (90% CL, a<sup>*</sup> = 0.99)</td><td>10<sup>3</sup></td><td>10 years</td></tr></tbody></table>
exclude an astrophysical origin of the EDGES signal. Nevertheless, some limits can
still be established on theoretical scenarios from the absence of a strong 21 cm signal.

In general, during that epoch, any additional sources, e.g. the dark matter
annihilation or decay, will heat up the intergalactic medium (IGM) and increase the
kinetic temperature that eventually determines the strength of the 21 cm signal.
Therefore, observations of the EDGES experiment can, in principle, impose
constraints on the properties of a dark matter scenario. It has been discussed that, in
the context of 21 cm cosmology, PBHs can lead to multiple distinct observable
effects: (i) The heating of the IGM due to their Hawking evaporation; (ii) The
modification of the primordial power spectrum (and hence the halo mass function)
induced by Poisson noise due to PBHs being point sources; (iii) A uniform heating
and ionization of the IGM due to X-rays produced during accretion and (iv) A local
modification of the temperature and density of the ambient medium surrounding
isolated PBHs. Due to their Hawking evaporation, PBHs can affect the evolution of
the IGM, and therefore their mass fraction can be constrained by the global 21 cm
signal. Using this, the initial mass fraction of PBHs is constrained to be $\beta_{\text{PBH}} \sim 2 \times 10^{-30}$ for very light PBHs in the mass range $6 \times 10^{13} \text{ g} \lesssim m_{\text{PBH}} \lesssim 3 \times 10^{14} \text{ g}$ (Yang 2020). Considering PBH evaporation, in Mittal et al. (2022) the following limits have
been found: $f_{\text{PBH}} = 10^{-6.84}(m_{\text{PBH}}/10^{15} \text{ g})^{3.75}$ in absence of X-ray heating and $f_{\text{PBH}} \le 10^{-9.73}(m_{\text{PBH}}/10^{15} \text{ g})^{3.96}$ in presence of X-ray heating. The mass and spin are
fundamental properties of a black hole, and they can substantially affect the BH
evaporation rate. in Natwariya et al. (2021), the lower limits on masses of PBHs
allowed to entirely constitute the DM have been found to be $1.5 \times 10^{17} \text{ g}$,
$1.9 \times 10^{17} \text{ g}$, $3.9 \times 10^{17} \text{ g}$ and $6.7 \times 10^{17} \text{ g}$ for PBH spins 0, 0.5, 0.9 and 0.9999,
respectively.

### 7.1.9 Limits from DM particle production

If evaporating PBHs are at present times a source of the boosted light DM particles,
then the constraints on PBH mass and abundance can be obtained from the
XENON1T data (Calabrese et al. 2022a). Considering DM masses smaller than 1
MeV, for low DM-nucleon cross-sections, the constraints on $f_{\text{PBH}}$ are almost
inversely proportional to the cross-section. However, for larger cross-sections the
results might depend on the model of the DM propagation through Earth and
atmosphere. The constraints are valid for the mass range $0.5 \lesssim m_{\text{PBH}}/(10^{15} \text{ g}) \lesssim 8.0$:
the kinetic energy of particles emitted by more massive PBHs is too low, the less
massive PBHs have evaporated nearly completely. For the DM particles interacting
with electrons similar analyses are presented in Li and Liao (2022), Calabrese et al.
(2022b).

### 7.1.10 Discussion and limitations

It is important to notice that all the limits on the PBH abundance in the asteroid mass
range rely on the common assumption that PBHs evaporate through Hawking
radiation. If the process of black hole evaporation is well established theoretically, its
existence is not proven observationally. Furthermore, it is still unclear whether
quantum gravity effects would affect or not this process, see for example (Anchor-
doqui et al. 2022, 2023; Thoss et al. 2024), where for instance it is shown that the
existence of one mesoscopic extra-dimension leads to a longer life-time of the PBHs
and one can have an all-dark-matter interpretation in terms of PBHs in the mass range
$10^{14} \lesssim m_{\text{PBH}}/g \lesssim 10^{21}$. The efficiency of the evaporation through the different
channels can also depend on the underlying particle physics theory. It is therefore
important to keep in mind these possible limitations. Nevertheless, it is also worth
noticing that such very light PBHs, if they exist, are the ideal target to probe and test
the Hawking radiation hypothesis for the first time. It was also recently claimed that
the memory burden effect limits by a factor two the PBH evaporation time, which
may reopen some part of the allowed mass range, down to $10^6$ kg (Alexandre et al.
2024).

The existence and importance of the Hawking radiation are not the only source of
uncertainty on the abundance limits from PBH evaporation. The different source of
uncertainties have been recently reviewed and critically analyzed in Auffinger
(2022), especially coming from the instrument characteristics, the prediction of the
(extra)galactic photon flux, the statistical method of signal-to-data comparison and
the computation of the Hawking radiation rate. The constraints on PBHs are found to
vary by several orders of magnitude, depending on the hypothesis. Even with an
“ideal” experiment, PBHs can only be probed through Hawking radiation at masses
below $10^{20}$ g, which is an intrinsic limitation of the phenomena.

In order to obtain those results, study the model dependence of the limits related to
Hawking radiation and include all these uncertainties in the data analysis, the
BlackHawk (Arbey and Auffinger 2019) and Isatis (Auffinger 2022) codes
have been developed, allowing to compute the evaporation primary and secondary
spectra of stable or long-lived particles, for any black hole distribution. We envisage
to interface the PBH numerical tool developed in the context of the present review
paper with these key tools.

7.2 Microlensing searches

One of the most stringent constraints on PBH as constituents of the dark matter in the
Universe comes from the microlensing amplification of light from a distant star as a
Massive Compact Halo Object (MACHO) moves across its line of sight.
Microlensing events from lenses of less than one solar mass towards the Large
Magellanic Cloud (LMC) where reported in the 1990s by several groups, but in the
subsequent decades only upper bounds were published, until SUBARU and later
OGLE reported events towards the galactic bulge and nearby Andromeda. The
constraints on the halo abundance of those objects depend very much on the mass,
velocity and spatial distribution in the halo. Most of these constraints are reported for
a standard Navarro–Frenk–White (NFW) halo with virial velocities and a spatially
uniform distribution of monochromatic (same mass) MACHOs, which put stringent
constraints in the range of masses from $10^{-9}$ to $30 M_{\odot}$. Above 10 $M_{\odot}$, it is essentially
impossible to derive bounds on MACHOs due to the fact that the amplification light
curve lasts more than a decade, while no microlensing survey has monitored continuously the light from stars for so long. Below $10^{-9} M_{\odot}$, the wavelength of visible light is comparable to the Schwarzschild radius of the compact object, and we enter the regime of wave optics, where diffraction makes the constraints essentially disappear. One could search for PBHs with $m_{\text{PBH}} < 10^{-9} M_{\odot}$ with light of shorter wavelength, like X-rays or gamma-rays, but sources that emit in those bands are typically not as abundant and stable as stars.

In the following subsections, we give a brief summary of the different microlensing constraints coming from different surveys.

### 7.2.1 Microlensing from stars in the Magellanic clouds

Microlensing observations of stars in the Large and Small Magellanic Clouds probe the fraction of the Galactic halo in MACHOs in a certain mass range (Paczynski 1986). The optical depth of the halo towards LMC and SMC is related to $f_{\text{PBH}}$ by $\tau_{\text{L}}^{\text{SMC}} = 1.4\tau_{\text{L}}^{\text{LMC}} = 6.6 \times 10^{-7} f_{\text{PBH}}$ for the standard halo model (Alcock et al. 2001). The MACHO project detected lenses with mass around $0.5 M_{\odot}$ but concluded that their halo contribution could be at most 10% (Hamadache et al. 2006), while the EROS project excluded $6 \times 10^{-8} M_{\odot} < m_{\text{PBH}} < 15 M_{\odot}$ from dominating the halo (Allsman et al. 2001; Tisserand et al. 2007). Since then further limits in the range $0.1 M_{\odot} < m_{\text{PBH}} < 20 M_{\odot}$ have come from the OGLE experiment (Wyrzykowski et al. 2011a, b). Very recently, even more stringent limits were claimed by OGLE that are also extended to heavier masses (Mróz et al. 2024a, b).

### 7.2.2 Microlensing from stars in the galactic bulge

Niikura et al. (2019a) have used data from a 5-year OGLE survey of the Galactic bulge to place much stronger limits in the range $10^{-6} M_{\odot} < m_{\text{PBH}} < 10^{-4} M_{\odot}$, although they also claim some positive detections. The galactic bulge is a dense region of stars towards the galactic center where microlensing events are likely to occur. In fact, the OGLE collaboration has detected numerous microlensing events towards the bulge (Wyrzykowski et al. 2016), searching specifically for dark objects that amplify the light of stars behind them. Moreover, these stars are sufficiently close to Earth that one can use parallax measurements by GAIA to break degeneracies between mass and distance which plague magnification microlensing events. This use of the so-called astrometric microlensing has allowed Wyrzykowski and Mandel (2020) to conclusively detect black holes in the mass gap, in the range $2 - 5 M_{\odot}$. Niikura et al. (2019a) also report microlensing events towards the galactic bulge in a lower mass range of Earth-mass primordial black holes. These are tantalizing hints of a multimodal population of primordial black holes in the halo of our galaxy, where non of the mass ranges are sufficient to comply 100% of all the dark matter halo.
7.2.3 Microlensing from stars in Andromeda

The Andromeda galaxy (M31) is too far away to resolve individual stars with a
typical dedicated (small, 2 m) telescope, as most microlensing surveys do. Therefore
they use a different technique called pixel lensing, which corresponds to single pixels
brightening up over the length of the survey with a characteristic Paczinsky
amplification curve due to a passing compact object. This technique was used by the
POINT-AGAPE collaboration (Calchi Novati et al. 2005), who reported six
microlensing events, which they argued could be interpreted as self-lensing of stars
and this allowed them to conclude that 20% of the Milky Way halo in the direction of
M31 could be in the form of 0.5 – 1.0 M<sub>☉</sub> primordial black holes.

More recently, it was reported by Niikura et al. (2019b) the observation, with the Hyper Supreme Cam (HSC) on Subaru Telescope, of a microlensing event towards Andromeda by a compact body with mass in the range 10<sup>-10</sup> – 10<sup>-5</sup> M<sub>☉</sub>. Give the detailed statistics of the observation, they claimed this could provide evidence for a non-zero component, at the level of few percent, for the dark matter halo to be in the form of primordial black holes.

7.2.4 Microlensing from quasars

Hawkins has claimed for many years that quasar microlensing data suggest the existence of PBH dark matter (Hawkins 1993). He originally argued for Jupiter-mass PBHs but later increased the mass estimate to 0.4*M<sub>☉</sub>* (Hawkins 2007). More recently, the detection of 24 microlensed quasars (Niikura et al. 2019a) suggests that up to 25% of galactic halos could be in PBHs with mass between 0.05 and 0.45*M<sub>☉</sub>*. These events could also be explained by intervening stars, but in several cases the stellar region of the lensing galaxy is not aligned with the quasar, which suggests a population of subsolar halo objects with *f<sub>PBH</sub>*≥0.01. Indeed, Hawkins has argued that the most plausible microlensers are PBHs, either in galactic halos or distributed along the lines of sight to the quasars (Hawkins 2020b).

7.2.5 Caustic crossing in giant arcs

Caustic crossing events in giant arcs are recently discovered phenomena that can be
used to probe the existence of dark compact objects such as PBHs (Oguri et al.
2018). These observations correspond to highly magnified images of individual stars,
inside giant arcs produced by the strong gravitational lensing of galaxy clusters. This
technique is complementary to microlensing of stars and quasars. Compared to
microlensing searches, the main difference comes from the existence of high
magnification regions, convergence and shear fields, due to the cluster potential.
Such environments are also required for quasar microlensing but here the source is
different and is a simple star, not a quasar.

Kelly et al. (2018) reported the observation of such an event, the MACS J1149
Lensed Star 1 (LS1), also known as *Icarus*. It corresponds to a faint transient near the
critical curve of the massive cluster MACS J1149.6+2223. It is interpreted as a
luminous star in a host galaxy at a redshift $z = 1.49$ which is highly magnified, by a factor of several thousands at the peak brightness, by a compact object that is very close to the critical curve of the foreground lens. The lens properties of this event seem to be consistent with a microlensing event induced by a star that contributes to the intra-cluster light but it can also be used to set constraints on the abundance of compact objects in this galaxy cluster. This leads to claimed limits, with different assumptions on the object transverse velocity, going down to $f_{PBH} < 0.1$ (Oguri et al. 2018). These limits are independent on the object mass in the range between $10^{-5} M_{\odot}$ and $10 M_{\odot}$ and become progressively less stringent at larger masses, up to $100 - 1000 M_{\odot}$ depending on the assumptions.

Caustic crossing in giant arcs therefore appears to be a promising technique to probe the existence of PBHs that is complementary to other microlensing searches.

### 7.2.6 Microlensing of supernovae

One can also constrain the abundance of PBHs with stellar masses using supernovae surveys. A certain number of supernova light curves are expected to be lensed by PBHs, exhibiting a magnification of order $\mu \sim 0.1$ if PBHs constitute all the dark matter (Zumalacarregui and Seljak 2018). Nevertheless, it has been debated whether the current JLA and Union data already exclude $f_{PBH} = 1$, especially in the case of broad PBH mass distributions (García-Bellido et al. 2018). The main discrepancy has to do with the assumed statistics and the treatment of outliers and other systematics, like the probability of amplification due to discrete sources along the line of sight (Boscá et al. 2022). Interestingly, the SNe limits are of the same order from $10 M_{\odot}$ up to $10^4 M_{\odot}$.

### 7.2.7 Femtolensing

Additional limits on the PBH abundance have been established based on the femtolensing of gamma-ray bursts. Compared to microlensing, the terminology of femtolensing is used to refer to the very small angular distance between lensed images and applies to objects of very small masses. The effect is probed by searching for interferometry patterns in the energy spectrum of lensed objects (Gould 1992). One must nevertheless take into account that the wavelength of the electromagnetic radiation can become comparable to the Schwarzschild radius of the PBH, which implies a lower bound on the mass range probed with this technique. In Barnacka et al. (2012), a search for femtolensing events has been performed on the FERMI observations of gamma-ray bursts for which the redshifts were known. The absence of femtolensing events has been used to set constraints on the dark matter fraction in PBHs of mass between $5 \times 10^{17} \text{ g}$ and $10^{20} \text{ g}$, with a limit down to $f_{PBH} \lesssim 0.04$.

### 7.2.8 Discussion and limitations

The robustness of microlensing limits is of very high importance and still controversial because this is the only observational technique to probe the abundance
of PBHs smaller than a solar mass, down to $10^{-10} M_{\odot}$. Another debated question is whether microlensing surveys of quasars or towards the galactic center only provide a limit on PBHs due to the lack of observations or also hint at their existence due to a series of observations.

The contribution of PBHs to the dark matter is also an interesting question. Deriving limits from microlensing surveys of stars in the Magellanic clouds depends on numerous astrophysical assumptions (velocity dispersion, galactic dark matter profile, etc...) and varying limits can be obtained with different hypothesis (Green 2016, 2017; Hawkins 2015, 2020b; Calcino et al. 2018; García-Bellido and Hawkins 2024). Wide mass distributions can also significantly alter the limits on their abundance, but without totally removing them. The clustering of PBHs was the main argument invoked to totally suppress the microlensing limits (García-Bellido and Clesse 2018; Carr et al. 2021b), even though it has been shown that this is not the case assuming initial Poisson distributed PBHs (Petač et al. 2022; Gorton and Green 2022). Furthermore, in the case of a large initial spatial clustering distribution beyond Poisson, it has been also recently shown that the combination of constraints from microlensing and Lyman-α forest completely rules out the parameter space for stellar-mass PBHs to be a dominant constituent of the dark matter (De Luca et al. 2022b).

7.3 Dynamical limits

7.3.1 PBH capture by neutron stars or white dwarfs

The asteroid-mass range, between $10^{-16} M_{\odot}$ and $10^{-10} M_{\odot}$ ($10^{17} - 10^{23}$ g), is the only one for which there is a general consensus that PBHs can contribute up to the entirety of the dark matter. Nevertheless, if this mass window remains open nowadays, this has not been always the case. Some constraints had been proposed, coming from observations of neutron stars and white dwarfs (WD) in globular clusters (Capela et al. 2013b, a; Pani and Loeb 2014; Defillon et al. 2014) that can capture asteroid-mass PBHs falling at the center of the star and rapidly swallowing its material until it becomes a black hole. Such constraints therefore arise from the combination of a dynamical and an accretion process. But those limits have been relieved due to large uncertainties. It was initially claimed that $f_{PBH}$ cannot exceed between $10^{-2}$ and 0.2, depending on the DM density in globular clusters. However, this approach is limited by our knowledge of the DM constitution of globular clusters. Using updated prescriptions, the limits from Capela et al. (2013a), Pani and Loeb (2014) have been then relieved in Montero-Camacho et al. (2019), re-opening the asteroid-mass range to full PBH-DM scenarios.

In Graham et al. (2015), it was proposed that the transit of PBHs through a white dwarf causes localized heating around the BH trajectory due to dynamical friction, eventually triggering a runaway thermonuclear fusion and causing the WD to explode as a supernova, even if the WD mass is below the Chandrasekhar limit. The constraints on PBHs were derived by using two classes of observations. One class comes from direct observations of WDs with known masses, whose existence places bounds on PBHs that are abundant enough that they would have transited the WD
with high probability. The second class arises from the measured rate of type Ia supernovae, that constrain PBHs with a low abundance and thus have a low probability of transit through a WD. The shape of the observed distribution of WD masses (with masses up to $1.25M_\odot$) excludes PBHs with masses $\sim 10^{19} - 10^{20}$g from being a dominant constituent of the local DM density, whereas the type Ia supernova rate disfavors BHs with masses between $10^{20} - 10^{22}$g. However, the latter is not a robust exclusion due to uncertainties in the methods used to determine the WD population with masses larger than $\sim 0.85M_\odot$ (necessary condition for a WD to explode due to a BH), and this mass could possibly change if the current simulation data used to derive it is refined. Measurements of WD binaries in GW observatories could be used in the future to further strengthen the BH bounds.

More recently, limits on PBH in the interesting asteroid-mass range between $10^{18}$ g and $10^{21}$ g based on the capture of PBHs by Sun-like stars in ultra-faint dwarf galaxies leading to its destruction, have been calculated in Esser and Tinyakov (2023). If one requires that no more than a fraction of stars $\xi = 0.5$ is destroyed, one obtains a limit on the PBH abundance of about $f_{PBH} \lesssim 0.3$ in the above-mentioned mass range. But the constraints strongly depend on the possible value of $\xi$ that is still fairly unknown. In Smirnov et al. (2024), it was suggested that the recent discovery of a population of faint supernovae (Ca-rich gap transients) may come from the explosion of white dwarfs, possibly due to PBH captures. Their unusual spatial distribution, predominantly at large distances from their presumed host galaxies, supports this hypothesis and the authors have shown that their spatial distribution matches well the distribution of dwarf spheroidal galaxies.

### 7.3.2 Ultra-faint dwarf galaxies

If PBHs above $10 M_\odot$ constitute an important fraction of the dark matter, their regular encounters in highly dark matter dominated ultra-faint dwarf galaxies (UFDGs), such as in Segues I or in Eridanus II, should have induced their dynamical heating, allowing them to reach a half-light radii larger than 10 parsecs. This argument also applies to star clusters in UFDGs, like for Eridanus II. Relatively stringent limits from UFDGs have been first established in Brandt (2016), for monochromatic mass functions. These have been extended to log-normal PBH mass functions, including the effect of a central intermediate-mass BH, in Green (2016), or in Li et al. (2017). All these works basically exclude $f_{PBH} = 1$ for $m_{PBH} > 10 M_\odot$. These limits are relatively solid because they rely on well-known Newtonian gravitational dynamics. Nevertheless, one should point out that these limits can be shifted if UFDGs host a central intermediate-mass black hole.

Another way to constrain PBHs from UFDGs is to consider the radial luminosity profile, given that stars must also be heated up by PBHs. In this case, Segue 1 provides the most stringent limit, excluding $f_{PBH} = 1$ for $m_{PBH} \sim 2 M_\odot$ (Koussiappas and Loeb 2017a; Stegmann et al. 2020). However, these limits rely on astrophysical assumptions about the stellar mass distribution and they could change for wide-mass PBH models, depending on the existence and mass of a central intermediate-mass black hole, etc. Finally, UFDGs also provide an argument in favor of solar-mass
PBHs constituting an important fraction of the dark matter, given that there is no
observation of such objects with a radius less than about 20 parsecs, whereas they
would still be above the magnitude limit of the current instruments.

7.3.3 Dark Matter profile of dwarf galaxies (core-cusp problem)

The regular PBH-PBH interactions in a dense environment should prevent the formation of cusps in the central region of dwarf galaxies, possibly explaining the core-cusp problem. As first pointed out in Clesse and García-Bellido (2018), the gravitational scattering cross-section for PBHs in the stellar-mass range are of order (0.1 − 1)cm²/g and so PBHs are in essence a well motivated self-interacting dark matter model. But contrary to particle dark matter, this does not need exotic physics but relies on the only well-known interaction that is not included in the standard model of particle physics, gravitation. The core-cusp problem in the context of PBHs was then further explored in Boldrini et al. (2020). They found that the mechanism works well for PBHs of mass between 20 and 100 M☉, if they have fPBH ∼ 0.01, but it would be less efficient and would allow larger abundances for PBHs of the order of the solar mass. Therefore, observations of low-mass dwarf galaxies provide at the same time a new limit on the PBH abundance and a possible hint for their existence. These limits imposed on the PBH abundance have been represented in Fig. 38.

7.3.4 Wide halo binaries

Wide stellar binaries can be destructed by the encounters with PBHs, therefore the existence of the undisrupted ones in the galactic halo also provides constraints on the PBH abundance (Carr and Sakellariadou 1999).

Using the analysis of Quinn et al. (2009), the following constraint at 2σ
confidence level has been found (Carr et al. 2021c): f<sub>PBH</sub> < 500 M<sub>☉</sub>/m<sub>PBH</sub> flattens
above M < 1000 M<sub>☉</sub> with f<sub>PBH</sub> ≲ 0.5. Using the results of (Monroy-Rodríguez and
Allen 2014), where the low-mass limit was reduced from 500 M<sub>☉</sub> to 21 − 78 M<sub>☉</sub>,
f<sub>PBH</sub> flattens at m<sub>PBH</sub> < 100 M<sub>☉</sub> with f<sub>PBH</sub> ≲ 0.1 (Carr et al. 2021c). In Tian et al.
(2020), by using a sample of 4351 halo wide binaries from the Gaia catalog, a break
at ≲ 0.1 pc in the separation distribution has been detected. Any break in the power
law for the separation distribution of wide halo binaries can be attributed to a
disruption of the binary as a result of encountering MACHOs (e.g. PBHs). Therefore,
this break implies m<sub>PBH</sub> > 10 M<sub>☉</sub>. These limits come from the expected disruption of
such wide binaries when they encounter PBHs. They therefore strongly depend on
the number density of PBHs, and they could vary if, for instance, a majority of PBHs
live in clusters.

7.3.5 Disruption of stellar streams

Stellar streams have been observed in the Milky Way halo and are likely formed by
the tidal stripping of progenitors. These would be perturbed by encounters with dark
matter subhalos, such as PBH clusters. Some methods have been proposed to
distinguish the origin of perturbers, in particular the PBH-DM hypothesis from standard particle CDM (Montanari and García-Bellido 2022). Stellar streams like GD-1 could probe the abundance of haloes down to $10^3 M_\odot$, providing a test of the existence of PBHs.

### 7.3.6 PBHs in the solar system

The dynamics of objects in the Oort cloud suggests the possible existence of a not yet observed planet, referred as Planet Nine. It could also be a dark compact object several times the mass of the Earth, in orbit around the Sun in the Oort cloud. It has been suggested that Planet Nine is a PBH captured by the Sun (Scholtz and Unwin 2020) a long time ago. It was later proposed to search for the existence of such a PBH in the Oort cloud by probing the flares from accretion of small Oort cloud objects with LSST (Siraj and Loeb 2020), which could result in improved limits on PBH abundances, or by probing the associated Hawking radiation (Arbey and Auffinger 2020).

Alternatively, it has been proposed to probe asteroid-mass PBHs in orbit around the Earth by searching for their gravitational signatures in gravimeter networks (Namigata et al. 2022). But the obtained limits on the PBH abundance with the current gravimeter sensitivity are not yet relevant. Near-Earth PBHs could also leave imprints in gravimeters or in GNSS orbit products. One way to detect a PBH passing through the Earth would be to observe the induced seismic signatures (Luo et al. 2012).

All this nevertheless remains highly speculative and it is not clear, for instance, if a PBH captured by the Sun (or by the Earth) would remain on a stable orbit, or if its orbit would likely pass near the Earth. However, it must also be pointed out that discovering a PBH inside the Solar System would have tremendous implications and could even, in principle, allow in-situ investigations of a black hole (Witten 2020).

### 7.3.7 Discussion and limitations

Again, there are a series of limitations and uncertainties that restrict the range of applicability of the dynamical limits. In particular, let us mention that if most PBHs are regrouped in clusters, the limits from wide halo binaries would be strongly affected and possibly suppressed. These clusters could even be the ultra-faint-dwarf galaxies, for which observations show a minimum radius around a few tens of parsecs. Dynamical heating of UFDGs by PBHs can therefore be seen both as a limit and as a clue for the existence of such PBHs. For UFDGs, other sources of uncertainties come from the stellar population model, the exact DM profile and the existence of a central intermediate-mass black hole that typically makes limits less stringent. This illustrates that the border between observational limits and evidence is blurred and connected, and above all strongly model dependent. A similar conclusion applies to the asteroid-mass region, for which uncertainties on the DM content of globular clusters and the details of the tidal capture of small PBHs by neutrons stars or white dwarfs still make this range viable for a full PBH-DM scenario.
7.4 Accretion limits

7.4.1 Cosmic microwave background (CMB) anisotropies

Primordial black holes which experience a phase of accretion at the present epoch are
subject to several constraints depending on their characteristic environment.

CMB anisotropies are impacted by accreting PBHs around the time of the
recombination and until the completion of the reionization, leading to stringent limits
on the abundance of PBHs heavier than around $10-100 M_{\odot}$. The main idea behind
these constraints is the following: in the early Universe, PBHs accrete primordial gas
and then convert a fraction of the accreted mass to radiation. This results in the
injection of energy into the primordial plasma, which affects its thermal and
ionization histories, leading to CMB distortions in its frequency spectrum, as well as
signatures in the CMB temperature and polarization angular power spectra.

When accretion comes into play, the first aspect to consider is its geometry, i.e. spherical accretion or accretion disk formation. The former is obtained when the characteristic angular momentum of the accreted gas at the Bondi radius (distance from the center at which the escape velocity equals the sound speed) is smaller than the angular momentum at the ISCO (Innermost Stable Circular Orbit) radius. The second aspect that must be taken into account is the local feedback. Indeed, radiation emanating from the accreting PBH may heat and ionize the accreting gas, influencing the radiative output. Moreover, as discussed in Ali-Haïmoud and Kamionkowski (2017), one can assume that (1) In the outermost region of the accretion flow, the ionization fraction is approximately equal to the background value. Radiative feedback is neglected, and Compton drag and cooling are included in the calculation. When the temperature of the gas reaches ~ $10^4$ K, the gas gets collisionally ionized. This first limiting case is refered to as collisional ionization; (2) In the innermost region, the radiation from the PBH photoionizes the gas, which becomes fully ionized and adiabatically compressed. This second limiting case is referred to as photoionization. Furthermore, the accreting PBHs are moving with respect to the ambient gas with some relative velocity, which does not allow for a perfectly spherically-symmetric accretion. This relative velocity is comprised of (1) A Gaussian linear contribution on large scales $v_L$ (derived in Ali-Haïmoud and Kamionkowski 2017); (2) A small-scale contribution $v_{NL}$ due to non-linear clustering of PBHs (not examined in Ali-Haïmoud and Kamionkowski 2017).

Another key quantity is the radiative efficiency factor $\epsilon$, which is strongly related to the geometry of the accretion flow and the accretion rate. It usually takes the value $\epsilon \sim 0.1$ for moderate or low disk accretion rate, that corresponds to the value inferred from BH observations. In the case of spherical accretion, one gets smaller values of the order of $10^{-5}$. Finally, to compute the CMB bound, it is crucial to evaluate the total energy injection rate, by taking into consideration that not all radiation is equally effective. To compute the impact on the CMB, one has to quantify what amount of injected energy is actually deposited into the medium by using the *energy deposition functions*, and the last step is to include the spectrum of the radiation emitted through BH accretion. Let's note that each accretion scenario has a
distinctive energy injection history, that could in principle be distinguished in the CMB angular power spectrum (Poulin et al. 2017).

For PBHs more massive than ~1 $M_\odot$, strong constraints were derived by Ricotti et al. (2008), while in Ali-Haïmoud and Kamionkowski (2017), CMB limits on the PBH abundance were re-examined in detail, based on previous works. In particular, in order to set the most conservative constraints on the PBH abundance, it is useful to compute the minimum PBH luminosity that is physically plausible and the accretion rate, taking into account the Compton drag and Compton cooling by CMB photons. Assuming a quasi-spherical accretion flow, a steady-state flow and a Newtonian approximation, the bounds derived show that in the collisional ionization limit, CMB anisotropy measurements by *Planck* exclude PBHs with masses $M \gtrsim 10^2 M_\odot$ from being the dominant DM component. In the photoionization limit, this threshold is lowered to ~10$M_\odot$. These constraints are weaker than those derived in Ricotti et al. (2008).

As a side note, it is shown in Ali-Haïmoud and Kamionkowski (2017) that local Compton heating by the radiation produced by the accreting PBH can be safely neglected, i.e. one can neglect local thermal feedback. It is also demonstrated that it is not self-consistent to assume that ionizations proceed exclusively through photoionizations or collisional ionization, i.e. the level of ionization feedback is expected to lie between the two limiting cases. Finally, it was also derived that global heating of the plasma due to accreting PBHs does not leave any (observable) imprint on CMB spectral distortions.

In Poulin et al. (2017), it is argued that the spherical accretion approximation is likely not valid and that an *accretion disk* should form during the dark ages, between recombination and reionization. More precisely, all plausible estimates suggest that a disk forms soon after recombination, mainly due to the fact that, already before recombination, stellar-mass PBHs are in a non-linear regime (i.e. they are clustered in halos of bound objects, e.g. PBH binaries or big clumps containing thousands of PBHs). The main criterion necessary for disk formation considered in Poulin et al. (2017) is that the angular momentum of the material at accretion distance is sufficient to keep the matter in Keplerian rotation at a distance $\gg 3r_s$, where $r_s$ is the Schwarzschild radius. To accumulate angular momentum, the accreted material must have significant velocity or density differences. Based on the angular momentum of the accreted gas, one can derive the condition required for a disk to form (Poulin et al. 2017). For instance, the disk formation criterion is likely satisfied if we consider an ideal, free-streaming homogeneous gas moving at a bulk motion with velocity of order $v_L$ without dispersion, or when PBHs constitute an important fraction of the DM with $M \gtrsim M_\odot$.

The conservative constraints obtained in Poulin et al. (2017) (at 95% C.L. using data from *Planck* high-l TT TE EE+lensing) are represented in Fig. 38. One notices that the presence of disks improves the CMB constraints on PBH by at least two orders of magnitude, so that PBHs with masses $M \gtrsim 2M_\odot$ (with monochromatic mass distribution) are excluded from accounting for the totality of the DM. The constraints derived for the monochromatic case were also extended to a log-normal mass
distribution, and become more restrictive for a broader distribution (Poulin et al. 2017).

In Serpico et al. (2020), they extended previous calculation (Poulin et al. 2017) to account for the accretion of dominant, non-PBH DM particles onto PBH. The increased gravitational potential felt by the baryons enhances baryonic accretion and the PBH luminosity constrained by the CMB. The accreting DM halos around each PBH are treated via a toy model and numerical simulations, assuming that the ionization of gas in the dark ages due to accretion onto PBHs is homogeneous.

The toy model considered in Serpico et al. (2020) is used to determine the upper limit to the mass growth of the PBH via DM accretion. The most optimistic scenario for PBH growth is for a cold DM with no dispersion. The time evolution of the radius of a mass-shell around a PBH is found to be $M_{halo} \simeq (3000/1+z)M_{PBH}$, which is expected to break down at very late times. This conjecture is then verified via cosmological N-body simulations, where in order to model the halo growth, it is required that the DM halo is composed of many DM particles when it becomes comparable to the PBH mass. An important result is that when the halo is compact, the accreting baryons see a BH whose effective mass is the sum of the PBH and the DM halo masses. If the halo is large, only a fraction of its mass contributes to the accretion. This is the model used in Serpico et al. (2020) to compute the impact of PBH accretion onto the CMB.

The CMB constraints of Serpico et al. (2020) are obtained by taking into account the effect of accretion within DM halos. Supposing that PBHs are accreting at Eddington luminosity (luminosity at which accretion is balanced by radiation pressure in a spherical system), one gets the limit $f_{PBH} < 2.9 \times 10^{-9}$. In absence of DM halos the constraints at 95% C.L. are stronger (by a factor ~4) than bounds derived in Poulin et al. (2017). Accounting for the DM halo leads to improved bounds for $f_{PBH} \lesssim 10^{-2}$ and $f_{PBH} \lesssim 0.2$ for the disk and spherical accretion cases, respectively. Eventually, the bounds flatten when $M \gtrsim 10^4 M_\odot$ as a consequence of the accretion that reaches the Eddington limit. CMB constraints exceed the ones coming from the non-observation of mergers by LVK, BBN, the non-observation of the SGWB, etc., for masses $M \gtrsim (20-50) M_\odot$. The CMB bounds remain dominant until $10^{3.5} M_\odot$, when they come close to the BBN ones. Specifically, it is excluded that PBH constitute the totality of the DM over the stellar mass range, and around $10 M_\odot$ no more than $f_{PBH} \sim 0.1$ is allowed (Serpico et al. 2020).

Finally, a possible relaxation of the limits from CMB anisotropies has recently been suggested in Facchinetti et al. (2023), Piga et al. (2022), Agius et al. (2024). The work of Facchinetti et al. (2023) relies on a more realistic accretion model based on hydrodynamical simulations and conservative assumptions for the emission efficiency. This leads to limits that are up to 2 orders of magnitude less stringent than previously estimated, between $10$ and $10^4 M_\odot$, which reopens the possibility that PBHs might explain at the same time (at least a fraction of) the dark matter, some of the LVK binary BH mergers, and the existence of supermassive BHs. The work of Piga et al. (2022) takes a step towards the development of a more realistic PBH accretion by accounting the contribution of outflows, for various accretion geometries, ionization models and mass distributions in absence and in presence
of mechanical feedback and non-thermal emissions due to the outflows. More conservative limits were more recently obtained from a critical re-investigation of all the possible effects (Agius et al. 2024). In general, it seems that PBH accretion is rather complex physical process that is hard to assess, which induces large uncertainties in all the accretion-based limits.

### 7.4.2 21 cm signal

There are several limits based on the EDGES detection of the absorption signal of the global redshifted 21 cm line in the range $15 < z < 22$, see however the previous discussion about the exclusion of an astrophysical origin of the EDGES signal by SARAS3. The Poisson distribution of PBHs can lead to isocurvature perturbations and consequently to early structure formation. The DM annihilation is enhanced in early formed dense halos and modifies the ionization and temperature evolution of baryons. Using the redshift dependence of the EDGES signal (Tashiro and Kadota 2021a), the following constraint on PBH abundance has been found: $f_{\text{PBH}} < 10^{-3.4}$ for $f_{\text{ann}}\langle\sigma v\rangle/m_\chi = 3 \times 10^{-28} \text{ cm}^3/\text{s/GeV}$, where $m_\chi$ denotes the DM mass.

PBHs can accrete the surrounding baryonic matter as well as DM particles from the surrounding halo. This results in the emission of high energy photons that are injected into the IGM, leading to changes in its thermal and ionization history via ionization, excitation and heating. For the PBHs accreting only baryons (Hektor et al. 2018), the following approximation for the constraint on the PBHs fraction has been found: $f_{\text{PBH}} \le C(\beta)(0.15/f_{\text{E}})(\lambda/0.01)^{-(1+\beta)}(M_{\text{PBH}}/10M_{\odot})^{-(1+\beta)}$, where $C(\beta) = 0.00015 + 0.00051\beta + 0.0091\beta^2$ if $\beta \le 0.37$ and $C(\beta) = 0.019\beta^{2.5}$ if $\beta > 0.37$ (here $\lambda$ is the mass accretion parameter, $T_k = 8K$ the maximal gas temperatures allowed, and $f_{\text{E}} < 1$ the effective energy absorption factor, assuming a monochromatic mass function). Including also DM particles accretion, in Yang (2021) for a differential 21 cm brightness temperature $\delta T_{21} \lesssim -100 \text{ mK}$ the limit $f_{\text{PBH}} \lesssim 2.6 \times 10^{-5}$ has been obtained for $m_{\text{PBH}} = 10M_{\odot}$ and $f_{\text{PBH}} \lesssim 2.6 \times 10^{-6}$ for $m_{\text{PBH}} = 10^4 M_{\odot}$. These limits are weakened by a factor of order 3 for $\delta T_{21} \lesssim -50 \text{ mK}$.

In Mena et al. (2019), it has been shown that future experiments like SKA and HERA could potentially improve the existing CMB bounds by more than an order of magnitude. The prospects for constraining PBHs with the future 21 cm forest observations have been discussed in Villanueva-Domingo and Ichiki (2022).

In Hasinger (2020) the contribution of accretion on PBHs to the observed low-frequency cosmic radio background (CRB), and thus to the EDGES signal, has been estimated. It has been found that for radio-quiet PBHs the contribution to the CRB constitutes at 1.4 GHz a fraction of 0.1% of the observed synchrotron radio background, mostly accumulated at $z \gtrsim 20$. This increases the depth of the 21 cm absorption line by only about 30%. However, some fraction of the radio-loud PBHs (e.g. 5% with 1000 times higher fluxes), as observed in the AGN population, will easily provide the 5% excess high-redshift radio background flux necessary to explain the EDGES observation.
7.4.3 Cosmic X-ray, infrared and radio backgrounds

Right after the first detection of gravitational-waves from a binary black hole merger, it has been proposed that a PBH origin would be supported by the spatial correlations between the source-subtracted infrared and X-ray backgrounds (Kashlinsky 2016), which have been measured at a 5σ level (Kashlinsky et al. 2018, 2005) and are unexplained by known sources. These correlations can be induced by DM halos seeded at higher redshift than in the standard CDM scenario by the Poisson fluctuations in the PBH distribution, if one has $f_{\text{PBH}}m_{\text{PBH}}/M_{\odot} \sim \mathcal{O}(1)$. In turn, matter accretion onto PBHs leads to X-ray emission while star formation in these halos leads to a spatially correlated infrared radiation. The possibility of this scenario was confirmed and explored in more details in Hasinger (2020), Cappelluti et al. (2022).

Recently, it has also been investigated if the observed excess of the cosmic X-ray and radio cosmic background by Chandra and ARCADE 2, suggesting an undiscovered population of emitters, could be explained by a high-redshift population of PBHs (Ziparo et al. 2022) distributed in DM halos and in the intergalactic medium. It was shown that the emission should dominantly come from small DM halos, with masses of order or smaller than $10^6 M_{\odot}$ at high redshifts $z \gtrsim 6$. Interestingly, it was also found that even if PBHs account for the excess in the radio background, it cannot at the same time explain the X-ray cosmic background. This leads to limits on the PBH abundance, in particular in the $1-100 M_{\odot}$ mass range, which are stronger than the microlensing limits or the bounds from ultra-faint dwarf galaxies, with $f_{\text{PBH}} \lesssim 9 \times 10^{-3} M_{\odot}/m_{\text{PBH}}$.

7.4.4 Limits from accretion at present epoch

A population of PBHs with masses around $\mathcal{O}(10 - 100) M_{\odot}$ in the Milky Way would be responsible for radio and X-ray emission if they experience a phase of accretion from the interstellar gas. By comparing the predicted emission with observation data in the VLA radio and Chandra X-ray catalogs, one can set a constraint on the fraction of dark matter in the form of PBHs to be of the order of $f_{\text{PBH}} \lesssim (10^{-2} - 10^{-3})$ in the relevant mass range, depending on the assumed PBH mass function (Gaggero et al. 2017; Manshanden et al. 2019). This constraint is shown in Fig. 38 and denoted “Xr”.

Second, the abundance of PBHs in the stellar or intermediate mass ranges, which interact with the interstellar medium emitting significant fluxes of X-ray photons, may be constrained by using the observed number density of X-ray binaries. In particular, by comparing the emitted radiation signals with the electromagnetic data, one can constrain the PBH abundance to be smaller than $f_{\text{PBH}} \lesssim (10^{-2} - 10^{-3})$ in the considered mass range (Inoue and Kusenko 2017). This constraint is labelled “XrB” in Fig. 38.

Finally, a cosmology-independent constraint can be set on the PBH abundance based on the absence of gas heating in the interstellar medium due to PBH interactions. In particular, photon emission from accreting PBHs, dynamical friction, winds and jets emission from accretion disks may be constrained using data from the
Leo T dwarf galaxy, setting a bound on PBHs with masses in the range ($10^{-2} - 10^6)M_{\odot}$ up to the percent level (Lu et al. 2021; Takhistov et al. 2022).

### 7.4.5 Discussion and limitations

Accretion limits, in particular limits from CMB temperature anisotropies, are claimed to set the most stringent limits on the abundance of PBHs between about $10M_{\odot}$ and $10^4 M_{\odot}$, i.e. in the range that is particularly relevant for GW observations with ground-based and space-based detectors. Their strengths nevertheless depend on the details of the accretion (e.g. disk or spherical accretion). Given the possible complexity of the accretion process, it is possible that those limits significantly change if different assumptions are chosen, for instance for the velocity-dependence of the accretion process. One should also notice that the Poisson-induced clustering of PBHs were not taken into account in the calculation of these limits. Other sources of uncertainties come in the computation of the injected energy and the resulting changes in the evolution of the ionized fraction at high redshfit. Limits from the X-ray and radio emission in the galactic center have their own uncertainties, e.g. related to the dark matter distribution. The limits that were presented were nevertheless considered as relatively conservative.

An important implication that was recently discovered for the limits from CMB temperature anisotropies is that the existence of a particle dark matter component in addition to PBHs reinforce those limits on the abundance of PBHs, in such a way that they may become incompatible with an explanation of GW observations.

## 7.5 Indirect constraints from density fluctuations

### 7.5.1 CMB distortions

Enhanced primordial density perturbations producing PBHs will also generate CMB distortions by dumping energy into the primordial plasma in the wave-number range $1 \lesssim k/\text{Mpc}^{-1} \lesssim 10^{4.5}$. The amount of distortion can be roughly estimated as $\langle\mu\rangle \sim O(1) P_{\zeta} \Delta N$, where $\Delta N$ is the e-fold number duration of the enhanced perturbations, which is equivalent to $\ln[k_{\text{max}}/k_{\text{min}}]$. Assuming $\Delta N \sim O(1)$, one can probe the distortion level from the power spectrum, and current bounds are $\langle\mu\rangle \lesssim 10^{-4}$ (Kohri et al. 2014; Nakama et al. 2018; Inomata and Nakama 2019; Byrnes et al. 2019; Kalaja et al. 2019; Gow et al. 2021; García-Bellido et al. 2017).

The PBH production is extremely sensitive to the tail of the statistical properties, i. e. probability distribution, of the primordial perturbations as discussed in Sect. 3. Hence, for highly efficient PBHs producing non-Gaussian perturbations, the distortion signal can be smaller and evade the current bounds. However, future experiments will reach remarkable sensitivities to probe both Gaussian and non-Gaussian primordial signals (Nakama et al. 2018; García-Bellido et al. 2017).

If scalar induced stochastic GW background is combined with CMB distortions, the probed PBH seed range will be $(0.1 - 10^{13})M_{\odot}$ (Ünal et al. 2021) (equivalently
<figure><img src="image_45.png" /><figcaption>**Fig. 40** Probing power of future PTA (SKA) and CMB distortion (PIXIE-like) experiments on primordial fluctuations for about 7 decades in wavenumbers and 13 decades in masses. Results are shown for inflationary perturbations that obey distinct probability distributions: Gaussian (red horizontal dashed line), chi-sqr ($\chi^2$, black horizontal dashed line), cubic-Gaussian ($G^3$, blue horizontal dashed line), detectable if</figcaption></figure>

the probed wavenumber range will be $1 \lesssim k/\text{Mpc}^{-1} \lesssim 10^7$). These probes will also reach conclusive results about the intriguing possibility that PBHs could be the seeds of SMBHs, and with the detection of billion-solar-mass SMBHs at redshift higher than 6 (Schneider et al. 2023), this has become an even more pressing question. Since CMB distortions and scalar induced GWs are global cosmic signals associated with the formation/production of PBHs, the conclusions will be robust to modifications in: (i) Astrophysical evolution, such as merger and accretion history; (ii) The amplitude and statistical properties of the primordial perturbations (Gaussian versus non-Gaussian); (iii) The clustering effects. The combined probes are shown in Fig. 40, where non-detection of both signals will robustly constrain the PBH abundance to the negligible amount of $f_{\text{PBH}} < 10^{-10}$ (Ünal et al. 2021).

Bosonic degrees of freedom scattering with spinning BHs can get energy from the rotational energy of BHs, and this phenomenon is called as superradiance. All BHs with enough spin can experience this process, but especially light PBHs, in the mass range, $10^{-33} - 1 M_{\odot}$, can experience superradiance even for tiny spin values, $10^{-5} - 10^{-2}$ (Ünal 2023). The instability rate of black hole is different for distinct types of perturbation, that is, scalar, vector, and tensor. Light PBHs can probe over 33 decades of bosons (scalar, vector and spin-2) in the mass range mentioned above, $10^{-12} - 10^{21}$ eV.<fnref n="15" /> This superradiance process boson mass range matches with the

<footnote n="15">¹⁵ Bosons could be beyond Standard model particles such as axions and dark photons, also it could be Standard Model photon with effective mass due to plasma interactions. If spinning light PBHs are observed, this implies non-existence of superradiance for the given BHs and boson mass range, therefore this observation can imply either there are no such particles or their self/external interactions prevent them from superradiance. This results in bounds on the axion decay constant (inversely related to the strength of self-interaction) and its energy density (Ünal 2023).</footnote>
plasma mass of photons, namely this allows light PBHs to amplify the photon production via superradiance and dump energy to CMB plasma, and if CMB cannot thermalize to have black-body spectrum this leads to CMB distortions, giving rise to bounds on PBHs depending on their masses (Pani and Loeb 2013).

## 7.5.2 Lyman-$\alpha$ forest

The constraints on PBHs can be obtained from the high-resolution high-redshift Lyman-$\alpha$ forest data since the Poisson fluctuations in the PBH number density induce a small-scale power enhancement departing from the standard CDM prediction. This enhancement can be interpreted as an isocurvature perturbation with the scale-invariant spectrum and isocurvature-to-adiabatic amplitude ratio for the monochromatic distribution: $f_{\text{iso}} = \sqrt{\frac{k_*^3 m_{\text{PBH}}/\text{PBH}}{2\pi^2 \Omega_{\text{CDM}} \rho_{\text{cr}} \frac{1}{A_s}}$ (here $k_* = 0.05$/Mpc is the pivot scale) (Murgia et al. 2019). Using MIKE and HIRES/KECK samples of quasar spectra and a grid of hydrodynamic simulations exploring different values of astrophysical parameters, in Murgia et al. (2019) the marginalized upper limit $f_{\text{PBH}} m_{\text{PBH}} \sim 60 M_{\odot}$ at $2\sigma$ on the PBH mass have been found for a monochromatic distribution. This has been obtained with a Gaussian prior on the reionization redshift, preventing the posterior distribution to peak on high redshift values that are disfavoured from the CMB and IGM observations.

## 7.5.3 El Gordo cluster and enhanced Halo Mass Functions

One of the most fascinating connections between PBHs and large scale structure (LSS) arises when one considers the effect of quantum diffusion in single field inflation with quasi-inflection points (Ezquiaga et al. 2020) as a natural way to generate large exponential tails in the density contrast PDF a fact that enhances the probability of collapse of small-scale fluctuations upon re-entry leading in this way to PBH formation. These highly non-Gaussian tails in the PDF can also be responsible at larger (galactic and cluster) scales for the enhancement of the halo mass function at high redshift (Ezquiaga et al. 2022), possibly explaining the increased abundance of massive galaxies seen by the James Webb Space Telescope at redshifs $z \sim 9-20$ (Finkelstein et al. 2022b), as well as significantly more massive clusters, like El Gordo (Asencio et al. 2020) at $z \sim 1$, than expected by the standard model $\Lambda$CDM. It remains to be shown whether such large deviation from the $\Lambda$CDM predictions for the halo mass function are compatible with LSS constraints. In the near future we may have deep imaging up to redshift $z \sim 20$ with both infrared and radio surveys like SKA, to test the enhanced exponential tails that could have given rise to PBHs in the early universe.

## 7.5.4 High-$z$ galaxies from JWST

The James Webb Space Telescope (JWST) imaging program via the Cosmic Evolution Early Release Science (CEERS) survey recently reported a population of massive galaxy candidates at redshift $z \ge 8$ with stellar masses of the order of