# Primordial black holes and their gravitational-wave signatures

Eleni Bagui, Sebastien Clesse, Valerio de Luca, Jose María Ezquiaga, Gabriele Franciolini, Juan García-Bellido, Cristian Joana, Rajeev Kumar Jain, Sachiko Kuroyanagi, Ilia Musco, et al.

► To cite this version:

Eleni Bagui, Sebastien Clesse, Valerio de Luca, Jose María Ezquiaga, Gabriele Franciolini, et al. Primordial black holes and their gravitational-wave signatures. Living Reviews in Relativity, 2025, 28 (1), pp.1. 10.1007/s41114-024-00053-w. hal-04284774

HAL Id: hal-04284774

https://hal.science/hal-04284774v1

Submitted on 10 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Primordial black holes and their gravitational-wave signatures

Eleni Bagui¹ · Sébastien Clesse¹ · Valerio De Luca² · Jose María Ezquiaga³ ·

Gabriele Franciolini⁴ · Juan García-Bellido⁵ · Cristian Joana⁶,⁷ ·

Rajeev Kumar Jain⁸ · Sachiko Kuroyanagi⁵,⁹ · Ilia Musco¹⁰,¹¹ ·

Theodoros Papanikolaou¹²,¹³,¹⁴ · Alvise Raccanelli⁴,¹⁵,¹⁶,¹⁷ ·

Sébastien Renaux-Petel¹⁸ · Antonio Riotto¹⁹ · Ester Ruiz Morales⁵,²⁰ ·

Marco Scalisi²¹ · Olga Sergijenko²²,²³,²⁴,²⁵ · Caner Ünal²⁶,²⁷,²⁸ ·

Vincent Vennin²⁹ · David Wands³⁰ · For the LISA Cosmology Working Group

Received: 22 December 2023 / Accepted: 5 November 2024

© The Author(s) 2024

## Abstract

In the recent years, primordial black holes (PBHs) have emerged as one of the most interesting and hotly debated topics in cosmology. Among other possibilities, PBHs could explain both some of the signals from binary black hole mergers observed in gravitational-wave detectors and an important component of the dark matter in the Universe. Significant progress has been achieved both on the theory side and from the point of view of observations, including new models and more accurate calculations of PBH formation, evolution, clustering, merger rates, as well as new astrophysical and cosmological probes. In this work, we review, analyze and combine the latest developments in order to perform end-to-end calculations of the various gravitational-wave signatures of PBHs. Different ways to distinguish PBHs from stellar black holes are emphasized. Finally, we discuss their detectability with LISA, the first planned gravitational-wave observatory in space.

**Keywords** Black holes · Primordial black holes · Cosmology · Gravitational waves · LISA

## List of symbols

$\zeta$ Curvature fluctuation

$\delta$ Density fluctuation

$P_{\zeta}$ Primordial curvature power spectrum (dimensionless)

$P_{\delta}$ Primordial density power spectrum (dimensionless)

$n_s$ Scalar spectral index of the power spectrum

<footnote>Extended author information available on the last page of the article</footnote>
<table>
  <tbody>
    <tr>
      <td>2.1</td>
      <td>Single-field inflationary models</td>
      <td>12</td>
    </tr>
    <tr>
      <td>2.1.1</td>
      <td>The basic idea and the slow-roll approximation</td>
      <td>12</td>
    </tr>
    <tr>
      <td>2.1.2</td>
      <td>Inflection point potentials</td>
      <td>14</td>
    </tr>
    <tr>
      <td>2.1.3</td>
      <td>Example model: critical Higgs inflation</td>
      <td>15</td>
    </tr>
    <tr>
      <td>2.1.4</td>
      <td>Reverse engineering approach</td>
      <td>18</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>Multi-field inflationary models</td>
      <td>19</td>
    </tr>
    <tr>
      <td>2.2.1</td>
      <td>Hybrid inflation</td>
      <td>20</td>
    </tr>
    <tr>
      <td>2.2.2</td>
      <td>Turns in multi-field inflation</td>
      <td>21</td>
    </tr>
    <tr>
      <td>2.2.3</td>
      <td>Two-stage models</td>
      <td>24</td>
    </tr>
    <tr>
      <td>2.2.4</td>
      <td>Axion-gauge scenario</td>
      <td>24</td>
    </tr>
    <tr>
      <td>2.2.5</td>
      <td>False-vacuum models</td>
      <td>26</td>
    </tr>
    <tr>
      <td>2.3</td>
      <td>Non-Gaussian models</td>
      <td>27</td>
    </tr>
    <tr>
      <td>2.3.1</td>
      <td>Motivations</td>
      <td>27</td>
    </tr>
    <tr>
      <td>2.3.2</td>
      <td>Quantum diffusion</td>
      <td>28</td>
    </tr>
    <tr>
      <td>2.3.3</td>
      <td>Stochastic spectator field</td>
      <td>30</td>
    </tr>
    <tr>
      <td>2.4</td>
      <td>Curvaton</td>
      <td>32</td>
    </tr>
    <tr>
      <td>2.5</td>
      <td>Preheating</td>
      <td>32</td>
    </tr>
    <tr>
      <td>2.6</td>
      <td>Early matter era</td>
      <td>36</td>
    </tr>
    <tr>
      <td>2.7</td>
      <td>Phase transitions</td>
      <td>36</td>
    </tr>
    <tr>
      <td>2.8</td>
      <td>Domain walls and cosmic strings</td>
      <td>38</td>
    </tr>
    <tr>
      <td>2.9</td>
      <td>Primordial magnetic fields</td>
      <td>38</td>
    </tr>
    <tr>
      <td>2.10</td>
      <td>Summary</td>
      <td>39</td>
    </tr>
    <tr>
      <td>3</td>
      <td>PBH formation, mass function and clustering</td>
      <td>39</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>From curvature perturbations to density contrast</td>
      <td>40</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Perturbation amplitude</td>
      <td>41</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Threshold for PBH formation</td>
      <td>43</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>Gravitational collapse of primordial curvature perturbations – standard formalism</td>
      <td>46</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>The inevitable non-Gaussianity of the primordial black hole abundance</td>
      <td>48</td>
    </tr>
    <tr>
      <td>3.6</td>
      <td>Collapse of some initial non-Gaussian perturbations</td>
      <td>49</td>
    </tr>
    <tr>
      <td>3.7</td>
      <td>Non-perturbative abundance from non-Gaussian perturbations</td>
      <td>49</td>
    </tr>
    <tr>
      <td>3.8</td>
      <td>Thermal history</td>
      <td>51</td>
    </tr>
    <tr>
      <td>3.8.1</td>
      <td>The equation of state and the number of relativistic degrees of freedom</td>
      <td>51</td>
    </tr>
    <tr>
      <td>3.8.2</td>
      <td>Thermal effects on PBH formation and abundance</td>
      <td>52</td>
    </tr>
    <tr>
      <td>3.9</td>
      <td>Evolution of the mass function through accretion</td>
      <td>57</td>
    </tr>
    <tr>
      <td>3.9.1</td>
      <td>Accretion onto isolated PBHs</td>
      <td>58</td>
    </tr>
    <tr>
      <td>3.9.2</td>
      <td>Accretion onto PBH binaries</td>
      <td>60</td>
    </tr>
    <tr>
      <td>3.10</td>
      <td>Spin distribution</td>
      <td>60</td>
    </tr>
    <tr>
      <td>3.11</td>
      <td>Comments on the clustering from Poisson initial conditions</td>
      <td>63</td>
    </tr>
    <tr>
      <td>3.12</td>
      <td>Summary</td>
      <td>65</td>
    </tr>
    <tr>
      <td>4</td>
      <td>PBH merging and encounter rates</td>
      <td>66</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Early binaries</td>
      <td>66</td>
    </tr>
    <tr>
      <td>4.1.1</td>
      <td>General rate formula</td>
      <td>66</td>
    </tr>
    <tr>
      <td>4.1.2</td>
      <td>Suppression from close PBHs and matter inhomogeneities</td>
      <td>67</td>
    </tr>
    <tr>
      <td>4.1.3</td>
      <td>Suppression from early Poisson-induced clustering</td>
      <td>70</td>
    </tr>
    <tr>
      <td>4.1.4</td>
      <td>Merging rates of perturbed binaries</td>
      <td>70</td>
    </tr>
    <tr>
      <td>4.1.5</td>
      <td>Link with observations of compact binary coalescences</td>
      <td>71</td>
    </tr>
    <tr>
      <td>4.1.6</td>
      <td>Limitations</td>
      <td>71</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Late binaries</td>
      <td>73</td>
    </tr>
    <tr>
      <td>4.2.1</td>
      <td>General formula</td>
      <td>73</td>
    </tr>
    <tr>
      <td>4.2.2</td>
      <td>Merging rates of binaries induced by three-body interactions</td>
      <td>74</td>
    </tr>
    <tr>
      <td>4.2.3</td>
      <td>Link with observations of compact binary coalescences</td>
      <td>75</td>
    </tr>
    <tr>
      <td>4.2.4</td>
      <td>Limitations</td>
      <td>76</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Hyperbolic encounters</td>
      <td>76</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>How to distinguish PBH from astrophysical binaries</td>
      <td>78</td>
    </tr>
    <tr>
      <td>4.5</td>
      <td>Summary</td>
      <td>80</td>
    </tr>
  </tbody>
</table>
<table><tbody><tr><td>5</td><td>Stochastic backgrounds</td><td>81</td></tr><tr><td>5.1</td><td>From second-order curvature fluctuations</td><td>81</td></tr><tr><td>5.2</td><td>Gauge invariance of the SGWB spectrum</td><td>86</td></tr><tr><td>5.3</td><td>Impact of primordial non-Gaussianity on scalar induced SGWB</td><td>87</td></tr><tr><td>5.4</td><td>Induced SGWB anisotropies</td><td>91</td></tr><tr><td>5.5</td><td>From PBH poisson fluctuations</td><td>92</td></tr><tr><td>5.6</td><td>From PBH mergers</td><td>95</td></tr><tr><td>5.6.1</td><td>Formulation</td><td>95</td></tr><tr><td>5.6.2</td><td>Early binaries</td><td>96</td></tr><tr><td>5.6.3</td><td>Late binaries in clusters</td><td>98</td></tr><tr><td>5.7</td><td>SGWB duty cycle</td><td>100</td></tr><tr><td>5.8</td><td>From close encounters</td><td>100</td></tr><tr><td>5.9</td><td>Summary</td><td>102</td></tr><tr><td>6</td><td>Gravitational-waves and large scale structure correlations</td><td>103</td></tr><tr><td>6.1</td><td>Introduction</td><td>103</td></tr><tr><td>6.2</td><td>GW-galaxy cross-correlations for resolved events</td><td>104</td></tr><tr><td>6.3</td><td>GW-galaxy cross-correlations for the stochastic background</td><td>105</td></tr><tr><td>6.4</td><td>GW × LSS forecasts</td><td>106</td></tr><tr><td>6.5</td><td>Summary</td><td>110</td></tr><tr><td>7</td><td>Current limits</td><td>111</td></tr><tr><td>7.1</td><td>PBH evaporation</td><td>113</td></tr><tr><td>7.1.1</td><td>Big-bang nucleosynthesis (BBN)</td><td>113</td></tr><tr><td>7.1.2</td><td>Gamma-rays and neutrinos from Hawking radiation</td><td>114</td></tr><tr><td>7.1.3</td><td>Positron annihilation in the Galactic Center</td><td>115</td></tr><tr><td>7.1.4</td><td>e<sup>±</sup> observations by Voyager 1</td><td>115</td></tr><tr><td>7.1.5</td><td>Ultra-faint dwarf galaxies</td><td>115</td></tr><tr><td>7.1.6</td><td>CMB limits on Dark radiation</td><td>116</td></tr><tr><td>7.1.7</td><td>Neutrinos</td><td>116</td></tr><tr><td>7.1.8</td><td>21 cm signal</td><td>117</td></tr><tr><td>7.1.9</td><td>Limits from DM particle production</td><td>118</td></tr><tr><td>7.1.10</td><td>Discussion and limitations</td><td>118</td></tr><tr><td>7.2</td><td>Microlensing searches</td><td>119</td></tr><tr><td>7.2.1</td><td>Microlensing from stars in the Magellanic clouds</td><td>120</td></tr><tr><td>7.2.2</td><td>Microlensing from stars in the galactic bulge</td><td>120</td></tr><tr><td>7.2.3</td><td>Microlensing from stars in Andromeda</td><td>121</td></tr><tr><td>7.2.4</td><td>Microlensing from quasars</td><td>121</td></tr><tr><td>7.2.5</td><td>Caustic crossing in giant arcs</td><td>121</td></tr><tr><td>7.2.6</td><td>Microlensing of supernovae</td><td>122</td></tr><tr><td>7.2.7</td><td>Femtolensing</td><td>122</td></tr><tr><td>7.2.8</td><td>Discussion and limitations</td><td>122</td></tr><tr><td>7.3</td><td>Dynamical limits</td><td>123</td></tr><tr><td>7.3.1</td><td>PBH capture by neutron stars or white dwarfs</td><td>123</td></tr><tr><td>7.3.2</td><td>Ultra-faint dwarf galaxies</td><td>124</td></tr><tr><td>7.3.3</td><td>Dark Matter profile of dwarf galaxies (core-cusp problem)</td><td>125</td></tr><tr><td>7.3.4</td><td>Wide halo binaries</td><td>125</td></tr><tr><td>7.3.5</td><td>Disruption of stellar streams</td><td>125</td></tr><tr><td>7.3.6</td><td>PBHs in the solar system</td><td>126</td></tr><tr><td>7.3.7</td><td>Discussion and limitations</td><td>126</td></tr><tr><td>7.4</td><td>Accretion limits</td><td>127</td></tr><tr><td>7.4.1</td><td>Cosmic microwave background (CMB) anisotropies</td><td>127</td></tr><tr><td>7.4.2</td><td>21 cm signal</td><td>130</td></tr><tr><td>7.4.3</td><td>Cosmic X-ray, infrared and radio backgrounds</td><td>131</td></tr><tr><td>7.4.4</td><td>Limits from accretion at present epoch</td><td>131</td></tr><tr><td>7.4.5</td><td>Discussion and limitations</td><td>132</td></tr><tr><td>7.5</td><td>Indirect constraints from density fluctuations</td><td>132</td></tr></tbody></table>
<table><tbody><tr><td>7.5.1</td><td>CMB distortions</td><td>132</td></tr><tr><td>7.5.2</td><td>Lyman-α forest</td><td>134</td></tr><tr><td>7.5.3</td><td>El Gordo cluster and enhanced Halo Mass Functions</td><td>134</td></tr><tr><td>7.5.4</td><td>High-z galaxies from JWST</td><td>134</td></tr><tr><td>7.5.5</td><td>Ultra-compact mini haloes</td><td>135</td></tr><tr><td>7.5.6</td><td>Discussion and limitations</td><td>137</td></tr><tr><td>7.6</td><td>Gravitational waves</td><td>137</td></tr><tr><td>7.6.1</td><td>LVK GW mergers</td><td>137</td></tr><tr><td>7.6.2</td><td>LVK sub-solar searches</td><td>140</td></tr><tr><td>7.6.3</td><td>LVK SGWB</td><td>140</td></tr><tr><td>7.6.4</td><td>PTAs and NANOGrav</td><td>141</td></tr><tr><td>7.6.5</td><td>Continuous and high-frequency GWs from planetary-mass PBHs</td><td>144</td></tr><tr><td>7.6.6</td><td>Discussion and limitations</td><td>144</td></tr><tr><td>7.7</td><td>Summary</td><td>145</td></tr><tr><td>8</td><td>Detectability with LISA</td><td>146</td></tr><tr><td>8.1</td><td>SGWB from second order curvature fluctuations</td><td>146</td></tr><tr><td>8.2</td><td>SGWB from ultralight evaporated PBHs</td><td>147</td></tr><tr><td>8.3</td><td>SGWB from early binaries</td><td>147</td></tr><tr><td>8.4</td><td>SGWB from late binaries</td><td>148</td></tr><tr><td>8.5</td><td>Intermediate-mass binary mergers</td><td>148</td></tr><tr><td>8.6</td><td>Extreme mass ratio inspirals</td><td>148</td></tr><tr><td>8.7</td><td>High-redshift binary mergers</td><td>149</td></tr><tr><td>8.8</td><td>Combination with LVK observations</td><td>149</td></tr><tr><td>8.9</td><td>Quasi-monochromatic continuous waves from subsolar PBHs</td><td>150</td></tr><tr><td>8.10</td><td>Observing near-Earth asteroid-mass primordial black holes</td><td>151</td></tr><tr><td>8.11</td><td>Summary</td><td>152</td></tr><tr><td>9</td><td>Conclusions</td><td>153</td></tr><tr><td></td><td>References</td><td>156</td></tr></tbody></table>

# 1 Introduction

## 1.1 A brief history of primordial black holes (PBHs)

The idea that black holes could have been formed in the early Universe dates back to the late 1960s with the pioneering work of Zel'dovich and Novikov (1967) and to the early 1970s with the work of Hawking (1971). This triggered Hawking's famous discovery of black hole evaporation (Hawking 1975). Subsequently, Carr, Hawking's PhD student at the time, continued to investigate PBHs (Carr and Hawking 1974; Carr 1975). Already in Hawking (1971), Carr and Hawking (1974), Chapline (1975), it was suggested that such primordial black holes (PBHs) could contribute to the suspected Dark Matter (DM) in the Universe or to the seeds of supermassive black holes. The first formation scenarios in the context of inflation were proposed in the 1990s (Dolgov and Silk 1993; Carr and Lidsey 1993; Carr et al. 1994; Ivanov et al. 1994; García-Bellido et al. 1996; Kim and Lee 1996; Kawasaki et al. 1998; Green and Liddle 1997; Ivanov 1997; Yokoyama 1998; Kotok and Naselsky 1998), but these usually led to (evaporating) PBHs of very small mass, except for García-Bellido et al. (1996) who predicted solar mass PBHs that could account for the MACHO microlensing events observed towards the Magellanic clouds (Aubourg et al. 1993; Alcock et al. 1997). Other mechanisms were also proposed, e.g., based
on phase transitions (Jedamzik 1997; Niemeyer and Jedamzik 1998), an early matter era (Khlopov and Polnarev 1980; Polnarev and Khlopov 1985; Green et al. 1997), scalar field instabilities (Khlopov et al. 1985), collapse of topological defects (Polnarev and Zembowicz 1991), modified gravity (Barrow and Carr 1996; Kawai and Kim 2021; Zhang 2022), string theory (Cicoli et al. 2018; Nanopoulos et al. 2020; Mavromatos et al. 2022), bouncing cosmological scenarios (Copeland et al. 1998; Carr and Coley 2011; Quintin and Brandenberger 2016; Chen et al. 2017; Banerjee et al. 2022; Papanikolaou et al. 2024a) as well as quantum gravity setups (Papanikolaou 2023a). In the late 1990s, stellar-mass PBHs were seriously considered as a dark matter candidate, following the possible detection in the MACHO survey (Aubourg et al. 1993; Alcock et al. 1997) of several microlensing events towards the Magellanic clouds. However, the EROS (Tisserand et al. 2007) and OGLE (Wyrzinski et al. 2010, 2011a, b; Calchi Novati et al. 2013) surveys later set more stringent limits on the PBH abundance, and at the same time, very stringent constraints from cosmic microwave background (CMB) observations were claimed in Ricotti et al. (2008).

Despite some pioneering numerical studies (Novikov and Polnarev 1978, 1980; Bicknell and Henriksen 1979), only more recently it has become possible to fully understand the mechanism of PBH formation with detailed spherically-symmetric numerical simulations (Shibata and Sasaki 1999; Musco et al. 2005; Polnarev and Musco 2007), showing that superhorizon cosmological perturbations would collapse to PBHs after re-entering the cosmological horizon, if their amplitude $\delta$ is larger than a certain threshold value $\delta_c$. This quantity, measured at horizon crossing, was initially estimated with a simplified Jeans length argument in Newtonian gravity (Carr 1975), giving $\delta_c \sim c_s^2$, where $c_s^2 = 1/3$ is the sound speed of the cosmological radiation fluid measured in units of the speed of light. More recently this argument has been generalized within general relativity using a three-zone model (Harada et al. 2013; Papanikolaou 2022b) which gives $\delta_c \simeq 0.4$ for a radiation-dominated Universe. Numerical simulations showed that the mechanism of critical collapse characterizing perfect fluids (Neilsen and Choptuik 2000) arises naturally in the context of PBH formation (Yokoyama 1998; Niemeyer and Jedamzik 1998; Jedamzik and Niemeyer 1999; Musco et al. 2005). In particular, the collapse is characterized by a relativistic wind which progressively separates the collapsing perturbation from the expanding background, approaching a self-similar critical solution characterizing the threshold (Musco et al. 2009; Musco and Miller 2013). Finally, a very recent detailed study has found a clear relation between the value of the threshold, $\delta_c$, and the initial configuration of the curvature (or energy density) profile used to set the initial condition of the numerical simulations, with $0.4 \le \delta_c \le 2/3$, where the shape is identified by a single parameter (Musco 2019; Escrivà et al. 2020). This range is reduced to $0.4 \le \delta_c \lesssim 0.6$ when the initial perturbations are computed, using peak theory (Bardeen et al. 1986), from the primordial power spectrum of cosmological perturbations (Musco et al. 2021).

Since 2016, the real game-changer that has rekindled the idea that PBHs may exist and constitute a significant fraction to the total dark matter (Bird et al. 2016; Clesse and García-Bellido 2017b; Sasaki et al. 2016) has been the first gravitational-wave
(GW) detection from a black hole merger by LIGO/Virgo (Abbott et al. 2016). Initially, the merger rate of early PBH binaries seemed to only allow a small fraction of the dark matter to be made of PBHs (Sasaki et al. 2016), but N-body simulations later showed that these rates are in fact significantly suppressed (Raidal et al. 2019; Trashorras et al. 2021). Nowadays, the importance of the different PBH binary formation channels, the possible abundance of PBHs, their viable mass function, are a subject of intense activity and are hotly debated (for recent reviews, see e.g. Carr and Kuhnel 2020; Carr et al. 2021c). Recent analysis based on the rate and mass (Clesse and García-Bellido 2018; Raidal et al. 2019; Carr et al. 2021b; Hall et al. 2020; Jedamzik 2020, 2021; Clesse and García-Bellido 2022; Hütsi et al. 2021c; Escrivà et al. 2023), or spin (Fernandez and Profumo 2019; De Luca et al. 2020e; García-Bellido et al. 2021c), or both (De Luca et al. 2020d, 2021c; Wong et al. 2021; Franciolini 2021; Franciolini et al. 2022a, b, c; Franciolini and Pani 2022) distributions of compact binary coalescences observed in the first, second (Abbott et al. 2019b) and third (Abbott et al. 2021a, 2023a) observing runs of the LIGO/Virgo/Kagra collaboration (LVK), suggest that some black holes may be primordial. Some are based on a Bayesian approach (Clesse and García-Bellido 2018; Raidal et al. 2019; Hall et al. 2020; Hütsi et al. 2021c; Wong et al. 2021; De Luca et al. 2021a, c; Kritos et al. 2021; Franciolini et al. 2022a, e; Franciolini 2021; Fernandez and Profumo 2019; García-Bellido et al. 2021c) that will ultimately allow us to compare statistically the different possible origins and scenarios.

Most of population analysis studies have concluded, for the considered scenarios, that only a fraction of dark matter of order of $10^{-3}$ made of PBHs in the LIGO/Virgo black hole mass range is allowed. While there is a rather general consensus among the community that PBHs of more than ten solar masses cannot explain the dark matter, the exact constraints on their abundance from GW obsevations at the solar-mass scale – motivated by the expected softening of the collapse threshold at the QCD transition – is still debated. On the one hand, the merger rates inferred from LIGO/Virgo are higher and still uncertain. On the other hand, the theoretical predictions for PBH merger rates are quite model dependent. There are also several sources of uncertainties in the rate prescriptions for broad PBH mass distributions, notably the importance of binary disruptions in the early Universe. Therefore, due to the broad variety of PBHs and astrophysical models as well as the large theoretical and observational uncertainties, it is still premature to firmly affirm (or deny) a PBH origin of GW observations or to relate them to a particular value of the PBH dark matter fraction. In this context, the next generation of GW detectors like Einstein Telescope (Maggiore et al. 2020), Cosmic Explorer (Reitze et al. 2019) and the Laser Interferometer Space Antenna (LISA) (Amaro-Seoane et al. 2017; Colpi et al. 2024) will play a crucial role. Also, the development of ultra-high frequency GW detectors may additionally complement the search for GW signatures of asteroid-mass PBHs (Herman et al. 2023; Franciolini et al. 2022d) (see also e.g. Aggarwal et al. 2021; Ringwald et al. 2021; Berlin et al. 2022a, b; Domcke et al. 2022).

Important progress has also been made on the limits on the allowed PBH abundance set by various types of observations. Shortly after the first GW detection, the CMB limits were re-analysed and found to be much less stringent than initially
thought (Ali-Haïmoud and Kamionkowski 2017), thereby re-opening the stellar mass region. At the same time, the degree of validity of the microlensing limits has been questioned (Hawkins 2015; Clesse and García-Bellido 2015; García-Bellido and Clesse 2018; Green 2017; Calcino et al. 2018; Carr et al. 2021b) for more realistic scenarios than monochromatic and homogeneous distributions of PBHs. New probes, like ultra-faint-dwarf galaxies (Brandt 2016; Green 2016; Li et al. 2017) and X-ray sources towards the galactic center (Gaggero et al. 2017), have revived the mass region between 1 and $100 M_\odot$, while the constraining power of other probes like high-cadence microlensing (Niikura et al. 2019b), neutron stars (Capela et al. 2013a) and white dwarves (Capela et al. 2013b), have been reduced, which has reopened the asteroid-mass region (Smyth et al. 2020). Finally, a series of recent microlensing events from OGLE (Niikura et al. 2019a; Wyrzykowski and Mandel 2020), HSC (Niikura et al. 2019b), as well as quasars (Hawkins 1993, 2020a, b, 2022), may hint at the existence of an important population of planetary-mass and solar-mass dark compact objects. The unexplained spatial correlations between infrared and X-ray backgrounds at high redshift (Kashlinsky et al. 2005) could also be due to an important population of stellar-mass PBHs (Kashlinsky 2016; Kashlinsky et al. 2018; Cappelluti et al. 2022).

Furthermore, since PBHs are formed by the collapse of large density fluctuations, there is an associated stochastic GW background (SGWB) sourced by these perturbations at second-order. It has been calculated (Ando et al. 2018; García-Bellido et al. 2017) that if BHs detected by LVK have a primordial origin, there is an inevitable accompanying SGWB peaking around pulsar-timing-array frequencies. In September 2020, NANOGrav claimed the possible detection of a SGWB at nanohertz frequency using pulsar timing arrays (Arzoumanian et al. 2020), which may have been sourced by primordial density perturbations which could be the origin of stellar-mass PBH formation (Vaskonen and Veermäe 2021; De Luca et al. 2021g; Kohri and Terada 2021). Remarkably, such claim was confirmed by more recent pulsar timing array data released in 2023 by the NANOGrav (Agazie et al. 2023b, c), EPTA (in combination with InPTA) (Antoniadis et al. 2023a, b, 2024), PPTA (Reardon et al. 2023a, b; Zic et al. 2023) and CPTA (Xu et al. 2023) collaborations, who found evidence for a Hellings and Downs angular correlation, typical of an homogeneous spin-2 GW background and consistent with the quadrupolar nature of GWs in GR (Hellings and Downs 1983).

On the theory side, a plethora of new models have been proposed in the last four years. Most of them are still subject to fine-tuning issues related to the amplitude of fluctuations required to form PBHs with a significant abundance. Nevertheless, one should explore a new class of models that do not require a modification of the primordial power spectrum on the cosmological scale. That is, one in which the PBH formation is related to an exponential, rather than Gaussian, tail of the distribution of primordial perturbations on small scales (Ezquiaga et al. 2020; Carr et al. 2021a), which do not affect the observations on CMB scales. In these models, a large population of PBHs could arise much more naturally. PBH masses of order of the solar mass and around $30 M_\odot$ could also arise naturally due to the thermal history of the Universe and the change in the equation of state at the QCD epoch (Jedamzik 1997; Niemeyer and Jedamzik 1998; Byrnes et al. 2018; Carr et al. 2021b;
Franciolini et al. 2022e; Escrivà et al. 2023), which should strongly boost PBH
formation. Finally, besides explaining the dark matter or a significant fraction of it,
PBHs could also be related to baryogenesis and different mechanisms have been
proposed (Dolgov et al. 2000; Baumann et al. 2007b; García-Bellido et al. 2021a;
Carr et al. 2021a; García-Bellido 2019; De Luca et al. 2021b).

Any firm detection would open a new window on the physics at play in the very early Universe and a possible way to solve various long-standing astrophysical and cosmological puzzles (Clesse and García-Bellido 2018; Carr et al. 2021b). Ultimately, the best way to distinguish PBHs from stellar black holes would be to detect sub-solar-mass black holes. Recently, four sub-solar black hole triggers have been found in the second observing run of LVK (Phukon et al. 2021; Morras et al. 2023), and three in the third observing run (Abbott et al. 2023c). These are not statistically significant enough to claim a firm detection, but if confirmed, these would strongly support a PBH origin of some LVK black holes and make a significant contribution to the dark matter in the Universe.

## 1.2 Why a(nother) review on PBHs?

There already exist a few recent and relatively complete review articles dedicated to PBHs, see e.g. Sasaki et al. (2018), García-Bellido (2017), Carr and Kuhnel (2020), Carr et al. (2021c, 2024). Earlier reviews that pre-date the first GW detection also exist, see e.g. Khlopov (2010), Carr et al. (2010). All of them are written by only a few authors. It has therefore been quite challenging for them to include detailed discussions and analysis of the broad range of topics covered by PBHs. Indeed, these include early Universe cosmology, numerical relativity, dark matter, astrophysics, celestial mechanics, GW astrophysics, etc. Those reviews therefore summarize some of the results obtained in these topics but do not include model specific A-to-Z calculations of PBH related observables. Very often, some of the most recent developments like the precise determination of the critical threshold density leading to PBH formation, which is highly model-dependent and varies with the evolution of the equation of state of the Universe, or the rate suppression/boost linked to Poisson clustering, are not always considered when models are compared to observations.

The goal of this paper is to review these most recent developments and discuss
how they impact the GW signatures of PBHs. Compared to other recent reviews, we
aim at integrating together the most recent and accurate models, e.g., of PBH
formation, evolution, clustering, merger rates, in order to compute the GW signatures
of PBHs, and to explore their detectability. For this purpose, we have developed in
parallel a numerical code that implements all those models and recent scenarios. The
code will be released soon in a separate publication, but it has been already used in
this paper in order to produce the key figures. The task of including these
increasingly complex models has been eased by the specific expertise of many
authors, who have actively contributed to these developments, covering a broad
range of topics.

Finally, another original aspect of this work is the fact that we plan to have a *living*
review format, with bi-yearly updates of the present manuscript, in order to track the
rapid progress in this very active field of research and provide at all times an up-to-
date reference publication. In future, this review could be expanded to include new sections on the detectability by PTAs and ground-based GW detectors, or on the cosmological signatures of PBH scenarios.

All of this encourages us to go beyond the present state, and provide a new kind of review on PBHs.

## 1.3 Probing PBHs with LISA

This project is hosted by the cosmology working group of the LISA consortium. Therefore our discussion has been focused on the GW signatures of PBHs in the particular context of the LISA mission.

PBHs cover a very wide range of masses and therefore frequency ranges of GWs. In particular, LISA will be sensitive to a broad band around the millihertz GW frequencies (García-Bellido et al. 2017; Cai et al. 2019a; Bartolo et al. 2019b; Ünal 2019), which is complementary to the ground-based GW detectors (Maggiore et al. 2020; Reitze et al. 2019) and electromagnetic probes (Kashlinsky et al. 2019). One of the positive prospects for LISA is the measurement of a SGWB spectrum which, if proven to have a cosmological origin, could provide some information on the physics of the early universe. One channel for example involves GWs from PBH formation in the radiation era. While it would be much more involved to robustly claim a connection to PBH formation, the search for a SGWB from PBH formation could allow to constrain their existence in different mass ranges, for example in the interesting asteroid-mass range around $10^{-12} M_{\odot}$ (Saito and Yokoyama 2009; García-Bellido et al. 2017; Cai et al. 2019a; Bartolo et al. 2019b; Ünal 2019). Moreover, the coalescence of the heavy seeds of supermassive black holes (SMBH) in the late Universe could leave clear signatures of their primordial origin, or even detect individual events at high redshift which could not have arisen from astrophysical mechanisms.

The great sensitivity of LISA at mHz frequencies (Baibhav et al. 2021; Amaro-Seoane et al. 2023; Arun et al. 2022; Karnesis et al. 2024) opens the possibility to detect the mergers of $10^3 - 10^4 M_{\odot}$ PBHs all the way to $z \simeq 100$. Furthermore, the isotropic SGWB from the coalescence of PBHs since recombination should have an amplitude and spectral shape that will make it easily detectable by LISA. In this review, we will describe the different features of known PBH scenarios that can be probed with LISA.

LISA will not only probe the GWs from individual PBH mergers and we will review numerous other types of signatures of PBHs that can be probed with LISA. Among them, the GW background from second order curvature fluctuations will probe the very interesting asteroid-mass range. Other backgrounds include the one from ultralight, evaporated PBHs, the ones from early and late PBH binaries as well as from hyperbolic encounters. GW from individual sources could also be observed by LISA, for instance from intermediate-mass PBH mergers, including at high-redshifts, as above mentioned, or from inspirals with extreme mass ratios which are expected for extended PBH mass distributions. LISA will also have a key role if data are combined to LVK or other detector observations at higher frequencies, for long-
duration signals including mergers with a subsolar component. Subsolar-mass PBHs could also produce detectable quasi-monochromatic continuous waves in the LISA frequency range. Finally, a more exotic possibility is to use LISA as a detector of near-Earth asteroid-mass PBHs.

## 1.4 Outline of this review

This review is organized as follows. In Sect. 2, we present the different classes of theoretical PBH formation scenarios, with a focus on recent inflation models but including a discussion of other production mechanisms such as (p)reheating, curvatons, phase transitions, topological defects and primordial magnetic fields. Among inflationary scenarios, we especially focus on recent proposals in which the usual Gaussian distribution of primordial curvature fluctuations is modified in a more complex way than can be described in terms of conventional non-Gaussian deviation parameters ($f_{NL}$, etc). Ultimately this could be one possible way to resolve the long-standing fine-tuning issue for primordial fluctuations as the origin of PBHs (Cole et al. 2023b).

Section 3 is dedicated to PBH formation and the computation of the PBH mass distribution. We start with the standard formation formalism but then we summarize and include in our calculations the results of the most recent and accurate studies based on simulations in numerical relativity, such as non-linear and non-Gaussian effects, the impact of the shape of the spectrum and the curvature/overdensity profiles, the changes in the equation of state linked to the thermal history, the mass and spin evolution after formation, etc.

Once the PBH mass distribution is known, the next important step is the calculation of PBH merger rates, covered in Sect. 4. PBH mergers are due to two dominant channels: early binaries that are formed before matter-radiation equality, and late binaries formed in PBH clusters. We use the latest prescriptions, e.g., based on N-body simulations, in order to derive merger rates for different mass scales and for some representative models, including for very low mass ratios and intermediate-mass BH mergers. Finally, we consider the case of hyperbolic encounters in PBH clusters that can also lead to a large number of GW burst events.

In Sect. 5, we review the different sources of stochastic backgrounds related to PBHs. One particularly relevant source of SGWB comes from the curvature fluctuations that source GWs at second order in perturbation theory. Other SGWBs are due to early or late PBH binaries as well as due to PBH isocurvature fluctuations.

In Sect. 6, we review the possible ways to use the cross-correlations between gravitational wave observations of individual mergers or the stochastic background, and large-scale structure (LSS), in order to constrain the existence and abundance of PBHs.

The current limits on the abundance of PBHs from various probes, including GWs, are reviewed and discussed in Sect. 7. In particular, we aim to present them in a model-dependent way. We emphasize and comment on the underlying hypotheses and discuss how these limits change according to the different classes of models.

Finally, the detectability of various GW signatures from PBHs with LISA is discussed in Sect. 8, including individual merger events with extreme mass ratios or
intermediate-masses, hyperbolic encounters, stochastic backgrounds, continuous waves, high-redshift signals, signal combination with other ground-based detectors, etc.

At the end of each section, we provide an augmented discussion of recent results, the current limitations and pave the way for future analyses by listing some interesting perspectives. Our conclusions are presented in Sect. 9.

## 2 Theoretical models

In this section, we provide an overview of the principal mechanisms that can lead to large curvature fluctuations and PBH formation. We start with inflationary scenarios (see Özsoy and Tasinato 2023 for a recent review), either in single-field (quasi-inflection point) or multi-field realisations (hybrid inflation, turns in field space, gauge field interactions, etc.). Where possible, we provide some general formulae for the calculation of the primordial power spectrum. Then we consider recent models relying on quantum diffusion of the inflaton or of a stochastic spectator field, generically leading to non-Gaussian distributions of the curvature fluctuations. The third part of this section is dedicated to other models of PBH formation, due to curvaton fields, preheating, phase transitions, early matter era, cosmic strings, domain walls and primordial magnetic fields. When PBH form due to primordial overdensities, one typically requires order one fluctuations to form PBHs due to their gravitational collapse, when they re-enter inside the Hubble horizon. One can therefore relate the PBH mass $m_{PBH}$, being of the order of the mass in the Hubble horizon, to the fluctuation wavenumber $k$, the formation redshift and the Universe temperature at that time, as sketched in Fig. 1.

### 2.1 Single-field inflationary models

In the following, we provide a review of the main production mechanisms of primordial black holes in the context of single-field inflationary models.

#### 2.1.1 The basic idea and the slow-roll approximation

The inflationary epoch provides an ideal setting for producing PBHs in the very early Universe, with single-field inflation being the minimal framework. While a nearly scale-invariant spectrum of primordial curvature perturbations is remarkably consistent with the observations of CMB temperature anisotropies at large scales, the spectrum on smaller scales remains largely unconstrained. When small-scale perturbations with sufficiently large overdensities re-enter the expanding horizon during the radiation-dominated epoch, their gradients induce a gravitational collapse that cannot be overcome even by the radiation pressure of the expanding plasma, thereby producing PBHs with a mass of the order of the horizon mass at the time of re-entry of a given wavenumber, $k$. The abundance of PBHs crucially depends on the collapse process and the nature of primordial curvature perturbations, which are discussed later. If PBHs start evaporating due to Hawking radiation right after their
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Relation between the PBH mass (here being set equal to the mass in the Hubble horizon), the wave-number $k$ associated to the primordial overdensity at the origin of the PBH formation, the Universe temperature and the redshift when this mode re-enters inside the Hubble radius. Image reproduced from Clesse et al. (2018) Image reproduced with permission</figcaption></figure>

production, one finds that PBHs with mass $m_{\text{PBH}} \gtrsim 10^{-18} M_{\odot}$ will survive until today and dominate the expansion of the Universe at matter-radiation equality unless their initial abundance is strongly suppressed.

Since single-field slow-roll inflationary models typically predict a nearly scale-invariant power spectrum of curvature perturbations on all scales, most slow-roll single-field inflation models do not produce PBHs, because a significant abundance of PBHs requires a substantial growth of the power spectrum on small scales. Nevertheless, a large positive running of the scalar primordial power spectrum can lead to PBH formation (Carr et al. 1994). Another possibility is to have a transition in the slope of the potential, with the scalar field slowly rolling towards a region of almost constant potential, as one gets for instance in the original hybrid inflation model along the valley. But the difficulty is then to end inflation in a limited number of e-folds while, in order to produce a large power, the potential is extremely flat. Another possibility for producing PBHs in single-field inflation is through the violation of the usual slow-roll conditions (Ivanov et al. 1994; García-Bellido and Ruiz Morales 2017; Motohashi and Hu 2017; Inomata et al. 2021b, 2022, 2023), although one should then pay a special attention to perturbative control (Inomata

<figure><img src="image_4.png" /><figcaption>**Fig. 2** A schematic plot of the inflationary potential with an inflection point. Image adapted from García-Bellido (2017) Image reproduced with permission</figcaption></figure>
et al. 2023). Moreover, not all single-field inflationary models can produce the relevant abundance of PBHs as it crucially depends on the shape of the potential and various model parameters. This typically leads to a high degree of fine-tuning of model parameters.

The power spectrum of curvature perturbations in single-field models is given at first order and in the slow-roll approximation, by

$$P_{\zeta}(k) = \frac{H_*^2}{\pi m_{\text{Pl}}^2 \epsilon_{1*}} \left[ 1 - 2(C+2)\epsilon_{1*} + C\epsilon_{2*} - (2\epsilon_{1*} + \epsilon_{2*}) \ln\left(\frac{k}{k_*}\right) \right], \quad (1)$$

where $C = 0.578 + 2 \ln 2 - 2$, $\epsilon_1 \equiv -\dot{H}/H^2$ is the first slow-roll Hubble-flow parameter, $\epsilon_2 \equiv d \ln \epsilon_1 / dN$ is the second Hubble-flow parameter, $m_{\text{Pl}}$ is the Planck mass (and $M_{\text{Pl}} = \sqrt{8\pi} m_{\text{Pl}}$ denotes the reduced Planck mass) and a star denotes a quantity evaluated at the time the comoving wavenumber, $k_*$, crosses the Hubble radius during inflation when $k_* = a_* H_*$. Therefore, one can naively expect that an enhancement of $P_\zeta$ can occur if there is a dynamical phase during the inflationary expansion in which $\epsilon_1$ becomes much smaller than unity. Such an epoch is often called an *ultra slow-roll* (USR) phase which can be achieved by means of a plateau region in the inflationary potential, which slows down the inflaton field more rapidly than in slow roll, before ending inflation. An efficient and interesting way to produce this plateau is with the introduction of an inflection point in the potential. However, one must ensure that the inflaton does not stay too long at the inflection point otherwise all the inflationary fluctuations that had successfully imprinted the metric perturbations on large scales to explain the CMB anisotropies will be expanded away. A schematic plot of the potential with an inflection point is shown in Fig. 2.

### 2.1.2 Inflection point potentials

The simplest way to achieve an inflection or a near-inflection point in the potential is to use a suitable polynomial potential or ratio of polynomials so as to obtain the desired nearly scale invariant spectrum on CMB scales (García-Bellido and Ruiz Morales 2017). A very generic feature of such a potential is that although the inflaton starts in a slow-roll regime, it dynamically enters into an USR stage close to the inflection point for a short duration of e-folds that induces an enhancement in the spectrum. Nevertheless, it is possible to evade the USR regime in quasi-inflection point potentials (Ezquiaga et al. 2018). However, one must always be careful that this new feature in the potential does not affect the large-scale dynamics which is already strongly constrained by CMB and large-scale structure (LSS). These types of potentials arise naturally in models of Higgs inflation (Bezrukov and Shaposhnikov 2008; Ezquiaga et al. 2018), MSSM inflation (Allahverdi et al. 2006, 2007), accidental inflation (Linde and Westphal 2008), string inflation (Cicoli et al. 2018; Özsoy et al. 2018), α-attractor models (Dalianis et al. 2019; Iacconi et al. 2022; Kallosh and Linde 2022), punctuated inflation (Jain et al. 2009, 2010), non-canonical inflation (Kamenshchik et al. 2019; Lin et al. 2020; Papanikolaou et al. 2023a), as an effective description of multifield inflation (Geller et al. 2022) and in more general inflection point inflation (Hotchkiss et al. 2011). However, to produce sufficient
PBHs in these models, one crucial difference is that, rather than imposing the initial conditions on the inflaton evolution close to the inflection point, one must assume that the field starts in the slow-roll regime at larger values of the potential well above the inflection point, and then slow-rolls down towards the minimum of the potential, crossing the inflection point, typically 30 to 40 e-folds before the end of inflation. Such a dynamical evolution must ensure that large-scale dynamics are consistent with the CMB constraints, PBHs are produced at some intermediate smaller scales during the USR regime and inflation naturally ends when the field finally rolls down towards the true minimum of the potential. In this context, note that there has been an intense discussion recently on loop corrections induced on CMB scales by enhanced scalar perturbations on PBH scales (Kristiano and Yokoyama 2024a; Riotto 2023a, b; Kristiano and Yokoyama 2024b; Firouzjahi 2023; Firouzjahi and Riotto 2024; Franciolini et al. 2024a; Tasinato 2023; Cheng et al. 2024). Hints for significant one-loop contribution, to the point of undermining the consistency of such set-ups, have appeared in computations including a subset of vertices as derived from the interaction Hamiltonian. However, it has been recently shown that these effects disappear due to cancellations once all relevant interaction terms, at least at cubic order, are taken into account (Fumagalli 2023), and this subject will likely deserve further studies. Note also that enhanced fluctuations can truly lead, through loop effects, to an increase of power on (intermediate) larger scales due to infrared rescattering (Fumagalli et al. 2024), in single-field models and beyond, with potentially important consequences for the PBH population.

2.1.3 Example model: critical Higgs inflation

Critical Higgs Inflation (CHI) is probably among the best physically-motivated single-field models of inflation which can give rise to a quasi-inflection point (Ezquiaga et al. 2018) and a high plateau in the matter power spectrum. The action of the non-minimally coupled Higgs field $\phi$ is given by

$$S = \int d^4 x \sqrt{g} \left[ \left( \frac{M_{\mathrm{Pl}}^2}{2} + \frac{\xi(\phi)}{2} \phi^2 \right) R - \frac{1}{2} (\partial \phi)^2 - \frac{1}{4} \lambda(\phi) \phi^4 \right] \quad (2)$$

where the running of Higgs couplings with the renormalization group scale $\mu$ can be parametrized as

$$\lambda(\phi) = \lambda_0 + b_{\lambda} \ln^2(\phi/\mu), \quad (3)$$

$$\xi(\phi) = \xi_0 + b_{\xi} \ln(\phi/\mu), \quad (4)$$

around the critical point $\phi = \mu$.

After standard metric and scalar field redefinitions,

$$g_{\mu\nu} \rightarrow (1 + M_{\mathrm{Pl}}^{-2} \xi(\phi) \phi^2) g_{\mu\nu}, \quad (5)$$
<figure><img src="image_5.png" /><figcaption>Fig. 3 Power spectrum $P_{R=\xi}(k)$ for the single-field Critical Higgs Inflation model, with an inflection point at $N \approx 36$ satisfying the Planck 2018 constraints. Image adapted from García-Bellido and Ruiz Morales (2017) Image reproduced with permission</figcaption></figure>

$$ \phi \to \varphi = \int d\phi \frac{\sqrt{1 + M_{\mathrm{Pl}}^{-2}\xi(\phi)\phi^2\left(1 + 6\left(\xi(\phi) + \frac{1}{2}\phi\xi'(\phi)\right)^2\right)}}{1 + M_{\mathrm{Pl}}^{-2}\xi(\phi)\phi^2}, \quad (6) $$

the effective inflationary potential becomes

$$ V(x) = \frac{V_0 (1 + a \ln^2 x) x^4}{(1 + c (1 + b \ln x) x^2)^2}, \quad (7) $$

with $V_0 = \lambda_0 \mu^4 / 4$, $a = b_\lambda / \lambda_0$, $b = b_\xi / \xi_0$ and $c = \xi_0 M_{\mathrm{Pl}}^{-2} \mu^2$. The potential has a flat plateau at large values of the field $x = \phi / \mu$, see Fig. 2, where $V(x \gg x_c) \simeq V_0 \frac{a}{(b_c)^2} = \frac{M_{\mathrm{Pl}}^4}{4} \frac{b_\xi}{b_\xi^2} \ll M_{\mathrm{Pl}}^4$. The potential also has a short secondary plateau around the critical point, $\phi_c = \mu$, where the inflaton-Higgs slows down and induces a large peak in the curvature power spectrum. This second plateau is induced by a quasi-inflection point at $x = x_c$, where $V'(x_c) \simeq 0$, $V''(x_c) \simeq 0$. As a consequence, the number of e-folds has a sharp jump, $\Delta N$, at that point, plus a slow rise towards larger field values, corresponding to CMB scales. This behavior is very similar to the one discussed in García-Bellido and Ruiz Morales (2017).

In order to have large PBH production and good agreement with the CMB constraints, the allowed range of CHI parameters can be found to be $\lambda_0 \approx (0.01 - 8) \times 10^{-7}$, $\xi_0 \approx (0.5 - 15)$, $x^2 \mu^2 \approx (0.05 - 1.2)$, $b_\lambda \approx (0.008 - 4) \times 10^{-6}$ and $b_\xi \approx (1 - 18)$, for the USR phase to occur $N \in (30, 35)$ e-folds before the end of inflation. These values ensure the right amplitude and tilt of CMB fluctuations as measured by Planck (2018), and would also generate sufficiently large non-Gaussian tails of curvature fluctuations to produce PBH at the QCD epoch, corresponding to masses or order a solar mass, see Fig. 3. The question arises whether these values, corresponding to the model parameters at the critical scale $\mu$,
are consistent with the values of the SM parameters at the EW scale. Given the latest
values of the top quark mass $m_{\text{top}}$ and strong coupling $\alpha_s$ (Sirunyan et al. 2020; Aad
et al. 2019), the values of $\lambda_0$ and $b_\lambda$ that we consider for the Higgs quartic coupling,
are consistent, within $2\sigma$, with the existence of a critical point $\lambda(\mu) \simeq \beta_\lambda(\mu) = 0$
around scales $\mu \simeq 10^{17} - 10^{18}$ GeV, where $\beta_\lambda$ is the beta function of the RGE of the
Higgs self-coupling in the SM.

Taking into account the running of both the Higgs self-coupling $\lambda$ and the non-minimal coupling to gravity $\xi$, there are regions of parameter space allowed by the Standard Model for which the inflaton-Higgs potential acquires a second plateau at smaller scales, around the critical point. This plateau gives a super-slow-roll evolution of the Higgs, inducing a high peak in the curvature power spectrum which is very broad. When those fluctuations reenter the horizon during the radiation era, they collapse to form PBHs with masses in the range 0.1 to $100 M_{\odot}$.

In single-field models with a near-inflection point, the inflaton slows down right
before the end of inflation, creating a stronger backreaction and a quick growth in
curvature fluctuations, giving rise to a significant increase in the power spectrum at
scales much smaller than those probed by CMB and LSS observations. Since the
spectral amplitude is essentially inversely proportional to the parameter $\epsilon_1$, whose
exact evolution can deviate significantly from the slow-roll approximation (García-
Bellido and Ruiz Morales 2017), the power spectrum has to be calculated by
integrating numerically the evolution of the inflaton field $\phi$ following the exact
(beyond slow-roll) equations. In Fig. 3, we show a typical $P(k)$ produced in quasi-
inflection-point models like CHI (Ezquiaga et al. 2018). One can notice the
difference between the exact power spectrum and that obtained in the slow-roll
approximation. The power spectrum can be parametrized as a double step with
different amplitudes and different tilts:

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Examples of primordial power spectra of curvature fluctuations, leading to (planetary-mass) PBH formation: power-law, broken power-law, Gaussian and log-normal models (solid lines), and particular examples of multifield and axion-gauge models (dashed lines). All amplitudes have been tuned to produce PBHs accounting for roughly 10% of all dark matter Image reproduced with permission</figcaption></figure>
$$ P_{\zeta}(k) = \begin{cases} A_1 \left(\frac{k}{k_1}\right)^{n_{s1}-1} & \text{for } k < k_c \\ A_2 \left(\frac{k}{k_2}\right)^{n_{s2}-1} & \text{for } k > k_c \end{cases} \qquad (8) $$

where $A_i$, $n_{si}$, $k_i$ are two different amplitudes, spectral tilts and arbitrary pivot scales for the power spectrum, before and after the cut at $k_c$. It is also possible to introduce some small amount of running of the tilt. This power spectrum neglects quantum diffusion, which induces non-Gaussian contributions (Biagetti et al. 2018; Ezquiaga and García-Bellido 2018) and exponential tails (Ezquiaga et al. 2020), which can be responsible for a large probability of PBH formation even with a relatively small power spectrum amplitude, as will be discussed in subsequent sections. Figure 4 illustrates examples of additional common shapes for the curvature power spectrum utilised in several models discussed below.

### 2.1.4 Reverse engineering approach

In single-field inflationary models, the commonly adopted strategy is to start from some scalar potential that features an approximate stationary inflection point, as shown in Fig. 2; the latter is controlled by various free parameters, which are fine-tuned to guarantee the desired enhancement of the power spectrum of curvature perturbations (Inomata et al. 2017; García-Bellido and Ruiz Morales 2017; Ballesteros and Taoso 2018; Hertzberg and Yamada 2018; Kannike et al. 2017; Dalianis et al. 2019; Inomata et al. 2018; Cheong et al. 2021; Ballesteros et al. 2020c; Iacconi et al. 2022; Kawai and Kim 2021). Following Franciolini and Urbano (2022), a different, and much more powerful, approach can be entertained. In a reverse engineering procedure, the model building can be moved away from the scalar potential—that from this perspective is the output of the analysis instead of the starting point—and focused on the underlying inflationary dynamics (see Ragavendra et al. 2021; Tasinato 2021; Ng and Wu 2021; Karam et al. 2023 for a similar viewpoint).

The inflationary background dynamics can be captured by specifying a minimal set of inputs. Consider the Hubble parameters $\epsilon_1$ and $\eta \equiv \epsilon_1 - \epsilon_2/2$, and the number of e-folds $N$, defined by $dN = Hdt$. Franciolini and Urbano (2022) devised the following analytical ansatz for the time-evolution of $\eta$

$$ 2\eta(N) = \left[ \eta_1 - \eta_{II} + (\eta_{II} - \eta_1) \tanh\left(\frac{N - N_1}{\delta N_1}\right) \right] + \left[ \eta_{II} + \eta_{III} + (\eta_{III} - \eta_{II}) \tanh\left(\frac{N - N_{II}}{\delta N_{II}}\right) \right] + \left[ \eta_{IV} - \eta_{III} + (\eta_{IV} - \eta_{III}) \tanh\left(\frac{N - N_{III}}{\delta N_{III}}\right) \right]. \qquad (9) $$

The inflationary dynamics is therefore divided in different stages: (i) A first stage where $\eta_1$ is constant and negative (necessarily taken to be small to reproduce the
conventional slow-roll dynamics) so that $\epsilon_1$ increases relatively gently; (ii) A USR phase characterized by $\eta_{II} > (3 + \epsilon_1)/2 \simeq 3/2$, where negative friction makes the parameter $\epsilon_1$ decrease abruptly (iii) An optional intermediate stage where $\epsilon_1$ remains constant, achieved if $\eta_{III} \simeq 0$, thus generating a second slow roll phase (i.e. an enhanced plateau in the curvature power spectrum; De Luca et al. 2021g; Franciolini and Urbano 2022); (iv) A final phase characterized by $\eta_{IV} < 0$, that brings $\epsilon_1$ back to $O(1)$ values for inflation to end. The parameters $\delta N_{I,II,III}$ control the sharpness of the transitions.

As shown in Franciolini and Urbano (2022), Franciolini et al. (2022e), it is possible to choose the parameters controlling the ansatz of Eq. (9) in such a way that the power spectrum of curvature perturbations is compatible with CMB constraints at large scales as well as featuring an enhanced peak or plateau at scales beyond the FIRAS constraints (i.e. for $k \gtrsim 10^4$/Mpc), where a sizeable PBH abundance can be obtained, corresponding to the possible formation of masses smaller than $O(10^4) M_\odot$.

Once the Hubble parameters are known, one can compute the inflationary potential by means of

$$ \frac{V(N)}{V(N_{\text{ref}})} = \exp\left\{-2 \int_{N_{\text{ref}}}^{N} dN' \left[\frac{\epsilon_1(3-\eta)}{3-\epsilon_1}\right]\right\}, \quad \phi(N) = \phi(N_{\text{ref}}) \pm \int_{N_{\text{ref}}}^{N} dN' \sqrt{2\epsilon_1}, \qquad (10) $$

where $N_{\text{ref}}$ denotes an arbitrary reference e-fold time, and where in the second equation we consider the minus sign having in mind a large-field model in which the field value decreases as inflation proceeds. Combining $V(N)$ and $\phi(N)$, we reconstruct the profile $V(\phi)$ of the inflationary potential in field space (Byrnes et al. 2019). Equation (10) shows the convenience of modelling the inflationary dynamics directly at the level of $\eta$ instead of $V(\phi)$. This is because the Hubble parameters enter in the exponent in the definition of $V(N)$, and thus allow for a much finer control on power spectral features when performing the reverse engineering procedure. This approach facilitates model building, even though it does not provide an interpretation of the reconstructed potential in terms of high-energy theories of inflation.

## 2.2 Multi-field inflationary models

As seen above, a common challenge of PBH formation within single-field inflationary models is keeping control over the two phases that lead respectively to the seeds for the CMB at large cosmic scales and for PBHs at smaller scales.

Multi-field inflationary scenarios can help separate and better control these two phases of perturbation generation, while also being very natural candidates from the UV perspective. Furthermore, violation of the standard slow-roll conditions can be efficiently realized in a context with more degrees of freedom (Clesse and García-Bellido 2015).
### 2.2.1 Hybrid inflation

In a setup with two scalar fields, the most studied scenario is perhaps *hybrid inflation* (Linde 1994; Copeland et al. 1994). In this class of models, inflation takes place along an almost flat valley in the direction of an initial slowly rolling scalar field, $\phi$. Inflation ends with a waterfall instability when $\phi$ crosses a critical point, $\phi_c$, due to the presence of a second field, $\psi$. The original hybrid potential can be written as

$$V(\phi, \psi) = \Lambda^4 \left[ \left( 1 - \frac{\psi^2}{M^2} + \frac{\phi^2}{\mu^2} + 2 \frac{\phi^2 \psi^2}{\phi_c^2 M^2} \right) \right]. \quad (11)$$

The general shape of the hybrid potential is shown in Fig. 5.

When $\mu M \ll M_{Pl}^2$, the waterfall phase transition is nearly instantaneous and lasts less than one e-fold. The global minima of the potential are located at $\phi = 0$ and $\psi = \pm M$. Inflation can take place along the valley $\psi = 0$ with the effective potential

$$V(\phi, \psi) = \Lambda^4 \left( 1 + \frac{\phi^2}{\mu^2} \right) \quad (12)$$

that most often produces a *blue-tilted* power spectrum for the original potential (a *red-tilted* spectrum can nevertheless be obtained for some specific parameter values; Clesse and Garcia-Bellido 2015). Other hybrid realisations, like F-term or D-term inflation, can produce other shapes for the effective potential, leading to a *red-tilted* spectrum, compatible with the Planck observations (see e.g. Braglia et al. 2023b, which realises this scenario via hybrid $\alpha$-attractors models with scalar potentials with no tachyonic directions and, therefore, bounded from below). PBHs of tiny mass can be produced during a fast waterfall phase (Garcia-Bellido et al. 1996), which would have evaporated today, eventually leading to the reheating of the Universe.

<figure><img src="image_7.png" /><figcaption>**Fig. 5** Representation of a typical hybrid inflation potential, $V$, with a possible trajectory in two-field space (dotted line). CMB perturbations are created along the valley at $\psi = 0$, during a first phase of inflation. Curvature perturbations suitable for production of PBHs are generated in a second flat part of the potential (red dotted line), when the mass-squared of the field $\psi$ changes sign. Image reproduced with permission</figcaption></figure>
In the other regime $\mu M \gg M_{\text{PI}}^2$, the waterfall phase is slow, lasting $N \gg 60$ e-folds (Clesse 2011), and it is then crucial to consider the exact two-field dynamics to calculate the shape of the primordial power spectrum, which is however not compatible with Planck CMB measurements (Clesse et al. 2014).

Hybrid inflation becomes interesting for the production of massive PBHs in the transient regime, when $\mu M \sim M_{\text{PI}}^2$, producing a significant but not too large number of e-folds, $10 \lesssim N \lesssim 40$. During the mild waterfall phase, the multi-field perturbation dynamics leads to a strong, exponential power spectrum enhancement (that can be sharper than the typical maximal $\propto k^4$ increase in single-field models; Byrnes et al. 2019; Cole et al. 2024) followed by a relatively sharp decrease, when going from large to small scales. The position and amplitude of the resulting broad peak is driven by a single combination of the hybrid potential parameters, $\Pi \equiv M \sqrt{\phi_c \mu_1 / M_{\text{PI}}^2}$ where $\mu_1$ comes from a linear expansion of the potential around the critical point $V \propto (\phi - \phi_c)/\mu_1$. Numerical calculations of the primordial power spectrum of curvature fluctuations using either the exact multi-field perturbation dynamics or the $\delta N$ formalism (where the primordial curvature perturbations are identified with the variation of the number of e-folds over its mean, see Sect. 2.3.2 for more details), were first done in Clesse and García-Bellido (2015). The resulting power spectrum is well approximated by a log-normal in mass spectrum,

$$P_{\zeta}(N_k) \approx P_{\text{p}} \exp \left[ - \frac{(N_k - N_{\text{p}})^2}{2\sigma_{\text{p}}^2} \right], \quad (13)$$

of maximal amplitude $P_{\text{p}}$ and width $\sigma_{\text{p}}$, where $N_k$ and $N_{\text{p}}$ denote the e-folds times associated to the Hubble exit of the scale $k$ and of the scale $p$ associated to the peak location in the power spectrum, respectively. Some examples of power spectra are shown in Fig. 6. Finally, let us emphasize that the exact amplitude of the power spectrum also depends on the quantum stochastic dynamics of the auxiliary field $\psi$ close to the critical instability point. In principle, this should lead to slightly different power spectrum amplitudes in different regions of the Universe, and therefore different PBH abundances. A natural outcome of the hybrid inflation scenario would thus be that PBHs are formed in clusters whose mass and size depends on these quantum stochastic dynamics. This can be also somehow used to reduce the fine-tuning linked to the required power spectrum amplitude, since many realisations would have been generically produced within the Universe.

### 2.2.2 Turns in multi-field inflation

Embeddings of inflation in high energy physics motivate scenarios in which multiple scalar fields participate in the inflationary dynamics. An ubiquitous phenomenon in this context is the one of turns in field space, corresponding to a bending of the inflationary trajectory, or more accurately to its deviation from a geodesic of the field
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Power spectrum of curvature perturbations for hybrid inflation parameters values $M = 0.1 M_{\text{Pl}}$, $\mu_1 = 3 \times 10^5 M_{\text{Pl}}$ and $\phi_c = 0.125 M_{\text{Pl}}$ (red), $0.1 M_{\text{Pl}}$ (blue), $0.075 M_{\text{Pl}}$ (green), and $0.05 M_{\text{Pl}}$ (cyan). Those parameters correspond respectively to $\Pi^2 = 375/300/225/150$. The power spectrum is degenerate for lower values of $M$, $\phi$ and larger values of $\mu_1$, keeping the combination $\Pi^2$ constant. For larger values of $M$, $\phi_c$ the degeneracy is broken: power spectra in orange and brown are obtained respectively for $M = \phi_c = M_{\text{Pl}}$ and $\mu_1 = 300 M_{\text{Pl}}/225 M_{\text{Pl}}$. Dashed lines assume $\psi_c = \psi_0$ whereas solid lines are obtained after averaging over 200 power spectra obtained from initial conditions on $\psi_c$ distributed according to a Gaussian of width $\psi_0$. The power spectra corresponding to these realizations are plotted in dashed light gray for illustration. The $\Lambda$ parameter has been fixed so that the spectrum amplitude on CMB anisotropy scales is in agreement with Planck data. The parameter $\mu_2 = 10 M_{\text{Pl}}$ so that the scalar spectral index on those scales is given by $n_s = 0.96$. Image reproduced with permission from Clesse and García-Bellido (2015), copyright by APS Image reproduced with permission</figcaption></figure>

space,<fnref n="1" /> quantified by a dimensionless parameter $\eta_{\perp}(t)$ measuring the acceleration of the trajectory perpendicular to its direction (Groot Nibbelink and van Tent 2000, 2002). When such a turn is sufficiently strong, with $\eta_{\perp} \gtrsim 1$ (see Renaux-Petel 2022 for a short overview), it leads to a transient tachyonic instability of the scalar field fluctuations before Hubble crossing (Cremonini et al. 2011; Brown 2018; Garcia-Saenz et al. 2018), similar in essence to the one arising in axion gauge-field inflation (see Sect. 2.2.4). This gives rise to an exponential enhancement of the power spectrum on those scales, compared to the ones that cross the Hubble radius before and after the turn, and thus constitutes a natural mechanism to seed PBH formation (Palma et al. 2020; Fumagalli et al. 2023; Braglia et al. 2020; Iacconi et al. 2022). Interestingly, owing to its inherently multi-field origin, the growth rate of the power spectrum can overcome the bound deduced in single-field setups (Byrnes et al. 2019; Carrilho et al. 2019; Özsoy and Tasinato 2020; Cole et al. 2024), i.e., the peak of the power spectrum can be sharper. An important parameter is the duration of the turn in e-folds $\delta \equiv \delta N$, which leads to qualitatively different behaviours for broad ($\delta \gtrsim 1$) and sharp turns ($\delta \ll 1$). While broad turns lead to a featureless bump of the power spectrum, sharp turns lead to a localised peak of the envelope of the power

<footnote n="1">¹ In the context of string theory compactifications, it has recently been shown (Freigang et al. 2023) that deviations from geodesic trajectories must be negligible near the boundary of the moduli space. Large turns can therefore be engineered just in the moduli space bulk.</footnote>
spectrum, modulated by order-one rapid oscillations, as described in simple cases by
the analytical power spectrum (assuming a constant value of η⊥ during the
turn) (Palma et al. 2020; Fumagalli et al. 2021a)

$$
\frac{\mathcal{P}_{\zeta}(k)}{\mathcal{P}_0} = \frac{e^{2\sqrt{(2-\kappa)\kappa}\eta_{\perp}\delta}}{2(2-\kappa)\kappa} \times \sin^2\left[e^{-\delta/2}\kappa\eta_{\perp} + \arctan\left(\frac{\kappa}{\sqrt{(2-\kappa)\kappa}}\right)\right]. \quad (14)
$$

Here, $P_0$ denotes the amplitude of the power spectrum in the absence of the transient instability, $\kappa = k/k_f$ is the dimensionless wavenumber in units of the maximally enhanced one, with $k_f$ the scale crossing the Hubble radius at the time of the sharp turn, and the formula is valid for $\kappa \le 2$ (see Fumagalli et al. 2021a for generalisations). The oscillations there, which are to a very good approximation periodic in $k$, are characteristic of a sharp feature, and go hand in hand with a boosted power spectrum generated by an event localised in time (Fumagalli et al. 2021a). The mass function of PBHs is highly sensitive to the tail of the probability density function (PDF) of the smoothed density contrast, which is not yet known for the class of models described here. However, tentative results for the mass function indicate that the oscillatory patterns are washed out when assuming Gaussian statistics of the primordial curvature fluctuations (Fumagalli et al. 2023), and hence further investigation is required to assess the observational consequences for PBHs of the oscillatory features. What is known, however, is that the SGWB generated by scalar fluctuations at horizon re-entry does keep an imprint of the oscillations of the scalar power spectrum: they manifest as oscillations in the frequency profile of the SGWB, providing a robust probe of small-scale features that would be inaccessible otherwise (Fumagalli et al. 2021a, b; Braglia et al. 2021a; Witkowski et al. 2022; Fumagalli et al. 2022b).

If one uses the standard criterion that the power spectrum has to reach values
$\mathcal{P}_{\zeta} \sim 10^{-2}$ to generate a substantial amount of PBHs, it is interesting that the required
enhancement from the baseline value $\mathcal{P}_0 \sim 10^{-9}$ necessitates a total angle $\Delta\theta = \eta_{\perp}\delta$
swept by the trajectory during the turn exceeding the value $\pi$ obtainable in flat field
space, and hence requires the inflationary field space to be curved (Palma et al.
2020). Such large values of the power spectrum are not without consequences
however, and considerations of backreaction and perturbative control are important
to take into account to delineate the space of theoretically viable models (Fumagalli
et al. 2021a). Indeed, the criterion for perturbative control parametrically depends on
the large parameter $\eta_{\perp}$ in a way that imposes the bound $\eta_{\perp}^6 \mathcal{P}_{\zeta} \lesssim 1$ (or $\eta_{\perp}^4 \mathcal{P}_{\zeta} \lesssim 1$) for
broad (Fumagalli et al. 2019; Bjorkmo et al. 2019) (or sharp) turns (Fumagalli et al.
2021a). Thus, an important open question is whether the inevitably present non-
Gaussianities of the curvature fluctuations in these models, enhanced in flattened
configurations (Garcia-Saenz and Renaux-Petel 2018; Fumagalli et al. 2019; Ferreira
2020), are such that a substantial amount of PBHs can be generated with a value of
the power spectrum less than $10^{-2}$ while ensuring perturbative control.
### 2.2.3 Two-stage models

There are several concrete examples of PBH formation resulting from a turn in the field space during inflation that separates two stages of inflation. The first stage produces the primordial power spectrum on CMB scales while the second stage is responsible for the enhancement of power on smaller scales, leading to PBH formation. Oscillations are expected at the transition, as explained in the previous Sect. Pi et al. (2018) have proposed a model combining Starobinsky inflation with a scalaron field that becomes massive at the end of the first stage of inflation and a curvature field that plays the role of the inflaton in a second stage. Kawasaki et al. (2016) proposed a scenario with a second phase of double inflation in a supergravity framework. One should notice that non-Gaussianity is typically small in these models, since most of the field trajectories are effectively single-field. Other advantages of these models are that one does not especially need to tune the position of a peak in the power spectrum, since it can be nearly scale invariant, as generically predicted by slow-roll inflation, and that the transition between the two regimes can be very sharp. Due to the features at the transition in the power spectrum, such models can also produce multi-modal distributions of black holes, as proposed in Carr and Kuhnel (2019).

### 2.2.4 Axion-gauge scenario

PBHs can also arise in models in which rolling axions interact with gauge fields, enhancing primordial density perturbations (Linde et al. 2013; Bugaev and Klimai 2014; Domcke et al. 2017; García-Bellido et al. 2016, 2017; Özsoy and Lalak 2021; Özsoy 2021). In natural inflation (Freese et al. 1990; Adams et al. 1993), the couplings of the inflaton to matter respect a shift symmetry, $\phi \to \phi + \text{constant}$. Therefore, they do not provide by themselves quantum corrections to the inflaton potential (that arise from a breaking of the shift symmetry, as for instance from instanton effects), which can therefore be technically naturally small and flat. Also, UV-completed theories, such as string theory or supergravity models of inflation, contain a large number of pseudo-scalar particles that could be identified with the axion inflaton, or some other dynamically relevant axion. The shift-symmetric coupling between the axion and the gauge field is

$$ \mathcal{L} \supset \frac{a}{f} F\tilde{F}, \qquad (15) $$

where a is the pseudoscalar axion, F is the gauge field strength, $\tilde{F}$ its dual, and f is the axion decay constant. Shift-symmetry is satisfied due to the fact that $FF$ is a total derivative and the constant term becomes a topological invariant. This interaction amplifies one helicity of the vector field around horizon crossing. The amplified gauge quanta might affect the evolution of the axion field (Anber and Sorbo 2010) or, even if produced in a smaller amount, source scalar density perturbations (Barnaby et al. 2011) and parity violating tensor modes (Sorbo 2011). These effects are exponentially sensitive to the parameter $\xi \equiv \frac{a}{2fH}$, and visible effects are obtained for
$f \simeq 10^{-2} M_{\text{PI}}$, namely for an axion decay constant parametrically equal to the string scale. On one hand, this enhances the possibility that models of string cosmology have a rich phenomenology. On the other hand, no effect is expected for models of a single axion inflaton, that are required to have a super-Planckian excursion, unless some specific mechanism allows excursions much greater than the scale of inflation, such is the case for monodromy (Silverstein and Westphal 2008) or for aligned axion inflation (Kim et al. 2005).

Scalar perturbations are produced via inverse decay $A + A \rightarrow a$ (where A denotes the gauge field), and therefore, they obey a non-Gaussian $\chi^2$ distribution. The amplified gauge fields also produce GWs via an analogous $2 \rightarrow 1$ process. We note that these primordial GWs, which are produced during inflation, form a different stochastic GW background with respect to the secondary GWs that are sourced by density perturbations in the process of horizon re-entry in radiation or matter-dominated eras.

The power spectrum of density perturbations can be parametrized as (Namba et al. 2016; Peloso et al. 2016; García-Bellido et al. 2016)

$$P_{\zeta}(k) = P_{\text{vacuum}}(k) + P_{\text{sourced}}(k), \qquad (16)$$

where $P_{\text{vacuum}} \simeq 2 \times 10^{-9} (k/k_{\text{cmb}})^{n_s-1}$ is the standard nearly scale invariant vacuum signal (with $k_{\text{cmb}}$ being some pivot wavenumber at CMB scales). As already mentioned, the sourced signal is exponentially sensitive to the parameter $\xi \equiv \dot{a}/2fH$, and it is therefore highly sensitive to the speed of the axion, and, ultimately, to its potential. For axion inflation with a single monomial potential the signal is typically blue, and one simply needs to require that the production is sufficiently small so that no stable PBH is produced via this mechanism (Linde et al. 2013). This does not need to be the case for a more complicated potential, or in the case in which the axion is not identified with the inflaton (García-Bellido et al. 2016, 2017). In this second case, one can imagine a situation in which the axion rolls significantly only for a limited number of e-folds (Namba et al. 2016; Peloso et al. 2016) and the sourced spectrum can be described as a nearly log-normal bump, peaked at the scales that were exiting the horizon while $a$ was rolling. The infra-red part of the spectrum scales as $P \propto k^3$ due to causality, while, close to the peak, the spectrum is well parametrized by (Namba et al. 2016)

$$P_{\text{sourced}}(k) = A(\xi_*, W) \exp \left[ - \frac{\ln^2(k/k_*)}{2\sigma^2(\xi_*, W)} \right], \qquad (17)$$

where the amplitude $A$, the central position $k_*$, and the width $\sigma$ of the peak are functions of the maximum value of the particle production parameter $\xi$ and of the amount of e-folds $W$ during which the axion is significantly rolling. An analogous parametrization also holds for the amount of sourced GWs (Namba et al. 2016).

The sourced scalar perturbations are highly non-Gaussian (Barnaby et al. 2011). For a nearly scale-invariant signal (as is the case if the axion is the inflaton), their bispectrum has a shape which is very close to the equilateral one. For the case of short-duration production leading to Eq. (17), the bispectrum is a maximum when the
magnitude of the three momenta is parametrically equal to $k_*$, and its scale-dependence can also be described accurately via a log-normal shape (Namba et al. 2016). The fact that the sourced scalar perturbations obey a $\chi^2$ statistics makes PBH production much more significant compared to Gaussian perturbations with the same power. Therefore, an equal amount of PBHs is produced with a much smaller power spectrum than in the Gaussian case. This results in a significantly smaller amount of stochastic GWs produced at second order by the scalar perturbations, with a different prospect for detection at LISA and PTA-SKA, see for example Fig. 7. This is discussed in more detail in Sect. 5.1.

### 2.2.5 False-vacuum models

It is also possible to form PBHs in both single and multi-field scenarios if some spatial regions end up in a false-vacuum. These regions can collapse into black holes before the field can pass the potential barrier through quantum tunneling and evolve towards the true vacuum. Several examples have been proposed.

For instance, the Higgs field is suspected to be metastable. If it is a spectator field during inflation, some regions may have ended in this region of the potential, thereby leading to the formation of black holes (Passaglia et al. 2020). Nevertheless, their abundance cannot be significant enough to contribute to the dark matter. Another possibility is a tiny step in the inflaton potential that can basically be overcome by most of field trajectories, but not all (Cai et al. 2022b). This can radically change the

<figure><img src="image_9.png" /><figcaption>**Fig. 7** The stochastic GWs produced by axion inflation at four main scales of interferometers for which we show the approximate sensitivities, nHz (PTAs and SKA), µHz (Gaia and Theia), mHz (LISA and BBO) and Hz (LVK, Einstein Telescope, and Cosmic Explorer). The enhanced density perturbations may produce PBHs which are a significant fraction of dark matter for $1-100 M_\odot$ and the totality in the $10^{-14}-10^{-11} M_\odot$ mass range and, remarkably, the enhanced perturbations leave inevitable GW backgrounds at the most sensitive regimes of GW detectors: for the first range of BHs this corresponds to PTA frequencies (and possibly future SKA) and LISA frequencies for the second range. The PBHs formed from fluctuations peaking at Hz scales will be so light that they are expected to be part of thermal history through Hawking radiation. The figure is updated (with slight modifications) from García-Bellido et al. (2017), Ünal (2019), García-Bellido et al. (2021b) Image reproduced with permission</figcaption></figure>
tail of the Gaussian fluctuations, even if the conventional non-linearity parameters ($f_{NL}$, $g_{NL}$, ...) remain small. The regions with large curvature fluctuations later form PBHs.

It has also been proposed that a light spectator field during inflation, having an
asymmetric polynomial potential, experiences stochastic fluctuations that displace
the field from the global minimum of its potential and populate a false vacuum
state (Maeso et al. 2022). Lattice simulations have been used to show how when the
Hubble radius reaches the false-vacuum bubble size, the potential energy is
transferred to the kinetic and gradient energies of the bubble wall, such that it begins
to contract and possibly collapse into a black hole. Finally, multi-field realisations
have been proposed, for instance in Garriga et al. (2016), Deng and Vilenkin (2017).

## 2.3 Non-Gaussian models

### 2.3.1 Motivations

In Sect. 2.1, we discussed how PBHs could be generated in single-field models due to an enhancement of the power spectrum caused, for example, by a second plateau in the potential. However, as we will discuss in more detail in Sect. 3, what really matters to determine the PBH formation is the full profile of the PDF of the primordial curvature perturbations, in particular the probability to be above a given threshold. This implies that PBH abundance is largely controlled by the tails of the PDF, which at the same time are highly sensitive to any non-Gaussian corrections.

Non-Gaussian corrections can either enhance or suppress the production of PBHs
depending on their effect on the tail of the PDF. One way of accounting for these
effects is to compute the different statistical moments of the distribution. Gaussian
profiles are fully characterized by the second moment (the power spectrum), while
general PDFs can have contributions from any of the higher moments. For illustrative
purposes, we plot in Fig. 8 the effect of the third and fourth moment as determined by
the skewness and kurtosis, respectively. As it can be clearly seen, a positive kurtosis
has a large impact in the tail of the PDF. Interestingly, one can analytically compute
the effect of non-zero higher order moments in the abundance of PBHs (Franciolini
et al. 2018; Matsubara and Sasaki 2022).

Non-Gaussian corrections are known to be relevant in multi-field models of
inflation. For instance, models of axion inflation in which the gauge field sources the
curvature perturbations (García-Bellido et al. 2016) display $\chi^2$ statistics for the
curvature perturbation. This has implications both for the PBH abundance and for the
GW signature, as studied in detail in García-Bellido et al. (2017). In particular, a
strong increase of the PBH abundance associated with a greater tail of the PDF does
not correspond to a strong increase of the associated GW production. Therefore, for
equal PBH abundance, a smaller GW signal is produced with respect to the Gaussian
case.

Moreover, these non-Gaussian corrections can also enhance the abundance of
PBHs in single-field models due to the dynamics of the ultra-slow roll phase (Biagetti
et al. 2021; Cai et al. 2022a) and the effect of quantum diffusion beyond slow-
<figure><img src="image_10.png" /><figcaption>**Fig. 8** Illustration of the effect of a positive skewness, $\bar{\kappa}_3$, and a positive (excess) kurtosis, $\bar{\kappa}_4$, for the probability density function $P(x)$ (left) and its logarithm (right) in comparison with a Gaussian distribution ($\bar{\kappa}_N > 2 = 0$). The dominant effect in the tail is given by a positive kurtosis. Image reproduced with permission from Ezquiaga and García-Bellido (2018), copyright by IOP and SISSA Image reproduced with permission</figcaption></figure>

roll (Firouzjahi et al. 2019; Biagetti et al. 2018; Ezquiaga and García-Bellido 2018; Pattison et al. 2021).<fnref n="2" /> In fact, in general terms, quantum diffusion introduces an exponential tail in the PDF of curvature perturbations (Pattison et al. 2017; Ezquiaga et al. 2020). For this reason, we discuss next the effect of quantum diffusion in the production of PBHs.

### 2.3.2 Quantum diffusion

Primordial curvature perturbations are expected to be well-described by a quasi-Gaussian distribution when they are small and close to the maximum of their probability distribution. This can be modelled using conventional cosmological perturbation theory where the free field fluctuations during inflation describe a Gaussian distribution and higher-order interactions lead to non-Gaussian corrections. However, PBHs result from the collapse of large density perturbations, which may be far from the peak of the distribution. In this regime, we require a non-perturbative approach, such as the stochastic-δN formalism (Enqvist et al. 2008; Fujita et al. 2013; Vennin and Starobinsky 2015).

Stochastic inflation (Starobinsky 1982) splits the inflaton field into small-scale quantum fluctuations and a large-scale effectively classical field above a coarse-graining scale (typically around the Hubble length). The large-scale field follows the non-linear evolution of a (locally) homogeneous and isotropic cosmology (Salopek and Bond 1990; Wands et al. 2000; Rigopoulos and Shellard 2003), except that small-scale modes are constantly swept up from the quantum vacuum state, crossing into the coarse-grained field, leading to quantum diffusion $\langle \Delta\varphi^2 \rangle \approx (H/2\pi)^2$ per Hubble volume per Hubble time. In regimes where PBHs are formed, the inflaton must undergo large quantum fluctuations (for instance as the result of travelling

<footnote n="2">² A violation of slow-roll is necessary for non-Gaussian effects to be sufficiently large to change the PBH abundance in single-field models (Passaglia et al. 2019).</footnote>
through a very flat potential), which implies that stochastic corrections are important
and must be taken into account (Pattison et al. 2017; Ezquiaga and García-Bellido
2018; Biagetti et al. 2018; Ezquiaga et al. 2020).

In order to derive the predicted abundance of PBHs in stochastic inflation, we
make use of the δN formalism (Starobinsky 1982, 1985; Sasaki and Stewart 1996;
Wands et al. 2000; Lyth et al. 2005), in which the primordial curvature perturbation is
identified with the perturbation in the number of e-folds $N = \int H dt$ with respect to
the average, $\zeta = N - \langle N \rangle$ (where $N$ denotes the duration of inflation starting from
an initial field configuration, measured in number of e-folds $N$). One thus has to
extract the distribution $P_{\phi}(N)$ of the number of e-folds of expansion $N$ from a given
inflaton field value $\phi$ to the end of inflation. This is the program of the stochastic δN
formalism, in which it can be shown (Vennin and Starobinsky 2015; Pattison et al.
2017; Ezquiaga et al. 2020; Ando and Vennin 2021) that the Fokker-Planck equation
driving the stochastic evolution of the coarse-grained field with respect to the
integrated expansion during inflation, $\partial P_N(\phi)/\partial N = \mathcal{L}_{\text{FP}}(\phi) \cdot P_N(\phi)$, can be recast
in terms of the adjoint Fokker-Planck operator acting on the duration distribution
$P_{\phi}(N)$, such that $\partial P_{\phi}(N)/\partial N = \mathcal{L}_{\text{FP}}^{\dagger}(\phi) \cdot P_{\phi}(N)$.

By solving this equation in different setups, one can show that quantum diffusion
leads unavoidably to non-Gaussian tails (Pattison et al. 2017; Ezquiaga et al. 2020).
In particular, if the classical evolution gives a finite number of e-folds $N_c$ from a
given field value to the end of inflation, then the large-$N$ limit (where $N \gg N_c$) is
inevitably dominated by quantum diffusion since it requires many quantum kicks to
evade the classical drift towards the end of inflation. Quite generally, the probability
distribution approaches an exponential tail, $P_{\phi}(N) = \sum_n a_n e^{-\Lambda_n N} \propto e^{-\Lambda_0 N}$ as
$N \to \infty$ (Pattison et al. 2017; Ezquiaga et al. 2020; Jackson et al. 2022), which
cannot be described by a quasi-Gaussian distribution (see Fig. 9). Thus predictions
for the abundance of PBHs, which will be described by the integrated probability
distribution function above a given threshold, may differ by several orders of
magnitude from results obtained assuming a Gaussian distribution (Pattison et al.
2017; Ezquiaga et al. 2020; Biagetti et al. 2021; Kitajima et al. 2021; De Luca and
Riotto 2022; Gow et al. 2023).

Let us also stress that these exponential tails cannot be properly described by
usual, perturbative parametrisations of non-Gaussian statistics (such as those based
on computing the few first moments of the distribution and the non-linearity
parameters $f_{\text{NL}}$, $g_{\text{NL}}$, etc.), which can only account for polynomial modulations of
Gaussian tails. A non-perturbative approach, such as the one presented here, is
therefore necessary.

It is also worth mentioning that the details of the tail (such as the value of $a_0$ and $\Lambda_0$) are determined by the specifics of the inflationary model under consideration. For instance, in the simplest setup where inflation is realised with a single scalar field, if the potential energy $V$ is constant across a field region of size $\Delta\phi$, then (Pattison et al. 2017) $\Lambda_0 = V/(96M_{\text{Pl}}^2\Delta\phi^2)$ and $\Lambda_n = (2n+1)^2\Lambda_0$. Wider “quantum wells” (i.e. larger values of $\Delta\phi$) thus lead to smaller $\Lambda_0$ and to wider tails. This result is valid if the initial velocity of the inflaton into the quantum well is vanishing. Otherwise, a transient period of “ultra slow roll” takes place until the inflaton relaxes to the slow-
<figure><img src="image_11.png" /><figcaption>**Fig. 9** Probability density functions of the curvature perturbations generated by an inflection-point inflationary potential (see inset). Solid lines correspond to the full distribution functions computed by mean of the stochastic-δN formalism, where different colours correspond to different locations in the potential v(φ) where the scale under consideration emerges from the Hubble radius and where Δφ<sub>well</sub> denotes the field range for which stochastic effects dominate over the classical dynamics. The dotted lines correspond to the standard result, which provides a good Gaussian approximation for the maximum of the distribution but that however fails to describe the exponential tails, where PBHs are nonetheless produced. Image reproduced with permission from Ezquiaga et al. (2020), copyright by IOP and SISSA Image reproduced with permission</figcaption></figure>

roll attractor, and on top of the eigenvalues Λ_n given above, a second set of decay rates Λ'_n = Λ_n + 3 also arises (Pattison et al. 2021). When the initial velocity is large compared to the velocity that would be required to cross the well classically, the amplitudes a'_n associated with those eigenvalues are larger than those for the rates Λ_n, hence the tails are more suppressed. Another way to suppress the tails is to add a tilt in the inflationary potential, V = 24π²M_Pl⁴v₀(1 + αφ/M_Pl), for which the decay rates are all shifted by a constant value (Ezquiaga et al. 2020), i.e. Λ'_n = Λ_n + α²/(4v₀). Unless α² is fine tuned to values less than the potential energy in Planckian units, the exponential tails are heavily suppressed. As a last example, let us mention the case where the potential features an inflection point, close to which one can approximate V = 24π²M_Pl⁴v₀[1 + β(φ/M_Pl)³]. In such models, if slow roll is not violated when approaching the inflection point (which depends on the details of the inflationary potential away from that point), one finds (Ezquiaga et al. 2020) Λ_n = (3/2)²/³π²(v₀α)¹/³(n + 1/2)², such that primordial black holes are always overproduced.

These examples show that the abundance of primordial black holes is ultimately related to the details of the inflationary potential where quantum diffusion plays an important role. PBHs would thus provide access to regions of the inflationary potential that can otherwise not be probed with other cosmological measurements such as the cosmic microwave background.

2.3.3 Stochastic spectator field

Another recent scenario (Carr et al. 2021a) takes the advantage of the dynamics of a
stochastic spectator field to produce PBHs, while the primordial power spectrum
remains at the CMB level almost everywhere in our Universe, except in subdominant regions where PBHs can form. If PBHs are linked to the generation of the baryon asymmetry of the Universe, thereby connecting (1) The PBHs and baryon abundance, (2) The baryon-to-photon ratio to the PBH density at the time of their formation, the stochastic and quantum nature of the spectator field fluctuations allows one to invoke an anthropic selection mechanism to resolve the long-standing parameter fine-tuning issues related to PBH formation (García-Bellido et al. 2021a; Carr et al. 2021a; García-Bellido 2019).

In this scenario, the stochastic spectator field, denoted $\psi$, experiences quantum fluctuations during inflation that lead it to explore a wide range of the potential, including its slow-roll part, in different Hubble patches. But the field does not have any impact on the inflationary dynamics, thus being a spectator field. Because the number of patches is huge, there necessarily exist some in which the spectator field $\psi$, at the time of Hubble exit, has the required mean value for subsequent quantum fluctuations (but sub-Hubble at that time) to later induce large curvature fluctuations. Our observable Universe would correspond to one of these patches.

After inflation, the spectator field and its fluctuations within these patches remain frozen during the radiation era until its potential energy starts to dominate the density of the Universe, well after inflation but also possibly well before the time of PBH formation. At this point, in the regions where the field is in the slow-roll part of the potential, a very short second inflationary phase (for at most a few e-folds) takes place, whereas in the rest of the patch the field quickly rolls down its potential without inflating. Such an extra expansion corresponds to the production of local non-linear curvature fluctuations that are still super-horizon at that time. Only later they re-enter the horizon and collapse into PBHs. But in the rest of the Universe the curvature fluctuations are statistically Gaussian and behave as expected in standard slow-roll inflation, unaffected by the spectator field.

The probability to form a PBH in a given region (that can be connected to the function $\beta$ introduced in the next section) has been calculated as (Carr et al. 2021a)

$$P_{\text{PBH}} = \sqrt{\frac{2}{\pi}} \frac{\Delta\psi^{\text{sr}}}{\sqrt{H_N^2/4\pi^2 + \langle\delta\psi^2\rangle_{N-1}}} \times \exp\left[-\frac{\langle\psi\rangle^2}{2(H_N^2/4\pi^2 + \langle\delta\psi^2\rangle_{N-1})}\right], \quad (18)$$

where $\Delta\psi^{\text{sr}}$ is the size of the slow-roll region where the field can induce $O(1)$ curvature fluctuations. The subscript $N$ is used to indicate that a quantity is evaluated at the e-fold time at which the scale exited the horizon during inflation. It can be noticed that the mathematical form of $P_{\text{PBH}}$ is quite similar to the one of $\beta$ (see next section) in more standard models of Gaussian inflationary perturbations. As a consequence, the shape of the PBH mass function is also very similar and driven by the evolution of $H$ during inflation. There is also no need to precisely fix the potential parameters in this model, the only condition being that the potential is of plateau or small-field type and that it can dominate the density of the Universe prior the QCD epoch (if one is interested to get stellar-mass or heavier PBHs). More precise calculations of the probability distribution of curvature fluctuations and some embedding in high-energy frameworks, including with the Brought-Englert-Higgs field,
were recently presented and discussed in Stamou and Clesse (2024a), Stamou and Clesse (2024b).

## 2.4 Curvaton

A different possibility to generate PBHs, while keeping the generation mechanisms for curvature perturbations at different scales decoupled, is represented by a simple modification of the original *curvaton scenario* (in the original curvaton model (Lyth and Wands 2002; Moroi and Takahashi 2001), one field, the inflaton, is responsible only for sustaining the inflationary background evolution while the other field, the curvaton, produces the primordial curvature perturbations). In this setup with two fields, the primordial curvature perturbations on CMB scales are produced by the inflaton, which acts very similarly to the standard single-field scenario, while the curvaton field becomes responsible for perturbations on smaller scales and thus formation of PBHs.

The power spectrum of curvature perturbations will be given as the sum of the contributions of the inflaton and of the curvaton

$$ \mathcal{P}(k) = \mathcal{P}_{\text{infl}}(k) + \mathcal{P}_{\text{curv}}(k). \qquad (19) $$

where $\mathcal{P}_{\text{curv}}(k)$ can be described by a lognormal function such as in Eq. (17).

At large scales (reasonably $k \lesssim 1 \text{Mpc}^{-1}$), the first term will be dominating, with its magnitude being determined by CMB normalization, that is $\mathcal{P}_{\text{infl}}(k) \simeq 2 \times 10^{-9}$. Moving to smaller scales (larger $k$), the second contribution, given by the curvaton field, increases and becomes dominant.

So, as long as the curvature perturbation on large cosmological scales is dominated by inflaton fluctuations with an almost scale-invariant spectrum, the curvaton field is not required to be light during inflation and its fluctuations may have a steep blue spectrum. Curvaton fluctuations could then give the dominant contribution to the primordial density perturbations on small scales after the curvaton decays. The curvaton may have a steep blue spectrum either due to interactions with the inflaton or other fields evolving during inflation (Yokoyama 1997), a non-trivial kinetic term (Pi and Sasaki 2023) or in an axion-like model, where the curvaton is identified with the phase of a complex field whose modulus decreases rapidly during inflation (Kawasaki et al. 2013).

An important characteristic of any curvaton model is that the fluctuations originate from non-adiabatic field fluctuations during inflation, a fact which gives rise to a non-Gaussian distribution in the primordial curvature perturbations, including a finite non-Gaussianity in the squeezed limit. Although this local-type non-Gaussianity is typically small ($f_{\text{NL}} \sim 1$), it can nonetheless have a significant effect on the abundance of PBHs, as discussed in Sect. 2.3.

## 2.5 Preheating

Another way through which one can produce PBHs is the preheating instability. During the preheating period after inflation, when the inflaton oscillates coherently at
its ground state and decays to other degrees of freedom, it is widely argued that resonant amplification of quantum field fluctuations, responsible for particle production, take place (Kofman et al. 1994, 1997). These amplified quantum fluctuations should be accompanied by a resonant amplification of the scalar metric fluctuations (usually quoted as metric preheating, see Finelli and Brandenberger 1999; Bassett et al. 1999b; Jedamzik and Sigl 2000; Bassett and Viniegra 2000), responsible for gravitational fluctuations in the curvature, since the two are coupled through Einstein's equations. One then expects that the amplified metric perturbations give rise to large curvature/matter perturbations, which in their turn collapse and form PBHs.

Historically, PBHs emanating from the preheating instability were speculated in the context of multi-field inflation and in particular in the context of two-field chaotic inflation (Bassett et al. 1999a; Green and Malik 2001; Bassett and Tsujikawa 2001; Suyama et al. 2005; and more recently in Torres-Lomas and Urena-LAlpez 2013; Torres-Lomas et al. 2014), since in this case the parametric amplification of entropy (isocurvature) fluctuations can source the parametric amplification of the adiabatic (curvature) fluctuations in the regime of broad resonance. Thus, these amplified curvature fluctuations could break the scale-invariance of the primordial power spectrum and give birth to PBHs through non-linear gravitational fluctuations inducing large density contrasts (Kou et al. 2021; Joana 2022).

It is predicted as well (Jedamzik et al. 2010; see also Easther et al. 2011) that in the context of single-field inflation, there is a pronounced resonant instability structure in the narrow regime where amplified metric perturbations can induce the production of PBHs as recently studied in Martin et al. (2020b), Auclair and Vennin (2021). At this point, one should point out that the narrow resonant structure of metric preheating was shown in Martin et al. (2020a) to be immune to the decay of the inflaton into a radiation fluid, ensuring in this way the transition to the Hot Big Bang phase of the Universe.

The case of metric preheating Most of the studies in which PBHs are produced during preheating treat them in the context of multi-field and especially in the context of two-field inflation. This makes it difficult to predict in an analytic way the matter power spectrum responsible for the PBH formation, which is constructed consequently numerically. On the contrary, in the context of single-field inflation, one can in principle extract in a direct way the matter power spectrum responsible for the PBH production. In what follows we emphasize on a specific scenario studied in Martin et al. (2020b), Auclair and Vennin (2021), in which PBHs are produced during the preheating phase due to the emergence of a resonant instability structure concerning the equation of motion of the scalar metric perturbations that are enhanced on small scales (metric preheating), which exit the Hubble radius close to the end of inflation. Specifically, the modes concerned satisfy the relation

$$aH < k < a\sqrt{3}Hm, \qquad (20)$$

where k is the comoving scale concerned, a is the scale factor, H is the Hubble parameter and m is the mass parameter of the inflationary potential. For these modes, as shown in Martin et al. (2020b), the density contrast scales linearly with the scale
factor, i.e. $\delta_k \sim a$, manifesting in this way an effective “dust” behavior during the instability phase in which PBHs can be easily produced.

Considering only the modes entering the instability band from above, namely the left hand side of the inequality (20), the matter power spectrum responsible for the PBH formation at the end of the preheating instability phase can be written analytically as (Martin et al. 2020b)

$$ \mathcal{P}_{\delta}\left(t_{\Gamma}, \frac{k}{k_{\text{end}}}\right)=\left(\frac{k}{k_{\text{end}}}\right)^{4}\left(\frac{\rho_{\text {inf }}}{\rho_{\Gamma}}\right)^{2 / 3}\left(\frac{6}{5}\right)^{2} \mathcal{P}_{\zeta, \text {end}}\left(\frac{k}{k_{\text {end }}}\right), \qquad (21) $$

where $k_{\text{end}}$ is the comoving scale exiting the Hubble radius at the end of inflation, $\rho_{\text{inf}}$ is the energy scale at the end of inflation, $\rho_{\Gamma}$ is the energy scale at the end of the preheating instability phase and $\mathcal{P}_{\zeta, \text{end}}(\frac{k}{k_{\text{end}}})$ is the curvature power spectrum at the end of inflation which in the slow-roll approximation reads as (Schwarz et al. 2001; Gong and Stewart 2001)

$$ \mathcal{P}_{\zeta, \text {end}}\left(\frac{k}{k_{\text {end }}}\right)=\frac{H_{*}^{2}(k)}{8 \pi^{2} M_{\text {Pl}}^{2} \epsilon_{1*}(k)}\left[1+\left(\frac{k}{k_{\text {end }}}\right)^{2}\right] \quad (22) \\ \qquad \qquad \qquad \qquad \qquad \qquad \left[1-2\left(C+1\right) \epsilon_{1*}(k)-C \epsilon_{2*}(k)\right] \text { for } k<k_{\text {end }} $$

The functions $H_*(k)$, $\epsilon_{1*}(k)$ and $\epsilon_{2*}(k)$ denote respectively the values of $H$, $\epsilon_1$ and $\epsilon_2$ at the time when the mode $k$ exits the Hubble radius during inflation. The parameter $C \simeq -0.7296$ is a numerical constant. It is important to mention here that Eq. (22) is the curvature power spectrum at NLO within the slow-roll approximation in contrast with Eq. (1) which gives $\mathcal{P}_{\zeta}(k)$ at leading-order (LO). This is why we obtain the extra $1 + (k/k_{\text{end}})^2$ term, which captures the first NLO correction coming from the decaying mode of $\zeta$ and is actually significant for modes $k$ close to $k_{\text{end}}$ as the ones considered here. Note also that Eq. (1) is expanded around a pivot scale $k_*$ whereas Eq. (22) is valid for all values of $k$ within the slow-roll approximation.

The relevant range of the density contrast at the end of the preheating phase is

$$ \left(\frac{\rho_{\text{inf}}}{\rho_{\Gamma}}\right)^{1/3} \left(\frac{k}{k_{\text{end}}}\right)^{2} \left(\frac{3\pi}{2}\right)^{2/3} \left[\left(\frac{k}{k_{\text{end}}}\right)^{3} \sqrt{\frac{\rho_{\text{inf}}}{\rho_{\Gamma}}} - 1\right]^{-2/3} < \delta_k(t_{\Gamma}) < 1, \qquad (23) $$

while the relevant parameter space $(\rho_{\Gamma}, \rho_{\text{inf}})$ so that we do not have PBH overproduction during preheating reads as

$$ \frac{\rho_{\Gamma}}{M_{\text{Pl}}^{4}} < \frac{4}{125 \sqrt{3} \pi^{5}} \left(\frac{\rho_{\text{inf}}}{M_{\text{Pl}}^{4}}\right)^{5 / 2} \qquad (24) $$

Below, in Fig. 10 we see the behavior of $\mathcal{P}_{\delta}(t_{\Gamma}, \frac{k}{k_{\text{end}}})$ [see Eq. (21)] by fixing $\rho_{\text{inf}}^{1/4} = 10^{14}$ GeV and $\rho_{\Gamma}^{1/4} = 10^7$ GeV. At this point, it is important to stress out that exactly at $k = k_{\text{end}}$, Eq. (22) is not trustful since there the slow-roll approximation breaks down.
The corresponding mass of the PBHs produced is related to the mass inside the Hubble radius at the time when a scale k re-enters the horizon, namely

$$m_{\text{PBH}} = \gamma M_{\text{H}} = \gamma \frac{(3M_{\text{pl}}^2)^{3/2}}{\sqrt{\rho_{\text{inf}}}} \left(\frac{k}{k_{\text{end}}}\right)^{-3}, \quad (25)$$

where $\gamma$ is an efficiency parameter of order one depending on the details of the gravitational collapse. Here, it is important to notice that in Eq. (25), one clearly sees the dependence of the PBH mass on the energy scale at the end of inflation and as a consequence on the underlying inflationary model. We need to stress here that Eq. (25) gives the PBH mass as a function of the comoving scale k for PBHs forming during preheating, when the Universe behaves effectively as pressureless dust. This gives rise to a k scaling of the PBH mass of the form of $k^{-3}$ in contrast with the $k^{-2}$ scaling valid for a PBH forming during a radiation-dominated era [See the discussion in Sect. 3].

Regarding the characteristic PBH mass produced during preheating one can choose $k/k_{\text{end}} = 0.1$, and using a relative high energy scale at the end of inflation, $\rho_{\text{inf}}^{1/4} = 10^{14}$ GeV, from Eq. (25) one obtains $m_{\text{PBH}} \simeq 10^{-27} M_{\odot}$. This mass, around which the matter power spectrum peaks, is very small compared to the mass range LISA can detect through the GW portal of PBH mergers. However, as pointed out recently in Auclair and Vennin (2021), in which the PBH mass function is computed more accurately with the use of the excursion-set formalism, the mass of PBHs produced from the metric preheating instability, at which the PBH mass function peaks, is higher than the PBH mass quoted above and can be as high as the solar mass, lying in this way within the mass range probed by LISA.

<figure><img src="image_12.png" /><figcaption>Fig. 10 The matter power spectrum at the end of the preheating instability phase. The energy scales at the end of inflation and at the end of the preheating instability phase are fixed as $\rho_{\text{inf}}^{1/4} = 10^{14}$ GeV and $\rho_{\text{pl}}^{1/4} = 10^{7}$ GeV Image reproduced with permission</figcaption></figure>
At this point, one should also stress that ultralight PBHs, such as the ones produced from the preheating instability, can be potentially detected by LISA through the stochastic background of GWs induced at second order from PBH isocurvature perturbations, as argued in Papanikolaou et al. (2021). For these GWs, the peak frequency could lie within the LISA band despite the smallness of the PBH mass.

## 2.6 Early matter era

An early prolonged period of matter domination – after the end of inflation – could arise due to the domination of moduli fields which are a natural prediction of any supergravity or superstring-inspired theory of the early Universe. Since the curvature threshold for the PBH collapse depends on the equation of state of the Universe, it might be easier to obtain the same abundance of PBH during the moduli dominated epoch from a peak in the inflationary power spectrum less strong than the typical $O(10^{-2})$ required in the standard scenario of radiation domination. Several mechanisms have been proposed to lead to such a production during an early matter dominated era, see e.g. Khlopov and Polnarev (1980), Khlopov et al. (1985), Green et al. (1997), Nayak and Singh (2012), Harada et al. (2016), Carr et al. (2017, 2018), Kokubu et al. (2018), Ballesteros et al. (2020b), de Jong et al. (2022), De Luca et al. (2022a). In such scenarios, the mass and fraction of PBHs produced from a given inflationary potential depend on whether the moduli dominated epoch is immediate or delayed after the end of inflation. Recently, it has also been pointed out that even if PBHs are formed during radiation domination, the presence of a reheating epoch can still crucially affect the PBH abundance (Cai et al. 2018; Bhaumik and Jain 2020).

In models of inflation based on moduli stabilisation in string theory compactifications, for example, the inflaton may be identified with a modulus field which is only gravitationally coupled to the visible sector, leading to a slow-reheating phase, where the inflaton decay is via Planck-suppressed operators, and hence a prolonged duration of matter-domination after inflation (Green et al. 1997; Ballesteros et al. 2020b). Ballesteros et al. (2020b) finds that PBH formation is most effective for reheating temperatures $\lesssim 10^6$ GeV, in which case PBHs can form a significant fraction of the dark matter density today when primordial curvature fluctuations from inflation are of order $10^{-4}$ in a matter-dominated era, rather than $\sim 10^{-2}$ usually required in a radiation-dominated era.

## 2.7 Phase transitions

The formation and abundance of PBHs can also be affected by phase transitions in the early Universe. During the radiation era, the formation of PBHs may be facilitated if pre-existing super-horizon density fluctuations enter the cosmological horizon during a first order phase transition which proceeds in approximate equilibrium (Khlopov et al. 1985; Jedamzik and Niemeyer 1999). During this phase transition the speed of sound tends to zero and the pressure response of the fluid
vanishes and does not counter-balance the gravitational collapse of horizon-sized
primordial over-densities. Detailed numerical simulations using general-relativistic
hydrodynamics show that the overdensity threshold required for PBH formation by
fluctuations that enter the cosmological horizon around the phase transition is smaller
than the ordinary value for PBH formation during the radiation epoch. This would
produce an increase in the PBH mass spectrum around the corresponding horizon
mass in which the phase transition occurs.

It has been argued that the QCD phase transition could produce an enhancement of the population of PBH's around $1M_{\odot}$ (Jedamzik 1997), but lattice QCD calculations have shown that the QCD phase transition in the Standard Model is not first order but a continuous cross-over (Bhattacharya et al. 2014). Nevertheless, it is possible to collapse a significant amount of PBHs at the QCD transition if there are large fluctuations and radiation pressure drops as the relativistic quark and gluon degrees of freedom disappear from the plasma (Byrnes et al. 2018; Carr et al. 2021b; Jedamzik 2021).

Alternative scenarios have been proposed, where the dynamics of the QCD phase
transition incorporates extra degrees of freedom, like lepton number viola-
tion (Bödeker et al. 2021) and solitosynthesis (García-Bellido et al. 2021b). In other
models, the QCD transition could be first order if the number of light quarks were
larger than three at the confinement. Models in which the dynamics of a light scalar
field have been tailored to suppress quark masses before the QCD transition, while
recovering the measured values after the transition, have been developed in Davoudi-
asl (2019).

Another possibility is that the collision of bubbles formed during a first order
phase transition in the early Universe could create either remnant regions dominated
by the false vacuum or a sufficiently large local concentration of the energy stored in
the walls, or even other more complex processes, to cause gravitational collapse and
the formation of a primordial black hole (Hawking et al. 1982; Liu et al. 2022; Jung
and Okui 2021; Baker et al. 2021b; Flores et al. 2024; Lewicki et al. 2023). This
process could be affected by spatial clustering properties of bubble generated in high-
energy phase transitions (Pirvu et al. 2022; De Luca et al. 2021c). Recent
models (Baker et al. 2021c) consider the possibility that there may be particles whose
mass increases significantly during the phase transition, suppressing the transmission
of the corresponding particles through the advancing bubble walls. This effect can
build a sufficiently large overdensity in front of the walls that collapse into a
primordial black hole. The primordial black hole density and mass distribution
depend on the model parameters. Other models (Kawana and Xie 2022) consider
fermion species that gain big masses in the true vacuum, so that the corresponding
particles get trapped in the false vacuum as they do not have sufficient energy to
penetrate the bubble wall. After the first order phase transition, the fermions are
compressed into the false vacuum remnants to form Fermi-balls that finally collapse
to PBH due to the Yukawa attractive force.

Higher-order phase transitions periods in the early Universe may also lead to the
formation of PBHs. For example, domain walls can be generated during a non-
equilibrium second order phase transition in the vacuum state of a scalar field with a
flat direction during inflation (Rubin et al. 2000). The background de Sitter
fluctuations of such an effectively massless scalar field could provide non-equilibrium redefinition of correlation length and give rise to the islands of one vacuum in the sea of another one. After the phase transition takes place in the Friedman-Robertson-Walker (FRW) epoch, the two vacua are separated by a wall. Some of these closed and sufficiently large walls could accumulate enough energy to collapse and form a massive black hole. The mass spectrum of the PBHs which can be created in such a way depends on the details of the potential of the scalar field that has the flat direction during inflation and triggers the phase transition.

## 2.8 Domain walls and cosmic strings

Topological defects in the early Universe may also lead to the production of PBHs. Hawking (1989) was the first to show that PBHs could form from the collapse of cosmic string loops. He argued that loops that collapse by a factor $(G\mu)^{-1}$, where $\mu$ is the mass per unit length of the strings, will inevitably collapse to form black holes. PBH production can be used to place bounds on the cosmic string network (Hawking 1990; Polnarev and Zembowicz 1991; Caldwell and Gates 1993; Garriga and Sakellariadou 1993; Caldwell and Casper 1996; MacGibbon et al. 1998; Helfer et al. 2019; James-Turner et al. 2020) or cosmic strings nucleated during inflation (Garriga and Vilenkin 1993). More recently it has been argued that cosmic string cusps can collapse gravitationally into PBHs with a mass function that could extend up to stellar-masses (Jenkins and Sakellariadou 2020), although this claim has been disputed in Blanco-Pillado et al. (2021a).

Moreover, domain walls, which arise in theories with a broken discrete symmetry, can also form PBHs when such domain walls are smaller than a critical radius fall within the cosmological horizon and collapse due to their own tension (Ferrer et al. 2019) or are produced by tunnelling during inflation (Liu et al. 2020). In such scenarios, the mass function of PBHs in general has a spikelike structure. Moreover, in QCD axion models, PBHs can form either by the collapse of long-lived string-domain wall networks or from the collapse of closed domain walls (Ferrer et al. 2019; Ge 2020). The mass and abundance of PBHs formed in these scenarios crucially depends on the QCD axion mass. However, interestingly, various observational constraints on the PBH abundance in turn limit the QCD axion parameter space.

## 2.9 Primordial magnetic fields

Primordial magnetic fields generated in the early Universe are considered one of the feasible candidates to provide the required seeds for the observed large-scale intergalactic magnetic fields. In addition, such primordial fields also generally induce an anisotropic stress that can act as a source for the evolution of the super-horizon curvature perturbations, leading to the formation of PBHs with a broad possible range of masses. Since a large abundance of PBHs requires a large amplitude of density perturbations, the strength of these primordial magnetic fields on small scales should also be large enough for this mechanism to work (Saga et al. 2020). Based on
this idea, strong constraints on the amplitude of primordial magnetic fields have been
established, a few hundreds nano-Gauss on scales in the wave-number range
$10^2 \text{ Mpc}^{-1} \le k \le 10^{18} \text{ Mpc}^{-1}$.

Alternatively, it has been proposed that PBHs can act as seeds of cosmic magnetic
fields due to their accretion disks, either due to a *Biermann battery* (Safarzadeh et al.
2018; Papanikolaou and Gourgouliatos 2023b) mechanism or more exotic mecha-
nisms such as monopoles (Maldacena 2021; Araya et al. 2021) or Kerr–Newman
PBHs (Hooper et al. 2023). In particular, in Papanikolaou and Gourgouliatos (2023b)
it was shown that locally isothermal disks around supermassive PBHs with masses
$M > 10^{10} M_{\odot}$, can generate a seed primordial magnetic field of the order of $10^{-30} \text{g}$,
which can be later amplified by various dynamo instability processes and provide the
seed for the present day magnetic field in intergalactic scales. Interestingly enough,
this population of magnetised supermassive PBHs can induce a stochastic GW
background which can be used to set tight constraints on the supermassive PBH
abundances (see Papanikolaou and Gourgouliatos 2023a for more details).

**2.10 Summary**

In this section, we have explored multiple mechanisms described in the literature for
generating primordial black holes, from gravitational collapse of primordial
fluctuations and inflation to preheating, phase transitions, collisions of bubble walls
and cosmic strings. Each of these mechanisms has its own mass spectrum and
individual characteristics, which makes it possible to differentiate between them. The
most developed one is that of primordial fluctuations, where specific models exist
that can be put to test by present and future observations. Other mechanisms have
also been explored and predictive models exist, but are not yet sufficiently developed
to compare with observations. But we also emphasized the possible role of non-
Gaussian tails in the distribution of curvature fluctuations, which arise in models with
an ultra-slow-roll phase of inflation or in models with a stochastic spectator field.
This is an interesting, recent and active research area since this could soften or even
resolve the fine-tuning issues, which have recently been quantified in Cole et al.
(2023b), associated to the formation of PBHs with significant abundances. In the next
Sections, we will explore the dynamics of gravitational collapse and how it depends
on the thermal history of the Universe, with its impact on the PBH mass function,
their spin distribution and merging rates.

**3 PBH formation, mass function and clustering**

In this section we will concentrate on the concrete PBH formation mechanism based
on gravitational collapse of large density or curvature fluctuations, and derive
expressions for the mass function and clustering of PBH formed this way.
## 3.1 From curvature perturbations to density contrast

PBHs form after cosmological perturbations re-enter the cosmological horizon. Assuming spherical symmetry on superhorizon scales, the local region of the Universe describing such perturbations has the following asymptotic form of the metric

$$ds^2 = -dt^2 + a^2(t) \left[ \frac{dr^2}{1 - K(r)r^2} + r^2 d\Omega^2 \right] = -dt^2 + a^2(t) e^{2\zeta(\hat{r})} [d\hat{r}^2 + \hat{r}^2 d\Omega^2], \quad (26)$$

where $a(t)$ is the scale factor, while $K(r)$ and $\zeta(\hat{r})$ are the conserved comoving curvature perturbations, defined on a super-Hubble scale, converging to zero at infinity where the universe is taken to be unperturbed and spatially flat. The equivalence between the radial and the angular parts of these two forms of the metric gives

$$\left\{ \begin{array}{l} r = \hat{r} e^{\zeta(\hat{r})}, \\ \\ \frac{dr}{\sqrt{1 - K(r)r^2}} = e^{\zeta(\hat{r})} d\hat{r} \end{array} \right. \qquad (27)$$

and the difference between the two Lagrangian coordinates $r$ and $\hat{r}$ is related to the parametrisation of the comoving coordinate, fixed by the curvature perturbation chosen in the metric (i.e. $K(r)$ or $\zeta(\hat{r})$). From a geometrical point of view, the coordinate $\hat{r}$ considers the perturbed region as a local FLRW separated universe with the curvature perturbation $\zeta(\hat{r})$ modifying the local expansion, while the curvature profile $K(r)$ is defined with respect to the background FLRW solution ($K = 0$) and measures more directly the spatial geometry of space-time.

On superhorizon scales, where the curvature profile is time independent, these two are related as

$$K(r)r^2 = -\hat{r}\zeta'(\hat{r})[2 + \hat{r}\zeta'(\hat{r})]. \quad (28)$$

In this regime, we can use the gradient expansion approach (Shibata and Sasaki 1999; Tomita 1975; Salopek and Bond 1990; Polnarev and Musco 2007), based on expanding the time dependent variables, like energy density and velocity profiles, as power series of a small parameter $\epsilon \ll 1$ up to the first non zero order, where $\epsilon$ is conveniently identified with the ratio between the Hubble radius and the length scale of the perturbation. Although this approach reproduces the time evolution of the linear perturbation theory, it also allows having non linear curvature perturbations if the spacetime is sufficiently smooth on the scale of the perturbation (see Lyth et al. 2005). This is equivalent to saying that in this regime pressure gradients are negligible, not playing an important role in the perturbation evolution, which grows with the Universe expansion in a self similar way.
In this approximation, the energy density profile can be written as<fnref n="3" /> (Yoo et al. 2018; Musco 2019)

$$
\begin{align}
\frac{\delta \rho}{\rho_b} &\equiv \frac{\rho(r, t) - \rho_b(t)}{\rho_b(t)} = \frac{1}{a^2 H^2} \frac{3(1+w)}{5+3w} \frac{[K(r) r^3]'}{3r^2} \tag{29} \\
&= -\frac{1}{a^2 H^2} \frac{4(1+w)}{5+3w} e^{-5\zeta(\hat{r})/2} \nabla^2 e^{\zeta(\hat{r})/2}, \nonumber
\end{align}
$$

where $H(t) = \dot{a}(t)/a(t)$ is the Hubble parameter. The parameter $w$ is the coefficient of the equation of state $p = w\rho_b$, relating the total (isotropic) pressure $p$ to the total background energy density $\rho_b$.

In the linear regime of curvature perturbations, when $\zeta \ll 1$, this expression can be written as

$$
\frac{\delta \rho}{\rho_b} \approx - \frac{1}{a^2 H^2} \frac{2(1+w)}{5+3w} \nabla^2 \zeta(\hat{r}) = - \frac{k^2}{a^2 H^2} \frac{2(1+w)}{5+3w} \zeta_k, \quad (30)
$$

where the second equality is obtained from the Fourier transformation $\nabla^2\zeta = k^2\zeta_k$,
showing that in the linear regime there is a simple one to one mapping between the
real space, where the perturbations collapses, and the Fourier space where the power
spectrum of cosmological perturbations is defined. If the curvature perturbation $\zeta$ is a
Gaussian variable, the density contrast $\delta\rho/\rho_b$ has also a Gaussian distribution within
the linear regime described by Eq. (30). However, because the amplitude of the
threshold $\delta_c$ for PBHs is non linear, the linear regime does not give an accurate
description of the statistics of the density contrast, and the mapping between the
Fourier space and the real space requires a more elaborate approach. As we will see
later, this allows to compute the threshold $\delta_c$ from the shape of the power spectrum
$P_\zeta$.

For this reason, the definition of the threshold in terms of $\zeta_c$ computed from the
linear approximation given by Eq. (30), does not give an accurate description and
should be dropped in favour of $\delta_c$ which is clearly defined in the following
section. We anticipate that the nonlinear relation between the density contrast and the
comoving curvature perturbations causes the former to obey non-Gaussian statistics,
even in absence of primordial non-Gaussianity. Computations based on peak theory
and on threshold statistics show that this unavoidable non-Gaussianity makes the
production of PBHs more difficult. In case of a peaked power spectrum, this effect is
compensated by an increase of the amplitude by an order 2–3 factor (De Luca et al.
2019c; Young et al. 2019).

**3.2 Perturbation amplitude**

To define the amplitude of a cosmological perturbation, it is useful to introduce the
compaction function (Shibata and Sasaki 1999; Musco 2019), defined as

<footnote n="3">³ K'(r) denotes differentiation with respect to r while ζ'(r̂) and ∇²ζ(r̂) denote differentiation with respect to r̂.</footnote>
$$C \equiv 2 \frac{\delta M(r, t)}{R(r, t)}, \qquad (31)$$

where $R(r, t)$ is the areal radius and $\delta M(r, t)$ is the difference between the Misner–Sharp mass within a sphere of radius $R(r, t)$, and the background mass $M_b(r, t) = 4\pi\rho_b(r, t)R^3(r, t)/3$ within the same areal radius but calculated with respect to a spatially flat FLRW metric. In the superhorizon regime (i.e. $\epsilon \ll 1$) the compaction function is time independent, and directly related to the curvature profile as

$$C(r) = \frac{3(1+w)}{5+3w} K(r) r^2, \qquad (32)$$

where this expression can be written in terms of $\zeta(\hat{r})$ using Eq. (28).

The comoving length scale of the perturbation can be conveniently identified with the distance where the compaction function reaches its peak (which is a maximum for a positive perturbation and a minimum for a negative one) when $r = r_m$ (Musco 2019)(i.e. $C'(r_m) = 0$), characterized by a constraint relation in terms of the curvature profile

$$K(r_m) + \frac{r_m}{2} K'(r_m) = 0, \quad \text{or} \quad \zeta'(\hat{r}_m) + \hat{r}_m \zeta''(\hat{r}_m) = 0. \qquad (33)$$

Given the curvature profile, the parameter $\epsilon$ of the gradient expansion is therefore defined as

$$\epsilon \equiv \frac{R_H(t)}{R_b(r_m, t)} = \frac{1}{aHr_m} = \frac{1}{aH\hat{r}_m e^{\zeta(\hat{r}_m)}}, \qquad (34)$$

where $R_H = 1/H$ is the cosmological horizon and $R_b(r, t) = a(t)r$ is the background component of the areal radius. With these definitions, the expression written in Eq. (29) is valid with good accuracy for $\epsilon \ll 1$.

According to this, we can now consistently define the perturbation amplitude as being the mass excess of the energy density within the scale $r_m$ when $\epsilon = 1$ ($aH\hat{r}_m = 1$). This corresponds to defining the horizon crossing time $t_H$ with a linear extrapolation from the superhorizon regime. Although this is not very accurate, it provides a well defined criterion to measure and compare consistently the amplitude of different perturbations, understanding how the threshold is varying with respect to the initial curvature profiles.

The amplitude of the perturbation measured at $t_H$, which we refer to as $\delta_m \equiv \delta(r_m, t_H)$, is given by the excess of mass averaged over a spherical volume of radius $R_m$, defined as

$$\delta_m = \frac{4\pi}{V_{R_m}} \int_0^{R_m} \frac{\delta\rho}{\rho_b} R^2 dR = \frac{3}{r_m^3} \int_0^{r_m} \frac{\delta\rho}{\rho_b} r^2 dr, \qquad (35)$$

where $V_{R_m} = 4\pi R_m^3/3$. The second equality is obtained by neglecting the higher order terms in $\epsilon$, approximating $R_m \simeq a(t)r_m$, which allows to simply integrate over the comoving volume of radius $r_m$. Inserting the expression for $\delta\rho/\rho_b$ given by
Eq. (29) into Eq. (35), one obtains $\delta_m = C(r_m)$, and with simple calculation seen in Musco (2019), one obtains the fundamental relation

$$ \delta_m = 3 \frac{\delta \rho}{\rho_b} (r_m, t_H), \qquad (36) $$

showing that at $r_m$ the perturbation amplitude is not affected by the location where this is measured.

### 3.3 Threshold for PBH formation

PBHs form when the perturbation amplitude $\delta_m > \delta_c$, where the value of the threshold $\delta_c$ depends on the shape of the energy density profile and the equation of state. For a radiation-dominated Universe ($p = \frac{1}{3}\rho$), it has been found that $2/5 \le \delta_c \le 2/3$, with the shape of the cosmological perturbation simply characterised by one dimensionless parameter corresponding to the width of the peak of the compaction function (Musco 2019; Escrivà et al. 2020)

$$ \alpha = - \frac{C''(r_m) r_m^2}{4C(r_m)}. \qquad (37) $$

The shape of the energy density contrast becomes peaked for $\alpha < 1$ corresponding to a broad profile of the compaction function, whereas the shape of the compaction function $C$ is more peaked for values of $\alpha > 1$, corresponding to broad profiles of the energy density contrast (Musco et al. 2021). See e.g. Fig. 2.1 of Franciolini (2021) for a graphical representation of this scaling. This allows to compute the threshold $\delta_c$ as a function of the shape parameter $\alpha$, with an analytic expression up to a few percent precision (Escrivà et al. 2020)

$$ \delta_c \approx \frac{4}{15} e^{-\frac{1}{2}} \frac{\alpha^{1-\frac{5}{2\alpha}}}{\Gamma(\frac{5}{2\alpha}) - \Gamma(\frac{5}{2\alpha}, \frac{1}{\alpha})}, \qquad (38) $$

where $\Gamma$ identifies the special Gamma-functions. This is consistent with the numerical analysis made in Musco (2019), where it was shown that the effects on the threshold of additional parameters are negligible, because the collapse in spherical symmetry is basically characterized by the configuration of the region within the forming apparent horizon ($r \le r_m$), whose initial configuration is fully described by the shape parameter $\alpha$.

The mass of a PBH is determined by the amplitude $\delta$ of the perturbation with respect to the corresponding threshold $\delta_c$, according to the scaling law of critical collapse

$$ m_{\text{PBH}} = \mathcal{K}(\delta - \delta_c)^{\eta} M_H \qquad (39) $$

where the critical exponent $\eta$ depends only on the equation of state ($\eta \simeq 0.36$ for a radiation-dominated Universe), while $\mathcal{K}$, as for $\delta_c$, depends on the initial configuration of the energy density profile, roughly varying between 1 and 10. For the standard scenario of a radiation-dominated Universe, expression (39) describes the
mass with good accuracy when $(\delta - \delta_c) \lesssim 10^{-2}$, corresponding to $m_{\text{PBH}} \approx \gamma M_{\text{H}}$, where $\gamma$, varying roughly between 0.1 and 1, is an efficiency parameter encapsulating the details of the gravitational collapse.

Although non linear cosmological density perturbations are described by a non Gaussian random field, when the Universe is still radiation-dominated one can compute the threshold $\delta_c$ from the shape of the Gaussian inflationary power spectrum following a simple prescription. The algorithm, divided into a few simple steps, accounts for both the non linear effects associated with the relation between the Gaussian curvature perturbation and the density contrast as well as for those arising at horizon crossing. While a more refined description of the various steps will be found in Musco et al. (2021), we give here a synthetic overview of the prescription one needs to follow.

1. The power spectrum of the curvature perturbation: take the primordial power spectrum $P_{\zeta}$ of the Gaussian curvature perturbation and compute, on superhorizon scales, its convolution with the transfer function $T(k, \eta)$

$$P_{\zeta}(k, \eta) = \frac{2\pi^2}{k^3} P_{\zeta}(k) T^2(k, \eta). \quad (40)$$

2. The comoving length scale $\hat{r}_m$ of the perturbation is related to the characteristic scale $k_*$ of the power spectrum $P_{\zeta}$. Compute the value of $k_*$ by solving the following integral equation

$$\int dk k^2 \left[ (k^2 \hat{r}_m^2 - 1) \frac{\sin(k \hat{r}_m)}{k \hat{r}_m} + \cos(k \hat{r}_m) \right] P_{\zeta}(k, \eta) = 0. \quad (41)$$

3. The shape parameter: compute the corresponding shape parameter $\alpha$ of the collapsing perturbation, including the correction from the non linear effects, by solving the following equation

$$F(\alpha)[1 + F(\alpha)]\alpha = 2\alpha_G \quad (42)$$

with

$$\alpha_G = -\frac{1}{4} \left[ 1 + \hat{r}_m \frac{\int dk k^4 \cos(k \hat{r}_m) P_{\zeta}(k, \eta)}{\int dk k^3 \sin(k \hat{r}_m) P_{\zeta}(k, \eta)} \right], \\ F(\alpha) = \sqrt{1 - \frac{2}{5} e^{-\frac{1}{\alpha}} \frac{\alpha^{1-\frac{5}{2\alpha}}}{\Gamma(\frac{5}{2\alpha}) - \Gamma(\frac{5}{2\alpha}, \frac{1}{\alpha})}} \quad (43)$$

A numerical fit of the shape parameter $\alpha$ as a function of $\alpha_G$, which represents the solution of Eq. (42) with a percent accuracy, is given by (Franciolini 2021)

$$\alpha \approx \begin{cases} 1.758 \alpha_G^{2.335} + 1.912 \alpha_G & 0.1 \lesssim \alpha \lesssim 4.5, \\ 4 \alpha_G^2 + 3.930 \alpha_G & \alpha \gtrsim 8. \end{cases} \quad (44)$$

In the intermediate region where $4.5 < \alpha < 8$, no simple power law fit with a percent accuracy was found, and one needs to solve Eq. (42) numerically.
4. The threshold $\delta_c$: compute the threshold as function of $\alpha$, fitting the numerical simulations.

* At superhorizon scales making a linear extrapolation at horizon crossing ($aH_{r_m} = 1$).

$$ \delta_c \approx \begin{cases} \alpha^{0.047} - 0.50 & 0.1 \lesssim \alpha \lesssim 7 \\ \alpha^{0.035} - 0.475 & 7 \lesssim \alpha \lesssim 13 \\ \alpha^{0.026} - 0.45 & 13 \lesssim \alpha \lesssim 30 \end{cases} \qquad (45) $$

* At horizon crossing taking into account also the non linear effects.

$$ \delta_c \approx \begin{cases} \alpha^{0.125} - 0.05 & 0.1 \lesssim \alpha \lesssim 3 \\ \alpha^{0.06} + 0.025 & 3 \lesssim \alpha \lesssim 8 \\ 1.15 & \alpha \gtrsim 8 \end{cases} \qquad (46) $$

This prescription is using the fact that, if $P_ζ$ follows a Gaussian distribution, $\Phi_L \equiv -\hat{r}_ζ'(\hat{r})$ is also a Gaussian variable, and one can write $\delta_m$ as

$$ \delta_c = \frac{4}{3} \Phi_{c,L} \left( 1 - \frac{1}{2} \Phi_{c,L} \right) \qquad (47) $$

where $\Phi_{m,L} = \Phi_L(\hat{r}_m)$, and $\Phi$ is such that $0.37 \lesssim \Phi_{c,L} \le 1$.

In Fig. 11, we show the threshold density contrast as a function of the shape parameter (left panel). Also, in the right panel, we report the critical collapse parameter $K$ defined in Eq. (39). We recall that the critical collapse exponent $\eta$ is independent from the shape of the collapsing overdensity peak (i.e. independent from $\alpha$) and it is fixed to be $\eta = 0.358$ in a perfect radiation fluid. In the following, we will discuss how these parameters are modified when thermal effects during the evolution of the universe, causing a significant departure from the perfect radiation fluid approximation.

Before concluding this Section, we mention two relevant points. First, it has been recently observed that corrections from the non-linear radiation transfer function and the determination of the true physical horizon crossing tend to decrease the PBH abundance (De Luca et al. 2023b). Secondly, a recent study (Musco and Papanikolaou 2022) investigated the formation of PBHs with an anisotropic relativistic perfect fluid, using a covariant form of the EoS in terms of pressure and energy density gradients. This indicates that the value of the threshold $\delta_c$ is increasing with the amplitude of the anisotropy, if this is not large, while further studies are necessary for a nonlinear modification of the EoS.
<figure><img src="image_13.png" /><figcaption>Fig. 11 Threshold $\delta_c$ and critical collapse parameter $\mathcal{K}$ as a function of the shape parameter $\alpha$, assuming a universe dominated by perfect radiation (i.e. $w = 1/3$). Figure adapted from Musco et al. (2021) Image reproduced with permission</figcaption></figure>

## 3.4 Gravitational collapse of primordial curvature perturbations – standard formalism

Although PBHs can be generated by many different mechanisms, we will focus here on the collapse of cosmological perturbations. In this case, the fraction $\beta_f$ of PBHs formed in the early Universe is determined by the probability that a given primordial curvature fluctuation $\zeta$, characterized by an amplitude $\delta$ of the density contrast, is above a certain threshold $\delta_c$. The fraction of PBHs at formation is then usually given by

$$ \beta_f(m_{\text{PBH}}) = \int_{\delta_c}^{\infty} P(m_{\text{PBH}}, \delta) d\delta, \qquad (48) $$

and the abundance of PBHs is sensitive to the PDF $P(m_{\text{PBH}}, \delta)$ as to the value of the threshold $\delta_c$. In particular, the PDF is characterized by the physics in the early Universe while the threshold $\delta_c$ depends on the conditions at the time of formation. In particular, $\delta_c$ has been actively studied in the past years (Musco et al. 2005; Polnarev and Musco 2007; Kopp et al. 2011; Harada et al. 2013; Young et al. 2014), with recent analysis on the dependence of the threshold on the shape of the power spectrum (Yoo et al. 2018; Germani and Musco 2019; Musco 2019), on the non-linear relation between the curvature perturbations and the density contrast (Young et al. 2019; De Luca et al. 2019c) as well as on the equation of state (EoS) (Harada et al. 2013; Escrivà et al. 2021; Papanikolaou 2022b) with a particular focus on the QCD phase transition (Jedamzik 1997; Byrnes et al. 2018).

Note also that in Eq. (48), we have expressed the PDF in terms of the mass $m_{\text{PBH}}$ of the PBHs formed instead of the number of e-folds because in first approximation the mass of PBHs is related to the size of the causal horizon collapsing, depending on the time of formation (See Eq. 39).

Given the fraction of PBHs $\beta_f(M)$, one can compute the contribution of PBHs to the energy density of the Universe and determine the fraction of DM they represent today. In this sense it is convenient to compute this quantity at the time of matter-radiation equality
$$
\Omega_{\text{PBH}}^{\text{eq}} = \int_{M_{\text{ev}}}^{M_{\text{eq}}} \beta_{\text{eq}}(m_{\text{PBH}}) \, \text{d} \ln m_{\text{PBH}} \quad (49)
$$

where $M_{\text{ev}}$ is the lower bound due to Hawking evaporation, i.e. PBHs with $M < M_{\text{ev}}$
would already been evaporated by now. We assume that the fraction of PBH at the
time of matter-radiation equality has grown because of the cosmic expansion fol-
lowed after the time of formation, which gives

$$
\beta_{\text{eq}}(m_{\text{PBH}}) = e^{(N_{\text{eq}} - N_f)} \beta_f(m_{\text{PBH}}). \quad (50)
$$

Large curvature fluctuations can be produced by very different means in the early
Universe. If the curvature perturbations were Gaussian, then the PDF would be
described only by its variance $\sigma^2$ and the fraction of PBHs measured at the time of
formation is usually computed as

$$
\beta_f(m_{\text{PBH}}) = 2 \int_{\delta_c}^{\infty} \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{\delta^2}{2\sigma^2}} d\delta = 1 - \operatorname{erf}\left(\frac{\delta_c}{\sqrt{2}\sigma}\right) = \operatorname{erfc}\left(\frac{\delta_c}{\sqrt{2}\sigma}\right), \quad (51)
$$

where the factor 2 outside the integral is sometimes introduced when the abundance
is computed using the excursion set method (see, e.g., Moradinezhad Dizgah et al.
2019; De Luca et al. 2020f). This expression, however, does not take into account
that when $\delta$ is larger then a certain value $\delta_{\max}$ (originally estimated $\delta_{\max} \sim 1$ by Carr
1975), the perturbation forms a separate closed universe, topologically disconnected.
A more accurate version can nevertheless be obtained re-normalizing the previous
expression as

$$
\beta_f(m_{\text{PBH}}) = \frac{\int_{\delta_c}^{\delta_{\max}} \frac{1}{\sqrt{2\pi}\sigma^2} e^{-\frac{\delta^2}{2\sigma^2}} d\delta}{\int_{0}^{\delta_{\max}} \frac{1}{\sqrt{2\pi}\sigma^2} e^{-\frac{\delta^2}{2\sigma^2}} d\delta} = 1 - \frac{\operatorname{erf}\left(\frac{\delta_c}{\sqrt{2}\sigma}\right)}{\operatorname{erf}\left(\frac{\delta_{\max}}{\sqrt{2}\sigma}\right)} = \frac{\operatorname{erfc}\left(\frac{\delta_c}{\sqrt{2}\sigma}\right) - \operatorname{erfc}\left(\frac{\delta_{\max}}{\sqrt{2}\sigma}\right)}{1 - \operatorname{erfc}\left(\frac{\delta_{\max}}{\sqrt{2}\sigma}\right)}. \quad (52)
$$

PBHs are produced through the tail of the PDF, ie. the area left under the PDF curve where $\delta$ is larger than the critical threshold $\delta_c$. Therefore, PBH abundance is exponentially sensitive to the threshold and the statistical properties of the primordial perturbations. Even tiny modifications in the tail of the distribution can change the PBH abundance by many orders of magnitude.

The variance of the field of density perturbations $\sigma$, according to the Gaussian distribution of $\delta$, consistent with the linear approximation of the relation between the energy density contrast and the curvature perturbation (see next subsection), is given by

$$
\sigma^2 = \langle \delta^2 \rangle = \int_0^\infty \frac{dk}{k} \mathcal{P}_\delta(k, r) = \frac{16}{81} \int_0^\infty \frac{dk}{k} (kr)^4 \tilde{W}^2(k, r) T^2(k, r) \mathcal{P}_\zeta(k), \quad (53)
$$

where $\mathcal{P}_\delta(k, r)$ and $\mathcal{P}_\zeta(k)$ are the density and the curvature power spectra, while
$\tilde{W}(k, r)$ is the Fourier transform of the top-hat smoothing function and $T(k, r)$ is the linear transfer function

$$ \tilde{W}(k, r) = 3 \left[ \frac{\sin(kr) - kr\cos(kr)}{(kr)^3} \right] \qquad (54) $$

$$ T(k, r) = 3 \left[ \frac{\sin\left(\frac{kr}{\sqrt{3}}\right) - \frac{kr}{\sqrt{3}}\cos\left(\frac{kr}{\sqrt{3}}\right)}{\left(\frac{kr}{\sqrt{3}}\right)^3} \right]. \qquad (55) $$

The quantities $\delta_c$ and $\sigma$, giving the critical threshold $v_c \equiv \delta_c/\sigma$, are typically computed at super horizon scale, when the curvature perturbation $\zeta$ for adiabatic perturbations is time independent, making a linear extrapolation at horizon crossing time, which is approximately the time when PBHs are formed (we will discuss this later in more details). All this shows that a larger amplitude of the power spectrum, combined with a shape decreasing the corresponding value of $v_c$, could increase significantly the fraction of PBHs because the abundance of PBHs is exponentially sensitive to the value of $\delta_c$.

The linear transfer function<fnref n="4" /> ensures that the variance would always converge, because the window function by itself, chosen consistently with the top-hat window function in real space used in the definition of $\delta_c$, does not guarantee convergence if the power spectrum considered is sufficiently broad. An alternative approach is to use a Gaussian window function, which provides always convergence of the variance, but this is contaminating the value of $\delta_c$, introducing an error in the final computation of the abundance (Young 2019).

## 3.5 The inevitable non-Gaussianity of the primordial black hole abundance

The non-linear relation between the curvature perturbation and the density contrast in (29) induces an inescapable non-Gaussianity of the density contrast, even assuming a Gaussian primordial initial condition (Musco 2019; Young et al. 2019; De Luca et al. 2019c). This can be captured in the non-linear expression (Young et al. 2019)

$$ \delta_m = \delta_l - \frac{3}{8}\delta_l^2, \qquad \delta_l = -\frac{4}{3}r_m\zeta'(r_m). \qquad (56) $$

It highlights two key points. First of all, the probability of forming PBHs does not depend on the comoving curvature perturbation itself $\zeta$, but on its derivative $\zeta'$. This is expected, given that on superhorizon scales one can always add or subtract to the comoving curvature perturbation a constant by a coordinate transformation and this may not influence any physical result. Secondly, thanks to the conservation of the probability

<footnote n="4">⁴ De Luca et al. (2023b) has recently investigated the role of non-linear corrections, coming from the non-linear radiation transfer function and the determination of the true physical horizon crossing, on the PBH abundance, showing that, while the critical threshold is larger than what routinely assumed, the variance is unlikely to be changed with respect to the linear one.</footnote>
$$P(\delta_1) = P[\delta_m(\delta_1)] \left| \frac{d\delta_m}{d\delta_1} \right|, \qquad (57)$$

one can compute the abundance of PBHs simply using the Gaussian probability of $\delta_1$, integrating it from the critical amplitude for $\delta_1$

$$\delta_{l,c} = \frac{4}{3} \left( 1 - \sqrt{1 - \frac{3}{2} \delta_c} \right), \qquad (58)$$

in terms of the critical amplitude $\delta_c$.

## 3.6 Collapse of some initial non-Gaussian perturbations

It has been observed that nearly all models discussed in Sect. 2 generate a considerable amount of primordial non-Gaussianity. We briefly summarize the $\beta$ function in case of large local non-Gaussianity, which turns the PDF into a $\chi^2$ distribution, explicitly $\delta_{\chi^2} = G^2 - \langle G^2 \rangle$; and also in case of large cubic interactions, which turn the PDF into a cubic-Gaussian distribution, explicitly $\delta_{G^3} = G^3$, with G a field obeying Gaussian statistics (Lyth 2012; Linde et al. 2013; Young and Byrnes 2013; Bugaev and Klimai 2014).

$$\begin{aligned} \beta_{f, \chi^2}(m_{\text{PBH}}) &= \text{erfc} \left( \sqrt{\frac{1}{2} + \frac{\delta_c}{\sqrt{2}\sigma}} \right), \\ \beta_{f, G^3}(m_{\text{PBH}}) &= \text{erfc} \left( \left( \frac{\delta_c}{\sqrt{\frac{8}{15}\sigma}} \right)^{1/3} \right). \end{aligned} \qquad (59)$$

We also note that the threshold density $\delta_c$ depends slightly on the non-Gaussianity (Kehagias et al. 2019; Kitajima et al. 2021; Escrivà et al. 2022).

## 3.7 Non-perturbative abundance from non-Gaussian perturbations

Ferrante et al. (2023b) and Gow et al. (2023) addressed the problem of computing the impact of a general form of local primordial NG on PBH abundance beyond the perturbative approach. In particular, let us focus on the functional form

$$\zeta(\vec{x}) = F(\zeta_G(\vec{x})), \qquad (60)$$

where F is a generic non-linear function of the Gaussian component $\zeta_G$. There is a large number of relevant cases in the literature of PBH formation where the NGs can be modelled as in Eq. (60). We summarize few relevant examples in Table 1 (adapted from Ferrante et al. 2023b).

It is important to stress here that going beyond the perturbative approach becomes a necessity at least when dealing with broad spectra of curvature perturbations. This is because, as shown in Ferrante et al. (2023b), in these cases the probability
distribution function of $\zeta_G$ is evaluated beyond the radius of convergence of the series expansion, leading to inaccurate results if the series is truncated at finite order. Notice that in non-attractor models the sign of non-Gaussianity is always positive (Firouzjahi and Riotto 2023), thus enhancing the PBH abundance.

Using the relation between the curvature perturbation and the compaction function in the gradient expansion (29), it is possible write

$$ C(r) = -2\Phi r \zeta'(r) \left[1 + \frac{r}{2}\zeta'(r)\right] = C_1(r) - \frac{1}{4\Phi} C_1(r)^2, \quad C_1(r) \equiv -2\Phi r \zeta'(r), \qquad (61) $$

where $C_1(r)$ defines the so-called linear component of the compaction function and we introduced the prefactor $\Phi$ that accounts for potential thermal effects (see more in depth discussion below Eq. (70)). Substituting in the previous expression the relation $\zeta = F(\zeta_G)$ (see Eq. (60)), one obtains that the linear component of the compaction function takes the form

$$ C_1(r) = -2\Phi r \zeta'_G(r) \frac{dF}{d\zeta_G} = C_G(r) \frac{dF}{d\zeta_G} \qquad (62) $$

with $C_G(r) \equiv -2\Phi r \zeta'_G(r)$. Consequently, the compaction function reads

$$ C(r) = C_G(r) \frac{dF}{d\zeta_G} - \frac{1}{4\Phi} C_G^2(r) \left(\frac{dF}{d\zeta_G}\right)^2 . \qquad (63) $$

Crucially, $C(r)$ depends on both the Gaussian linear component $C_G$ and the Gaussian curvature perturbation $\zeta_G$. Both these random variables are Gaussian; $\zeta_G$ is Gaussian by definition while $C_G$ is defined by means of the derivative of the Gaussian variable $\zeta_G$.

Adopting threshold statistics on the compaction function, one can express the non-Gaussian mass fraction as an integration of the joined PDF in the domain which results in over-threshold perturbations, i.e.

<table><caption>**Table 1** Some examples of PBH formation scenarios characterised by local non-Gaussianity of the form defined in Eq. (60). The function $X(r_{\text{dec}}, \zeta_G)$ given in the second row is defined in (Sasaki et al. 2006)</caption><tr><td>Power-series expansion (Bugaev and Klimai 2013; Nakama et al. 2017; Byrnes et al. 2012; Young and Byrnes 2013; Yoo et al. 2018; Kawasaki and Nakatsuka 2019; Yoo et al. 2019; Riccardi et al. 2021; Taoso and Urbano 2021; Meng et al. 2022; Escrivà et al. 2022)</td><td>$\zeta = \zeta_G + \frac{3}{5}f_{NL}\zeta_G^2 + \frac{9}{25}g_{NL}\zeta_G^3 + \dots$</td></tr><tr><td>Curvaton (Sasaki et al. 2006; Pi and Sasaki 2023)</td><td>$\zeta = \log[X(r_{\text{dec}}, \zeta_G)]$</td></tr><tr><td>USR (Atal et al. 2019)</td><td>$\zeta = -(\frac{6}{5}f_{NL})^{-1}\log(1 - \frac{6}{5}f_{NL}\zeta_G)$</td></tr><tr><td>USR with an upward step (Cai et al. 2022a)</td><td>$\zeta = -\frac{2}{|h|}[\sqrt{1 - |h|\zeta_G} - 1]$</td></tr></table>
$$
\beta_{\text{NG}} = \int_{D} \mathcal{K}(C - C_{\text{th}})^{\gamma} \mathbf{P}_{G}(C_{G}, \zeta_{G}) \mathrm{d}C_{G} \mathrm{d}\zeta_{G}, \qquad (64)
$$

$$
\mathcal{D} = \{C_G, \zeta_G \in \mathbb{R} : C(C_G, \zeta_G) > C_{\text{th}} \wedge C_1(C_G, \zeta_G) < 2\Phi\}, \quad (65)
$$

where the multivariate normal distribution of ($C_G, \zeta_G$) can be written as

$$
\mathbf{P}_G(C_G, \zeta_G) = \frac{1}{(2\pi)\sigma_c \sigma_r \sqrt{1 - \gamma_{cr}^2}} \exp\left(-\frac{\zeta_G^2}{2\sigma_r^2}\right) \exp\left[-\frac{1}{2(1 - \gamma_{cr}^2)}\left(\frac{C_G}{\sigma_c} - \frac{\gamma_{cr}\zeta_G}{\sigma_r}\right)^2\right]. \quad (66)
$$

To shorten the notation, we have followed Young (2022) and defined the element of
the covariance matrix as

$$
\langle C_G C_G \rangle = \sigma_c^2 = \frac{4\Phi^2}{9} \int_0^\infty \frac{dk}{k} (kr_m)^4 \tilde{W}^2(k, r_m) T^2(k, r_m) P_\zeta(k), \quad (67)
$$

$$
\langle C_G \zeta_G \rangle = \sigma_{cr}^2 = \frac{2\Phi}{3} \int_0^\infty \frac{dk}{k} (kr_m)^2 \tilde{W}(k, r_m) W_s(k, r_m) T^2(k, r_m) P_\zeta(k), \quad (68)
$$

$$
\langle \zeta_G \zeta_G \rangle = \sigma_r^2 = \int_0^\infty \frac{dk}{k} W_s^2(k, r_m) T^2(k, r_m) P_\zeta(k), \quad (69)
$$

where $W_s(k,r) = \sin(kr)/kr$, $\tilde{W}(k,R)$ and $T(k,\tau)$ are given in Eqs. (55) and (54),
while $\gamma_{cr} \equiv \sigma_{cr}^2/\sigma_c\sigma_r$.

3.8 Thermal history

3.8.1 The equation of state and the number of relativistic degrees of freedom

The reheating at the end of inflation should have filled the Universe with radiation. In the absence of extensions beyond the Standard Model (SM) of particle physics, the Universe remains dominated by relativistic particles with an energy density decreasing as the fourth power of the temperature as the Universe expands and cools down. The number of relativistic degrees of freedom remains constant ($g_* = 106.75$) until around 200 GeV, when the temperature of the Universe falls to the mass thresholds of SM particles.

As shown in Fig. 12 (left panel), the first particle to become non-relativistic is the top quark at $T \simeq m_t = 172$ GeV, followed by the Higgs boson at 125 GeV, and the Z and W bosons at 92 and 81 GeV, respectively. These particles become non-relativistic at nearly the same time and this induces a significant drop in the number of relativistic degrees of freedom down to $g_* = 86.75$. There are further changes at the b and c quark and $\tau$-lepton thresholds but these are too small to appear in Fig. 12. Thereafter, $g_*$ remains approximately constant until the QCD transition at around 200 MeV, when protons and neutrons condense out of the free light quarks and gluons. The number of relativistic degrees of freedom then falls abruptly to
<figure><img src="image_14.png" /><figcaption>**Fig. 12** Relativistic degrees of freedom $g_*$ (Left panel) and equation-of-state parameter $w$ (Right panel), both as a function of temperature $T$ (in MeV). The grey vertical lines correspond to the masses of the electron, pion, proton/neutron, W, Z bosons and top quark, respectively. The grey dashed horizontal lines indicate values of $g_*=100$ and $w=1/3$, respectively. Image reproduced with permission</figcaption></figure>

$g_* = 17.25$. A little later the pions become non-relativistic and then the muons, giving $g_* = 10.75$. After that, $g_*$ remains constant until $e^+e^-$ annihilation and neutrino decoupling at around 1 MeV, when it finally drops to $g_* = 3.36$.

Whenever the number of relativistic degrees of freedom suddenly drops, it changes the effective equation of state parameter $w$. As shown in Fig. 12 (right panel), there are thus three main periods in the thermal history of the Universe when $w$ decreases. After each period, $w$ resumes its relativistic value of 1/3 but each sudden drop modifies the probability of gravitational collapse of any large curvature fluctuations present at that time.

### 3.8.2 Thermal effects on PBH formation and abundance

Because the threshold $\delta_c$ is a measure of the pressure gradients working against the gravitational collapse, it depends on the equation-of-state, and since the PBH abundance depends exponentially on the value of the threshold, the thermal history should have left imprints in the PBH mass distribution. In particular, assuming a nearly scale-invariant primordial power spectrum, one finds that the PBH mass distribution is broad with features in the form of a high peak at the solar-mass scale, a bump around $30M_\odot$, which would have implications for LVK observations, as well as one bump at planetary-masses from the electroweak scale and a final one for intermediate-mass PBHs (Byrnes et al. 2018; Carr et al. 2021b).

In order to gain an intuition on the relevance of these effects for on the mass distribution derived from (nearly) scale invariant power spectra, we show in Fig. 13 two representative models of PBH formation from power spectra with spectral index $n_s = 0.97$ and $n_s = 1$. These models are later used in this review for computing the PBH merger rates and SGWB. However, to obtain Fig. 13 the effects of the variation of the equation-of-state and of the pressure gradients *during* the PBH formation process, discussed later, were neglected.
<figure><img src="image_15.png" /><figcaption>**Fig. 13** PBH density fraction at formation $\beta^{\text{form}}$ (left panel) and the corresponding PBH mass function $f_{\text{PBH}}$ today (right panel), neglecting the effects of PBH growth by accretion and hierarchical mergers, for two models with a power-law primordial power spectrum and including the effects of thermal history: Model 1 from Carr et al. (2021b), Clesse and García-Bellido (2022) with spectral index $n_s = 0.97$; Model 2 from De Luca et al. (2021g), Byrnes et al. (2018) with $n_s = 1$, and a cut-off mass of $10^{-14} M_{\odot}$. The transition between the large-scale and small-scale power spectrum is fixed at $k = 10^3 \text{ Mpc}^{-1}$. The power spectrum amplitude is normalized such that both models produce an integrated PBH fraction $f_{\text{PBH}} = 1$, i.e. PBH constitute the totality of Dark Matter. A value of $\gamma = 0.8$ (ratio between the PBH mass and the Hubble horizon mass at formation) was assumed. Image reproduced with permission from Auclair et al. (2023), copyright by the author(s). Image reproduced with permission</figcaption></figure>

More recently, new numerical simulations of PBH formation at the QCD epoch have been performed (Escrivà et al. 2023; Musco et al. 2024), including the time and radial variations of the equation of state. These modify significantly the PBH mass function characterised by a lower but broader QCD-induced peak. The physical reason of this broadening is the time duration of the formation process, during which the changes in the equation-of-state favor the collapse. As a result, the threshold and the resulting PBH mass distribution do not only depend on the fluctuation mass at horizon crossing, but also on the precise evolution of the equation of state, as well as the details of the curvature fluctuation profile.

The thermal history, affecting the equation of state of the fluid dominating the energy density at the time of PBH collapse, impacts many aspects of the PBH formation. Following Franciolini et al. (2022e), Musco et al. (2024), here we discuss how the threshold is varying according to the equation of state, focusing on the modification of the relation between the density contrast and the curvature perturbation when the equation of state is varying with time. This has an impact also on the PBH masses with a modified scaling behaviour of the critical collapse.

When the equation-of-state is not constant, Eq. (29) is modulated by an overall factor $\Phi(t)$, which we can define in the long wavelength approximation as (Franciolini et al. 2022e)

$$ \frac{\delta \rho}{\rho_b} (r, t) = -\frac{4}{3} \Phi \left( \frac{1}{aH} \right)^2 e^{-5\zeta(r)/2} \nabla^2 e^{\zeta(r)/2}. \qquad (70) $$

$\Phi(t)$ can be computed by solving the equation (Polnarev and Musco 2007)
$$ \frac{1}{H} \frac{d\Phi(t)}{dt} + \frac{5 + 3w(t)}{2} \Phi(t) - \frac{3}{2}(1 + w(t)) = 0 \qquad (71) $$

integrated from past infinity to the time when the amplitude of the perturbation is computed. Assuming the standard models of the very early Universe, the initial condition for Eq. (71) is derived assuming a radiation dominated medium, with an equation-of-state $p = w\rho$ and $w = 1/3$. When a constant $w(t) = \bar{w}$ characterises the fluid dominating the energy budget of the Universe, we have $d\Phi(t)/dt = 0$ and one obtains

$$ \bar{\Phi} = \frac{3(1 + \bar{w})}{(5 + 3\bar{w})}, \qquad (72) $$

yielding $\bar{\Phi} = 2/3$ for a radiation fluid with $\bar{w} = 1/3$. Such a prefactor is standardly reported in Eq. (29). The behavior of $\Phi$ across the QCD phase transition, obtained by solving Eq. (71), differs from the average $\bar{\Phi}$, particularly in the region where $w$ and $c_s^2$ are quickly varying with respect to $M_H$ (left panel of Fig. 14). The resulting evolution of $\Phi$ as a function of horizon crossing mass is shown in the middle panel of Fig. 14.

Overall, the QCD phase transition introduces an additional degree of freedom into the problem, which is the characteristic scale of the horizon crossing of the cosmological perturbation. In Fig. 15, obtained from the computations done in Musco et al. (2024), we summarise the modified behaviour of the threshold $\delta_c$ for different values of the shape parameter $\alpha$ defined in (37). The different lines with a color varying between red, for smaller values of $\alpha$ and blue, for larger values, shows as the threshold is varying during the QCD epoch, as function of the cosmological horizon mass $M_H$. This is in agreement with the results obtained in Escrivà et al. (2023), apart from a few percent difference due to the different set up of the initial conditions.

Musco et al. (2024) have computed also the modified critical scaling behaviour defined in Eq. (39), which is modified as

$$ m_{\text{PBH}} = K(M_H)[\delta - \delta_c(M_H)]^{\gamma(M_H)} M_H, \qquad (73) $$

where $\delta_c(M_H)$, $\gamma(M_H)$ and $K(M_H)$ are fitted to numerical simulations performed in Musco et al. (2024) as a function of $M_H$, i.e. when the perturbation is crossing the cosmological horizon, defined as

$$ M_H \approx 17 M_{\odot} \left( \frac{g_*}{10.75} \right)^{-1/6} \left( \frac{k/\kappa}{10^6 \text{Mpc}^{-1}} \right)^{-2}, \qquad (74) $$

where $g_*$ is the number of degrees of freedom of relativistic particles and $\kappa \equiv r_m k$ relates the spectral wavenumber with the horizon crossing size of perturbations as derived in Sec. 3.3. The different lines shown in Fig. 16 with a color varying between red, for smaller values of $M_H$, and blue for larger values, shows how the scaling law is modified by the characteristic scale of the problem. Notice that an exact power-law critical behaviour is only obtained close enough to the density threshold
<figure><img src="image_16.png" /><figcaption>**Fig. 14** Left panel: the equation-of-state parameter $w = p/\rho$ (red) and squared speed of sound (blue) as functions of the cosmological horizon mass $M_H$. Right panel: Evolution of the equation-of-state dependent parameter $\Phi$, relating the density contrast to the curvature perturbation as functions of the cosmological horizon mass $M_H$. Image reproduced with permission from Franciolini et al. (2022e), copyright by APS Image reproduced with permission</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 15** Relative variation of the threshold compared to what is obtained assuming perfect radiation as a function of the horizon crossing time (parametrised here with $M_H$) induced by the QCD thermal effects. The color code indicates the different values of $\log_{10}(\alpha)$ as indicated by the bar on top of the frame. Figure adapted from Musco et al. (2024) Image reproduced with permission</figcaption></figure>

($\delta - \delta_c \lesssim 10^{-5}$), where the PBH masses are significantly smaller than the cosmological horizon mass, not able to affect significantly the collapse, while for larger values the equation-of-state during the QCD epoch induces further modifications. Franciolini et al. (2022e) fit the relation between the PBH and horizon mass using the power-law template (73) in the range of $\delta$ which most contributes to the abundance, i. e. ($\delta - \delta_c \in [10^{-5}, 2 \times 10^{-2}]$), and finds that deviations from the functional form used in Eq. (73) would only induce a small correction which we can neglect. The resulting values of $K(M_H)$ and $\gamma(M_H)$ used here are shown in the right plot of Fig. 16. For reference, we also show the critical collapse during the radiation dominated epoch of the early Universe, indicated with a black dashed line. A general
trend is observed: for $M_H \lesssim 3M_\odot$, there is a tendency to generate heavier PBHs, while the opposite is found when $M_H \gtrsim 3M_\odot$. This can be seen in the left panel of Fig. 16, where orange (light blue) lines fall above (below) the dashed black line indicating the result for a radiation perfect fluid.

In order to compute the full mass distribution, including both the thermal effects as well as non-linearities (56), it is convenient to invert the relation between horizon and PBH mass through the critical collapse relation (73), focusing only on the so called type-I branch for which $\delta_l < \Phi$, see Musco et al. (2021) for more details, as

$$ \delta_l = 2\Phi(1 - \sqrt{\Lambda}), \quad \text{where} \quad \Lambda = 1 - \frac{\delta_c}{\Phi} - \frac{1}{\Phi} \left( \frac{m_{\text{PBH}}}{\kappa M_H} \right)^{1/7}. \qquad (75) $$

Finally, one can arrive at a compact expression for the computation of the mass distribution which is

$$ \psi(m_{\text{PBH}}) \equiv \frac{1}{\Omega_{\text{PBH}}} \frac{d\Omega_{\text{PBH}}}{dm_{\text{PBH}}} = \frac{8}{3\pi \Omega_{\text{PBH}} m_{\text{PBH}}} \int \frac{dM_H}{M_H} \left(\frac{M_{\text{eq}}}{M_H}\right)^{1/2} \left(\frac{\sigma_1}{aH\sigma_0}\right)^3 \\ \times \frac{\Phi^3 \kappa}{\gamma \sigma_0^4} \left(\frac{m_{\text{PBH}}}{\kappa M_H}\right)^{\frac{1+\gamma}{\gamma}} \frac{(1 - \sqrt{\Lambda})^3}{\Lambda^{1/2}} \exp\left[-\frac{2\Phi^2}{\sigma_0^2}(1 - \sqrt{\Lambda})^2\right], \qquad (76) $$

and the integration range of $M_H$ is subject to the condition $\Lambda > 0$ (because we require $\delta > \delta_c$). An extension of this formula that includes both thermal effects and non-Gaussianities of $\zeta$, as discussed in Sect. 3.7, is presented in Ferrante et al. (2023b). The quantities $\kappa(M_H)$, $\gamma(M_H)$, $\Phi(M_H)$, $\delta_c(M_H)$, and $\sigma_l(M_H)$ are left within the integration over the horizon mass scale, as they all explicitly depend on $M_H$ when thermal effects are included. The variances can be computed integrating the curvature power spectrum as

<figure><img src="image_18.png" /><figcaption>Fig. 16 Left panel: PBH mass $m_{\text{PBH}}$ plotted as a function of $\delta - \delta_c$ computed at the cosmological horizon crossing (see Musco et al. 2024 for more details). The behavior for a radiation dominated medium is plotted with a black dashed line. Right panel: The values of the power law coefficients in Eq. (73) found by fitting the results of numerical simulations shown in the left panel. Image adapted from Franciolini et al. (2022e) Image reproduced with permission</figcaption></figure>
$$
\sigma_i^2(r_m) = \frac{4}{9} \Phi^2 \int_0^\infty \frac{dk}{k} (kr_m)^4 \tilde{W}^2(k, r_m) T^2(k, r_m) k^{2i} P_\zeta(k), \quad (77)
$$

where $i = 0, 1$ and the window and transfer functions are defined in Eqs. (54) and (55), respectively. If a nearly scale invariant spectrum of curvature perturbations is assumed, the smoothing scale $r_m$ is related with power spectral modes as $r_m k \equiv \kappa = 4.49$, as discussed in Sect. 3.3. This relation is strictly valid for a shape parameter $\alpha = 3$ (Musco et al. 2021). Notice that in the low mass limit, that is when we consider masses much below the smallest $M_H$, one can find that the mass distribution (76) scales as

$$
\psi(m_{\text{PBH}}) \propto (m_{\text{PBH}})^{1/\gamma}, \quad (78)
$$

where $1/\gamma \approx 2.8$ if one assumes the energy density of the Universe behaving as a relativistic fluid with $w = 1/3$, which gives $\gamma \approx 0.36$ (Niemeyer and Jedamzik 1998).

We conclude this section with a note of caution. While it is true that the QCD era
produces a natural boost of formation around the solar mass, the physics of collapse
alone is not predictive. This is because all the model dependence is retained in the
spectrum and in the statistics of curvature perturbations. Without strong assumptions
on the tilt and running of the spectrum, for instance, no prediction on the height of
the QCD peak at the solar mass compared to other scales can be obtained. This is
reflected in Fig. 17, where different mass distributions are derived, assuming a
power-law curvature power spectrum with varying tilt $n_s$ (not necessarily related to
the one at large CMB scales). The different ways one can introduce the curvature or
density profiles can explain some of the differences observed between Musco et al.
(2024) and Escrivà et al. (2023), as well as the still low degree of accuracy to which
we know the evolution of the equation-of-state from QCD lattice simulations.
Overall, the PBH formation remains an intrinsically non-linear and fully relativistic
process, which makes the physics quite complex, as the recent developments have
shown. The improvement of our understanding of PBH formation has been very
important and will probably be pushed further ahead in the next years, with
significant implications for the PBH mass distribution expected in the various
proposed models.

3.9 Evolution of the mass function through accretion

The primordial PBH mass function would be preserved up to the present epoch unless phenomena affecting the PBH masses modify it during the cosmological history. In particular, it was shown that PBH mass accretion may be relevant for the PBH evolution (De Luca et al. 2020e). In this Subsection and in the following, we will describe the impact on the mass and spin of PBHs by a phase of accretion (see De Luca and Bellomo 2023 for a recent review).
<figure><img src="image_19.png" /><figcaption>Fig. 17 Mass function obtained with a few choices of the curvature power spectrum. This plot assumes $f_{\text{PBH}} = 10^{-3}$, the minimum horizon mass to be $\lesssim 10^{-2.5} M_{\odot}$, the largest mass $M_{\text{H}}^{\text{max}} = 10^{2.8} M_{\odot}$ and a variable tilt $n_s$. The black dashed line reports the lognormal mass distribution found as the best fit in the analysis of Franciolini et al. (2022a). Image reproduced with permission from Franciolini et al. (2022e), copyright by APS Image reproduced with permission</figcaption></figure>

### 3.9.1 Accretion onto isolated PBHs

PBHs can accrete mass if particles from the surrounding environment fall into the gravitational potential well generated by the compact object and eventually get eaten by the PBH. It can be shown that the strongest impact on the PBH mass is coming from baryonic accretion. In case PBHs are not accounting for the entirety of the DM, there exists an additional DM fluid in the Universe. Such a secondary DM component however would only act as a catalyst of baryonic accretion while it marginally contributes to the PBH mass growth (Ricotti et al. 2008).

During the cosmological history, one can model the mass variation using the Bondi-Hoyle accretion rate (Bondi 1952; Hoyle and Lyttleton 1939, 1940a, b; Ricotti 2007; Ricotti et al. 2008)

$$ \dot{m}_{\text{PBH}} = 4\pi\lambda m_{\text{p}} n_{\text{gas}} v_{\text{eff}}^{-3} m_{\text{PBH}}^2, \qquad (79) $$

where the PBH effective velocity $v_{\text{eff}} = \sqrt{v_{\text{rel}}^2 + c_s^2}$ is given in terms of the relative velocity $v_{\text{rel}}$ with respect to the surrounding gas of hydrogen with mass $m_{\text{p}}$, sound speed $c_s$ and number density $n_{\text{gas}}$. We denote with the Bondi radius $r_b$ the region of space under the influence of the PBH potential well which is also determining the characteristic scale of the gas cloud accreting onto the PBHs. The parameter $\lambda$, whose explicit expression is summarized for example in Appendix B of De Luca et al. (2020e), takes into account the effects of gas viscosity, Hubble expansion and the coupling of the CMB radiation to the gas through Compton scattering (Ricotti 2007). Notice that the aforementioned treatment is based on newtonian mechanics, and general-relativistic effects may lead to a significant increase in the mass accretion rate, depending on the specific environmental settings, see for example Cruz-Osorio and Rezzolla (2020), Cruz-Osorio et al. (2021).
The presence of observational constraints for PBHs with masses larger than $O(M_\odot)$ suggest that they can comprise only a fraction of the DM in the Universe, and therefore an additional DM halo has to be considered when modelling gas accretion (Mack et al. 2007; Adamek et al. 2019). The presence of a DM halo catalyses the accretion of gas onto the PBHs, with a halo characteristic size and mass which grows with time as (Ricotti et al. 2008)

$$r_{\text{halo}} = 0.019 \text{ pc} \left(\frac{M_{\text{halo}}}{M_{\odot}}\right)^{1/3} \left(\frac{1+z}{1000}\right)^{-1}, \quad M_{\text{halo}} = 3m_{\text{PBH}} \left(\frac{1+z}{1000}\right)^{-1} \quad (80)$$

and enhances the accretion rate (Ricotti 2007). In the limit in which the characteristic halo radius is smaller than the PBH Bondi radius, accretion occurs onto a PBH with effective point mass $M_{\text{halo}}$, while for comparable sizes or in the opposite limit, the proper contribution is captured by corrections to the accretion parameter $\lambda$, see Appendix B of De Luca et al. (2020e) and references therein for additional details.

For PBHs with masses larger than $O(10)M_\odot$, accretion can reach super-Eddington values before the reionization epoch and play an important role in both the mass and spin evolution of PBHs (see Ricotti 2007; De Luca et al. 2020e). The increase of the PBH masses due to accretion also affects their mass distribution function in a non-linear fashion (García-Bellido 2018; De Luca et al. 2019a). Indeed, for an initial mass function $\psi(m_{\text{PBH}}^1, z_1)$ at formation redshift $z_1$, its evolution is governed by the identity (De Luca et al. 2020e, c)

$$\psi(m_{\text{PBH}}(m_{\text{PBH}}^1, z), z) \, dm_{\text{PBH}} = \psi(m_{\text{PBH}}^1, z_1) \, dm_{\text{PBH}}^1, \quad (81)$$

where $m_{\text{PBH}}(m_{\text{PBH}}^1, z)$ is the final mass at redshift $z$ for a PBH with mass $m_{\text{PBH}}^1$ at redshift $z_1$. The main effect of accretion on the mass distribution is to make the latter broader at high masses, producing a high-mass tail that can be orders of magnitude above its corresponding value at formation (De Luca et al. 2020c).

Accretion may be also particularly relevant for the formation of SMBHs with masses $m_{\text{BH}} \gtrsim 10^9 M_\odot$, from PBHs lighter than $10^4 M_\odot$ formed at high redshift (Serpico et al. 2020). In particular, assuming Gaussian initial conditions and PBH abundance $f_{\text{PBH}} \lesssim 10^{-9}$ in order to be compatible with the strong observational constraints by CMB spectral distortions (Kohri et al. 2014; Carr and Silk 2018; Nakama et al. 2018; Ünal et al. 2021), epochs of Eddington accretion onto PBHs with masses $\sim 200 M_\odot$ may be responsible for the formation of SMBHs at redshift $z \gtrsim 6$ (Serpico et al. 2020), whose merging will possibly be detected by LISA.

At the onset of structure formation and reionization epoch, the increase of the PBH characteristic velocities and the gas speed of sound reduce the accretion efficiency (Ali-Haïmoud et al. 2017; Hasinger 2020; Hütsi et al. 2019). One can adopt an agnostic view and consider a cut-off redshift $z_{\text{cut-off}}$ below which accretion is negligible (De Luca et al. 2020e, d). Its value is relatively unconstrained due to the large uncertainties in modelling accretion at relatively small redshift, such as X-ray pre-heating (Oh and Haiman 2003), details of the structure formation and feedback effects (Park and Ricotti 2011, 2012, 2013). In particular, the emission of X-rays from accreting PBHs may heat the plasma locally (local type) or produce Strömgren
spheres of ionised H-II which may impact the gas sound speed and thus the PBH accretion rate (global type) (Ricotti et al. 2008; Ali-Haïmoud and Kamionkowski 2017), while the generation of outflows may reduce the accretion efficiency depending on the outflow’s velocity and direction (mechanical type) (Bosch-Ramon and Bellomo 2020). Similarly, the nonspherical nature of the accretion flow could impact on the efficiency of the process. These effects make the full accretion process more complicated, and dedicated works beyond the Bondi-Hoyle formalism are needed to draw firmer conclusions on its final impact onto PBHs.

### 3.9.2 Accretion onto PBH binaries

When accretion occurs onto binary PBHs, one has to take into account both global accretion processes (i.e., of the binary as a whole) and local accretion processes (i.e., onto the individual components of the binary). In particular, as the typical binary separation is smaller than the binary Bondi radius, the infall of gas is driven by the total mass $M_{tot}$ of the binary system. In this configuration, both PBHs experience accretion from the gas with an enhanced density, and their individual accretion rates are modulated by their masses and orbital velocities (De Luca et al. 2020d)

$$ \dot{m}_1 = \dot{M}_{tot} \frac{1}{\sqrt{2(1+q)}}, \qquad \dot{m}_2 = \dot{M}_{tot} \frac{\sqrt{q}}{\sqrt{2(1+q)}} \qquad (82) $$

where the accretion rate $\dot{M}_{tot}$ is computed using Eq. (79) adopting the binary system parameters. For binary systems, one always finds that the relative accretion rate for the secondary component is larger than the one for the primary component, resulting in a growth of the mass ratio of the binary towards unity according to the equation (De Luca et al. 2020d)

$$ \dot{q} = q \left( \frac{\dot{m}_2}{m_2} - \frac{\dot{m}_1}{m_1} \right) > 0. \qquad (83) $$

## 3.10 Spin distribution

In this subsection, we discuss the theoretical expectation for the PBH initial spin and how it evolves when accretion is effective.

The collapse of density perturbations generating a PBH is expected to be a rare event due to the fact that their cosmological abundance cannot exceed that of the DM, and applying the peak theory formalism one finds that the high (and rare) peaks in the density contrast giving rise to PBHs are primarily spherical (Bardeen et al. 1986). At first order in perturbation theory, however, small asymmetries introduce torques induced by the surrounding matter perturbations, which are ultimately responsible for the generation of a small angular momentum before collapse. The small time of collapse makes however the action of these torques limited in time. At formation, the mean PBH dimensionless Kerr parameter $\chi \equiv J/m_{PBH}^2$ is estimated to be (De Luca et al. 2019a)
$$
\chi_i \sim 10^{-2} \sqrt{1 - \gamma_{\text{PS}}^2}, \qquad (84)
$$

in terms of the power spectrum shape parameter γ<sub>PS</sub>, defined as the ratio between the cross and auto-correlations of the density contrast and its spatial derivatives (close to unity for very narrow power spectra). The initial PBH spin is therefore expected to be below the percent level (see also Mirbabayi et al. 2020; Harada et al. 2021).<sup>5</sup> The distribution of the dimensionless Kerr parameter at formation takes the form as found in (De Luca et al. 2019a)

$$
P(\chi_i) d\chi_i = u \exp[-2.37 - 4.12 \ln(u\chi_i) - 1.53 \ln^2(u\chi_i) - 0.13 \ln^3(u\chi_i)] d\chi_i, \quad (85)
$$

where for simplicity we defined

$$
u \equiv \left( \frac{5}{2^{7/2}} \frac{\gamma_{\text{PS}}^6 v}{\sigma_{\delta} \sqrt{1 - \gamma_{\text{PS}}^2}} \right), \qquad (86)
$$

in terms of the rescaled critical peak amplitude v ≡ δ_c / σ_δ and the density contrast variance σ_δ.

The initial PBH spin will be retained until the time of merger unless an efficient
phase of accretion underwent during the PBH cosmological evolution. Indeed, since
the accretion rate and the geometry of the accretion flow are intertwined, for
accretion rates slightly above the super-Eddington limit, the angular momentum
carried by the infalling accreting gas may lead to the formation of a geometrically
thin accretion disk along the equatorial plane of the PBH (Shakura and Sunyaev
1973; Novikov and Thorne 1973; Ricotti et al. 2008).<fnref n="6" /> Apart from efficient accretion,
the formation of a disk happens if the typical gas velocity is larger than the Keplerian
velocity near the PBH, which translates into a condition on the minimum PBH
mass (De Luca et al. 2020e)

$$
m_{\text{PBH}} \gtrsim 6 \cdot 10^2 M_{\odot} D^{1.17} \xi^{4.33}(z) \frac{(1+z/1000)^{3.35}}{\left[1 + 0.031(1+z/1000)^{-1.72}\right]^{0.68}}, \quad (87)
$$

in terms of the constant factor $D \sim \mathcal{O}(1-10)$, which takes into account relativistic effects, and the parameter $\xi(z) = \max[1, \langle v_{\text{eff}} \rangle/c_s]$, which describes the effect of the PBH motion with respect to the gas.

In such a configuration, mass accretion is accompanied by an increase of the PBH
spin perpendicularly to the disk plane, whose growth rate can be described following
a geodesic model for circular disk motion as (Bardeen et al. 1972; De Luca et al.

<footnote n="5">⁵ The formation of PBHs in non-standard scenarios, like during an early matter-dominated epoch following inflation (Harada et al. 2017) or from the collapse of Q-balls (Cotner and Kusenko 2017), may lead to higher values of the initial spin.</footnote>

<footnote n="6">⁶ For accretion rates slightly smaller than the Eddington limit, accretion is non-spherical and an advection-dominated accretion flow (ADAF) may form (Narayan and Yi 1995). In the opposite case, the accretion luminosity might be so strong that the disk “puffs up” and becomes thicker, resulting in less efficient angular momentum transfer.</footnote>
2020e) (see also Bardeen et al. 1972; Thorne 1974; Brito et al. 2015; Volonteri et al. 2005)

$$ \dot{\chi} = (\mathcal{F}(\chi) - 2\chi) \frac{\dot{m}_{\text{PBH}}}{m_{\text{PBH}}}, \qquad (88) $$

in terms of the combination $\mathcal{F}(\chi) \equiv L(m_{\text{PBH}}, J)/m_{\text{PBH}}E(m_{\text{PBH}}, J)$, which is a function of the energy and angular momentum per unit mass (Bardeen et al. 1972)

$$ E(m_{\text{PBH}}, J) = \sqrt{1 - 2 \frac{m_{\text{PBH}}}{3 r_{\text{ISCO}}} \quad \text{and} \quad L(m_{\text{PBH}}, J) = \frac{2 m_{\text{PBH}}}{3 \sqrt{3}} \left( 1 + 2 \sqrt{3 \frac{r_{\text{ISCO}}}{m_{\text{PBH}}} - 2} \right)}, \quad (89) $$

depending on the ISCO radius

$$ r_{\text{ISCO}}(m_{\text{PBH}}, J) = m_{\text{PBH}} \left[ 3 + Z_2 - \sqrt{(3 - Z_1)(3 + Z_1 + 2Z_2)} \right], \quad (90) $$

with $Z_1 = 1 + (1 - \chi^2)^{1/3} [(1 + \chi)^{1/3} + (1 - \chi)^{1/3}]$ and $Z_2 = \sqrt{3\chi^2 + Z_1^2}$. The spin evolution continues within the relevant accretion timescale until it reaches the extremal limit $\chi_{\text{max}} = 0.998$ dictated by radiation effects (Thorne 1974). An analytical fit of the spin as a function of mass induced by accretion onto PBH binaries can be found in Franciolini et al. (2022b). This fit is based on the accretion model presented in De Luca et al. (2020e, 2020d) which is valid in the mass range currently observed by LVKC experiment. Detailed modelling of accretion on PBH with masses above around $\mathcal{O}(10^2) M_{\odot}$, particularly relevant for PBH mergers observable by LISA, is still lacking in the literature, due to the difficulty in modelling accretion with the inclusion of feedback effects at such large masses.

If instead accretion occurs onto a PBH binary, the fact that total mass of the binary drives accretion leads to a much more efficient angular momentum transfer on each PBH (De Luca et al. 2020d, e). In this case, the accretion onto the binary components is never spherical and disks can form easier around the PBHs. Moreover, since the lighter component of the binary experiences a stronger accretion rate, this results into a higher spin with respect to the one of the primary and heavier component (De Luca et al. 2020d).

An important parameter measurable through GW observations is the effective spin

$$ \chi_{\text{eff}} \equiv \frac{\chi_1 \cos \alpha_1 + q \chi_2 \cos \alpha_2}{1 + q}, \qquad (91) $$

which is a function of the mass ratio $q$, of both BH spin magnitudes $\chi_j$ ($j = 1, 2$), and of their orientation with respect to the orbital angular momentum, parametrized by the tilt angles $\alpha_j$. In Fig. 18, we show the expected distribution of $\chi_{\text{eff}}$ by averaging over the spin angles, as a function of PBH masses in binaries for various choices of $Z_{\text{cut-off}}$.

An alternative scenario for inducing spin on PBH occurs in the context of dense PBH clusters, where close hyperbolic encounters (García-Bellido and Nesseris 2017, 2018) may spin up both black holes (Nelson et al. 2019; Jaraba and García-Bellido 2021) due to the fundamental frame-dragging effect of general relativity, analogous to magnetic induction in electromagnetism, in this case due to a current of
matter as one black hole scatters off another. In the case of highly asymmetric binaries with large mass ratios, the induced spin can reach up to $\chi \simeq 0.8$ for the most massive one (Jaraba and García-Bellido 2021). Nevertheless, the distribution of spins of PBH in dense clusters is still characterized by a Gaussian around zero-spin with a width of order $\sigma \simeq 0.2$ (García-Bellido et al. 2021c).

## 3.11 Comments on the clustering from Poisson initial conditions

In the absence of non-Gaussian initial conditions, PBHs are Poisson distributed, as shown in Desjacques and Riotto (2018). In such a case, the PBH clustering evolution has been studied with a cosmological N-body simulations in Inman and Ali-Haïmoud (2019) up to redshift of $\mathcal{O}(10^2)$ and confirmed analytically in De Luca et al. (2020a) assuming a monochromatic PBH mass distribution. In this Subsection we briefly summarise the main conclusions.

To characterize the PBH two-point correlation function as a function of the comoving separation $x = |\vec{x}|$, we introduce the overdensity of discrete PBH centers at position $\vec{x}_i$ with respect to the total background DM energy density,

$$ \frac{\delta \rho_{\text{PBH}}(\vec{x}, z)}{f_{\text{PBH}} \bar{\rho}_{\text{DM}}} = \frac{1}{\bar{n}_{\text{PBH}}} \sum_{i} \delta_{\text{d}}(\vec{x} - \vec{x}_{i}(z)) - 1, \quad (92) $$

where $\delta_d(\vec{x})$ is the three-dimensional Dirac distribution, and

$$ \bar{n}_{\text{PBH}} \approx 3.2 f_{\text{PBH}} \left( \frac{20 M_{\odot}/h}{m_{\text{PBH}}} \right) (h/\text{kpc})^3 \quad (93) $$

is the average number density of PBHs per comoving volume. Here, $i$ runs over the positions of PBHs. The two-point correlation function of this discrete point process takes the general form

<figure><img src="image_20.png" /><figcaption>Fig. 18 Left: Predicted primary ($\chi_1$) and secondary ($\chi_2$) spins as a function of primary mass $m_1$ and mass ratio $q$ for various values of $z_{\text{cut-off}}$ (indicated by colors specified in the right panel). Right: Predicted distribution of $\chi_{\text{eff}}$ as a function of PBH mass $m_1$ (assuming equal mass binaries) for three choices of $z_{\text{cut-off}}$. Image reproduced with permission from Franciolini et al. (2022b), copyright by APS Image reproduced with permission</figcaption></figure>
$$ < \frac{\delta \rho_{\text{PBH}}(\vec{x}, z) \delta \rho_{\text{PBH}}(0, z)}{\bar{\rho}_{\text{DM}}} > = \frac{f_{\text{PBH}}^2}{\bar{h}_{\text{PBH}}} \delta_d(\vec{x}) + \xi(x, z). \quad (94) $$

This expression emphasizes that $\xi(x, z)$ is the so-called reduced PBH correlation function and, thus, is distinct from the additive Poisson noise proportional to the Dirac delta. As mentioned before, initially the reduced correlation function is negligible. The corresponding PBH power spectrum

$$ \Delta^2(k, z) = \frac{k^3}{2\pi^2} \int d^3x e^{i\vec{k}\cdot\vec{x}} < \frac{\delta\rho_{\text{PBH}}(\vec{x}, z) \delta\rho_{\text{PBH}}(0, z)}{\bar{\rho}_{\text{DM}} \bar{\rho}_{\text{DM}}} >, \quad (95) $$

is conveniently defined relative to the total cold dark matter average density.

After a linear growth, the PBH power spectrum enters in the so called “Quasi-Linear” regime during which

$$ \Delta_{\text{QL}}^2(k) \approx 0.04 f_{\text{PBH}}^{3/4} \left( \frac{20 M_{\odot}/h}{m_{\text{PBH}}} \right)^{-3/4} \left[ 1 + 26 f_{\text{PBH}} \left( \frac{100}{1+z} \right) \right]^{3/2} \left( \frac{k}{h/\text{kpc}} \right)^{9/4} . \quad (96) $$

When the power spectrum becomes larger than around 200, it enters in the “Non-Linear” regime during which

$$ \Delta_{\text{NL}}^2(k) \approx 0.2 f_{\text{PBH}}^{3/5} \left( \frac{m_{\text{PBH}}}{20 M_{\odot}/h} \right)^{3/5} \left[ 1 + 26 f_{\text{PBH}} \left( \frac{100}{1+z} \right) \right]^{6/5} \left( \frac{k}{h/\text{kpc}} \right)^{9/5} . \quad (97) $$

PBHs form halos whose mass distribution may be described by the Press-Schechter theory (Press and Schechter 1974) with an initial Poisson power spectrum. The resulting number density of PBH halos with mass between M and (M + dM) reads

$$ \frac{dn(M, z)}{dM} = \frac{\bar{\rho}_{\text{PBH}}}{\sqrt{\pi}} \left( \frac{M}{M_*(z)} \right)^{1/2} \frac{e^{-M/M_*(z)}}{M^2}, \quad (98) $$

where $\bar{\rho}_{\text{PBH}}$ is the average PBH energy density and (Hütsi et al. 2019)

$$ M_*(z) = N_*(z) \cdot m_{\text{PBH}} \approx f_{\text{PBH}}^2 \left( \frac{2600}{1+z} \right)^2 m_{\text{PBH}} \quad (99) $$

is the typical mass of halos collapsing at redshift z. This prediction agrees with the results of cosmological N-body simulations, as shown in Fig. 19.

If we restrict ourselves to small scales, both members of a PBH pair are almost certainly drawn from the same PBH halo. In this limit, if the PBH density profile is $\rho_{\text{PBH}}(x) \sim x^{-\epsilon_{\text{PBH}}}$, then the two-point correlation function must behaves like $\sim x^{-2\epsilon_{\text{PBH}}+3}$ (Peebles 1974; McClelland and Silk 1977) as it is proportional to the square density profile. Imposing $(-2\epsilon_{\text{PBH}} + 3) = -9/5$, we infer that the PBH density profile should satisfy
$$ \rho_{\text{PBH}}(x) \sim x^{-12/5}, \qquad (100) $$

which is confirmed by numerical simulations presented in Raidal et al. (2019). Notice that since clustering is hierarchical, each halo has a certain survival time and, therefore, a given probability to be absorbed by a bigger halo formed at a later redshift. As shown in De Luca et al. (2020a) the evaporation time (due to dynamical relaxation) of PBH halos is typically larger than their survival time, which implies that PBH halos are stable against evaporation. We will discuss later on the impact of clustering on the merger rate of PBH binaries.

## 3.12 Summary

In this section, we did not only review the standard formalism to compute the PBH mass function from the curvature or density power spectrum, but also the most recent developments on the computation of the critical overdensity threshold and of the relation between curvature and density fluctuations, which depend on the perturbation statistics (Gaussian, non-Gaussian), the primordial power spectrum and the evolution of the equation-of-state through different phases like the QCD transition. The impact of these effects is important on the PBH abundance and mass function. Even if it is possible to rescale the power spectrum amplitude in order to obtain the required PBH abundance when one only focuses on a single observable, this is not the case anymore when one tries to correlate multiple observations (e.g. PBH abundance and induced GWs, as discussed in the following sections). In particular, the effects of the QCD transition may have greatly imprinted on the mass function for stellar-mass PBHs with direct consequences for GW observations that will be discussed in the next Sections. We also went beyond the common assumption that the

<figure><img src="image_21.png" /><figcaption>Fig. 19 The abundance of halos $N_{\text{HL}}$ containing a given number of PBHs, $N_{\text{PBH/HL}}$, i.e. halo mass function. Solid lines report the results of the N-body simulations while dashed lines are indicates the theoretical prediction assuming Poisson statistics. Image reproduced with permission from Inman and Ali-Haïmoud (2019), copyright by APS Image reproduced with permission</figcaption></figure>
PBH mass and spin distributions do not change much with time, by considering effects of accretion on isolated PBHs and PBH binaries. Finally we emphasized the importance of the inevitable clustering of PBHs due to their inevitable Poisson fluctuations.

# 4 PBH merging and encounter rates

In this section, we calculate the expected PBH merging rates and their evolution with redshift for two different binary formation channels: before matter-radiation equality and due to tidal capture in clusters. For this purpose, we use a general mass distribution $f(m_{\text{PBH}})$ that must be specified by the underlying theoretical model and assumptions about PBH formation, as explained in the previous sections. For the two channels, we also discuss some effects linked to PBH clustering induced by the inevitable Poisson fluctuations in their spatial distribution. This can cause a merging rate suppression for early binaries due to their tidal disruption by early-forming clusters and a merger rate boost for late-binaries due to an enhanced clustering compared to common halo mass functions in $\Lambda$-CDM cosmology. We also review the expected rates of disrupted binaries. Finally, we study the rate of hyperbolic encounters in PBH clusters that can induce GW bursts, possibly observable in Earth-based and space-based detectors.

## 4.1 Early binaries

### 4.1.1 General rate formula

PBH binaries may have formed in the early Universe before matter-radiation equality. Because of their random spatial distribution at formation, two PBHs can form sufficiently close to each other for their dynamics to decouple from the Universe expansion. They form a binary instead of directly merging, because of the gravitational influence of one or several PBHs nearby. Possibly, the binary is sufficiently stable and it takes of the order of the age of the Universe for the two black holes to merge.

Under the hypothesis that there is no local non-Gaussianity cross-correlating large and small scales, PBHs do not form in clusters, but follow a Poissonian distribution (Desjacques and Riotto 2018; Ali-Haïmoud 2018; Ballesteros et al. 2018; Moradinezhad Dizgah et al. 2019; Inman and Ali-Haïmoud 2019). The merging rates (per unit logarithmic mass of the two binary components) of such early binaries (denoted EB) at a time $t$ were first estimated in (Ioka et al. 1998) and are given by (Sasaki et al. 2016; Ali-Haïmoud et al. 2017; Kocsis et al. 2018; Raidal et al. 2019; De Luca et al. 2020d; Hütsi et al. 2021c)
$$
\begin{equation}
\begin{aligned}
R_{\text{EB}} = & \frac{1.6 \times 10^6}{\text{Gpc}^3 \text{yr}} \times f_{\text{sup}}(m_1, m_2, f_{\text{PBH}}) f_{\text{PBH}}^{53/37} f(m_1) f(m_2) \\
& \times \left(\frac{t}{t_0}\right)^{-34/37} \left(\frac{m_1 + m_2}{M_{\odot}}\right)^{-32/37} \left[\frac{m_1 m_2}{(m_1 + m_2)^2}\right]^{-34/37},
\end{aligned}
\tag{101}
\end{equation}
$$

where $m_1$ and $m_2$ are the two black hole masses and $t_0$ is the age of the Universe. Compared to the predictions of (Sasaki et al. 2016), one has to include a suppression factor $f_{\text{sup}}$ to take into account several physical effects that have been initially neglected but were revealed by N-body simulations of (Raidal et al. 2019).

Some analytical prescriptions have been proposed, e.g. in (Raidal et al. 2019; Hütsi et al. 2021c), to calculate the suppression factor $f_{\text{sup}}(f_{\text{PBH}}, m_1, m_2)$. It can be written as the product of $S_1(f_{\text{PBH}}, m_1, m_2)$ and $S_2(f_{\text{PBH}})$ that take into account the rate suppression due to nearby PBHs or matter fluctuations and PBH clusters seeded by Poisson fluctuations, respectively. These analytical prescriptions have been compared with N-body simulations, but only in the case of monochromatic and log-normal PBH mass distributions. Therefore, one must be cautious when applying them to broad mass functions, even if they exhibit a peak, like the one from the QCD transition.

4.1.2 Suppression from close PBHs and matter inhomogeneities

The first suppression factor we consider ($S_1$) takes into account the binary disruption by either matter fluctuations, with a (rescaled) variance $\sigma_m^2$, usually around 0.005, or by the number $\bar{N}$ of PBHs within a sphere centered on the binary of radius corresponding to the maximum comoving distance at which a nearby PBH would fall onto the binary to disrupt it, before matter-radiation equality. An analytical prescription for $S_1$ has been proposed in (Raidal et al. 2019; Hütsi et al. 2021c), and is given by

$$
S_1 \approx 1.42 \left[ \frac{\left\langle \left( \frac{m_{\mathrm{PBH}}^2}{\langle m_{\mathrm{PBH}} \rangle^2} \right)^2 \right\rangle}{\bar{N} + C} + \frac{\sigma_m^2}{f_{\mathrm{PBH}}^2} \right]^{-21/74} E^{-\bar{N}}, \quad (102)
$$

with

$$
\bar{N} = \frac{m_1 + m_2}{\langle m_{\text{PBH}} \rangle} \frac{f_{\text{PBH}}}{f_{\text{PBH}} + \sigma_M}. \quad (103)
$$

Note that in Eqs. (102) and (103), the mean PBH mass and their variance are
calculated from the mass function through

$$
\langle m_{\text{PBH}} \rangle = \frac{\int m_{\text{PBH}} \, \mathrm{d}n_{\text{PBH}}}{n_{\text{PBH}}} = \left[ \int \frac{f(m_{\text{PBH}})}{m_{\text{PBH}}} \, \mathrm{d} \ln m_{\text{PBH}} \right]^{-1} \quad (104)
$$