# Curbing PBHs with PTAs

A.J. Iovino  0000-0002-8531-5962<EMAIL_ADDRESS>Dipartimento
di Fisica, “Sapienza” Università di Roma, Piazzale Aldo Moro 5, 00185, Roma,
Italy  Istituto Nazionale di Fisica Nucleare, sezione di Roma, Piazzale Aldo
Moro 5, 00185, Roma, Italy  Department of Theoretical Physics and
Gravitational Wave Science Center,
24 quai E. Ansermet, CH-1211 Geneva 4, Switzerland  G. Perna
0000-0002-7364-1904<EMAIL_ADDRESS>Dipartimento di Fisica e
Astronomia “Galileo Galilei”, Università degli Studi di Padova, Via Marzolo 8,
I-35131, Padova, Italy  INFN, Sezione di Padova, Via Marzolo 8, I-35131,
Padova, Italy  Department of Theoretical Physics and Gravitational Wave
Science Center,
24 quai E. Ansermet, CH-1211 Geneva 4, Switzerland  A. Riotto
0000-0001-6948-0856<EMAIL_ADDRESS>Department of Theoretical Physics
and Gravitational Wave Science Center,
24 quai E. Ansermet, CH-1211 Geneva 4, Switzerland  H. Veermäe
0000-0003-1845-1355<EMAIL_ADDRESS>Keemilise ja Bioloogilise Füüsika
Instituut, Rävala pst. 10, 10143 Tallinn, Estonia

###### Abstract

Sizeable primordial curvature perturbations needed to seed a population of
primordial black holes (PBHs) will be accompanied by a scalar-induced
gravitational wave signal that can be detectable by pulsar timing arrays
(PTA). We derive conservative bounds on the amplitude of the scalar power
spectrum at the PTA frequencies and estimate the implied constraints on the
PBH abundance. We show that only a small fraction of dark matter can consist
of stellar mass PBHs when the abundance is calculated using threshold
statistics. The strength and the shape of the constraint depend on the shape
of the power spectrum and the nature of the non-Gaussianities. We find that
constraints on the PBH abundance arise in the mass range
$0.1-10^{3}\,M_{\odot}$, with the sub-solar mass range being constrained only
for narrow curvature power spectra. These constraints are softened when
positive non-Gaussianity is introduced and can be eliminated when $f_{\rm
NL}\gtrsim 5$. On the other hand, if the PBH abundance is computed via the
theory of peaks, the PTA constraints on PBHs are significantly relaxed,
signalling once more the theoretical uncertainties in assessing the PBH
abundance. We further discuss how strong positive non-Gaussianites can allow
for heavy PBHs to potentially seed supermassive BHs.

##  I Introduction

Primordial black holes (PBHs) Zel’dovich and Novikov (1967); Hawking (1971);
Carr and Hawking (1974) are currently in the spotlight as they may solve
several open questions in astrophysics and cosmology. Different mechanisms and
models have been proposed as viable for PBH production and, depending on the
formation scenario considered and on the models for PBH formations, a variety
of mass functions have been predicted (see Ref. Bagui  _et al._ (2023) for a
recent review). Throughout this work, we will assume the standard formation
scenario which suggests that PBHs form out of the gravitational collapse after
the horizon re-entry of large over-densities in the primordial density
contrast field Ivanov  _et al._ (1994); Ivanov (1998); Blinnikov  _et al._
(2016). In such a scenario, a non-negligible PBH abundance requires an
enhanced primordial curvature power spectrum at scales smaller than the scale
of the Cosmic Microwave Background (CMB), at which the spectral amplitude is
around $10^{-9}$ Akrami  _et al._ (2020a).

In addition to producing PBHs, the same primordial curvature perturbations
generate tensor perturbations, also known as scalar-induced Gravitational
Waves (SIGW). As first shown by Refs. Matarrese  _et al._ (1994); Acquaviva
_et al._ (2003); Mollerach  _et al._ (2004); Ananda  _et al._ (2007); Baumann
_et al._ (2007), when going to second order, scalar and tensor perturbations
are coupled and the former provide a source for Gravitational Waves (GW) (see
also Refs. Espinosa  _et al._ (2018); Domènech (2021)). The statistical
properties of the curvature fluctuations can leave an imprint on this signal.
Indeed, the presence of primordial non-Gaussianity (NG) Gangui  _et al._
(1994); Bartolo  _et al._ (2002); Maldacena (2003); Bartolo  _et al._ (2004a,
b); Celoria and Matarrese (2020) can modify the spectrum inducing specific
features depending on the amplitude of the sourcing scalar power spectrum
Garcia-Saenz  _et al._ (2023); Unal (2019); Adshead  _et al._ (2021); Abe  _et
al._ (2023); Perna  _et al._ (2024). Moreover, the specific features of the
spectra and the NG modifications, as well as the amplitude of the induced SIGW
background, depend not only on the shape of the curvature power spectrum
considered, but also on the specific scales involved (see, e.g., Cai  _et al._
(2019a, b)). Hence detectable SIGWs offer a potential way to probe early
universe cosmology.

The recent Pulsar Timing Array (PTA) data release by the NANOGrav Agazie  _et
al._ (2023a, b), EPTA (in combination with InPTA) Antoniadis  _et al._ (2023a,
b, c), PPTA Reardon  _et al._ (2023a); Zic  _et al._ (2023); Reardon  _et al._
(2023b) and CPTA Xu  _et al._ (2023) collaborations shows evidence of a
Hellings-Downs pattern in the angular correlations which is characteristic of
GWs. The most stringent constraints and largest statistical evidence come from
the NANOGrav 15-year data (NANOGrav15). Even though it is currently not
possible to pinpoint the sources of this signal, this discovery is of extreme
importance since it opens a further window into the early universe and has an
impact on the GW search at the nHz scale (see, e.g., Ellis  _et al._ (2024);
Figueroa  _et al._ (2024)).

Since PTAs are sensitive to frequencies of the order of the nHz, they are
associated with the production of PBHs in the stellar mass range Chen  _et
al._ (2020); De Luca  _et al._ (2021); Vaskonen and Veermäe (2021); Kohri and
Terada (2021); Domènech and Pi (2022); Dandoy  _et al._ (2023); Franciolini
_et al._ (2023a). Crucially, requiring that the produced SIGW background does
not exceed the one registered by the NANOGrav15, strongly limits the maximum
amplitude of the curvature power spectrum, and thereby the abundance of PBHs
in the corresponding mass range. Additionally, if the power spectrum is
enhanced at scales larger than accessible by PTAs, PBHs can serve as seeds of
supermassive black holes (SMBHs) Duechting (2004); Bernal  _et al._ (2018);
Serpico  _et al._ (2020). Nevertheless, in this case, the scales related to
the formation process are strongly constrained by CMB spectral distortion
analysis of the FIRAS collaboration Fixsen  _et al._ (1996); Chluba  _et al._
(2012a, b); Chluba and Grin (2013); Bianchini and Fabbian (2022) and, without
NGs, the allowed PBH abundance is too small to provide a primordial origin for
SMBHs seeds.

In this work constraints on PBHs arising from NANOGrav15 are considered in
detail with an emphasis on local primordial NGs and assuming different shapes
for the primordial power spectrum. After having reviewed the origin of NGs in
the standard formation scenario for PBHs in Sec. II, in Sec. III we discuss
different formalisms for the computation of the PBH abundance and the role of
NGs. In Sec. IV we summarise the computation of the SIGW background. Then, in
Sec. V, we derive PTA constraints on the power spectrum and the abundance of
PBHs, and comment on how large positive NGs affect the generation of PBHs as
seeds for SMBHs. We conclude in Sec. VI.

##  II Non-Gaussian primordial curvature perturbations

Generally, different models for PBH production require an enhancement of the
curvature power spectrum that depends on the model considered. Typical shapes
widely used in literature are the _log-normal_ (LN)

$\mathcal{P}_{\zeta}(k)=\frac{{A}}{\sqrt{2\pi}\Delta}\exp\left(-\frac{1}{2\Delta^{2}}\ln^{2}\left(\frac{k}{k_{*}}\right)\right)\,,$
(1)

where $k_{*}$ denotes the peak frequency, or the _broken power-law_ (BPL)

$\mathcal{P}_{\zeta}(k)={A}\,\frac{\left(\alpha+\beta\right)^{\lambda}}{\left[\beta\left(k/k_{*}\right)^{-\alpha/\lambda}+\alpha\left(k/k_{*}\right)^{\beta/\lambda}\right]^{\lambda}},$
(2)

where the parameters $\alpha,\beta>0$ describe the growth and decay of the
spectrum around the peak and $\lambda$ the width of the peak. Typically,
$\alpha\approx 4$ Byrnes  _et al._ (2019); Karam  _et al._ (2023a) and
$\beta\lesssim 4$. The cases considered in this work are reported in Tab. 1.

Case |  PS shape |  $f_{\rm NL}$
---|---|---
LN |  LN $(\Delta=0.5)$ |  $[0,10]$
BPL1 |  BPL $(\alpha=4,\beta=3,\lambda=1)$ |  $[0,10]$
BPL2 |  BPL $(\alpha=4,\beta=0.5,\lambda=1)$ |  $[0,10]$
Table 1: List of benchmark cases considered in this work.

The standard formation scenario assumes that PBHs arise from the gravitational
collapse of significant over-densities in the primordial density contrast
field Ivanov  _et al._ (1994); Ivanov (1998). Given the random nature of this
field, the calculation of PBH abundance is inherently statistical and requires
precise knowledge of the probability density function (PDF) of the density
fluctuations.

Limiting the analysis to the assumption that the PDF of density fluctuations
$\delta$ follows the Gaussian statistics is theoretically flawed for two
reasons. First, density fluctuations in the primordial radiation field
originate from curvature perturbations $\zeta$, previously stretched on super-
horizon scales during inflation, after their horizon re-entry. In the long-
wavelength approximation, the equation which relates curvature perturbations
to density fluctuations is intrinsically non-linear Harada  _et al._ (2015)

$\displaystyle\delta(\vec{x},t)\\!=\\!-\frac{2}{3}\frac{\Phi}{(aH)^{2}}e^{-2\zeta(\vec{x})}\bigg{[}\nabla^{2}\zeta(\vec{x})\\!+\\!\frac{1}{2}\partial_{i}\zeta(\vec{x})\partial_{i}\zeta(\vec{x})\bigg{]},\\!\\!$
(3)

where $a$ denotes the scale factor, $H$ the Hubble rate, $\Phi$ is related to
the equation of state parameter ${w_{\rm u}}$ of the universe. For ${w_{\rm
u}}$ constant, $\Phi=3(1+{w_{\rm u}})/(5+3{w_{\rm u}})$ Polnarev and Musco
(2007). We dropped the explicit $\vec{x}$ and $t$ dependence for brevity. For
this reason, even under the assumption that curvature perturbations follow
exact Gaussian statistics, density fluctuations inherit an unavoidable amount
of NG from non-linear (NL) corrections De Luca  _et al._ (2019); Young  _et
al._ (2019); Germani and Sheth (2020). Second, there is no guarantee that
$\zeta$ is a Gaussian field – we refer to such cases as primordial NGs. Local
primordial NGs can be described by a relation

$\zeta(\vec{x})=F(\zeta_{\rm G}(\vec{x}))$ (4)

between $\zeta$ and its Gaussian counterpart $\zeta_{G}$, which depends on the
mechanism that generates the enhanced power spectrum at small scales (see the
related discussion in Ref. Ferrante  _et al._ (2023a)). Such NGs are thus
generically independent of large-scale NGs constrained by CMB data (e.g.
Akrami  _et al._ (2020b)). Copious PBH production via critical collapse of
large overdensities is achievable in a wide range of scenarios including
single-field inflation with an inflation point in the inflaton’s potential
(usually called USR models) Garcia-Bellido and Ruiz Morales (2017); Pi  _et
al._ (2018); Kannike  _et al._ (2017); Ballesteros  _et al._ (2020); Inomata
_et al._ (2017); Iacconi  _et al._ (2022); Kawai and Kim (2021); Bhaumik and
Jain (2020); Cheong  _et al._ (2021); Inomata  _et al._ (2018); Dalianis  _et
al._ (2019); Motohashi  _et al._ (2020); Hertzberg and Yamada (2018);
Ballesteros and Taoso (2018); Karam  _et al._ (2023a); Rasanen and Tomberg
(2019); Balaji  _et al._ (2022); Frolovsky and Ketov (2023); Dimopoulos
(2017); Germani and Prokopec (2017); Choudhury and Mazumdar (2014); Ragavendra
and Sriramkumar (2023); Cheng  _et al._ (2022); Franciolini  _et al._ (2023b);
Karam  _et al._ (2023b); Mishra  _et al._ (2023); Cole  _et al._ (2023); Karam
_et al._ (2023b); Frosina and Urbano (2023); Franciolini and Urbano (2022);
Choudhury  _et al._ (2024); Wang  _et al._ (2024); Stamou (2021, 2024);
Heydari and Karami (2022a, b, 2024) and models with spectator field, i.e., the
curvaton Enqvist and Sloth (2002); Lyth and Wands (2002); Sloth (2003); Lyth
_et al._ (2003); Dimopoulos  _et al._ (2003); Kohri  _et al._ (2013); Kawasaki
_et al._ (2013a, b); Carr  _et al._ (2017); Ando  _et al._ (2018a, b); Chen
and Cai (2019); Liu and Prokopec (2021); Pi and Sasaki (2023); Cai  _et al._
(2021); Liu (2023); Chen  _et al._ (2023); Torrado  _et al._ (2018); Chen  _et
al._ (2023); Wilkins and Cable (2024); Inomata  _et al._ (2024).

To be as model-independent as possible we follow the common approach and
describe $\zeta$ by the usual power series expansion

$\zeta=\zeta_{\rm G}+\frac{3}{5}f_{\rm NL}(\zeta_{\rm G}^{2}-\langle\zeta_{\rm
G}^{2}\rangle)+\frac{9}{25}g_{\rm NL}\zeta_{\rm G}^{3}+\dots\,,$ (5)

where $\zeta_{\rm G}$ obeys the Gaussian statistics while the coefficients
$f_{\rm NL}$, $g_{\rm NL}$ encode deviations from the Gaussian limit. In
general, the coefficients can depend on the scale of the perturbation. We omit
this possibility here.

##  III PBH formation

A PBH is formed when a sufficiently large primordial overdensity collapses
after horizon re-entry. The mass of the resulting PBH follows a critical
scaling law

$M_{\rm PBH}=\mathcal{K}M_{k}(\mathcal{C}-\mathcal{C}_{\rm th})^{\gamma}\,,$
(6)

where

$M_{k}\approx 17M_{\odot}\left(\frac{k}{10^{6}{\rm
Mpc}^{-1}}\right)^{-2}\,\left(\frac{g_{*}}{10.75}\right)^{-1/6}\,$ (7)

is the mass contained within a Hubble horizon corresponding to a comoving
scale $k$. The threshold $\mathcal{C}_{\rm th}$ is estimated following Ref.
Musco  _et al._ (2021) and depends on the shape of the power spectrum. We also
account for the time dependence of $\mathcal{C}_{\rm th}$, the critical
exponent $\gamma$ and the prefactor $\mathcal{K}$ during the QCD phase
transition Musco  _et al._ (2024).

The fraction of energy density $\beta_{k}(M_{\rm PBH}){\rm d}\ln M_{\rm PBH}$
collapsing into PBHs can be estimated as

$\\!\\!\\!\beta_{k}(M_{\rm PBH})\\!=\\!\int_{\mathcal{C}_{\rm
th}}^{\infty}\\!{\rm d}\mathcal{C}\frac{M_{\rm
PBH}(\mathcal{C})}{M_{k}}P_{k}(\mathcal{C})\delta\left[\ln\frac{M_{\rm
PBH}}{M_{\rm PBH}(\mathcal{C})}\right]\\!,$ (8)

where $P_{k}(\mathcal{C})$ denotes the probability a BH will form in the
Hubble patch. The PBH mass function can be obtained directly from the collapse
fraction:

$\displaystyle\\!\frac{{\rm d}f_{\rm PBH}}{{\rm d}\ln M_{\rm PBH}}$
$\displaystyle\\!=\\!\frac{1}{\Omega_{\rm DM}}\int\frac{{\rm
d}M_{k}}{M_{k}}\,\beta_{k}(M_{\rm PBH})\left(\frac{M_{\rm
eq}}{M_{k}}\right)^{1/2}\\!\\!,$ (9)

where $M_{\rm eq}\approx 2.8\times 10^{17}\,\,M_{\odot}$ is the horizon mass
at the time of matter-radiation equality and $\Omega_{\rm DM}=0.12h^{-2}$ is
the cold dark matter density Akrami  _et al._ (2020a). To characterise the PBH
population, we will consider the PBH abundance and the mean PBH mass,

$\displaystyle f_{\rm PBH}$ $\displaystyle=\int\frac{{\rm d}M_{\rm
PBH}}{M_{\rm PBH}}\frac{{\rm d}f_{\rm PBH}}{{\rm d}\ln M_{\rm PBH}}\,,$ (10)
$\displaystyle\langle M_{\rm PBH}\rangle$ $\displaystyle=f_{\rm
PBH}\left(\int\frac{{\rm d}M_{\rm PBH}}{M_{\rm PBH}^{2}}\frac{{\rm d}f_{\rm
PBH}}{{\rm d}\ln M_{\rm PBH}}\right)^{-1}\,,$

where the mean PBH mass is computed with respect to the PBH number density
(see e.g. Andrés-Carcasona  _et al._ (2024)). This offers a clear basis for
comparing constraints on PBHs with extended mass functions.

The above description can be applied to estimate PBH abundance using both the
threshold statistics and peaks theory. However, the shape and the origin of
$P_{k}(\mathcal{C})$ in these two approaches are different.

#### Threshold statistics

dictates that the PBH formation probability can be estimated from the
statistics of the compaction function $\mathcal{C}$ Ferrante  _et al._
(2023a); Gow  _et al._ (2023), generically defined as twice the local mass
excess over the areal radius111Using the threshold statistics approach, which
relies on average compaction profiles Ianniccari  _et al._ (2024), corrections
to the horizon crossing and from the non-linear radiation transfer function
Franciolini  _et al._ (2023c); De Luca  _et al._ (2023) can affect the PBH
abundance. Since the amount of these corrections is still not well understood
we leave their inclusions for a future work.. The latter can be estimated in
scenarios with non-Gaussian curvature perturbations $\zeta$ in case they can
be characterized via an auxiliary Gaussian field $\zeta_{G}$ as outlined in
Sec. II. The compaction function
$\mathcal{C}=\mathcal{C}_{1}-\mathcal{C}_{1}^{2}/(4\Phi)$ can then be
constructed from $\mathcal{C}_{1}=\mathcal{C}_{\rm G}\,{\rm d}F/{\rm
d}\zeta_{\rm G}$ which can be expressed in terms of the Gaussian field
$\mathcal{C}_{\rm G}=-2\Phi\,r\,\zeta_{\rm G}^{\prime}$.

The PBH formation probability in a Hubble patch is then

$P_{k}(\mathcal{C})=\int_{\mathcal{D}}\delta(\mathcal{C}-\mathcal{C}(\mathcal{C}_{\rm
G},\zeta_{\rm G}))P_{{\rm G},k}(\mathcal{C}_{\rm G},\zeta_{\rm G}){\rm
d}\mathcal{C}_{\rm G}{\rm d}\zeta_{\rm G}\,,$ (11)

where the domain of integration
$\mathcal{D}=\left\\{\mathcal{C}(\mathcal{C}_{\rm G},\zeta_{\rm
G})>\mathcal{C}_{\rm th}~{}\land~{}\mathcal{C}_{1}(\mathcal{C}_{\rm
G},\zeta_{\rm G})<2\Phi\right\\}$. and the Gaussian components are distributed
as

$P_{\mathrm{G}}\left(\mathcal{C}_{\mathrm{G}},\zeta_{\mathrm{G}}\right)=\frac{e^{\left[-\frac{1}{2\left(1-\gamma_{cr}^{2}\right)}\left(\frac{\mathcal{C}_{\mathrm{G}}}{\sigma_{c}}-\frac{\gamma_{cr}\zeta_{\mathrm{G}}}{\sigma_{r}}\right)^{2}\\!-\\!\frac{\zeta_{\mathrm{G}}^{2}}{2\sigma_{r}^{2}}\right]}}{2\pi\sigma_{c}\sigma_{\tau}\sqrt{1-\gamma_{cr}^{2}}}.$
(12)

The correlators are given by

$\displaystyle\sigma_{c}^{2}=\frac{4\Phi^{2}}{9}\int_{0}^{\infty}\frac{{\rm
d}k}{k}\left(kr_{m}\right)^{4}W^{2}\left(k,r_{m}\right)P^{T}_{\zeta}\,,$ (13a)
$\displaystyle\sigma_{cr}^{2}=\frac{2\Phi}{3}\\!\int_{0}^{\infty}\\!\\!\frac{{\rm
d}k}{k}\\!\left(kr_{m}\right)^{2}\\!W\\!\\!\left(k,r_{m}\right)\\!W_{s}\\!\left(k,r_{m}\right)\\!P^{T}_{\zeta}\\!\,,$
(13b) $\displaystyle\sigma_{r}^{2}=\int_{0}^{\infty}\frac{{\rm
d}k}{k}W_{s}^{2}\left(k,r_{m}\right)P^{T}_{\zeta}\,,$ (13c)

with $P^{T}_{\zeta}=T^{2}\left(k,r_{m}\right)P_{\zeta}(k)$, and
$\gamma_{cr}\equiv\sigma_{cr}^{2}/\sigma_{c}\sigma_{\tau}$. We have defined
$W\left(k,r_{m}\right),$ $W_{s}\left(k,r_{m}\right)$ and
$T\left(k,r_{m}\right)$ as the top-hat window function, the spherical-shell
window function, and the radiation transfer function Young (2022).

The PBH mass function (9) can be expressed as

$\displaystyle\frac{{\rm d}f_{\rm PBH}}{{\rm d}\ln M_{\rm
PBH}}\\!=\\!\frac{1}{\Omega_{\rm DM}}\int_{M_{\rm H}^{\rm min}}\frac{dM_{\rm
H}}{M_{\rm H}}\left(\frac{M_{\rm eq}}{M_{\rm
H}}\right)^{1/2}\left(\frac{M_{\rm PBH}}{\mathcal{K}M_{\rm
H}}\right)^{\frac{1+\gamma}{\gamma}}$ $\displaystyle\times\frac{{\cal
K}}{\gamma\sqrt{\Lambda}}\int{\rm d}\zeta_{{\rm G},k}P_{\rm G}({\cal C}_{\rm
G}(M_{\rm PBH},\zeta_{\rm G}),\zeta_{\rm G})\left(\frac{dF}{d\zeta_{\rm
G}}\right)^{-1}\\!\\!\\!\\!,$ (14)

where $\Lambda=1-({\cal C}_{\rm th}-\left(M_{\rm
PBH}/(\mathcal{K}M_{H})\right)^{1/\gamma})/\Phi$.

#### Theory of peaks

It is well known that the Press-Schechter approach does not agree with the
theory of peaks (see, e.g., Refs. Green  _et al._ (2004); Young  _et al._
(2014); De Luca  _et al._ (2019)). A technical drawback of peaks theory is
that it is not clear how to include NGs in the computation of the abundance in
this approach using a generic functional form for the curvature perturbation
field as in Eq. (4). Hence, for comparison, we limit our analysis for the peak
theory approach to the case of negligible primordial NGs222The impact of
local-type NG on PBH abundance in peak theory is expected to be reduced,
especially for higher-order terms Young (2022) and for peaked power spectra.
Indeed, the compaction is volume-averaged over the scale of the perturbation
and depends on the curvature perturbation at its edge, rather than its
centre..

In this approach, the starting point is the number density of peaks with a
height in the interval $(\nu,\nu+{\rm d}\nu)$. When $\nu\gg 1$, it is given by
Bardeen  _et al._ (1986)

$\frac{{\rm d}\mathcal{N}}{{\rm
d}\nu}=\frac{1}{4\pi^{2}r_{m}^{3}}\left(\frac{\sigma_{cc}}{\sigma_{c}}\right)^{3}\nu^{3}\exp\left(-\frac{\nu^{2}}{2}\right)\,,$
(15)

where we introduced the rescaled peak height
$\nu\equiv\mathcal{C}_{1}/\sigma_{c}$ and the first rescaled moment of the
distribution Balaji  _et al._ (2023)

$\sigma_{cc}^{2}=\frac{4\Phi^{2}}{9}\int_{0}^{\infty}\frac{{\rm
d}k}{k}\left(kr_{m}\right)^{6}W^{2}\left(k,r_{m}\right)P^{T}_{\zeta}\,.$ (16)

The number of peaks within a Hubble volume is therefore

$\frac{{\rm d}N_{k}}{{\rm d}\nu}=\frac{4\pi}{3}r_{m}^{3}\mathcal{N}$ (17)

and, since $N_{k}\ll 1$, the probability of finding at least one peak within a
range $(\nu,\nu+{\rm d}\nu)$ within a Hubble volume is $P_{k}\approx{\rm
d}N_{k}$. The mass fraction of PBHs is then Yoo  _et al._ (2018, 2019); Gow
_et al._ (2021); Franciolini  _et al._ (2022)

$P_{k}(\mathcal{C})=\frac{1}{3\pi}\left(\frac{\sigma_{cc}}{\sigma_{c}}\right)^{3}\nu^{3}\exp\left(-\frac{\nu^{2}}{2}\right)$
(18)

and the domain of integration agrees with Eq. 8. Thus, peaks theory predicts
an additional $\nu^{3}$ factor to the collapse probability in the Gaussian
case.

Recasting the quantities in terms of $M_{k}$, as in (III), we can express the
PBH abundance as

$\displaystyle\frac{{\rm d}f_{\rm PBH}}{{\rm d}\ln M_{\rm PBH}}=$
$\displaystyle\frac{1}{\Omega_{\rm DM}}\int_{M_{\rm H}^{\rm min}}\frac{dM_{\rm
H}}{M_{\rm H}}\left(\frac{M_{\rm eq}}{M_{\rm
H}}\right)^{1/2}\left(\frac{M_{\rm PBH}}{\mathcal{K}M_{\rm
H}}\right)^{\frac{1+\gamma}{\gamma}}$ $\displaystyle\times\frac{{\cal
K}}{\gamma}\frac{\left(2\Phi\left(1-\sqrt{\Lambda}\right)\right)^{3}}{3\pi\sigma_{c}^{4}\Lambda^{1/2}}\left(\frac{\sigma_{cc}}{\sigma_{c}}\right)^{3}$
$\displaystyle\times\exp\left[{\frac{-2\Phi^{2}}{\sigma_{c}^{2}}\left(1-\sqrt{\Lambda}\right)^{2}}\right]\,.$
(19)

##  IV Scalar-induced gravitational waves

Primordial curvature fluctuations large enough to generate PBHs will also
induce a non-negligible SIGW background. At the second order in perturbation
theory, scalar modes can source tensor modes (for a review, see Ref. Domènech
(2021)). Starting from Einstein’s equations and omitting anisotropic stress,
the equations of motion for GWs in Fourier space are

$h^{\prime\prime}_{\lambda}(\textbf{k},\eta)+2\mathcal{H}h^{\prime}_{\lambda}(\textbf{k},\eta)+k^{2}h_{\lambda}(\textbf{k},\eta)=4\mathcal{S}_{\lambda}(\textbf{k},\eta)\,,$
(20)

with $\lambda$ indicating the polarization and the ′ being the derivative with
respect to the conformal time $\eta$. The source term on the right-hand side
reads Baumann  _et al._ (2007)

$\mathcal{S}_{\lambda}(\textbf{k},\eta)\\!=\\!\int\frac{{\rm
d}^{3}\textbf{q}}{(2\pi)^{\frac{3}{2}}}Q_{\lambda}(\textbf{k},\textbf{q})f(|\textbf{k}-\textbf{q}|,q,\eta)\zeta_{\textbf{q}}\zeta_{\textbf{k}-\textbf{q}}\,,$
(21)

where $Q$ is a projection factor and $f$ is a function that contains the
linear evolution of the Newtonian potential after horizon re-entry and the
transfer function mapping them to the curvature perturbation $\zeta$. Eq. (20)
can be solved using the Green’s functions approach

$h_{\lambda}(\textbf{k})\\!=\\!4\int\frac{{\rm
d}^{3}\textbf{q}}{(2\pi)^{\frac{3}{2}}}Q_{\lambda}(\textbf{k},\textbf{q})I(|\textbf{k}-\textbf{q}|,q,\eta)\zeta_{\textbf{q}}\zeta_{\textbf{k}-\textbf{q}}\,,$
(22)

where $I$ is the transfer function. To speed up our analysis, we assume
perfect radiation domination and do not account for the variation of sound
speed during the QCD era (see, for example, Hajkarim and Schaffner-Bielich
(2020); Abe  _et al._ (2021)) which also leads to specific imprints in the
low-frequency tail of any cosmological SGWB Franciolini  _et al._ (2024). On
top of that, cosmic expansion may additionally be affected by unknown physics
in the dark sector, which can, e.g., lead to a brief period of matter
domination of kination Ferreira and Joyce (1998); Pallis (2006); Redmond  _et
al._ (2018); Co  _et al._ (2022); Gouttenoire  _et al._ (2021); Chang and Cui
(2022). Explicit analytical expressions for the kernel in different epochs can
be found in Domènech (2020), after performing the oscillation average well
inside the horizon.

The GW dimensionless power spectrum $\mathcal{P}_{h}(k)$ is defined as

$\langle
h_{\lambda_{1}}(\mathbf{k}_{1})h_{\lambda_{2}}(\mathbf{k}_{2})\rangle=\delta^{3}(\mathbf{k}_{1}+\mathbf{k}_{2})\delta_{\lambda_{1}\lambda_{2}}\frac{2\pi^{2}}{k^{3}}\mathcal{P}_{h,\lambda_{1}}(k_{1})$
(23)

and hence

$\displaystyle\Omega_{\rm
GW}(k,\eta)=\frac{1}{48}\left(\frac{k}{a(\eta)H(\eta)}\right)\sum_{\lambda=+,\times}\overline{\mathcal{P}_{h,\lambda}(k)}$
(24)

at the production. In the Gaussian case, one obtains

$\displaystyle\Omega_{\rm GW}(k,\eta)|_{\rm G}=$ $\displaystyle\hskip
2.84544pt\frac{1}{12}\left(\frac{k}{a(\eta)H(\eta)}\right)^{2}\int_{0}^{\infty}{\rm{d}}t\int_{-1}^{1}{\rm{d}}s\hskip
2.84544pt\frac{1}{u^{2}v^{2}}$
$\displaystyle\times\overline{\tilde{J}^{2}(u,v,x)}\hskip
2.84544pt\mathcal{P}_{\zeta_{g}}\mathcal{P}_{\zeta_{g}}(vk)\mathcal{P}_{\zeta_{g}}(uk)\,,$
(25)

where $x=k\eta$, $t$ and $s$ are dimensionless variables introduced to ease
the numerical analysis, $\tilde{J}$ is related to the transfer function $I$
defined before and the overbar indicates an oscillation average after the
modes are well inside the horizon. Furthermore $u=(t+s+1)/2$ and
$v=(t-s+1)/2$. The spectrum observed today can be obtained using entropy
conservation

$h^{2}\Omega_{\rm GW}|_{k,0}=h^{2}\Omega_{\rm
rad,0}\left(\frac{g_{s,\eta_{0}}}{g_{s,\eta_{k}}}\right)^{\frac{4}{3}}\frac{g_{\rho,\eta_{k}}}{g_{\rho,\eta_{0}}}\Omega_{\rm
GW}(k,\eta)\,,$ (26)

where we used the energy ($g_{\rho}$) and entropy density ($g_{s}$) degrees of
freedom of Ref. Borsanyi  _et al._ (2016). According to Ref. Aghanim  _et al._
(2020), we consider $h^{2}\Omega_{\rm rad,0}=4.2\cdot 10^{-5}$.

###  IV.1 Effects of primordial NG on SIGW

As shown in Kohri and Terada (2018); Cai  _et al._ (2019a, b); Hajkarim and
Schaffner-Bielich (2020); Atal and Domènech (2021); Yuan and Huang (2021);
Domènech  _et al._ (2022); Garcia-Saenz  _et al._ (2023); Liu  _et al._
(2024); Yuan  _et al._ (2023); Li  _et al._ (2023); Unal (2019); Adshead  _et
al._ (2021); Abe  _et al._ (2023); Perna  _et al._ (2024) primordial NGs can
leave a non-negligible imprint on the SIGW power spectrum. Hence, considering
only the quadratic contribution $f_{\rm NL}$ in the local expansion in Eq.
(5)333We specify that, in principle, all the terms in the power series
expansion in Eq. (5) should be included. Being interested mainly in the
qualitative effect of NGs, we truncated the power series at the quadratic
contribution, $f_{\rm NL}$, setting all other terms to 0. These additional NG
corrections would further shift the results in Fig. 3., at the next-to-leading
order one obtains (see, e.g., Unal (2019); Adshead  _et al._ (2021))

$\displaystyle\Omega_{\rm GW}(k,\eta)|_{\rm t}$
$\displaystyle=\frac{1}{12\pi}\left(\frac{k}{a(\eta)H(\eta)}\right)^{2}f_{\rm
NL}^{2}$ $\displaystyle\times$
$\displaystyle\int_{0}^{\infty}{\rm{d}}t_{1}\int_{-1}^{1}{\rm{d}}s_{1}\int_{0}^{\infty}{\rm{d}}t_{2}\int_{-1}^{1}{\rm{d}}s_{2}\int_{0}^{2\pi}{\rm{d}}\varphi_{12}$
$\displaystyle\times$ $\displaystyle\cos
2\varphi_{12}\frac{u_{1}v_{1}}{(u_{2}v_{2})^{2}}\frac{1}{w_{a,12}^{3}}\overline{\tilde{J}(u_{1},v_{1},x)\tilde{J}(u_{2},v_{2},x)}$
$\displaystyle\times$
$\displaystyle\mathcal{P}_{\zeta_{g}}(v_{2}k)\mathcal{P}_{\zeta_{g}}(u_{2}k)\mathcal{P}_{\zeta_{g}}(w_{a,12}k)\,,$
(27) $\displaystyle\Omega_{\rm GW}(k,\eta)|_{\rm u}$ $\displaystyle=\hskip
2.84544pt\frac{1}{12\pi}\left(\frac{k}{a(\eta)H(\eta)}\right)^{2}f_{\rm
NL}^{2}$ $\displaystyle\times$
$\displaystyle\int_{0}^{\infty}{\rm{d}}t_{1}\int_{-1}^{1}{\rm{d}}s_{1}\int_{0}^{\infty}{\rm{d}}t_{2}\int_{-1}^{1}{\rm{d}}s_{2}\int_{0}^{2\pi}{\rm{d}}\varphi_{12}$
$\displaystyle\times$ $\displaystyle\cos
2\varphi_{12}\frac{u_{1}u_{2}}{(v_{1}v_{2})^{2}}\frac{1}{w_{b,12}^{3}}\overline{\tilde{J}(u_{1},v_{1},x)\tilde{J}(u_{2},v_{2},x)}$
$\displaystyle\times$
$\displaystyle\mathcal{P}_{\zeta_{g}}(v_{1}k)\mathcal{P}_{\zeta_{g}}(v_{2}k)\mathcal{P}_{\zeta_{g}}(w_{b,12}k)$
(28) $\displaystyle\Omega_{\rm GW}(k,\eta)|_{\rm hyb}$ $\displaystyle=\hskip
2.84544pt\frac{1}{12}\left(\frac{k}{a(\eta)H(\eta)}\right)^{2}f_{\rm NL}^{2}$
$\displaystyle\times$
$\displaystyle\int_{0}^{\infty}{\rm{d}}t_{1}\int_{-1}^{1}{\rm{d}}s_{1}\int_{0}^{\infty}{\rm{d}}t_{2}\int_{-1}^{1}{\rm{d}}s_{2}$
$\displaystyle\times$ $\displaystyle\frac{1}{(u_{1}v_{1}u_{2}v_{2})^{2}}\hskip
2.84544pt\overline{\tilde{J}^{2}(u_{1},v_{1},x)}$ $\displaystyle\times$
$\displaystyle\mathcal{P}_{\zeta_{g}}(u_{1}k)\mathcal{P}_{\zeta_{g}}(v_{2}v_{1}k)\mathcal{P}_{\zeta_{g}}(u_{2}v_{1}k)\,.$
(29)

As in the Gaussian case, the $u_{i}$ and $v_{i}$ are related similarly to the
integration variables $t_{i}$ and $s_{i}$.

The GW spectrum can then be expressed as

$\displaystyle\Omega_{{\rm GW}}|_{k,0}$ $\displaystyle=A^{2}\Omega_{\rm
GW}|_{{\rm G}}(k/k_{*},0)$ (30) $\displaystyle+A^{3}f_{\rm NL}^{2}\Omega_{\rm
GW}|_{{\rm NG}}(k/k_{*},0)\,,$

where $\Omega_{\rm GW}|_{{\rm NG}}$ is the sum of (27), (28) and (29) with
$A=1$, $f_{\rm NL}=1$ computed assuming $k_{*}=1$.444Factoring out $k_{*}$
dependence in this way is possible in pure radiation dominance, that is, when
QCD effects are neglected. In this way, given the spectral shape is fixed, the
GW spectrum can be straightforwardly computed for given $A$, $k_{*}$ and
$f_{\rm NL}$.

Figure 1: Left panel: Constraints from NANOGrav15 Agazie  _et al._ (2023a) on
the amplitude of a log-normal curvature power spectrum (1) with $\Delta=0.5$
assuming Gaussian fluctuations and changing the window function. The
horizontal lines correspond to $f_{\rm PBH}=1$ using threshold statistics
(solid) and theory of peaks (dashed), respectively, assuming NGs arising
solely from the non-linear corrections. Right panel: Inferred constraints on
PBH abundance.

##  V Constraints on PBHs

###  V.1 PTA observations

Consider the constraints on the SIGW signal from the current PTA observations.
Using from the NANOGrav 15 year dataset Agazie  _et al._ (2023a);
Collaboration (2023), we compute the observational upper bound at 95%
confidence level on $h^{2}\Omega^{\rm bound}_{\rm GW}$ in each of the first 14
frequency bins. The frequency of bin ”$i$” is given by $f_{i}=i/T$, where
$T=16.03\,{\rm yr}=(1.98\,{\rm nHz})^{-1}$ is the time of observation. We then
compare the observation with theoretical prediction $h^{2}\Omega_{\rm GW}$ in
Eq. (30) in each frequency bin. The constraint on $h^{2}\Omega_{\rm GW}$ in
any given bin implies an upper bound on $A$ for a given model of NGs and a
fixed shape of the curvature power spectrum. To combine the constraints of
individual bins, we consider the strongest constraint on $A$. The upper bound
on $A$ is then mapped to an upper bound on $f_{\rm PBH}$. More specifically,
for every fixed benchmark model in Table 1, we can map the pair $(A,k_{*})$
uniquely into the pair $(f_{\rm PBH},\langle M_{\rm PBH}\rangle)$.

Due to the finite observational time, the spectral resolution of the detector
is limited and can be estimated to be approximately $\Delta f=1/T=1.98\,{\rm
nHz}$. This will limit the PTA sensitivity to sharp features in the SIGW
spectrum. To account for this, the contribution to $i$-th bin can be estimated
as

$\Omega_{{\rm GW},i}=\int\frac{{\rm d}f}{f}\Omega_{\rm GW}^{\rm
th}(f)W(f-f_{i})\,$ (31)

where $W$ denotes a window function. To understand the effect of limited
spectral resolution, we will consider three different idealized cases

$\displaystyle W_{\rm D}(f)$ $\displaystyle\propto\delta\left(f\right)\,,$
(32) $\displaystyle W_{\rm S}(f)$ $\displaystyle\propto{\rm
sinc}^{2}\left(\pi\frac{f}{\Delta f}\right)\,,$ (33) $\displaystyle W_{\rm
TH}(f)$ $\displaystyle\propto\theta\left(f+\frac{\Delta
f}{2}\right)\theta\left(\frac{\Delta f}{2}-f\right)\,,$ (34)

which correspond to an infinitely accurate frequency resolution, a sharp
binning of frequencies and a time-domain top-hat filter on $h_{\rm c}(t)$ with
duration $T$, respectively. The window functions are normalized so that
$\int{\rm d}\ln f\,W(f-f_{i})=1$. These window functions represent different
limiting cases and can thus capture systematics related to the spectral
resolution not otherwise accessible in our simplified analysis. As it can be
seen from Fig. 1, the effect is expectedly larger at large scales and for
heavy PBHs as these are mostly related to the signal at small frequencies. For
instance, at $k\simeq 10^{6}\,{\rm{Mpc}}^{-1}$, $W_{\rm D}(f)$ gives a mildly
stronger constraint on $A$ when compared to the other choices, as it can
resolve sharp peaks in the SIGW spectrum. This is reflected in the constraints
on the PBH abundance, as shown on the right panel of Fig. 1. This behaviour is
reverted at even smaller $k_{*}$, and $W_{\rm D}(f)$ is less constraining than
the other choices. For the remainder of the analysis, we will consider the
Dirac window function $W_{\rm D}(f)$.

Figure 2: Left panel: Constraints from NANOGrav15 Agazie  _et al._ (2023a) on
the amplitude of a log-normal curvature power spectrum, Eq. (1), with
$\Delta=0.5$ (blue), and on a broken power-law curvature power spectrum, Eq.
(2), with $\alpha=4$, $\gamma=1$ for $\beta=3$ (pink) and $\beta=0.5$ (red),
assuming Gaussian fluctuations. The horizontal lines correspond to $f_{\rm
PBH}=1$ using threshold statistics (solid) and theory of peaks (dashed),
respectively, given the non-linearity-only scenario. Right panel: The implied
NANOGrav15 constraints on PBH abundance.

Figure 3: Left panels: Constraints from NANOGrav15 Agazie  _et al._ (2023a) on
the amplitude of a log-normal curvature power spectrum (1) with $\Delta=0.5$
(top panels) and on a broken power-law curvature power spectrum (2) with
$\alpha=4$, $\gamma=1$ for $\beta=3$ (middle panels) and $\beta=0.5$ (bottom
panels), assuming quadratic primordial NGs only with $|f_{\rm NL}|\in[0,10]$.
The nearly horizontal lines correspond to $f_{\rm PBH}=1$ using threshold
statistics with several $f_{NL}$ values. Right panels: Inferred constraints on
the PBH abundance using threshold statistics.

Another important theoretical uncertainty when estimating constraints on the
PBH abundance arises due to the choice of the prescription when computing the
PBH abundance. The right panels of Figs. 1 and 2 show the constraints on the
PBH abundance arising when the abundance is computed using threshold
statistics (solid lines) and peaks theory (dashed lines). In both figures,
primordial NGs are neglected ($f_{\rm NL}=0$) so the NGs arise only from the
non-linear relation between density and curvature perturbations. As one can
see, almost the entire parameter space accessible by LVK via GWs from PBH
binary mergers Hütsi  _et al._ (2021); Romero-Rodriguez  _et al._ (2022);
Franciolini  _et al._ (2022); Andrés-Carcasona  _et al._ (2024) is excluded
when threshold statistics is assumed. On the other hand, with the peaks
theory, only a small portion of the parameter space relevant to LVK is
constrained. Despite this significant theoretical uncertainty, PTA
observations can constrain the PBH abundance in both cases.

The dependence of the constraints on the shape of the power spectrum is
illustrated in Fig. 2. While relatively narrow spectra, such as the LN and the
BPL1 cases, produce similar results, broader spectra, represented by the BPL2
case are constrained also at larger scales due to the broadness of the
spectrum. We remark that the constraints at large scales $k_{*}\lesssim
5\times 10^{6}{\rm Mpc}^{-1}$ arise mostly from the first bin by NANOGrav15.
This is also the region in which the SIGW are most strongly constrained. In a
similar vein, the constraints on the PBH abundance in the threshold statistics
case are quite similar for all spectra except at large masses $\langle M_{\rm
PBH}\rangle$. In particular, the strongest constraints touch the dashed region
in Fig. 2 which corresponds to having less than a single PBH within the
current Hubble volume. This line is given by $f_{\rm PBH}4\pi\Omega_{\rm
DM}M_{\rm pl}^{2}/H_{0}<\langle M_{\rm PBH}\rangle$.

Fig. 3 shows the effect of primordial NGs. Firstly, we underline that the NG
corrections to SIGW depend on $f_{\rm NL}^{2}$ and are thus indifferent to the
sign of $f_{\rm NL}$. On the other hand, this latter strongly impacts $f_{\rm
PBH}$. Indeed, negative primordial NGs suppress the tail of the PDF with
respect to the Gaussian case, that is, they reduce the probability of large
fluctuations that cross the threshold. So, a larger amplitude would be needed
to obtain the same PBH abundance. Positive primordial NGs, on the other hand,
have the opposite effect as they make large fluctuations more likely.

As shown on the left panels of Fig. 3, for small $f_{\rm NL}$ we observe that
the corrections to the SIGW affect only amplitudes at the tails, both at
smaller and higher scales. A larger positive $f_{\rm NL}$, on the other hand,
further shifts the constant $f_{\rm PBH}$ curves downward. As expected, this
reduces the allowed amplitude $A$ for a given $k_{*}$. Conversely, for a fixed
PBH abundance, positive primordial NGs decrease the required amplitude $A$
with respect to the Gaussian case. Since the latter effect is larger than the
reduction of the SIGW strength, positive primordial NGs can completely remove
constraints on the PBH abundance.

The situation differs when considering negative NGs. Since negative NGs imply
an upward shift in the power spectrum amplitude in order to achieve the same
PBH abundance, the constraints become more stringent.

One might mistakenly assume that by increasing the negative value of the
coefficient $f_{\rm NL}$, the required power spectral amplitude for the same
PBH abundance would always rise, tightening the constraints indefinitely.
However, this is not the case. As shown in Ref. Franciolini  _et al._ (2023a),
the required amplitude actually reaches a maximum for $f_{\rm NL}\simeq-2$ and
then decreases, eventually approaching the amplitude of the Gaussian case for
very large negative values of $f_{\rm NL}$. Consequently, the constraints for
the case of $f_{\rm NL}=-2$ can be considered as the most stringent.

Nevertheless, this claim should be taken cum grano salis. Indeed, it is
important to note that the prescription outlined in Ref. Musco  _et al._
(2021) to compute the threshold for PBH collapse only accounts for NGs arising
from the non-linear relation between the density contrast and the curvature
perturbations. In principle, also primordial NGs beyond the quadratic
approximation should be taken into account when computing the threshold value.
Following Refs. Kehagias  _et al._ (2019); Escrivà  _et al._ (2022), it
appears that the effect on the threshold is small only for positive NGs and
for small negative NGs ($f_{\rm NL}\gtrsim-2$) and at most of the order of a
few percent. For more negative NGs, there are huge uncertainties on the exact
value of the threshold, which could potentially modify the constraints
introduced above555One alternative approach is to use the averaged value over
a sphere of radius equal to the location of the maximum of the compaction
function, $\bar{C}_{\rm th}=2/5$, since this value is fully non-perturbative
and independent of NGs Escrivà  _et al._ (2020); Kehagias  _et al._ (2024). As
discussed in Ref. Ianniccari  _et al._ (2024), this requires determining the
statistics of the curvature perturbation and computing the connected
cumulants, which is a complex task and is left for future work..

For completeness, we analyze a common scenario for PBH formation.

In _Ultra slow roll_ (USR) models of single-field inflation, the peak in
$\mathcal{P}_{\zeta}$ arises from a brief phase of ultra-slow-roll followed by
constant-roll inflation dual to it Atal and Germani (2019); Biagetti  _et al._
(2018); Karam  _et al._ (2023a). In this case, the NGs can be related to the
large $k$ spectral slope Atal  _et al._ (2019); Tomberg (2023),

$\zeta=-\frac{2}{\beta}\log\left(1-\frac{\beta}{2}\zeta_{\rm G}\right).$ (35)

In Fig. 4 we show a similar plot as in Fig. 2 but computing the abundance of
PBHs using the full non-gaussian relation reported in Eq.(35) for the two
broken power law cases reported in Tab. 1. We compute the abundance using
exclusively the method based on the threshold statistics on the compaction
function, and as we can see from the plot, with respect to the Gaussian
approximation the constraints are highly relaxed and, differently from the
Gaussian case, they are of the same order of magnitude. The reason is the
following: expanding Eq. (35) at second order, we get $f_{\rm NL}=5\beta/12$,
then using the reference points $\beta=3$ and $\beta=0.5$ shown in Fig. 4, one
finds $f_{\rm NL}=5/4$ and $f_{\rm NL}=5/24$. As a consequence, the narrower
power spectrum, which has stronger constraints in the Gaussian approximation,
due to larger primordial non-gaussianities, has constraints similar to the
broader case.

Moreover, with respect the impact on the PBH abundance, the NG corrections to
the SIGW are irrelevant.

Figure 4: The NANOGrav15 constraints on the PBH abundance using threshold
statistics assuming the two cases for the broken power law reported in Tab. 1
assuming the gaussian approximation (solid) and with the full non-gaussian
relation in Eq. (35).

The other constraints reported here are GW O3 Andrés-Carcasona  _et al._
(2024), EROS Tisserand  _et al._ (2007), OGLE Mroz  _et al._ (2024a, b), Seg1
Koushiappas and Loeb (2017), Planck Serpico  _et al._ (2020); Agius  _et al._
(2024); Facchinetti  _et al._ (2023), Eri II Brandt (2016), WB Monroy-
Rodríguez and Allen (2014), Ly$-\alpha$ Murgia  _et al._ (2019) and SNe
Zumalacarregui and Seljak (2018). These constraints are shown for
monochromatic mass functions.

Figure 5: Constraints on the amplitude of the power spectrum, assuming
negligible primordial NGs, from the FIRAS experiments, CMB, $\alpha$-Lyman and
NANOGrav15 for the same cases reported in Tab. 1.

###  V.2 $\mu$ distortions and PBH seeds for SMBHs

At the largest scales, the primordial power spectrum is strongly constrained
by the CMB observations Aghanim  _et al._ (2020) and the Lyman-$\alpha$ forest
data Bird  _et al._ (2011). The CMB observations strongly constrain the
curvature power spectrum at scales $10^{-4}\mathrm{Mpc}^{-1}\lesssim k\lesssim
1\mathrm{Mpc}^{-1}$. At redshifts $z\lesssim 10^{6}$, energy injections into
the primordial plasma cause persisting spectral distortions in the CMB. These
distortions are divided into chemical potential $\mu$-type distortions created
at early times and Compton $y$-type distortions created at $z\lesssim 5\times
10^{4}$. For a given curvature power spectrum $\mathcal{P}_{\zeta}(k)$ the
spectral distortions are Chluba  _et al._ (2012b); Chluba and Grin (2013)

$X=\int_{-\infty}^{\infty}\frac{\mathrm{d}k}{k}\mathcal{P}_{\zeta}(k)W_{X}(k)$
(36)

with $X=\mu,y$ where $k_{\min}=1$ $\mathrm{Mpc}^{-1}$ and the window functions
can be approximated by

$W_{\mu}(k)=2.2\left[e^{-\frac{(\hat{k}/1360)^{2}}{1+(\hat{k}/260)^{0.6}+\hat{k}/340}}-e^{-(\hat{k}/32)^{2}}\right],$
(37) $W_{y}(k)=0.4e^{-(\hat{k}/32)^{2}}$ (38)

with $\hat{k}=k/\left(1\mathrm{Mpc}^{-1}\right)$. The COBE/Firas observations
constrain the $\mu$ distortions as $\mu\leq 4.7\times 10^{-5}$ Bianchini and
Fabbian (2022) and $y\leq 1.5\times 10^{-5}$ at $95\%$ confidence level Fixsen
_et al._ (1996). The situation is summarised in Fig. 5. Including NGs modifies
only the tails of the constraints and not the most constrained region Sharma
_et al._ (2024).

Consequently, the amplitude of the power spectrum in the range of scales
related to the FIRAS experiment is constrained to be at most $A=10^{-5}$.

We focus on the threshold statistics approach and we compute the mass fraction
$\beta$ as a function of the parameters $f_{\rm NL}$ and $g_{\rm NL}$,
assuming the power series ansatz as in Eq. 5. The power series expansion holds
until we are in the perturbative regime, that is, the terms are ordered
hierarchically with the higher orders being typically smaller than the lower
ones. For a narrow power spectrum, this translates into a constraint for the
amplitude of the power spectrum, $(3/5)\left|f_{\mathrm{NL}}\right|A^{1/2}\ll
1$, and $(9/25)\left|g_{\mathrm{NL}}\right|A\ll 1$.

Such large positive NGs can lead to the production of a non-negligible
population of massive PBHs that can seed SMBHs.

Figure 6: Left panel: Computation of the mass fraction $\beta$ changing the NG
parameters $f_{\rm NL}$ and $g_{\rm NL}$ with two benchmark cases for the log-
normal power spectrum with an amplitude fixed to be $A=10^{-5}$. Right panel:
Computation of the PBH abundance $f_{\rm PBH}$ changing the scale of the peak
$k_{*}$ and the NG parameter $g_{\rm NL}$. We fix $\Delta=1$, $A=10^{-5}$ and
$f_{\rm NL}=1$.

Indeed, the origin of SMBHs presents a significant challenge in the field of
astrophysics. Although it is well known that they occupy the centres of most
galaxies Kormendy and Richstone (1995); Magorrian  _et al._ (1998); Richstone
_et al._ (1998), the processes leading to their formation remain unclear.

In the context of PBHs, even a minimal presence of massive PBHs within the
mass range $10^{3}-10^{6}\,M_{\odot}$ could potentially seed for SMBHs
Duechting (2004); Kohri  _et al._ (2014); Bernal  _et al._ (2018). As a naive
computation to estimate the necessary abundance of these seeds, we follow
Refs. Vaskonen and Veermäe (2021); Serpico  _et al._ (2020). Considering that
SMBHs make up approximately $\mathcal{O}(10^{-4})$ of the stellar mass in
their host galaxies Reines and Volonteri (2015), and that stars account for
roughly $\mathcal{O}(10^{-2})$ of the total cosmic matter content Fukugita
_et al._ (1998), we deduce that the overall SMBH density is about $10^{6}$
times less than the dark matter density.

The primordial seeds’ abundance can then be expressed as

$f_{\mathrm{PBH}}\sim\mathcal{O}(10^{-6})\times\left\langle
M_{\mathrm{PBH}}\right\rangle/M_{\rm SMBH}\,,$ (39)

with $\left\langle M_{\mathrm{PBH}}\right\rangle>10^{3}M_{\odot}$.
Nevertheless, the scales related to the formation of SMBHs are strongly
constrained by the analysis of the CMB spectral distortion by the FIRAS
collaboration Fixsen  _et al._ (1996); Chluba  _et al._ (2012a, b); Chluba and
Grin (2013); Bianchini and Fabbian (2022), and in the gaussian approximation,
the corresponding abundance of PBHs is too small to furnish a primordial
origin for Supermassive black holes seeds.

As shown in the left panel of Fig. 6, in agreement with the literature Ünal
_et al._ (2021); Nakama  _et al._ (2018); Hooper  _et al._ (2024); Byrnes  _et
al._ (2024), quadratic primordial NGs are not enough to produce a sizeable
amount of PBHs to be the seed of SMBHs. As a consequence, differently from
Ref. Byrnes  _et al._ (2024), we find that already a sizeable cubic parameter
$g_{\rm NL}$ is enough to overcome this hurdle. This discrepancy arises
because their analysis is subject to a few simplifications: the use of the
curvature perturbation instead of the density contrast to determine the
collapse and the adoption a different threshold that is independent of the
shape of the curvature power spectrum.

According to Fig. 6, assuming a typical SMBH mass of approximately $10^{8}$
$M_{\odot}$, condition (39), without the violation of the perturbative
criterium, holds across the range of PBH masses analyzed.

##  VI Conclusions and outlook

In this work, we have updated the constraints on the amplitude of the
curvature power spectrum, ensuring that the amplitude of the produced SIGWs is
not larger than the signal recently detected by the NANOGrav collaboration.

Considering two approaches to determine the PBH abundance, we have recast
constraints on the SIGW background into constraints on PBHs. We found that,
regardless of the specific shape of the power spectrum, when we determine the
abundance using threshold statistics on the compaction function, only a small
fraction of DM in the form of PBHs is available in the solar mass range. For
narrow spectra we find that peaked power spectrum in the range around $20-50$
$M_{\odot}$ cannot exist and significant constraints ($f_{\rm PBH}\lesssim
10^{-5}$) in the sub-solar mass range are present, while they are less
constraining for broader spectra.

This tension can be significantly relaxed by introducing positive NGs, even
small as those that naturally arise in USR models.

When we use the method based on peak theory, we find that these constraints
are of the same order of magnitude or less stringent than those already
existing.

We then discuss how, by including large primordial NGs ($g_{\rm NL}\simeq
10^{5}$), it is possible to produce seeds for SMBHs without violating the
perturbativity criterion. Despite this, a major limitation of our claim is
that currently, there are no models in the literature with such large positive
NGs. Indeed, in commons scenarios, such as curvaton and USR models,
considering the full NG relation, instead of a power series truncated at
second order, reduces the general impact of the NGs Ferrante  _et al._ (2023a,
b).

To improve future analysis, we not only need more data to better reconstruct
the signal, but it is also crucial to maximize the refinement of the PBH
abundance calculation, which is still currently subject to a series of
theoretical uncertainties.

###### Acknowledgements.

We thank G. Franciolini, A. Ianniccari and D. Perrone for useful discussions.
A.J.I. acknowledges the financial support provided under the “Progetti per
Avvio alla Ricerca Tipo 1”, protocol number AR1231886850F568. H.V. is
supported by the Estonian Research Council grants PSG869 and RVTT7 and the
Center of Excellence program TK202.

## References

  * Zel’dovich and Novikov (1967) Y. B. Zel’dovich and I. D. Novikov, Sov. Astron. 10, 602 (1967).
  * Hawking (1971) S. Hawking, Mon. Not. Roy. Astron. Soc. 152, 75 (1971).
  * Carr and Hawking (1974) B. J. Carr and S. W. Hawking, Mon. Not. Roy. Astron. Soc. 168, 399 (1974).
  * Bagui  _et al._ (2023) E. Bagui _et al._ (LISA Cosmology Working Group), (2023), arXiv:2310.19857 [astro-ph.CO] .
  * Ivanov  _et al._ (1994) P. Ivanov, P. Naselsky, and I. Novikov, Phys. Rev. D 50, 7173 (1994).
  * Ivanov (1998) P. Ivanov, Phys. Rev. D 57, 7145 (1998), arXiv:astro-ph/9708224 .
  * Blinnikov  _et al._ (2016) S. Blinnikov, A. Dolgov, N. K. Porayko, and K. Postnov, JCAP 11, 036 (2016), arXiv:1611.00541 [astro-ph.HE] .
  * Akrami  _et al._ (2020a) Y. Akrami _et al._ (Planck), Astron. Astrophys. 641, A10 (2020a), arXiv:1807.06211 [astro-ph.CO] .
  * Matarrese  _et al._ (1994) S. Matarrese, O. Pantano, and D. Saez, Phys. Rev. Lett. 72, 320 (1994), arXiv:astro-ph/9310036 .
  * Acquaviva  _et al._ (2003) V. Acquaviva, N. Bartolo, S. Matarrese, and A. Riotto, Nucl. Phys. B 667, 119 (2003), arXiv:astro-ph/0209156 .
  * Mollerach  _et al._ (2004) S. Mollerach, D. Harari, and S. Matarrese, Phys. Rev. D 69, 063002 (2004), arXiv:astro-ph/0310711 .
  * Ananda  _et al._ (2007) K. N. Ananda, C. Clarkson, and D. Wands, Phys. Rev. D 75, 123518 (2007), arXiv:gr-qc/0612013 .
  * Baumann  _et al._ (2007) D. Baumann, P. J. Steinhardt, K. Takahashi, and K. Ichiki, Phys. Rev. D 76, 084019 (2007), arXiv:hep-th/0703290 .
  * Espinosa  _et al._ (2018) J. R. Espinosa, D. Racco, and A. Riotto, JCAP 09, 012 (2018), arXiv:1804.07732 [hep-ph] .
  * Domènech (2021) G. Domènech, Universe 7, 398 (2021), arXiv:2109.01398 [gr-qc] .
  * Gangui  _et al._ (1994) A. Gangui, F. Lucchin, S. Matarrese, and S. Mollerach, Astrophys. J. 430, 447 (1994), arXiv:astro-ph/9312033 .
  * Bartolo  _et al._ (2002) N. Bartolo, S. Matarrese, and A. Riotto, Phys. Rev. D 65, 103505 (2002), arXiv:hep-ph/0112261 .
  * Maldacena (2003) J. M. Maldacena, JHEP 05, 013 (2003), arXiv:astro-ph/0210603 .
  * Bartolo  _et al._ (2004a) N. Bartolo, S. Matarrese, and A. Riotto, Phys. Rev. D 69, 043503 (2004a), arXiv:hep-ph/0309033 .
  * Bartolo  _et al._ (2004b) N. Bartolo, E. Komatsu, S. Matarrese, and A. Riotto, Phys. Rept. 402, 103 (2004b), arXiv:astro-ph/0406398 .
  * Celoria and Matarrese (2020) M. Celoria and S. Matarrese, Proc. Int. Sch. Phys. Fermi 200, 179 (2020), arXiv:1812.08197 [astro-ph.CO] .
  * Garcia-Saenz  _et al._ (2023) S. Garcia-Saenz, L. Pinol, S. Renaux-Petel, and D. Werth, JCAP 03, 057 (2023), arXiv:2207.14267 [astro-ph.CO] .
  * Unal (2019) C. Unal, Phys. Rev. D 99, 041301 (2019), arXiv:1811.09151 [astro-ph.CO] .
  * Adshead  _et al._ (2021) P. Adshead, K. D. Lozanov, and Z. J. Weiner, JCAP 10, 080 (2021), arXiv:2105.01659 [astro-ph.CO] .
  * Abe  _et al._ (2023) K. T. Abe, R. Inui, Y. Tada, and S. Yokoyama, JCAP 05, 044 (2023), arXiv:2209.13891 [astro-ph.CO] .
  * Perna  _et al._ (2024) G. Perna, C. Testini, A. Ricciardone, and S. Matarrese, (2024), arXiv:2403.06962 [astro-ph.CO] .
  * Cai  _et al._ (2019a) R.-g. Cai, S. Pi, and M. Sasaki, Phys. Rev. Lett. 122, 201101 (2019a), arXiv:1810.11000 [astro-ph.CO] .
  * Cai  _et al._ (2019b) R.-G. Cai, S. Pi, S.-J. Wang, and X.-Y. Yang, JCAP 10, 059 (2019b), arXiv:1907.06372 [astro-ph.CO] .
  * Agazie  _et al._ (2023a) G. Agazie _et al._ (NANOGrav), Astrophys. J. Lett. 951, L8 (2023a), arXiv:2306.16213 [astro-ph.HE] .
  * Agazie  _et al._ (2023b) G. Agazie _et al._ (NANOGrav), Astrophys. J. Lett. 951, L9 (2023b), arXiv:2306.16217 [astro-ph.HE] .
  * Antoniadis  _et al._ (2023a) J. Antoniadis _et al._ (EPTA), (2023a), arXiv:2306.16214 [astro-ph.HE] .
  * Antoniadis  _et al._ (2023b) J. Antoniadis _et al._ (EPTA),  (2023b), 10.1051/0004-6361/202346841, arXiv:2306.16224 [astro-ph.HE] .
  * Antoniadis  _et al._ (2023c) J. Antoniadis _et al._ (EPTA), (2023c), arXiv:2306.16227 [astro-ph.CO] .
  * Reardon  _et al._ (2023a) D. J. Reardon _et al._ , Astrophys. J. Lett. 951, L6 (2023a), arXiv:2306.16215 [astro-ph.HE] .
  * Zic  _et al._ (2023) A. Zic _et al._ , (2023), arXiv:2306.16230 [astro-ph.HE] .
  * Reardon  _et al._ (2023b) D. J. Reardon _et al._ , Astrophys. J. Lett. 951, L7 (2023b), arXiv:2306.16229 [astro-ph.HE] .
  * Xu  _et al._ (2023) H. Xu _et al._ , Res. Astron. Astrophys. 23, 075024 (2023), arXiv:2306.16216 [astro-ph.HE] .
  * Ellis  _et al._ (2024) J. Ellis, M. Fairbairn, G. Franciolini, G. Hütsi, A. Iovino, M. Lewicki, M. Raidal, J. Urrutia, V. Vaskonen, and H. Veermäe, Phys. Rev. D 109, 023522 (2024), arXiv:2308.08546 [astro-ph.CO] .
  * Figueroa  _et al._ (2024) D. G. Figueroa, M. Pieroni, A. Ricciardone, and P. Simakachorn, Phys. Rev. Lett. 132, 171002 (2024), arXiv:2307.02399 [astro-ph.CO] .
  * Chen  _et al._ (2020) Z.-C. Chen, C. Yuan, and Q.-G. Huang, Phys. Rev. Lett. 124, 251101 (2020), arXiv:1910.12239 [astro-ph.CO] .
  * De Luca  _et al._ (2021) V. De Luca, G. Franciolini, and A. Riotto, Phys. Rev. Lett. 126, 041303 (2021), arXiv:2009.08268 [astro-ph.CO] .
  * Vaskonen and Veermäe (2021) V. Vaskonen and H. Veermäe, Phys. Rev. Lett. 126, 051303 (2021), arXiv:2009.07832 [astro-ph.CO] .
  * Kohri and Terada (2021) K. Kohri and T. Terada, Phys. Lett. B 813, 136040 (2021), arXiv:2009.11853 [astro-ph.CO] .
  * Domènech and Pi (2022) G. Domènech and S. Pi, Sci. China Phys. Mech. Astron. 65, 230411 (2022), arXiv:2010.03976 [astro-ph.CO] .
  * Dandoy  _et al._ (2023) V. Dandoy, V. Domcke, and F. Rompineve, SciPost Phys. Core 6, 060 (2023), arXiv:2302.07901 [astro-ph.CO] .
  * Franciolini  _et al._ (2023a) G. Franciolini, A. Iovino, Junior., V. Vaskonen, and H. Veermae, Phys. Rev. Lett. 131, 201401 (2023a), arXiv:2306.17149 [astro-ph.CO] .
  * Duechting (2004) N. Duechting, Phys. Rev. D 70, 064015 (2004), arXiv:astro-ph/0406260 .
  * Bernal  _et al._ (2018) J. L. Bernal, A. Raccanelli, L. Verde, and J. Silk, JCAP 05, 017 (2018), [Erratum: JCAP 01, E01 (2020)], arXiv:1712.01311 [astro-ph.CO] .
  * Serpico  _et al._ (2020) P. D. Serpico, V. Poulin, D. Inman, and K. Kohri, Phys. Rev. Res. 2, 023204 (2020), arXiv:2002.10771 [astro-ph.CO] .
  * Fixsen  _et al._ (1996) D. J. Fixsen, E. S. Cheng, J. M. Gales, J. C. Mather, R. A. Shafer, and E. L. Wright, Astrophys. J. 473, 576 (1996), arXiv:astro-ph/9605054 .
  * Chluba  _et al._ (2012a) J. Chluba, R. Khatri, and R. A. Sunyaev, Mon. Not. Roy. Astron. Soc. 425, 1129 (2012a), arXiv:1202.0057 [astro-ph.CO] .
  * Chluba  _et al._ (2012b) J. Chluba, A. L. Erickcek, and I. Ben-Dayan, Astrophys. J. 758, 76 (2012b), arXiv:1203.2681 [astro-ph.CO] .
  * Chluba and Grin (2013) J. Chluba and D. Grin, Mon. Not. Roy. Astron. Soc. 434, 1619 (2013), arXiv:1304.4596 [astro-ph.CO] .
  * Bianchini and Fabbian (2022) F. Bianchini and G. Fabbian, Phys. Rev. D 106, 063527 (2022), arXiv:2206.02762 [astro-ph.CO] .
  * Byrnes  _et al._ (2019) C. T. Byrnes, P. S. Cole, and S. P. Patil, JCAP 06, 028 (2019), arXiv:1811.11158 [astro-ph.CO] .
  * Karam  _et al._ (2023a) A. Karam, N. Koivunen, E. Tomberg, V. Vaskonen, and H. Veermäe, JCAP 03, 013 (2023a), arXiv:2205.13540 [astro-ph.CO] .
  * Harada  _et al._ (2015) T. Harada, C.-M. Yoo, T. Nakama, and Y. Koga, Phys. Rev. D 91, 084057 (2015), arXiv:1503.03934 [gr-qc] .
  * Polnarev and Musco (2007) A. G. Polnarev and I. Musco, Class. Quant. Grav. 24, 1405 (2007), arXiv:gr-qc/0605122 .
  * De Luca  _et al._ (2019) V. De Luca, G. Franciolini, A. Kehagias, M. Peloso, A. Riotto, and C. Ünal, JCAP 07, 048 (2019), arXiv:1904.00970 [astro-ph.CO] .
  * Young  _et al._ (2019) S. Young, I. Musco, and C. T. Byrnes, JCAP 11, 012 (2019), arXiv:1904.00984 [astro-ph.CO] .
  * Germani and Sheth (2020) C. Germani and R. K. Sheth, Phys. Rev. D 101, 063520 (2020), arXiv:1912.07072 [astro-ph.CO] .
  * Ferrante  _et al._ (2023a) G. Ferrante, G. Franciolini, A. Iovino, Junior., and A. Urbano, Phys. Rev. D 107, 043520 (2023a), arXiv:2211.01728 [astro-ph.CO] .
  * Akrami  _et al._ (2020b) Y. Akrami _et al._ (Planck), Astron. Astrophys. 641, A9 (2020b), arXiv:1905.05697 [astro-ph.CO] .
  * Garcia-Bellido and Ruiz Morales (2017) J. Garcia-Bellido and E. Ruiz Morales, Phys. Dark Univ. 18, 47 (2017), arXiv:1702.03901 [astro-ph.CO] .
  * Pi  _et al._ (2018) S. Pi, Y.-l. Zhang, Q.-G. Huang, and M. Sasaki, JCAP 05, 042 (2018), arXiv:1712.09896 [astro-ph.CO] .
  * Kannike  _et al._ (2017) K. Kannike, L. Marzola, M. Raidal, and H. Veermäe, JCAP 09, 020 (2017), arXiv:1705.06225 [astro-ph.CO] .
  * Ballesteros  _et al._ (2020) G. Ballesteros, J. Rey, M. Taoso, and A. Urbano, JCAP 07, 025 (2020), arXiv:2001.08220 [astro-ph.CO] .
  * Inomata  _et al._ (2017) K. Inomata, M. Kawasaki, K. Mukaida, Y. Tada, and T. T. Yanagida, Phys. Rev. D 95, 123510 (2017), arXiv:1611.06130 [astro-ph.CO] .
  * Iacconi  _et al._ (2022) L. Iacconi, H. Assadullahi, M. Fasiello, and D. Wands, JCAP 06, 007 (2022), arXiv:2112.05092 [astro-ph.CO] .
  * Kawai and Kim (2021) S. Kawai and J. Kim, Phys. Rev. D 104, 083545 (2021), arXiv:2108.01340 [astro-ph.CO] .
  * Bhaumik and Jain (2020) N. Bhaumik and R. K. Jain, JCAP 01, 037 (2020), arXiv:1907.04125 [astro-ph.CO] .
  * Cheong  _et al._ (2021) D. Y. Cheong, S. M. Lee, and S. C. Park, JCAP 01, 032 (2021), arXiv:1912.12032 [hep-ph] .
  * Inomata  _et al._ (2018) K. Inomata, M. Kawasaki, K. Mukaida, and T. T. Yanagida, Phys. Rev. D 97, 043514 (2018), arXiv:1711.06129 [astro-ph.CO] .
  * Dalianis  _et al._ (2019) I. Dalianis, A. Kehagias, and G. Tringas, JCAP 01, 037 (2019), arXiv:1805.09483 [astro-ph.CO] .
  * Motohashi  _et al._ (2020) H. Motohashi, S. Mukohyama, and M. Oliosi, JCAP 03, 002 (2020), arXiv:1910.13235 [gr-qc] .
  * Hertzberg and Yamada (2018) M. P. Hertzberg and M. Yamada, Phys. Rev. D 97, 083509 (2018), arXiv:1712.09750 [astro-ph.CO] .
  * Ballesteros and Taoso (2018) G. Ballesteros and M. Taoso, Phys. Rev. D 97, 023501 (2018), arXiv:1709.05565 [hep-ph] .
  * Rasanen and Tomberg (2019) S. Rasanen and E. Tomberg, JCAP 01, 038 (2019), arXiv:1810.12608 [astro-ph.CO] .
  * Balaji  _et al._ (2022) S. Balaji, J. Silk, and Y.-P. Wu, JCAP 06, 008 (2022), arXiv:2202.00700 [astro-ph.CO] .
  * Frolovsky and Ketov (2023) D. Frolovsky and S. V. Ketov, Universe 9, 294 (2023), arXiv:2304.12558 [astro-ph.CO] .
  * Dimopoulos (2017) K. Dimopoulos, Phys. Lett. B 775, 262 (2017), arXiv:1707.05644 [hep-ph] .
  * Germani and Prokopec (2017) C. Germani and T. Prokopec, Phys. Dark Univ. 18, 6 (2017), arXiv:1706.04226 [astro-ph.CO] .
  * Choudhury and Mazumdar (2014) S. Choudhury and A. Mazumdar, Phys. Lett. B 733, 270 (2014), arXiv:1307.5119 [astro-ph.CO] .
  * Ragavendra and Sriramkumar (2023) H. V. Ragavendra and L. Sriramkumar, Galaxies 11, 34 (2023), arXiv:2301.08887 [astro-ph.CO] .
  * Cheng  _et al._ (2022) S.-L. Cheng, D.-S. Lee, and K.-W. Ng, Phys. Lett. B 827, 136956 (2022), arXiv:2106.09275 [astro-ph.CO] .
  * Franciolini  _et al._ (2023b) G. Franciolini, A. Iovino, Junior., M. Taoso, and A. Urbano, (2023b), arXiv:2305.03491 [astro-ph.CO] .
  * Karam  _et al._ (2023b) A. Karam, N. Koivunen, E. Tomberg, A. Racioppi, and H. Veermäe, JCAP 09, 002 (2023b), arXiv:2305.09630 [astro-ph.CO] .
  * Mishra  _et al._ (2023) S. S. Mishra, E. J. Copeland, and A. M. Green, JCAP 09, 005 (2023), arXiv:2303.17375 [astro-ph.CO] .
  * Cole  _et al._ (2023) P. S. Cole, A. D. Gow, C. T. Byrnes, and S. P. Patil, JCAP 08, 031 (2023), arXiv:2304.01997 [astro-ph.CO] .
  * Frosina and Urbano (2023) L. Frosina and A. Urbano, Phys. Rev. D 108, 103544 (2023), arXiv:2308.06915 [astro-ph.CO] .
  * Franciolini and Urbano (2022) G. Franciolini and A. Urbano, Phys. Rev. D 106, 123519 (2022), arXiv:2207.10056 [astro-ph.CO] .
  * Choudhury  _et al._ (2024) S. Choudhury, A. Karde, S. Panda, and M. Sami, (2024), arXiv:2401.10925 [astro-ph.CO] .
  * Wang  _et al._ (2024) X. Wang, Y.-l. Zhang, and M. Sasaki, (2024), arXiv:2404.02492 [astro-ph.CO] .
  * Stamou (2021) I. D. Stamou, Phys. Rev. D 103, 083512 (2021), arXiv:2104.08654 [hep-ph] .
  * Stamou (2024) I. Stamou, (2024), arXiv:2404.14321 [astro-ph.CO] .
  * Heydari and Karami (2022a) S. Heydari and K. Karami, Eur. Phys. J. C 82, 83 (2022a), arXiv:2107.10550 [gr-qc] .
  * Heydari and Karami (2022b) S. Heydari and K. Karami, JCAP 03, 033 (2022b), arXiv:2111.00494 [gr-qc] .
  * Heydari and Karami (2024) S. Heydari and K. Karami, JCAP 02, 047 (2024), arXiv:2309.01239 [astro-ph.CO] .
  * Enqvist and Sloth (2002) K. Enqvist and M. S. Sloth, Nucl. Phys. B 626, 395 (2002), arXiv:hep-ph/0109214 .
  * Lyth and Wands (2002) D. H. Lyth and D. Wands, Phys. Lett. B 524, 5 (2002), arXiv:hep-ph/0110002 .
  * Sloth (2003) M. S. Sloth, Nucl. Phys. B 656, 239 (2003), arXiv:hep-ph/0208241 .
  * Lyth  _et al._ (2003) D. H. Lyth, C. Ungarelli, and D. Wands, Phys. Rev. D 67, 023503 (2003), arXiv:astro-ph/0208055 .
  * Dimopoulos  _et al._ (2003) K. Dimopoulos, G. Lazarides, D. Lyth, and R. Ruiz de Austri, JHEP 05, 057 (2003), arXiv:hep-ph/0303154 .
  * Kohri  _et al._ (2013) K. Kohri, C.-M. Lin, and T. Matsuda, Phys. Rev. D 87, 103527 (2013), arXiv:1211.2371 [hep-ph] .
  * Kawasaki  _et al._ (2013a) M. Kawasaki, N. Kitajima, and T. T. Yanagida, Phys. Rev. D 87, 063519 (2013a), arXiv:1207.2550 [hep-ph] .
  * Kawasaki  _et al._ (2013b) M. Kawasaki, N. Kitajima, and S. Yokoyama, JCAP 08, 042 (2013b), arXiv:1305.4464 [astro-ph.CO] .
  * Carr  _et al._ (2017) B. Carr, T. Tenkanen, and V. Vaskonen, Phys. Rev. D 96, 063507 (2017), arXiv:1706.03746 [astro-ph.CO] .
  * Ando  _et al._ (2018a) K. Ando, K. Inomata, M. Kawasaki, K. Mukaida, and T. T. Yanagida, Phys. Rev. D 97, 123512 (2018a), arXiv:1711.08956 [astro-ph.CO] .
  * Ando  _et al._ (2018b) K. Ando, M. Kawasaki, and H. Nakatsuka, Phys. Rev. D 98, 083508 (2018b), arXiv:1805.07757 [astro-ph.CO] .
  * Chen and Cai (2019) C. Chen and Y.-F. Cai, JCAP 10, 068 (2019), arXiv:1908.03942 [astro-ph.CO] .
  * Liu and Prokopec (2021) L.-H. Liu and T. Prokopec, JCAP 06, 033 (2021), arXiv:2005.11069 [astro-ph.CO] .
  * Pi and Sasaki (2023) S. Pi and M. Sasaki, Phys. Rev. D 108, L101301 (2023), arXiv:2112.12680 [astro-ph.CO] .
  * Cai  _et al._ (2021) R.-G. Cai, C. Chen, and C. Fu, Phys. Rev. D 104, 083537 (2021), arXiv:2108.03422 [astro-ph.CO] .
  * Liu (2023) L.-H. Liu, Chin. Phys. C 47, 015105 (2023), arXiv:2107.07310 [astro-ph.CO] .
  * Chen  _et al._ (2023) C. Chen, A. Ghoshal, Z. Lalak, Y. Luo, and A. Naskar, JCAP 08, 041 (2023), arXiv:2305.12325 [astro-ph.CO] .
  * Torrado  _et al._ (2018) J. Torrado, C. T. Byrnes, R. J. Hardwick, V. Vennin, and D. Wands, Phys. Rev. D 98, 063525 (2018), arXiv:1712.05364 [astro-ph.CO] .
  * Wilkins and Cable (2024) A. Wilkins and A. Cable, JCAP 02, 026 (2024), arXiv:2306.09232 [astro-ph.CO] .
  * Inomata  _et al._ (2024) K. Inomata, M. Kawasaki, K. Mukaida, and T. T. Yanagida, Phys. Rev. D 109, 043508 (2024), arXiv:2309.11398 [astro-ph.CO] .
  * Musco  _et al._ (2021) I. Musco, V. De Luca, G. Franciolini, and A. Riotto, Phys. Rev. D 103, 063538 (2021), arXiv:2011.03014 [astro-ph.CO] .
  * Musco  _et al._ (2024) I. Musco, K. Jedamzik, and S. Young, Phys. Rev. D 109, 083506 (2024), arXiv:2303.07980 [astro-ph.CO] .
  * Andrés-Carcasona  _et al._ (2024) M. Andrés-Carcasona, A. J. Iovino, V. Vaskonen, H. Veermäe, M. Martínez, O. Pujolàs, and L. M. Mir, (2024), arXiv:2405.05732 [astro-ph.CO] .
  * Gow  _et al._ (2023) A. D. Gow, H. Assadullahi, J. H. P. Jackson, K. Koyama, V. Vennin, and D. Wands, EPL 142, 49001 (2023), arXiv:2211.08348 [astro-ph.CO] .
  * Ianniccari  _et al._ (2024) A. Ianniccari, A. J. Iovino, A. Kehagias, D. Perrone, and A. Riotto, (2024), arXiv:2402.11033 [astro-ph.CO] .
  * Franciolini  _et al._ (2023c) G. Franciolini, A. Ianniccari, A. Kehagias, D. Perrone, and A. Riotto, (2023c), arXiv:2311.03239 [astro-ph.CO] .
  * De Luca  _et al._ (2023) V. De Luca, A. Kehagias, and A. Riotto, Phys. Rev. D 108, 063531 (2023), arXiv:2307.13633 [astro-ph.CO] .
  * Young (2022) S. Young, JCAP 05, 037 (2022), arXiv:2201.13345 [astro-ph.CO] .
  * Green  _et al._ (2004) A. M. Green, A. R. Liddle, K. A. Malik, and M. Sasaki, Phys. Rev. D 70, 041502 (2004), arXiv:astro-ph/0403181 .
  * Young  _et al._ (2014) S. Young, C. T. Byrnes, and M. Sasaki, JCAP 07, 045 (2014), arXiv:1405.7023 [gr-qc] .
  * Bardeen  _et al._ (1986) J. M. Bardeen, J. R. Bond, N. Kaiser, and A. S. Szalay, Astrophys. J. 304, 15 (1986).
  * Balaji  _et al._ (2023) S. Balaji, G. Domènech, and G. Franciolini, JCAP 10, 041 (2023), arXiv:2307.08552 [gr-qc] .
  * Yoo  _et al._ (2018) C.-M. Yoo, T. Harada, J. Garriga, and K. Kohri, PTEP 2018, 123E01 (2018), arXiv:1805.03946 [astro-ph.CO] .
  * Yoo  _et al._ (2019) C.-M. Yoo, J.-O. Gong, and S. Yokoyama, JCAP 09, 033 (2019), arXiv:1906.06790 [astro-ph.CO] .
  * Gow  _et al._ (2021) A. D. Gow, C. T. Byrnes, P. S. Cole, and S. Young, JCAP 02, 002 (2021), arXiv:2008.03289 [astro-ph.CO] .
  * Franciolini  _et al._ (2022) G. Franciolini, I. Musco, P. Pani, and A. Urbano, Phys. Rev. D 106, 123526 (2022), arXiv:2209.05959 [astro-ph.CO] .
  * Hajkarim and Schaffner-Bielich (2020) F. Hajkarim and J. Schaffner-Bielich, Phys. Rev. D 101, 043522 (2020), arXiv:1910.12357 [hep-ph] .
  * Abe  _et al._ (2021) K. T. Abe, Y. Tada, and I. Ueda, JCAP 06, 048 (2021), arXiv:2010.06193 [astro-ph.CO] .
  * Franciolini  _et al._ (2024) G. Franciolini, D. Racco, and F. Rompineve, Phys. Rev. Lett. 132, 081001 (2024), arXiv:2306.17136 [astro-ph.CO] .
  * Ferreira and Joyce (1998) P. G. Ferreira and M. Joyce, Phys. Rev. D 58, 023503 (1998), arXiv:astro-ph/9711102 .
  * Pallis (2006) C. Pallis, Nucl. Phys. B 751, 129 (2006), arXiv:hep-ph/0510234 .
  * Redmond  _et al._ (2018) K. Redmond, A. Trezza, and A. L. Erickcek, Phys. Rev. D 98, 063504 (2018), arXiv:1807.01327 [astro-ph.CO] .
  * Co  _et al._ (2022) R. T. Co, D. Dunsky, N. Fernandez, A. Ghalsasi, L. J. Hall, K. Harigaya, and J. Shelton, JHEP 09, 116 (2022), arXiv:2108.09299 [hep-ph] .
  * Gouttenoire  _et al._ (2021) Y. Gouttenoire, G. Servant, and P. Simakachorn, (2021), arXiv:2111.01150 [hep-ph] .
  * Chang and Cui (2022) C.-F. Chang and Y. Cui, JHEP 03, 114 (2022), arXiv:2106.09746 [hep-ph] .
  * Domènech (2020) G. Domènech, Int. J. Mod. Phys. D 29, 2050028 (2020), arXiv:1912.05583 [gr-qc] .
  * Borsanyi  _et al._ (2016) S. Borsanyi _et al._ , Nature 539, 69 (2016), arXiv:1606.07494 [hep-lat] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Kohri and Terada (2018) K. Kohri and T. Terada, Phys. Rev. D 97, 123532 (2018), arXiv:1804.08577 [gr-qc] .
  * Atal and Domènech (2021) V. Atal and G. Domènech, JCAP 06, 001 (2021), [Erratum: JCAP 10, E01 (2023)], arXiv:2103.01056 [astro-ph.CO] .
  * Yuan and Huang (2021) C. Yuan and Q.-G. Huang, Phys. Lett. B 821, 136606 (2021), arXiv:2007.10686 [astro-ph.CO] .
  * Domènech  _et al._ (2022) G. Domènech, S. Passaglia, and S. Renaux-Petel, JCAP 03, 023 (2022), arXiv:2112.10163 [astro-ph.CO] .
  * Liu  _et al._ (2024) L. Liu, Z.-C. Chen, and Q.-G. Huang, Phys. Rev. D 109, L061301 (2024), arXiv:2307.01102 [astro-ph.CO] .
  * Yuan  _et al._ (2023) C. Yuan, D.-S. Meng, and Q.-G. Huang, JCAP 12, 036 (2023), arXiv:2308.07155 [astro-ph.CO] .
  * Li  _et al._ (2023) J.-P. Li, S. Wang, Z.-C. Zhao, and K. Kohri,  (2023), 10.1088/1475-7516/2024/06/039, arXiv:2309.07792 [astro-ph.CO] .
  * Collaboration (2023) T. N. Collaboration, “KDE Representations of the Gravitational Wave Background Free Spectra Present in the NANOGrav 15-Year Dataset,”  (2023).
  * Hütsi  _et al._ (2021) G. Hütsi, M. Raidal, V. Vaskonen, and H. Veermäe, JCAP 03, 068 (2021), arXiv:2012.02786 [astro-ph.CO] .
  * Romero-Rodriguez  _et al._ (2022) A. Romero-Rodriguez, M. Martinez, O. Pujolàs, M. Sakellariadou, and V. Vaskonen, Phys. Rev. Lett. 128, 051301 (2022), arXiv:2107.11660 [gr-qc] .
  * Kehagias  _et al._ (2019) A. Kehagias, I. Musco, and A. Riotto, JCAP 12, 029 (2019), arXiv:1906.07135 [astro-ph.CO] .
  * Escrivà  _et al._ (2022) A. Escrivà, Y. Tada, S. Yokoyama, and C.-M. Yoo, JCAP 05, 012 (2022), arXiv:2202.01028 [astro-ph.CO] .
  * Escrivà  _et al._ (2020) A. Escrivà, C. Germani, and R. K. Sheth, Phys. Rev. D 101, 044022 (2020), arXiv:1907.13311 [gr-qc] .
  * Kehagias  _et al._ (2024) A. Kehagias, D. Perrone, and A. Riotto, (2024), arXiv:2405.05208 [astro-ph.CO] .
  * Atal and Germani (2019) V. Atal and C. Germani, Phys. Dark Univ. 24, 100275 (2019), arXiv:1811.07857 [astro-ph.CO] .
  * Biagetti  _et al._ (2018) M. Biagetti, G. Franciolini, A. Kehagias, and A. Riotto, JCAP 07, 032 (2018), arXiv:1804.07124 [astro-ph.CO] .
  * Atal  _et al._ (2019) V. Atal, J. Garriga, and A. Marcos-Caballero, JCAP 09, 073 (2019), arXiv:1905.13202 [astro-ph.CO] .
  * Tomberg (2023) E. Tomberg, Phys. Rev. D 108, 043502 (2023), arXiv:2304.10903 [astro-ph.CO] .
  * Tisserand  _et al._ (2007) P. Tisserand _et al._ (EROS-2), Astron. Astrophys. 469, 387 (2007), arXiv:astro-ph/0607207 .
  * Mroz  _et al._ (2024a) P. Mroz _et al._ , (2024a), arXiv:2403.02386 [astro-ph.GA] .
  * Mroz  _et al._ (2024b) P. Mroz _et al._ , (2024b), arXiv:2403.02398 [astro-ph.GA] .
  * Koushiappas and Loeb (2017) S. M. Koushiappas and A. Loeb, Phys. Rev. Lett. 119, 041102 (2017), arXiv:1704.01668 [astro-ph.GA] .
  * Agius  _et al._ (2024) D. Agius, R. Essig, D. Gaggero, F. Scarcella, G. Suczewski, and M. Valli, (2024), arXiv:2403.18895 [hep-ph] .
  * Facchinetti  _et al._ (2023) G. Facchinetti, M. Lucca, and S. Clesse, Phys. Rev. D 107, 043537 (2023), arXiv:2212.07969 [astro-ph.CO] .
  * Brandt (2016) T. D. Brandt, Astrophys. J. Lett. 824, L31 (2016), arXiv:1605.03665 [astro-ph.GA] .
  * Monroy-Rodríguez and Allen (2014) M. A. Monroy-Rodríguez and C. Allen, Astrophys. J. 790, 159 (2014), arXiv:1406.5169 [astro-ph.GA] .
  * Murgia  _et al._ (2019) R. Murgia, G. Scelfo, M. Viel, and A. Raccanelli, Phys. Rev. Lett. 123, 071102 (2019), arXiv:1903.10509 [astro-ph.CO] .
  * Zumalacarregui and Seljak (2018) M. Zumalacarregui and U. Seljak, Phys. Rev. Lett. 121, 141101 (2018), arXiv:1712.02240 [astro-ph.CO] .
  * Bird  _et al._ (2011) S. Bird, H. V. Peiris, M. Viel, and L. Verde,  413, 1717 (2011), arXiv:1010.1519 [astro-ph.CO] .
  * Sharma  _et al._ (2024) D. Sharma, J. Lesgourgues, and C. T. Byrnes, (2024), arXiv:2404.18474 [astro-ph.CO] .
  * Kormendy and Richstone (1995) J. Kormendy and D. Richstone, Ann. Rev. Astron. Astrophys. 33, 581 (1995).
  * Magorrian  _et al._ (1998) J. Magorrian _et al._ , Astron. J. 115, 2285 (1998), arXiv:astro-ph/9708072 .
  * Richstone  _et al._ (1998) D. Richstone _et al._ , Nature 395, A14 (1998), arXiv:astro-ph/9810378 .
  * Kohri  _et al._ (2014) K. Kohri, T. Nakama, and T. Suyama, Phys. Rev. D 90, 083514 (2014), arXiv:1405.5999 [astro-ph.CO] .
  * Reines and Volonteri (2015) A. E. Reines and M. Volonteri, The Astrophysical Journal 813, 82 (2015).
  * Fukugita  _et al._ (1998) M. Fukugita, C. J. Hogan, and P. J. E. Peebles, Astrophys. J. 503, 518 (1998), arXiv:astro-ph/9712020 .
  * Ünal  _et al._ (2021) C. Ünal, E. D. Kovetz, and S. P. Patil, Phys. Rev. D 103, 063519 (2021), arXiv:2008.11184 [astro-ph.CO] .
  * Nakama  _et al._ (2018) T. Nakama, B. Carr, and J. Silk, Phys. Rev. D 97, 043525 (2018), arXiv:1710.06945 [astro-ph.CO] .
  * Hooper  _et al._ (2024) D. Hooper, A. Ireland, G. Krnjaic, and A. Stebbins, JCAP 04, 021 (2024), arXiv:2308.00756 [astro-ph.CO] .
  * Byrnes  _et al._ (2024) C. T. Byrnes, J. Lesgourgues, and D. Sharma, (2024), arXiv:2404.18475 [astro-ph.CO] .
  * Ferrante  _et al._ (2023b) G. Ferrante, G. Franciolini, A. Iovino, Junior., and A. Urbano, JCAP 06, 057 (2023b), arXiv:2305.13382 [astro-ph.CO] .