# Constraining the impact of standard model phase transitions on primordial
black holes

Xavier Pritchard and  0009-0007-6543-8563 Christian T. Byrnes
0000-0003-2583-6536

(October 2023)

###### Abstract

Phase transitions in the early universe lead to a reduction in the equation of
state of the primordial plasma. This exponentially enhances the formation rate
of primordial black holes. However, this sensitivity to the equation of state
is the same that primordial black hole abundances show to the primordial
curvature power spectrum amplitude. In this paper, we investigate peaked power
spectra and show the challenges associated with motivating populations of
primordial black holes with standard model enhancements. The parametrisation
of different power spectra plays an important role in this discussion. The
allowed parameter space consistent with a large QCD phase transition impact on
the primordial black hole abundance differs greatly. This is particularly
evident for wider spectra.

##  1 Introduction

The idea that compact objects in the early universe could have observational
consequences in the present day was first postulated over fifty years ago [1].
Shortly after this work, it was suggested that these objects could be black
holes, formed from the collapse of density perturbations [2, 3]. And so was
born the field of primordial black holes (PBHs).

Relatively early on in the field’s lifetime, it was known that PBHs were a
possible dark matter candidate [4] (see [5, 6] for reviews). More recently,
there have been several advancements in gravitational wave (GW) science, such
as the Laser Interferometer Gravitational-Wave Observatory (LIGO), Virgo and
the Kamioka Gravitational Wave Detector (KAGRA) [7, 8, 9]. These advancements
have led to new perspectives on PBHs, such as accounting for the detected GWs
from BH mergers [10, 11, 12]. Amongst several other factors, this has led to a
huge surge in the literature (see [13, 14, 15] for recent reviews of the
field).

PBH studies typically require an understanding of primordial curvature
perturbations. This information is encoded in the primordial curvature power
spectrum $\mathcal{P}_{\zeta}(k)$, assuming the perturbations follow a
Gaussian distribution (for non-Gaussian PBH studies see [16, 17, 18], for
example). In the context of PBHs, $\mathcal{P}$ζ needs to be enhanced by many
orders of magnitude from the amplitude measured on the large scales ($10^{-4}$
Mpc-1 $\lesssim$ $k$ $\lesssim$ 1 Mpc-1) of the cosmic microwave background
(CMB), where the power spectrum has been measured with extreme accuracy [19].

A popular model used to achieve such an enhancement, whilst matching CMB
observations, is ultra-slow-roll inflation [20, 21]. The level of fine-tuning
needed to achieve such a boost of power on small scales is significant in
several theoretically motivated single-field inflation models [22] (see also
[23, 24] for more promising multi-field studies). Regardless, we assume there
is a peak and study how different parametrisations (and thus shapes) of the
peak impact PBH abundances and constraints.

PBH studies not only require an understanding of curvature perturbations in
the early universe but also of the universe itself when these perturbations
re-enter the horizon. Typically, PBHs are thought to form when the universe is
radiation dominated (see [25, 26, 27] for examples of PBHs forming during an
early matter dominated epoch and [28] for a review). The pressure of said
radiation is dictated by the standard model (SM) of particle physics. As the
universe expands and cools, particles become nonrelativistic, and the pressure
of the radiation drops more rapidly than the energy density.

It has been understood for a long time that the reduction in pressure of the
primordial plasma leads to an exponential enhancement in the formation rate of
PBHs [29, 30]. Making use of this SM-related enhancement and assuming a flat
$\mathcal{P}_{\zeta}(k)$, it was argued that several cosmological questions
such as LIGO-Virgo-KAGRA observations and the origin of supermassive black
holes could be explained in a unified PBH scenario [31]. This scenario has
since been shown to be unviable [32], or at the very least it is finely-tuned.
One of the main goals of this paper is to extend this investigation to peaked
power spectrum. Thus, we will consider whether the SM-related enhancement in
abundance can motivate answers to the originally considered cosmic conundra
for more realistic power spectra.

A difficulty often faced when motivating PBHs at a given mass scale is the
sensitivity of PBH abundance to the amplitude of $\mathcal{P}$ζ, which is also
exponential. Hence, we expect PBH formation to take place only on scales very
close to where the primordial power spectrum peaks, regardless of whether this
coincides with a scale at which a SM particle decouples. If these scales don’t
coincide but you would like a population of PBHs motivated by the SM, you must
then have a sufficiently broad power spectrum. Thus, we investigate the
challenges associated with having such a broad peak.

The key challenge to generating PBHs on larger scales than the horizon scale
during the QCD transition are CMB spectral distortions, whose constraints on
$\mathcal{P}$ζ are about three orders of magnitude tighter than observational
PBH constraints [33]. Thus, power spectra peaked at these scales with broad
tails to small $k$ are ruled out. On smaller scales, pulsar timing array (PTA)
constraints from the possible induced stochastic GWs are comparably tight to
direct PBH constraints. This places a less tight constraint on the tails of
the primordial power spectrum towards smaller scales.

We find that power spectrum constraints are largely dependant on the
parametrisation you assume. For some parametrisations of $\mathcal{P}$ζ,
motivating a population of PBHs at a given SM-related scale through the
associated reduction in pressure is tightly constrained. This requires the
power spectrum to be peaked at this scale, which we argue is fine-tuned. For
other shapes of the power spectrum we find that the standard model is able to
have a large impact on primordial black hole abundances.

The paper is organised as follows. In Section 2 we review early universe phase
transitions, tracking both the pressure and energy density of the primordial
plasma and thus the equation of state and threshold for collapse. In Section 3
we describe the power spectrum parametrisations that we consider, before
explaining how we calculate our PBH abundances. In Section 4 we summarise the
current constraints on the primordial power spectrum which we then use to
constrain the parameters of our chosen power spectra. Throughout this section,
we pay special attention to the mass scales relevant to the SM. We conclude
our findings in Section 5. Finally, we discuss additional details and
extensions of our work in the appendix.

##  2 Relativistic degrees of freedom in the early universe

Reheating leaves the universe filled with a hot, dense plasma of relativistic
SM particles [34, 35] (see [36, 37, 38] for reviews). As the universe expands
the density and pressure of the plasma decrease at the same rate, along with
the temperature. As the temperature decreases below the mass threshold of a
given SM particle the particle decouples from the plasma, resulting in a drop
in relativistic degrees of freedom. This results in the ratio between pressure
$p$ and the total energy density $\rho$, or equation of state, becoming non-
constant and varying with temperature as

$\omega(T)\equiv\frac{p(T)}{\rho(T)}.$ (2.1)

To obtain our early universe equation of state we use the following formulas
for $\rho_{i}$ and $p_{i}$ for a given particle, $i$, in thermal equilibrium

$\displaystyle\rho_{i}$
$\displaystyle=\frac{g_{i}}{2\pi^{2}}\int_{m_{i}}^{\infty}\frac{\sqrt{E^{2}-m_{i}^{2}}}{\text{exp}(E/T)\pm
1}E^{2}\text{d}E,$ (2.2) $\displaystyle
p_{i}=\frac{g_{i}}{6\pi^{2}}\int_{m_{i}}^{\infty}\frac{(E^{2}-m_{i}^{2})^{3/2}}{\text{exp}(E/T)\pm
1}\text{d}E,$ (2.3)

where $m_{i}$ is the particle’s mass, $g_{i}$ is the particle’s degrees of
freedom, and the plus (minus) sign in the denominator applies to fermions
(bosons). This information can be found in [39], for instance. To obtain the
total pressure and density and hence $\omega(T)$ we sum over all particles.

To relate temperature to horizon mass, wavenumber and time we use the
following approximate relations [40, 41]

$\begin{split}M_{H}&\simeq 1.5\times
10^{5}\left(\frac{g_{\ast}}{10.75}\right)^{-1/2}\left(\frac{T}{\text{1MeV}}\right)^{-2}M_{\odot},\\\
&\simeq
17\left(\frac{g_{\ast}}{10.75}\right)^{-1/6}\left(\frac{k}{10^{6}\text{Mpc}^{-1}}\right)^{-2}M_{\odot},\\\
&\simeq 2\times 10^{5}\left(\frac{t}{1\text{s}}\right)M_{\odot},\end{split}$
(2.4)

where $g_{\ast}$ is the number of relativistic degrees of freedom as a
function of $T$ and $k$ respectively. We note the cancellation of $g_{\ast}$
contributions in the final equation.

###  2.1 Standard model phase transitions

The first drop in the equation of state in relation to SM particles occurs at
around 100 GeV (corresponding to $M_{H}\approx 10^{-5}$ $M_{\odot}$) when the
electroweak (EW) gauge symmetry SU(2)$\times$U(1) is spontaneously broken by
the vacuum expectation value of the Higgs field. The observed mass of the
Higgs boson [42] suggests that this symmetry breaking is a smooth crossover
rather than a first-order phase transition [43], thus this period is named the
electroweak crossover. During this crossover the top quark, Z, W, and Higgs
bosons each become non-relativistic. In Fig. 1 we track the contributions to
$\omega$ during the EW PT as an example.

Figure 1: Total reduction in the equation of state $\omega$ during the EW
crossover, as a function of the horizon mass of the universe. We also track
each relevant particle’s contribution to this reduction.

In minimal extensions of the standard model, the EW PT can be strongly first
order, for example in Next-to-Minimal supersymmetry models [44, 45, 46],
composite Higgs models [47, 48] and many more. This would result in violent
bubble collisions and the scalar order parameter and fluid of light particles
could produce GWs [49] within the sensitivity region of future GW detectors.

The most drastic drop in $\omega(T)$ occurs during the QCD PT when the
temperature of the universe is approximately 100 MeV(corresponding to
$M_{H}\approx 1$ $M_{\odot}$). This occurs when strong interactions confine
quarks into hadrons. Lattice QCD studies suggest that similar to the EW PT,
the QCD PT is also a crossover [50, 51]. Lattice studies are needed due to QCD
being highly non-perturbative. To provide a complete equation of state, we
thus interpolate the results of [52].

Shortly after the QCD PT pions and muons decouple from the plasma. Finally, at
around 1 MeV (corresponding to $M_{H}\approx 10^{6}$ $M_{\odot}$),
$e^{-}e^{+}$ annihilation and neutrino decoupling occur. It has recently been
argued that, in the context of PBH formation, the reduction in $\omega(T)$
during the $e^{-}e^{+}$ annihilation epoch doesn’t necessarily enhance
abundances [53, 38]. The argument is that, since neutrino interactions with
the rest of the plasma freeze out shortly before this PT, they are able to
free-stream out of overdense regions.

###  2.2 Threshold for collapse

The critical threshold for collapse, or $\delta_{c}$ as it is often written,
is a key quantity when calculating PBH abundances. This quantity defines the
minimum value a density perturbation amplitude must have at horizon entry for
a BH horizon to form and collapse. Generally, a black hole collapses if the
inwards gravitational force overcomes the outwards force of pressure
gradients. In our context, the density (pressure) of the primordial plasma
when a perturbation re-enters the horizon will set the gravitational
(pressure) forces. Thus, it is naturally a function of the equation of state
$\omega$. We note that the collapse threshold value also depends on the shape
of the energy density profile [54, 55, 56].

Using the equation of state, we calculate the threshold for collapse following
the method of [57, 58]. To do this, we track the equation of state at the
horizon entry time of the perturbation and the time at which the overdensity
begins to collapse. These times correspond to the minimum and maximum peak
formation times, which can be related to minimum and maximum peak masses using
Eq. 2.4. This can be intuitively understood as the horizon mass growing with
time.

It is also useful to track time-averaged and logarithmic time-averaged values
from the equation of state. These values are both similar, suggesting an
insensitivity to our method of averaging. It should also be stated that the
averaged values are likely to be closer to the true answer, due to the
equation of state varying between horizon entry and overdensity collapse. We
plot $\delta_{c}(M_{H})$ in Fig. 2, obtained using the time-averaging method.

Figure 2: The drop in the threshold for the collapse of PBHs, as a function of
the horizon mass of the universe. Here we take the radiation domination value
of the threshold as 0.453, which has been found numerically [57]. Clearly,
there are three departures from this value, corresponding to the three PTs
described in the previous subsection.

When calculating many quantities involved with PBHs there are typically many
layers of detail, both theoretically and numerically. For example, in the
above derivation of $\delta_{c}$ from $\omega$. As noted in [11], for a more
complete calculation we should include the effect of both the sound speed
squared and the gravitational field potential on the dynamics of PBH collapse.
In this example, we find that although these inclusions shift each drop in
$\delta_{c}$ to slightly higher masses, the shape remains similar. Throughout
the paper, we have tried to be careful in selecting methods that lead to
results that are insensitive to increases in complexity, which we exemplify in
Appendix C.

##  3 PBH mass distribution

The abundance of PBHs is usually stated in terms of the mass fraction of the
universe that collapses to form PBHs at the time of formation or the fraction
of CDM comprised of PBHs. These, in themselves, are interesting quantities to
study, with a plethora of techniques used in calculations throughout the
literature. Moreover, knowing the PBH abundance is pivotal to discussions
about GWs from black hole mergers, constraints on PBHs, and gravitational wave
backgrounds to name a few.

###  3.1 The curvature power spectrum

As previously mentioned, there are several challenges with having an extended,
nearly flat primordial power spectrum if we wish to motivate PBH abundance
enhancement due to the drops in $\omega(T)$ [32]. We aim to test these
conclusions for different shapes of power spectra, of which we focus on two
commonly used in the literature. The two we focus on are different
parametrisations of peaks with tuneable positions and widths, which offer a
more realistic situation than a flat peak with sharp cutoffs; generating a 6–7
orders of magnitude increase in power at a single scale is easier for a model
than creating the same increase in power over a large range of scales. The
first peaked parametrisation we consider is the lognormal power spectrum,
defined as

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

where $\Delta$ is the width of the power spectrum, $k_{p}$ is the peak scale
of the power spectrum and the normalisation is chosen such that the variance
with respect to $\zeta$ is equal to $A$. We recover the Dirac delta function
for $\Delta\rightarrow 0$.

Secondly, we focus on the broken power law power spectrum, defined as

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

where $\alpha$ describes the growth, $\tau$ describes the decay and $\lambda$
the width of the spectrum at the scale $k_{p}$. This parametrisation arises
naturally from USR models [16]. In both these examples, $A$ is the power
spectrum amplitude.

In the literature, it is common to see BPL shape parameters described with
different notation. We have chosen our notation so as not to confuse later
discussion, and have been consistent throughout the paper.

###  3.2 Variance of perturbations

Once we have chosen a curvature power spectrum we can define the smoothed
density power spectrum during radiation domination as

$\mathcal{P}_{\delta_{R}}(k)=\frac{16}{81}(kR)^{4}W^{2}(k,R)\mathcal{P}_{\zeta}(k),$
(3.3)

where $W(k,R)$ is a window function used to relate $\mathcal{P}_{\zeta}(k)$ in
Fourier space to the probability distribution function of $\delta$ in real
space [59, 60]. The 16/81 factor depends on the equation of state which is
non-constant during a given PT, making this equation approximate. We also
choose to omit a transfer function $T(k,R)$ because it’s order one at horizon
entry. It has, however, recently been argued that a non-linear transfer
function is needed to be consistent in threshold and variance calculations
[61]. Moreover, the renormalisation procedure [62] must also be applied on
super-horizon scales.

As discussed in [63], we must use the same smoothing function when calculating
the PBH formation criterion. In this work, we mainly considered a Gaussian
window function modified by a factor of 2 in the exponent [63] defined as

$W(k,R)=\text{exp}\left(-\frac{(kR)^{2}}{4}\right).$ (3.4)

We also consider the Fourier transform of the top-hat smoothing function,
given by

$\tilde{W}(k,R)=3\frac{\text{sin}(kR)-kR\,\text{cos}(kR)}{(kR)^{3}}.$ (3.5)

We can then use $\mathcal{P}_{\delta_{R}}(k)$ to calculate the momenta of the
distribution, given as

$\sigma^{2}_{i}(R)=\int_{0}^{\infty}\frac{\text{d}k}{k}k^{2i}\mathcal{P}_{\delta_{R}},$
(3.6)

where $i$ is a non-negative integer. We note here an important point regarding
both of the considered power spectra, the lognormal parametrisation has been
constructed such that the variance with respect to $\zeta$ is equal to the
power spectrum amplitude $A$. The broken power law parametrisation, on the
other hand, does not share this property. We find that, upon increasing
respective width parameters, the variance of perturbations with respect to the
density contrast (Eq. 3.6, with $i$=0) increases (decreases) for LN (BPL).
This will be an important distinction when discussing PBH constraints later in
this paper.

###  3.3 Press–Schechter formalism

The mass fraction of the universe that collapses to form PBHs at the time of
formation, $\beta$, is a possible way of describing PBH abundances. This can
loosely be thought of as the effect the curvature power spectrum and the
equation of state in the early universe have on each other. Typically, this is
calculated using one of two formalisms: Press–Schechter theory [64] or peaks
theory [65]. The focus of this subsection and the vast majority of the paper
is the former. We note that the power spectrum amplitude required to generate
a given PBH abundance is insensitive to the calculation technique [66].

In the context of PBHs, Press-Schechter theory says that density perturbations
with a peak value greater than the threshold value, $\delta_{c}$, will
collapse to form a black hole. The abundance is thus given as

$\beta=2\int_{\delta_{c}}^{\infty}\text{d}\delta P(\delta),$ (3.7)

where the factor of 2 is the Press-Schechter factor. If we assume that
perturbations follow a Gaussian distribution then the Gaussian probability
density, $P$($\delta$), is given as

$P(\delta)=\frac{1}{\sqrt{2\pi\sigma_{0}^{2}}}\text{exp}\left(-\frac{\delta^{2}}{2\sigma_{0}^{2}}\right).$
(3.8)

If we further assume that PBHs form with exactly the horizon mass we can then
approximate the abundance as

$\beta=\text{erfc}\left(\frac{\delta_{c}}{\sqrt{2\sigma_{0}^{2}}}\right),$
(3.9)

where erfc is the complementary error function. Finally, we can compute the
mass distribution of PBHs stated in terms of the fraction of CDM made up of
PBHs of a given mass $M$, given as

$f(M_{\text{PBH}})=\frac{1}{\Omega_{\text{CDM}}}\frac{\text{d}\Omega_{\text{PBH}}}{\text{dln}M},$
(3.10)

where $\Omega_{\text{PBH}}$ is the total abundance of PBHs relative to the
critical density and $\Omega_{\text{CDM}}=0.245$ [67]. This equation is
simplified if $\beta$ is constant when $\delta_{c}$ is constant, which is the
case if we assume that the density power spectrum at horizon-entry is scale-
invariant over a relevant but limited range of scales. If we further assume
that black holes form with mass $M$=$M_{H}$, we arrive at

$f(M_{\text{PBH}})=\frac{2.4}{f_{\text{PBH}}}\beta(M_{\text{PBH}})\left(\frac{M_{\text{PBH}}}{M_{\text{eq}}}\right)^{-1/2},$
(3.11)

where the factor 2.4 comes from 2(1+$\Omega_{\text{b}}/\Omega_{\text{CDM}})$,
with $\Omega_{\text{b}}=0.0456$ [67] and $M_{\text{eq}}$ is the horizon mass
at the time of matter-radiation equality, $M_{\text{eq}}\approx 2.8\times
10^{17}$ $M_{\odot}$.

In a more realistic calculation, we should take into consideration the effect
critical collapse has on the PBH mass spectrum. The previous calculation can
thus be extended to account for PBHs not forming with exactly the horizon
mass. It was found that the PBH mass follows a scaling relation [57],
described by the critical collapse equation

$M_{\text{PBH}}=\kappa M_{H}(\delta-\delta_{c})^{\gamma},$ (3.12)

where $M_{H}$ is the horizon mass at horizon re-entry, $\kappa\simeq 0.36$,
$\gamma\simeq 3.3$ and $\delta_{c}\simeq 0.453$ are found numerically assuming
radiation domination.

We arrive at the amended expression

$\beta=2\int_{\delta_{c}}^{\infty}\text{d}\delta\frac{M_{\text{PBH}}}{M_{H}}P(\delta)=2\int_{\delta_{c}}^{\infty}\text{d}\delta\kappa(\delta-\delta_{c})^{\gamma}P(\delta).$
(3.13)

Defining $\mu\equiv\frac{M_{\text{PBH}}}{\kappa M_{H}}$ and
$\frac{\text{d}\delta}{\text{d}\text{ln}M_{\text{PBH}}}=\frac{1}{\gamma}\mu^{1/\gamma}$,
we obtain the PBH mass distribution today

${}\begin{split}f(M_{\text{PBH}})=&\frac{1}{\Omega_{\text{CDM}}}\int_{-\infty}^{\infty}\text{d}\text{ln}M_{H}\frac{2}{\sqrt{2\pi\sigma_{0}^{2}(M_{H})}}\text{exp}\left[-\frac{(\mu^{1/\gamma}+\delta_{c}(M_{H}))^{2}}{2\sigma_{0}^{2}(M_{H})}\right]\\\
&\qquad\qquad\qquad\qquad\qquad\qquad\times\frac{M_{\text{PBH}}}{\gamma
M_{H}}\mu^{1/\gamma}\sqrt{\frac{M_{\text{eq}}}{M_{H}}}.\end{split}$ (3.14)

###  3.4 Phase transition impacts on PBH abundances

We would like to quantify a given PT’s impact on PBH abundances. To do this we
first calculate the scale corresponding to the peak of $f(M_{\text{PBH}})$
assuming radiation domination (which is not exactly the same scale as
$k_{p}$), for every point in a given power spectrum’s parameter space. We
label this scale $k_{f_{\rm peak}}$. We then consider the full equation of
state and calculate the ratio of $f(M_{\text{PBH}})$ at a given PT scale
versus $f(M_{\text{PBH}})$ at $k_{f_{\rm peak}}$. We choose to normalise our
power spectrum amplitude $A$ such that $f_{\text{PBH}}=1$, which gives us a
better understanding of the effect of varying each parameter. This value of
$A$ also depends on power spectrum parameters. In reality, this value of
$f_{\text{PBH}}$ is constrained at the QCD scale to be $\mathcal{O}$(10-3)
[68]. Notably, the ratio $f_{\rm PT}/f_{\rm peak}$ for lower values of
$f_{\text{PBH}}$ is similar, for a given choice of parameters.

If this ratio is unity, the value of $f(M_{\text{PBH}})$ is the same at
$k_{f_{\rm peak}}$ and $k_{PT}$. Of course, this will occur if the scales are
the same (i.e. $k_{f_{\rm peak}}=k_{PT}$). We argue that in this case, despite
a given PT increasing $f(M_{\text{PBH}})$, there would have to be some level
of fine-tuning in a given model to realise the scenario. Of interest are the
parameter values that lead to the ratio being greater than one, signifying a
significant change in the shape of our abundance curve which will now have its
highest value at the same scale as the PT. In these scenarios, we can say that
a given PT has impacted the abundance of PBHs.

Initially, we focus on a lognormal $\mathcal{P}_{\zeta}$ peaked at scales
around the QCD crossover. The biggest drop in relativistic degrees of freedom
occurs on these scales so we expect the biggest change in the shape of
$f(M_{\text{PBH}})$ to be around this QCD scale, upon making $\delta_{c}$ a
function of $M_{H}$. In Fig. 3, we exemplify our process of calculating
$f_{\rm PT}/f_{\rm peak}$, for lognormal parameters $k_{p}=2\times 10^{6}$ and
$\Delta=2.5$, recovering a value of $\approx 9.14$. This value is relatively
large, so we would say that the QCD PT has had a large impact on the PBH
abundance curve, which can be seen from this figure.

Interestingly, we see in the plot that the peak in the full
$f(M_{\text{PBH}})$ does not quite occur at the QCD scale. This property
becomes more evident the further $k_{p}$ gets from $k_{\text{QCD}}$. We
considered calculating the ratio at the peak scale in the full
$f(M_{\text{PBH}})$, instead of at the QCD scale. However, regardless of at
which which scale you evaluate $f(M_{\text{PBH}})$ in the numerator, the ratio
remains relatively similar. We simply point out that, as the difference in
$k_{p}$ and $k_{\text{QCD}}$ increases, our ratio values increasingly
underestimate this new ratio.

Figure 3: We exemplify our method of calculating the impact of the QCD PT on
the PBH abundance. We show the abundance curves assuming radiation domination
and considering the full equation of state, for a lognormal power spectrum
with $k_{p}=2\times 10^{6}$ and $\Delta=2.5$. The vertical gridlines show the
peak scale in the radiation domination abundance curve and the QCD scale. The
horizontal gridlines are the corresponding values of the full
$f(M_{\text{PBH}})$ evaluated at these scales. In this example, we recover a
ratio of $9.14$.

In Fig. 4, we show the value of the ratio, $f_{\rm PT}/f_{\rm peak}$, assuming
a lognormal peak. We can see from the figure that, for a given peak scale in
$\mathcal{P}_{\zeta}$, increasing the width will lead to an increase in the
ratio. This is to be expected, for wider spectra scales further from the peak
will have more power. As previously mentioned, the PBH abundance is very
sensitive to this increase. By widening our spectra we are thus enhancing the
impact of the drop in $\delta_{c}$ on the shape of our abundance curve.

As well as the QCD PT, we investigated the EW PT in this context. These
results are also plotted in Fig. 4 where we can see that none of the [$k_{p}$,
$\Delta$] parameter space leads to a ratio $f_{\text{EW}}$/$f_{\rm peak}$
greater than unity. This doesn’t come as much of a surprise, despite the
exponential sensitivity of abundances to $\delta_{c}$($\omega$) we only see a
drop of $\sim$ 1$\%$ in this parameter (as opposed to $\sim$ 10$\%$ during the
QCD PT). In [31], two planetary-mass microlensing events [69, 70] were argued
to possibly be explained by PBHs. The motivation for this was the softening in
$\omega$ at the EW scale producing a peak in their abundance curve. However,
our results suggest difficulty motivating PBHs at this precise scale, assuming
a lognormal power spectrum.

Figure 4: The ratio of the full $f(M_{\text{PBH}})$ evaluated at PT scales
versus the peak scale in the abundance curve assuming radiation domination. We
show the results assuming a lognormal power spectrum, over the QCD and EW PTs.
This gives us some intuition about the effect of a PT on the shape of
$\mathcal{P}$ζ as its parameters are varied. At the EW scale, we see the
maximum of this ratio equaling unity.

Despite having the ability to choose lognormal parameter values which lead to
a large enhancement on the impact of the QCD PT on PBH abundances, many
parameter values are ruled out by observational constraints. These will be
discussed in Section 4.

Lastly, we investigated the effect of the QCD PT on the PBH abundance curve
assuming a broken power law power spectrum. To make a more direct comparison
with Fig. 4, we focus on the width parameter of the broken power law,
$\lambda$. Our results are plotted in Fig. 5, where we have set shape
parameters $\alpha=4$ and $\tau=3$. Plotting up to $\lambda=40$, we reproduce
a similar plot to the lognormal, with some interesting distinctions. It would
be informative to investigate whether such a large value for $\lambda$ is
natural from a model-building perspective. One could also perform a log-
likelihood analysis on a relevant dataset as was done in [71] (although only
up to $\lambda=10$), to reveal whether this value of $\lambda$ was favorable.
However, these discussions are beyond the scope of this paper.

The figure exemplifies how important the cut-offs of a power spectrum are when
calculating $f_{\text{PT}}$/$f_{\rm peak}$. For $k_{p}\lesssim 6\times 10^{6}$
Mpc-1 (which is the scale at which $f_{\rm peak}$=$f_{\text{QCD}}$), $\tau$ is
the relevant quantity. For $k_{p}$ greater than this scale, $\alpha$ is the
relevant quantity. Decreasing the steepness of the slope by a value of 1, we
see a large increase in the ratio. We point out that this effect isn’t solely
due to the slope. As can also be seen in the lognormal case (which has the
same slope on either side of the peak), the ratio is typically larger for
lower values of $k_{p}$. This is due to pions ($\pi^{\pm,0}$) and muons
($\mu$) decoupling shortly after the QCD PT, leaving an additional drop in
$\delta_{c}$, which can be seen in Fig. 2.

Figure 5: The ratio of the full $f(M_{\text{PBH}})$ evaluated at the QCD scale
versus the peak scale in the abundance curve assuming radiation domination,
for a BPL power spectrum. In this plot, we vary the width parameter,
$\lambda$, whilst keeping $\alpha=4$ and $\tau=3$. We can see several
similarities to Fig. 4.

##  4 PBH constraints

Primordial black holes with initial masses greater than $\sim 10^{15}$ g, in
the standard cosmological picture, wouldn’t have evaporated completely by
today. Interestingly, in higher dimensions this evaporation bound is weakened
[72], leaving a larger range in which PBHs can account for the totality of
dark matter. In [73], the observational constraints on PBHs with initial
masses greater than 1015 g are discussed in detail. These constraints can be
categorised in two ways, direct and indirect observational constraints.
Indirect observational constraints will be discussed in the following
subsections. For direct constraints, we use Fig. 18 from [73] which shows the
combined constraints on the mass fraction $\beta$ with common parameters
normalised out. However, see [74] for hints of detections from a variety of
observations. We use equation (8) of [73] to relate this parameter to $\beta$
as follows

$\beta^{\prime}(M)\equiv\gamma_{\text{eff}}^{1/2}\left(\frac{g_{\ast
i}}{106.75}\right)^{-1/4}\left(\frac{h}{0.67}\right)^{-2}\beta(M),$ (4.1)

where $h$ and $g_{\ast i}$ are known up to 10 percent errors which are
unimportant to us. We set $\gamma_{\text{eff}}$ equal to unity, however, the
precise value is an ongoing debate [75].

We would like to use the constraints on the mass fraction to constrain
$\mathcal{P}_{\zeta}$. To do this we simply plug Eq. 3.9 into Eq. 4.1 and
rearrange for the power spectrum amplitude. This results in the amplitude
constraint plot being a function of power spectrum parameters. We account for
the non-linearities between the density contrast and curvature perturbation
[76] following the arguments of [66]. We thus multiply our constraint line by
a factor of 1.98. In Fig. 6 we compare the full non-linear calculation
(detailed in Appendix A) to the simpler approximation. We see excellent
agreement between the two, justifying the simplification.

Figure 6: The direct observational constraints of PBHs on the primordial power
spectrum. We compare the full calculation including non-linearities to the
approximation made in [66], which we use for the remainder of the paper. We
use a lognormal power spectrum with $\Delta$=1. We see excellent agreement
between this approximation and the more computationally expensive full
calculation.

As previously mentioned, it was found that the power spectrum amplitude
required to generate a given value of $f_{\text{PBH}}$ is insensitive to
whether we calculate abundances via Press–Schechter or peaks theory [66]. We
find it important to state that, if we instead fixed the power spectrum
amplitude and then calculated the abundance, we expect a 1–2 order of
magnitude difference for each technique [77, 78]. This is an interesting
distinction between the two techniques however, for the purposes of this
paper, we are interested in constraints on power spectrum amplitude. Thus, it
is more important to us that the power spectrum amplitude remains similar
regardless of abundance calculation technique. This means that, if we were to
use peaks theory, we would expect similar limits on the power spectrum from
direct PBH constraints.

###  4.1 CMB $\mu$-constraints

Spectral distortions quantify deviations from the CMB’s black-body temperature
distribution. If there were an increase in power on small scales the energy
stored in density perturbations would be dissipated through photon diffusion,
via Silk damping. On scales relevant to $\mu$-distortions (the most relevant
for our analysis), Bremsstrahlung and double Compton scattering are
inefficient due to the expansion of the universe. These inefficiencies result
in photon distributions on different scales being described by different
blackbodies [79].

The spectral distortions place one of the tightest current constraints on the
primordial curvature power spectrum. During the $\mu$-era, due to no detection
of spectral distortions, the power spectrum amplitude must be lower than 10-4,
meaning PBH formation at this scale is completely ruled out in the standard
picture of Gaussian curvature perturbations [33]. Interestingly, this scale
corresponds to PBHs of supermassive black hole masses and $e^{-}$$e^{+}$
annihilation. Adding non-Gaussianity with positive skewness is not enough to
stop $\mu$-constraints from being significantly tighter than the PBH
constraints unless very extreme and completely non-perturbative non-
Gaussianity is considered [80, 81, 82, 83, 84].

To quantify how much additional power is needed to cause observable spectral
distortions we must calculate the final $\mu$-distortions induced by scalar
perturbations, which can be approximated as [79]

$\mu\approx\int_{k_{\text{min}}}^{\infty}\frac{\text{d}k}{k}\mathcal{P}_{\zeta}(k)W_{\mu}(k),$
(4.2)

with $k$-space window functions of the form

$W_{\mu}(k)\approx
2.27\left[\text{exp}\left(-\left[\frac{\hat{k}}{1360}\right]^{2}\middle/\left[1+\left[\frac{\hat{k}}{260}\right]^{0.3}+\frac{\hat{k}}{340}\right]\right)-\text{exp}\left(-\left[\frac{\hat{k}}{32}\right]^{2}\right)\right],$
(4.3)

where $\hat{k}\equiv k$/1 Mpc-1 and $k_{\text{min}}$ $\simeq$ 1 Mpc-1. We can
compare the total induced $\mu$-distortions with observations, resulting in a
power spectrum-dependant function of scale $k$. The strongest current
observational constraint on $\mu$-distortions comes from the FIRAS instrument
on the COBE satellite. A recent reanalysis of FIRAS data suggests
$\mu$-distortions to be smaller than 4.7$\times$10-5 [85], which is 2 times
tighter than the original analysis [86].

###  4.2 Pulsar timing arrays

On smaller scales than $\mu$-distortion constraints, the primordial curvature
power spectrum is constrained by PTAs. An increase in power at a given scale
will result in large first-order scalar perturbations at that scale. In
cosmological perturbation theory at second order, scalar and tensor
perturbations are coupled meaning this increase in power can be responsible
for a stochastic gravitational wave background (see [87] for a review).

Recently, PTA collaborations NANOGrav [88, 89], EPTA [90, 91, 92], PPTA [93,
94, 95] and CPTA [96] have released data showing evidence for a stochastic
gravitational wave background. There is already much literature investigating
PTA results and $\mathcal{P}_{\zeta}$ [97, 98, 99, 100], as well as PBH
studies [101, 102, 103, 104, 105, 106, 107, 108]. In addition to scalar-
induced GWs there are many other cosmological explanations for this evidence
[109], such as phase transitions [110, 111, 112, 113, 114, 115, 116] and
topological defects [117, 118, 119, 120, 121, 122]. There are also
astrophysical explanations for the signal [123, 124, 125, 126, 127], such as
merging supermassive black holes. We assume that the stochastic gravitational
wave background is due to something other than scalar-induced GWs, allowing us
to place tighter constraints on the curvature power spectrum.

To calculate the scalar-induced GW spectrum from a given curvature power
spectrum during radiation domination, $\Omega_{\text{GW,RD}}$, we used [128].
We relate this to the GW spectrum today using the following formula [129]

$\Omega_{\text{GWB}}(f)h^{2}\approx 1.62\times 10^{-5}\Omega_{\text{GW,RD}}.$
(4.4)

Then, we find the limiting amplitude of the given power spectrum such that PTA
data may constrain the scalar-induced GW spectrum. We use the all-sky GW
strain 95$\%$ upper limit data from NANOGrav15 [130], which we convert to a
limit on primordial GWB energy density through the following equation

$\Omega_{\text{GWB}}(f)h^{2}=\frac{2\pi^{2}}{3H_{0}^{2}}f^{2}h_{c}^{2}(f),$
(4.5)

where $h\approx$ 0.67 [67] is the dimensionless Hubble parameter, $H_{0}$=100
$h$ km s-1 Mpc-1, and $h_{c}$ is the GW strain. We convert this to units of
wavelength via the following equation

$f=1.6\text{nHz}\left(\frac{k}{10^{6}\text{Mpc}^{-1}}\right).$ (4.6)

Finally, we arrive at our formula for constraining the power spectrum
amplitude using PTA data

$A_{\text{constraint}}=\text{Min}\left(\sqrt{\frac{\Omega_{\text{GWB}}(k)h^{2}}{\Omega_{\text{GWB,signal}}(k,k_{p},\Delta)h^{2}}}\right).$
(4.7)

###  4.3 Constraining $\mathcal{P}_{\zeta}$

####  4.3.1 Lognormal peak

As has been discussed, both direct and indirect observations can constrain the
primordial curvature power spectrum. In Fig. 7 we show an example plot for the
constraints on the lognormal power spectrum amplitude with $\Delta=1.1$. A
point above a given constraint line is ruled out by the respective constraint.
We also show the values of the power spectrum amplitude leading to
$f_{\text{PBH}}=1$ and $f_{\text{PBH}}=10^{-10}$. The small value of $f_{\rm
PBH}=10^{-10}$ is the order of magnitude constraint we would reach if most of
the dark matter consisted of WIMPs [131, 132, 133, 134, 135, 136, 137].

Figure 7: We show the indirect constraints on an example lognormal power
spectrum, with $\Delta$=1.1. We also show the amplitude required to produce a
PBH abundance $f_{\text{PBH}}$=1 and $f_{\text{PBH}}$=10-10 and the current
(direct) PBH constraint which lies in between these extreme values.

Using Fig. 7, we can identify the scales at which spectral distortion and PTA
observations place tighter constraints on the power spectrum amplitude than
direct PBH constraints. For the given values of $f_{\rm PBH}$ we can calculate
the values of $k_{p}$ that are ruled out by indirect constraints. We point out
that each constraint line and each $f_{\rm PBH}$ line will differ upon
variation of the width of the lognormal, $\Delta$. This can be seen from each
respective constraint equation, which all depend on $\mathcal{P}_{\zeta}$.
Using all this information we can make a parameter space constraint plot,
which we show in Fig. 8.

Figure 8: We show the constraints on the lognormal power spectrum from
$\mu$-distortions (blue line) and PTA (red line) for different scenarios. The
parameter points above these lines cannot produce PBH populations with $f_{\rm
PBH}$ values above those shown in the caption.

Notably, $\mu$-distortion constraints on power spectrum parameters remain
similar regardless of each of the scenarios considered. This shows an
insensitivity to $f_{\text{PBH}}$ that comes from the relatively small
variation in direct observational constraint lines, which is due to the
exponential sensitivity of $f_{\text{PBH}}$ to the power spectrum amplitude.
In other words, a relatively small change in $A$ will result in a relatively
big change in $f_{\text{PBH}}$. With regard to the PTA parameter space
constraints, we see more dependence on our PBH constraint line. This is
because PTA and PBH constraints are comparably tight, resulting in a
sensitivity to relatively small changes in both PTA and PBH constraint lines.
We also note the coincidence of scales between PTA observations and the QCD
scale ($\approx 3\times 10^{6}$ Mpc-1), which is evident from the plot.

We can see that a wider lognormal power spectrum leads to tighter
$\mu$-distortion and PTA constraints. In particular, if the power spectrum
peaks at scales $k_{p}\lesssim 10^{7}$ Mpc-1, then to evade all constraints
the width of the power spectrum is required to be $\Delta\lesssim 1$. In this
case, it is difficult to realise the QCD PT largely impacting the shape of
$f_{\text{PBH}}$ (as can be seen in Fig. 4). In [31], PBHs in this mass range
were argued to offer solutions to several ongoing issues [138, 139, 140].
Moreover, LIGO-Virgo-KAGRA detections of BHs of potentially primordial origin
have also been studied in detail, including the QCD impact [11, 12]. We have
shown that, if a lognormal $\mathcal{P}_{\zeta}$ doesn’t peak around the QCD
scale, then it must be broad for the PT to impact the overall abundance. We
have also shown that this scenario is tightly constrained. We reiterate that
this conclusion is drawn assuming the PTA signal is not due to scalar-induced
GWs and hence we are not trying to fit PTA data.

####  4.3.2 Broken power law peak

The broken power law power spectrum could be said to have several advantages
over the lognormal if we want to motivate PBHs from the drops in $\delta_{c}$.
Firstly, instead of falling logarithmically about the peak it has sharper cut-
offs. This results in scales around the peak having less of an effect on a
given constraint. Secondly, with two extra parameters used to describe the
BPL, we will have more freedom to control the overall power spectrum shape.

To produce the parameter space constraint plots for the BPL equivalent to Fig.
8, we had to fix two of ($\alpha$, $\tau$, $\lambda$). For example, in the
[$k_{p}$, $\alpha$] plot we fixed both $\tau$ and $\lambda$. The fixed values
we chose were based on [101], who found $\alpha=4$, $\tau=3$, and $\lambda=1$
to be a relatively conservative parameter combination. In other words, a small
change in each respective fixed parameter won’t affect our conclusions. The
results are shown in Fig. 9.

Figure 9: The constraints on BPL parameters ($Top$ $left$: $\alpha$, $Top$
$right$: $\tau$, $Bottom$: $\lambda$) and $k_{p}$ coming from
$\mu$-distortions and PTA, where shaded regions are ruled out. We have also
included $f_{\text{PBH}}=1$, 10-4 constraints. When varying each shape
parameter respectively, we set the remaining two according to $\alpha=4$,
$\tau=3$ and $\lambda=1$.

We can see from Fig. 9, that for shape parameters $\tau$ and $\gamma$
(describing the decay and width of the BPL power spectrum) the variation in
constraining power of $\mu$-distortion observations is not large. We note that
this is particularly evident for the $\tau$ parameter, which is to be expected
as this describes the tail of the distribution. For the growth parameter,
$\alpha$, we see that this is not the case. Decreasing $\alpha$ can be thought
of as flattening out the power spectrum on scales $k<k_{p}$. Thus, it is also
expected that decreasing $\alpha$ will lead to tighter $\mu$-distortion
constraints on power spectra peaked at smaller scales. This can be seen from
the plot, which suggests a difficulty motivating BPL spectra with
$\alpha\lesssim 1$.

Similar to the lognormal constraint plot, $\mu$-distortion constraints appear
to be insensitive to the value of $f_{\text{PBH}}$. On the other hand, we find
that if we take $f_{\text{PBH}}=10^{-10}$ (as was done in the lognormal case)
then PTA null-detections rule out none of the parameter space. We thus compare
constraints with the power spectrum amplitude leading to
$f_{\text{PBH}}=10^{-4}$. Using these constraints we consider the BPL
parameter space permitting multiple PT imprints on PBH abundance in Appendix
B. To make a more direct comparison between the lognormal we focus on the
$\lambda$ (width) parameter of the BPL for the remainder of this subsection.

As mentioned in Section 3.2, increasing the width parameter for both
considered power spectra affects the variance of perturbations with respect to
the density contrast differently. Direct PBH constraints depend on this
quantity (as can be seen from rearranging Eq. 4.1), making this an important
distinction between the two spectra. Notably, increasing the width of the BPL
(LN) results in more (less) tight direct constraints of on power spectrum
amplitude. For the case of a BPL power spectrum with $\lambda\gtrsim$ 6, we
find that PTA constraints are no longer comparable to direct observational
constraints. In other words, direct observations more tightly constrain the
power spectrum amplitude than PTA constraints. From Fig. 9, we can see the
other parameter space constraints don’t share this property, with PTA
constraints remaining relatively similar upon increasing $\alpha$ and $\tau$.

Looking at Fig. 5, we see that a large portion of interesting (as defined in
Section 3.4) parameter space is unconstrained by $\mu$-distortion and PTA
constraints. Thus, the shape of the power spectrum plays an important role if
we wish for the QCD PT to affect the overall shape of PBH abundance curves. We
have also shown that the constraining power of both direct and indirect
observations can differ greatly based on the assumed peaked power spectrum
shape. Upon variation of power spectrum parameters, the behaviour of the
variance, $\sigma_{0}^{2}$, plays an important role in this reasoning. We have
shown that these differences in behaviour become more apparent when increasing
the width of our respective power spectrum.

Recently, it was found that the constraints on the total abundance of PBHs in
the solar mass range were insensitive to the shape of the power spectrum [84].
Interestingly, we find our calculation of PTA constraints on the power
spectrum amplitude to be slightly less constraining than theirs, whilst our
PBH constraints are slightly tighter. Moreover, their constraints on the power
spectrum amplitude obtained through Press–Schecter and peaks theory differ
more than we would have expected. We have found that the shape of
$\mathcal{P}_{\zeta}$ can have a large impact on PBH abundances. Specifically,
we have shown that it is much more difficult to motivate PBHs due to
reductions in the equation of state if we assume a lognormal
$\mathcal{P}_{\zeta}$ over a BPL. In contrast to [84], we have focused on
wider spectra in our investigation where we find the contrast of constraints
on different parameterisations to be most evident.

##  5 Conclusions

In this paper, we have calculated constraints on two popular parametrisations
of peaked power spectra. Combining direct observational non-detections of
primordial black holes, recent pulsar timing array data, and $\mu$-distortions
we have put constraints on the parameters of each power spectrum. We have also
compared this to $\mathcal{P}_{\zeta}$ constraints assuming different values
of $f_{\text{PBH}}$. We find the lognormal power spectrum parameter space to
be more tightly constrained than that of the broken power law, offering less
flexibility in describing populations of primordial black holes. We believe
the additional freedom of the broken power law comes from having extra
parameters and sharper cut-offs about the peak.

We have paid particular interest to standard model phase transitions
throughout our analysis, investigating the impact that the corresponding drop
in $\omega$ (and thus threshold for collapse) has on $f(M_{\text{PBH}})$. This
was in part due to each phase transition scale corresponding to an interesting
cosmological scale such that if a population of primordial black holes were to
form around said scale, it could offer a solution to a cosmological unknown.
In the context of the lognormal power spectrum, we find it difficult to
motivate any of these scenarios.

For the broken power law we find that, upon increasing the width of the power
spectrum, the variance of perturbations with respect to the density contrast
decreases. This, in turn, leads to tighter direct PBH constraints on the power
spectrum amplitude. We find that for width parameter $\lambda\gtrsim 6$ the
pulsar timing array constraints are no longer competitive to the direct PBH
constraints on the power spectrum. This results in a large region of parameter
space in which the QCD phase transition has a large impact on primordial black
hole abundances, whilst evading pulsar timing array constraints. Thus,
different parameterisations of peaked power spectra can lead to vastly
different conclusions as to whether the shape of our abundance curves can be
impacted by the standard model-related drops in the equation of state.

Nonetheless, we found that for peaked power spectra it was difficult
motivating primordial black holes with the drop in equation of state
associated to the electroweak phase transition. We have also discussed the
difficulties related to power spectra peaked at scales around $e^{-}e^{+}$
annihilation. Moreover, capturing the impact of more than one phase transition
on primordial black hole abundances is not permitted for any of the parameter
space values we considered.

Lastly, we show in the appendices that our results are insensitive to several
extensions we could make to our calculations. Namely, non-linearities,
collapse dynamic effects, and constraints including the full equation of
state. This gives us more confidence in our conclusions above.

Acknowledgements

The authors thank S. Young and A. Gow for commenting on a draft of this paper
and K. Saikawa, S. Shirai, A. Aljazaeri, I. Aldecoa Tamayo, R. Munoz, D.
Gillies, R. Cereskaite, and G. Assant for helpful conversations. XP is
supported by an STFC studentship. CB is supported by STFC grants ST/X001040/1
and ST/X000796/1.

##  Appendix A Non-linearities

It has been known for some time that the density contrast and curvature
perturbation have a non-linear relationship [76]. In the super-horizon regime,
this can be expressed as [141]

$\delta(\hat{r})=-\frac{4(1+\omega)}{5+3\omega}\left(\frac{1}{aH}\right)^{2}\text{exp}\left(-\frac{5\zeta(r)}{2}\right)\nabla^{2}\text{exp}\left(\frac{\zeta(r)}{2}\right),$
(A.1)

where $\hat{r}$ is the radial coordinate in the comoving synchronous gauge,
and in spherical symmetry

$\nabla^{2}=\frac{\partial^{2}}{\partial
r^{2}}+\frac{2}{r}\frac{\partial}{\partial r}.$ (A.2)

The result of this relationship is that we expect the density contrast to be
non-Gaussian, regardless of whether $\zeta$ follows a Gaussian distribution.
However, when calculating PBH abundances we are more interested in the
smoothed density contrast. This can be written in terms of a linear, Gaussian
component $\delta_{l}$ as follows

$\delta_{m}=\delta_{l}-\frac{3}{8}\delta_{l}^{2}.$ (A.3)

The critical amplitude of this component can be calculated by inverting the
above, giving

$\delta_{c,l\pm}=\frac{4}{3}\left(1\pm\sqrt{\frac{2-3\delta_{c}}{2}}\right).$
(A.4)

We can further note that the number density of sufficiently rare peaks in a
comoving volume for a random Gaussian field is [65]

$\mathcal{N}=\left(\frac{\sigma_{1}^{3}}{4\pi^{2}\sigma_{0}^{3}}\right)\nu^{3}\text{exp}\left(-\frac{\nu^{2}}{2}\right),$
(A.5)

where $\nu$$\equiv\delta_{l}/\sigma_{0}$. Following [142], peaks theory is
then used to calculate $\beta$ such that the non-linear relationship between
$\delta$ and $\zeta$ is accounted for, given as

$\beta(M_{H})=\int^{\frac{4}{3\sigma_{0}}}_{\nu_{c-}}\text{d}\nu\frac{\kappa}{3\pi}(\nu\sigma_{0}-\frac{3}{8}(\nu\sigma_{0})^{2}-\delta_{c})^{\gamma}\left(\frac{\sigma_{1}}{aH\sigma_{0}}\right)^{3}\nu^{3}\text{exp}\left(-\frac{\nu^{2}}{2}\right),$
(A.6)

where $\nu_{c-}\equiv\delta_{c,l-}$/$\sigma_{0}$.

To account for non-linearities when calculating our PBH constraints, we follow
the arguments of [66]. They assume that each peak that collapses to form a PBH
has an amplitude close to the critical value. Using the Fourier transform of
the top-hat window function, which provides an analytic relation between
$\delta$ and $\zeta$, $\delta_{c}=0.55$ [143] which leads to the critical
amplitude of the Gaussian component $\delta_{c,l-}\approx 0.77$. Thus, they
conclude that peaks should have an amplitude 1.41 times larger if we are to
approximate non-linearities. The overall factor we times our PBH constraint
line is then 1.412=1.98, because $\mathcal{P}_{\zeta}\propto\delta_{c}^{2}$.

##  Appendix B A broad BPL example

As mentioned, the $e^{-}e^{+}$ annihilation scale coincides with tight
$\mu$-distortion constraints. Thus, when asking whether we can have a peaked
power spectrum that captures features of two early universe phase transitions
we are talking about the QCD and EW crossovers. The BPL offers a more
interesting answer to this question than the lognormal, which requires an
unreasonable width. As can be seen from Fig. 9 there is a large portion of
unconstrained [$k_{p}$, $\tau$] parameter space with small $\tau$. We can thus
realise a situation in which two PTs significantly impact our abundance
calculations: $k^{4}$ growth in power [144] (or steeper [145, 146, 147, 148]
to a peak at 105 Mpc-1$\lesssim k_{p}\lesssim k_{QCD}\approx 3\times 10^{6}$
Mpc-1 with $\tau\ll 1$.

In this scenario, we are able to capture the effects of the QCD PT on PBH
abundances, due to the position of the peak. Combined with the steep growth we
are also able to evade $\mu$ constraints. Lastly, the shape parameter $\tau\ll
1$ results in a long tail in our distribution, allowing us to also capture EW
effects on our abundance. This seems interesting until we consider
correlations between $\delta$ and $\zeta$. The correlation coefficient,
discussed in the appendix of [149], is defined as

$\gamma_{\delta\zeta}=\frac{\sigma_{\delta\zeta}^{2}}{\sigma_{\delta}\sigma_{\zeta}},$
(B.1)

with

$\begin{split}\sigma_{\delta\zeta}^{2}&=\frac{4}{9}\int^{\infty}_{0}\frac{\text{d}k}{k}(kR)^{2}\tilde{W}(k,R)\mathcal{P}_{\zeta}(k),\\\
\sigma_{\delta}^{2}&=\frac{16}{81}\int^{\infty}_{0}\frac{\text{d}k}{k}(kR)^{4}\tilde{W}^{2}(k,R)\mathcal{P}_{\zeta}(k),\\\
\sigma_{\zeta}^{2}&=\int^{\infty}_{0}\frac{\text{d}k}{k}\mathcal{P}_{\zeta}(k).\end{split}$
(B.2)

If $\gamma_{\delta\zeta}$ is close to one, then $\delta$ and $\zeta$ are
highly correlated. In other words, a peak in $\delta$ corresponds to a peak in
$\zeta$. This validates the high-peak limit assumption, and thus spherical
symmetry. On the contrary, if the correlation coefficient is close to zero
then perturbation theory based on $\zeta$ doesn’t work. In other words, whilst
small values of $\delta$ could give $f_{\text{PBH}}\approx 1$, this could
result in $\zeta$ greater than unity. This leads to our results being
unreliable, so we would ideally be using power spectrum parameters such that
$\gamma_{\delta\zeta}\approx 1$.

In Fig. 10 we plot the correlation coefficient for the $\tau$ parameter of the
BPL $\mathcal{P}_{\zeta}$. The figure shows that for $\tau\ll$1 we obtain
$\gamma_{\delta\zeta}\ll 1$. Thus, we are unable to draw any reliable
conclusions from our method regarding $\tau\ll 1$. A log-likelihood analysis
of EPTADR2 and NANOGrav15 data, however, has shown that this scenario is
largely disfavoured [71].

Figure 10: The correlation coefficient, $\gamma_{\delta\zeta}$ for the BPL
shape parameter $\tau$. This value quantifies the correlation between $\delta$
and $\zeta$. As $\tau$$\ll 1$, these quantities are highly uncorrelated, which
would lead our results to be unreliable.

We calculate $\gamma_{\delta\zeta}$ for the lognormal and BPL parameter values
leading to the $f_{\text{QCD}}$/$f_{\rm peak}$$\approx 9$, as discussed in
Section 3.4. For the BPL, the width parameter $\lambda=40$ leads to
$\gamma_{\delta\zeta}\approx 0.39$. For the lognormal, the width parameter
$\Delta=2.5$ leads to $\gamma_{\delta\zeta}\approx 0.3$. This shows that our
results leading to this large impact on abundance are more reliable if we
assume the shape of the power spectrum is a BPL instead of a lognormal.

##  Appendix C Extensions of our work

###  C.1 Threshold calculations

As mentioned in the main text, when calculating the threshold of collapse for
PBHs from the equation of state there are many extensions we can make. Our
method, for example, neglects the effects of sound speed and doesn’t compute
the gravitational potential correctly. Calculating the threshold including
these effects was done in [11]. The shape differs slightly from ours, however
the total relative drop is similar. We argue that for our discussions
involving the threshold the more important factor is the relative change and
not the shape. Thus, we are not concerned by this discrepancy and believe our
conclusions remain reliable despite being approximate. We compare these
thresholds in Fig. 11.

Figure 11: We plot three comparisons of the relative change in threshold for
collapse during the QCD PT. We can see that our threshold calculation leads to
a similar drop in threshold to [11]. We also plot results obtained from
fitting functions of energy and entropy densities, which leads to a slightly
reduced drop.

The third alternative available to us is to rewrite our equation of state in
terms of relativistic degrees of energy density, $g_{\ast\rho}$, and entropy
density, $g_{\ast s}$. These are related to energy and entropy via the
following equations

$g_{\ast\rho}(T)\equiv\frac{30\rho(T)}{\pi^{2}T^{4}},\qquad g_{\ast
s}(T)\equiv\frac{45s(T)}{2\pi^{2}T^{3}},$ (C.1)

where $s$ is the entropy density.

We can retrieve the pressure $p$ from the relationship between the energy
density and entropy density via $p$=$sT-\rho$. Then we are able to recast
$\omega(T)$ as

$\omega(T)\equiv\frac{p(T)}{\rho(T)}=\frac{4g_{\ast
s}(T)}{3g_{\ast\rho}(T)}-1.$ (C.2)

Using fitting functions for these densities, provided by [150], we can plot
the corresponding threshold for collapse. The lack of reduction in threshold
is noticeable, which we expect would lead to different results to ones we
found. In particular, to obtain the same impact from a given phase transition
on PBH abundance at a fixed $k_{p}$, we would require a wider power spectrum
which is more tightly constrained.

This is due to the authors specifically constructing smooth functions which
reproduce the features of $g_{\ast s}(T)$ and $g_{\ast\rho}(T)$. This does not
necessarily mean that fine features in some composite functions of these two
quantities can be reproduced. If we wanted to construct composite functions
such as $g_{\ast s}(T)/g_{\ast\rho}(T)$ we would require information on the
correlations between the two functions, else we have very large errors. We
conclude that, in the context of $\omega(T)$ and primordial black hole
abundances, it is better to follow the method outlined in the main text.

###  C.2 $\kappa$($M_{H}$), $\gamma(M_{H})$

Another extension of our work that could be made is regarding Eq. 3.14. As
discussed in [11], in a full calculation the quantities $\kappa$ and $\gamma$
are nonconstant during a PT. These numerical factors depend on the radial
profile of the perturbations being considered as well as the characteristic
scale of the horizon crossing of said perturbations. We plot the impact of the
$M_{H}$ dependence of these quantities on the PBH abundance in Fig. 12,
assuming a nearly scale-invariant curvature power spectrum.

Figure 12: The PBH mass function f($M_{\rm PBH}$) calculated varying the
$M_{H}$ dependence of the quantities $\kappa$ and $\gamma$. The dashed line
shows the resulting abundance with no PTs. The green line follows our typical
calculation.

As can be seen from the plot, the shape of each abundance curve remains
similar. There are, however, some interesting distinctions. In particular,
accounting for variation in $\gamma$ boosts PBH abundance by a factor of 2.
Secondly, varying $\kappa$ seems to shift the peak to slightly lower masses.
Once more, we believe our conclusions would not be altered by the inclusion of
the variation of these parameters in our calculations.

###  C.3 Constraint dependence on $M_{H}$

Lastly, we comment on the impact of a varying equation of state on our
constraints. For example, the scalar-induced GW spectrum we calculated assumes
the constant radiation domination value of $\omega$. As previously discussed,
this is not the case in the early universe. We show the PTA constraints for
different constant values of $\omega$ in Fig. 13. As can be seen, if $\omega$
decreases the constraints on $\mathcal{P}_{\zeta}$ are less severe. We can go
further toward the complete picture by using $\omega$($T$), however, this
makes the computation much more involved. From the figure we expect these
changes to be relatively small. We note that this has been done for a cut
power-law power spectrum [151].

As well as PTA constraints depending on the equation of state, direct PBH
constraints also depend on this quantity. This can be seen from Eq. 3.7, where
it is understood that the threshold depends on $\omega$. Moreover, the overall
factor we multiplied our PBH constraint line in Appendix A in order to
approximate non-linearities depends on both threshold and gravitational
potential. In Fig. 13, we also show the impact of including these additional
contributions to the PBH constraint line. Again, we see the impact of this
inclusion is relatively small.

Figure 13: This plot shows the ($Top$: PTA, $Bottom$: PBH) constraints on the
lognormal power spectrum, varying the equation of state parameter. As we can
see, a reduction in the equation of state results in a weaker constraint on
the power spectrum from PTA data, but tighter constraints from the null-
detection of PBH. In both cases, the change is relatively small.

## References

  * [1] Ya.. Zeldovich and I.. Novikov  “The Hypothesis of Cores Retarded during Expansion and the Hot Cosmological Model”  In _Soviet Astronomy_ 10, 1967, pp. 602
  * [2] Stephen Hawking  “Gravitationally Collapsed Objects of Very Low Mass”  In _Monthly Notices of the Royal Astronomical Society_ 152.1, 1971, pp. 75–78  DOI: 10.1093/mnras/152.1.75
  * [3] B.. Carr and S.. Hawking  “Black Holes in the Early Universe”  In _Monthly Notices of the Royal Astronomical Society_ 168.2, 1974, pp. 399–415  DOI: 10.1093/mnras/168.2.399
  * [4] G.. Chapline  “Cosmological effects of primordial black holes”  In _Nature_ 253.5489, 1975, pp. 251–252  DOI: 10.1038/253251a0
  * [5] Anne M Green and Bradley J Kavanagh  “Primordial black holes as a dark matter candidate”  In _Journal of Physics G: Nuclear and Particle Physics_ 48.4 IOP Publishing, 2021, pp. 043001  DOI: 10.1088/1361-6471/abc534
  * [6] Bernard Carr and Florian Kühnel  “Primordial black holes as dark matter candidates”  In _SciPost Physics Lecture Notes_ Stichting SciPost, 2022  DOI: 10.21468/scipostphyslectnotes.48
  * [7] B.. al. Abbott  “Binary Black Hole Mergers in the First Advanced LIGO Observing Run”  In _Physical Review X_ 6.4 American Physical Society (APS), 2016  DOI: 10.1103/physrevx.6.041015
  * [8] R. al. Abbott  “GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo during the Second Part of the Third Observing Run”  In _Physical Review X_ 13.4 American Physical Society (APS), 2023  DOI: 10.1103/physrevx.13.041039
  * [9] B.. Abbott  “Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA”  In _Living Rev. Rel._ 19, 2016, pp. 1  DOI: 10.1007/s41114-020-00026-9
  * [10] Simeon Bird et al.  “Did LIGO Detect Dark Matter?”  In _Physical Review Letters_ 116.20 American Physical Society (APS), 2016  DOI: 10.1103/physrevlett.116.201301
  * [11] Gabriele Franciolini, Ilia Musco, Paolo Pani and Alfredo Urbano  “From inflation to black hole mergers and back again: Gravitational-wave data-driven constraints on inflationary scenarios with a first-principle model of primordial black holes across the QCD epoch”  In _Physical Review D_ 106.12 American Physical Society (APS), 2022  DOI: 10.1103/physrevd.106.123526
  * [12] Albert Escrivà, Eleni Bagui and Sebastien Clesse  “Simulations of PBH formation at the QCD epoch and comparison with the GWTC-3 catalog”  In _Journal of Cosmology and Astroparticle Physics_ 2023.05 IOP Publishing, 2023, pp. 004  DOI: 10.1088/1475-7516/2023/05/004
  * [13] Eleni Bagui et al.  “Primordial black holes and their gravitational-wave signatures”, 2023  arXiv: https://arxiv.org/abs/2310.19857
  * [14] Albert Escrivà, Florian Kühnel and Yuichiro Tada  “Primordial black holes”  In _Black Holes in the Era of Gravitational-Wave Astronomy_ Elsevier, 2024, pp. 261–377  DOI: 10.1016/b978-0-32-395636-9.00012-8
  * [15] Chul-Moon Yoo  “The basics of primordial black hole formation and abundance estimation”, 2024  arXiv: https://arxiv.org/abs/2211.13512
  * [16] Vicente Atal and Cristiano Germani  “The role of non-gaussianities in Primordial Black Hole formation”  In _Phys. Dark Univ._ 24, 2019, pp. 100275  DOI: 10.1016/j.dark.2019.100275
  * [17] Andrew D. Gow et al.  “Non-perturbative non-Gaussianity and primordial black holes”  In _EPL_ 142.4, 2023, pp. 49001  DOI: 10.1209/0295-5075/acd417
  * [18] Shi Pi  “Non-Gaussianities in primordial black hole formation and induced gravitational waves”, 2024  arXiv: https://arxiv.org/abs/2404.06151
  * [19] Y. Akrami  “Planck 2018 results. X. Constraints on inflation”  In _Astron. Astrophys._ 641, 2020, pp. A10  DOI: 10.1051/0004-6361/201833887
  * [20] Hayato Motohashi and Wayne Hu  “Primordial black holes and slow-roll violation”  In _Phys. Rev. D_ 96 American Physical Society, 2017, pp. 063503  DOI: 10.1103/PhysRevD.96.063503
  * [21] Juan Garcia-Bellido and Ester Ruiz Morales  “Primordial black holes from single field models of inflation”  In _Physics of the Dark Universe_ 18, 2017, pp. 47–54  DOI: https://doi.org/10.1016/j.dark.2017.09.007
  * [22] Philippa S. Cole, Andrew D. Gow, Christian T. Byrnes and Subodh P. Patil  “Primordial black holes from single-field inflation: a fine-tuning audit”  In _JCAP_ 08, 2023, pp. 031  DOI: 10.1088/1475-7516/2023/08/031
  * [23] Wenzer Qin et al.  “Planck constraints and gravitational wave forecasts for primordial black hole dark matter seeded by multifield inflation”  In _Phys. Rev. D_ 108.4, 2023, pp. 043508  DOI: 10.1103/PhysRevD.108.043508
  * [24] Ioanna Stamou  “Mechanisms for producing Primordial Black Holes from Inflationary Models Beyond Fine-Tuning”, 2024  arXiv:2404.14321 [astro-ph.CO]
  * [25] Tomohiro Harada et al.  “PRIMORDIAL BLACK HOLE FORMATION IN THE MATTER-DOMINATED PHASE OF THE UNIVERSE”  In _The Astrophysical Journal_ 833.1 American Astronomical Society, 2016, pp. 61  DOI: 10.3847/1538-4357/833/1/61
  * [26] Tomohiro Harada, Chul-Moon Yoo, Kazunori Kohri and Ken-Ichi Nakao  “Spins of primordial black holes formed in the matter-dominated phase of the Universe”  In _Physical Review D_ 96.8 American Physical Society (APS), 2017  DOI: 10.1103/physrevd.96.083517
  * [27] Valerio De Luca et al.  “Primordial black holes in matter-dominated eras: The role of accretion”  In _Physics Letters B_ 832 Elsevier BV, 2022, pp. 137265  DOI: 10.1016/j.physletb.2022.137265
  * [28] Rouzbeh Allahverdi  “The First Three Seconds: a Review of Possible Expansion Histories of the Early Universe”, 2020  DOI: 10.21105/astro.2006.16182
  * [29] Christian Y. Cardall and George M. Fuller  “Semianalytic analysis of primordial black hole formation during a first order QCD phase transition”, 1998  arXiv:astro-ph/9801103
  * [30] K. Jedamzik and Jens C. Niemeyer  “Primordial black hole formation during first order phase transitions”  In _Phys. Rev. D_ 59, 1999, pp. 124014  DOI: 10.1103/PhysRevD.59.124014
  * [31] Bernard Carr, Sebastien Clesse, Juan Garcia-Bellido and Florian Kuhnel  “Cosmic Conundra Explained by Thermal History and Primordial Black Holes”  In _arXiv e-prints_ , 2019, pp. arXiv:1906.08217  DOI: 10.48550/arXiv.1906.08217
  * [32] Joaquim Iguaz Juan, Pasquale D. Serpico and Guillermo Franco Abellán  “The QCD phase transition behind a PBH origin of LIGO/Virgo events?”  In _JCAP_ 07.07, 2022, pp. 009  DOI: 10.1088/1475-7516/2022/07/009
  * [33] Jens Chluba, Adrienne L. Erickcek and Ido Ben-Dayan  “Probing the inflaton: Small-scale power spectrum constraints from measurements of the CMB energy spectrum”  In _Astrophys. J._ 758, 2012, pp. 76  DOI: 10.1088/0004-637X/758/2/76
  * [34] L.. Abbott, Edward Farhi and Mark B. Wise  “Particle production in the new inflationary cosmology”  In _Physics Letters B_ 117.1-2, 1982, pp. 29–33  DOI: 10.1016/0370-2693(82)90867-X
  * [35] Jennie H. Traschen and Robert H. Brandenberger  “Particle production during out-of-equilibrium phase transitions”  In _Phys. Rev. D_ 42 American Physical Society, 1990, pp. 2491–2504  DOI: 10.1103/PhysRevD.42.2491
  * [36] Bruce A. Bassett, Shinji Tsujikawa and David Wands  “Inflation dynamics and reheating”  In _Reviews of Modern Physics_ 78.2 American Physical Society (APS), 2006, pp. 537–589  DOI: 10.1103/revmodphys.78.537
  * [37] Rouzbeh Allahverdi, Robert Brandenberger, Francis-Yan Cyr-Racine and Anupam Mazumdar  “Reheating in Inflationary Cosmology: Theory and Applications”  In _Annual Review of Nuclear and Particle Science_ 60.1 Annual Reviews, 2010, pp. 27–51  DOI: 10.1146/annurev.nucl.012809.104511
  * [38] Karsten Jedamzik  “Primordial black hole formation during cosmic phase transitions”, 2024  arXiv:2406.11417 [astro-ph.CO]
  * [39] Lars Husdal  “On Effective Degrees of Freedom in the Early Universe”  In _Galaxies_ 4.4, 2016, pp. 78  DOI: 10.3390/galaxies4040078
  * [40] Tomohiro Nakama, Joseph Silk and Marc Kamionkowski  “Stochastic gravitational waves associated with the formation of primordial black holes”  In _Physical Review D_ 95.4 American Physical Society (APS), 2017  DOI: 10.1103/physrevd.95.043511
  * [41] B.. Carr, Kazunori Kohri, Yuuiti Sendouda and Jun’ichi Yokoyama  “New cosmological constraints on primordial black holes”  In _Phys. Rev. D_ 81, 2010, pp. 104019  DOI: 10.1103/PhysRevD.81.104019
  * [42] Georges Aad  “Combined Measurement of the Higgs Boson Mass in $pp$ Collisions at $\sqrt{s}=7$ and 8 TeV with the ATLAS and CMS Experiments”  In _Phys. Rev. Lett._ 114, 2015, pp. 191803  DOI: 10.1103/PhysRevLett.114.191803
  * [43] K. Kajantie, M. Laine, K. Rummukainen and Mikhail E. Shaposhnikov  “Is there a hot electroweak phase transition at $m_{H}\gtrsim m_{W}$?”  In _Phys. Rev. Lett._ 77, 1996, pp. 2887–2890  DOI: 10.1103/PhysRevLett.77.2887
  * [44] Stephan J Huber and Thomas Konstandin  “Production of gravitational waves in the nMSSM”  In _Journal of Cosmology and Astroparticle Physics_ 2008.05 IOP Publishing, 2008, pp. 017  DOI: 10.1088/1475-7516/2008/05/017
  * [45] Stephan J. Huber, Thomas Konstandin, Germano Nardini and Ingo Rues  “Detectable gravitational waves from very strong phase transitions in the general NMSSM”  In _Journal of Cosmology and Astroparticle Physics_ 2016.03 IOP Publishing, 2016, pp. 036–036  DOI: 10.1088/1475-7516/2016/03/036
  * [46] S.V. Demidov, D.S. Gorbunov and D.V. Kirpichnikov  “Gravitational waves from phase transition in split NMSSM”  In _Physics Letters B_ 779 Elsevier BV, 2018, pp. 191–194  DOI: 10.1016/j.physletb.2018.02.007
  * [47] Sebastian Bruggisser, Benedict Harling, Oleksii Matsedonskyi and Geraldine Servant  “Electroweak Phase Transition and Baryogenesis in Composite Higgs Models”, 2018  arXiv:1804.07314 [hep-ph]
  * [48] Ligong Bian, Yongcheng Wu and Ke-Pan Xie  “Electroweak phase transition with composite Higgs models: calculability, gravitational waves and collider searches”  In _Journal of High Energy Physics_ 2019.12 Springer ScienceBusiness Media LLC, 2019  DOI: 10.1007/jhep12(2019)028
  * [49] Mark Hindmarsh, Stephan J. Huber, Kari Rummukainen and David J. Weir  “Gravitational Waves from the Sound of a First Order Phase Transition”  In _Physical Review Letters_ 112.4 American Physical Society (APS), 2014  DOI: 10.1103/physrevlett.112.041301
  * [50] A. Bazavov et al.  “Equation of state in ($2+1$)-flavor QCD”  In _Phys. Rev. D_ 90 American Physical Society, 2014, pp. 094503  DOI: 10.1103/PhysRevD.90.094503
  * [51] Szabocls Borsanyi et al.  “Full result for the QCD equation of state with 2+1 flavors”  In _Phys. Lett. B_ 730, 2014, pp. 99–104  DOI: 10.1016/j.physletb.2014.01.007
  * [52] Sz. Borsanyi et al.  “Lattice QCD for Cosmology”, 2016  arXiv:1606.07494 [hep-lat]
  * [53] Ilia Musco, Karsten Jedamzik and Sam Young  “Primordial black hole formation during the QCD phase transition: threshold, mass distribution and abundance”, 2023  arXiv:2303.07980 [astro-ph.CO]
  * [54] Ilia Musco  “Threshold for primordial black holes: Dependence on the shape of the cosmological perturbations”  In _Physical Review D_ 100.12 American Physical Society (APS), 2019  DOI: 10.1103/physrevd.100.123524
  * [55] Albert Escrivà, Cristiano Germani and Ravi K. Sheth  “Universal threshold for primordial black hole formation”  In _Physical Review D_ 101.4 American Physical Society (APS), 2020  DOI: 10.1103/physrevd.101.044022
  * [56] Ilia Musco, Valerio De Luca, Gabriele Franciolini and Antonio Riotto  “Threshold for primordial black holes. II. A simple analytic prescription”  In _Physical Review D_ 103.6 American Physical Society (APS), 2021  DOI: 10.1103/physrevd.103.063538
  * [57] Ilia Musco and John C Miller  “Primordial black hole formation in the early universe: critical behaviour and self-similarity”  In _Classical and Quantum Gravity_ 30.14 IOP Publishing, 2013, pp. 145009  DOI: 10.1088/0264-9381/30/14/145009
  * [58] Christian T. Byrnes, Mark Hindmarsh, Sam Young and Michael R.. Hawkins  “Primordial black holes with an accurate QCD equation of state”  In _Journal of Cosmology and Astrophysics_ 2018.8, 2018, pp. 041  DOI: 10.1088/1475-7516/2018/08/041
  * [59] Kenta Ando, Keisuke Inomata and Masahiro Kawasaki  “Primordial black holes and uncertainties in the choice of the window function”  In _Phys. Rev. D_ 97.10, 2018, pp. 103528  DOI: 10.1103/PhysRevD.97.103528
  * [60] Alba Kalaja et al.  “From Primordial Black Holes Abundance to Primordial Curvature Power Spectrum (and back)”  In _JCAP_ 10, 2019, pp. 031  DOI: 10.1088/1475-7516/2019/10/031
  * [61] Valerio De Luca, Alex Kehagias and Antonio Riotto  “How well do we know the primordial black hole abundance: The crucial role of nonlinearities when approaching the horizon”  In _Phys. Rev. D_ 108.6, 2023, pp. 063531  DOI: 10.1103/PhysRevD.108.063531
  * [62] Gabriele Franciolini et al.  “Renormalized Primordial Black Holes”, 2023  arXiv: https://arxiv.org/abs/2311.03239
  * [63] Sam Young  “The primordial black hole formation criterion re-examined: Parametrisation, timing and the choice of window function”  In _International Journal of Modern Physics D_ 29.02 World Scientific Pub Co Pte Lt, 2019, pp. 2030002  DOI: 10.1142/s0218271820300025
  * [64] William H. Press and Paul Schechter  “Formation of Galaxies and Clusters of Galaxies by Self-Similar Gravitational Condensation”  In _Astrophysical Journal_ 187, 1974, pp. 425–438  DOI: 10.1086/152650
  * [65] J.. Bardeen, J.. Bond, N. Kaiser and A.. Szalay  “The Statistics of Peaks of Gaussian Random Fields”  In _Astrophysical Journal_ 304, 1986, pp. 15  DOI: 10.1086/164143
  * [66] Andrew D. Gow, Christian T. Byrnes, Philippa S. Cole and Sam Young  “The power spectrum on small scales: robust constraints and comparing PBH methodologies”  In _Journal of Cosmology and Astroparticle Physics_ 2021.02 IOP Publishing, 2021, pp. 002–002  DOI: 10.1088/1475-7516/2021/02/002
  * [67] N. al. Aghanim  “Planck2018 results: VI. Cosmological parameters”  In _Astronomy & Astrophysics_ 641 EDP Sciences, 2020, pp. A6  DOI: 10.1051/0004-6361/201833910
  * [68] M. Andrés-Carcasona et al.  “Constraints on primordial black holes from LIGO-Virgo-KAGRA O3 events”, 2024  arXiv:2405.05732 [astro-ph.CO]
  * [69] Hiroko Niikura et al.  “Constraints on Earth-mass primordial black holes from OGLE 5-year microlensing events”  In _Physical Review D_ 99.8 American Physical Society (APS), 2019  DOI: 10.1103/physrevd.99.083503
  * [70] Hiroko Niikura et al.  “Microlensing constraints on primordial black holes with Subaru/HSC Andromeda observations”  In _Nature Astronomy_ 3, 2019, pp. 524–534  DOI: 10.1038/s41550-019-0723-1
  * [71] Gabriele Franciolini, Antonio Junior Iovino, Ville Vaskonen and Hardi Veermäe  “Recent Gravitational Wave Observation by Pulsar Timing Arrays and Primordial Black Holes: The Importance of Non-Gaussianities”  In _Physical Review Letters_ 131.20 American Physical Society (APS), 2023  DOI: 10.1103/physrevlett.131.201401
  * [72] Raf Guedens, Dominic Clancy and Andrew R Liddle  “Primordial black holes in braneworld cosmologies: Formation, cosmological evolution and evaporation”  In _Phys. Rev. D_ 66, 2002, pp. 043513  DOI: 10.1103/PhysRevD.66.043513
  * [73] Bernard Carr, Kazunori Kohri, Yuuiti Sendouda and Jun’ichi Yokoyama  “Constraints on primordial black holes”  In _Reports on Progress in Physics_ 84.11 IOP Publishing, 2021, pp. 116902  DOI: 10.1088/1361-6633/ac1e31
  * [74] B.J. Carr et al.  “Observational evidence for primordial black holes: A positivist perspective”  In _Physics Reports_ 1054 Elsevier BV, 2024, pp. 1–68  DOI: 10.1016/j.physrep.2023.11.005
  * [75] Jonathan Kozaczuk, Tongyan Lin and Ethan Villarama  “Signals of primordial black holes at gravitational wave interferometers”  In _Physical Review D_ 105.12 American Physical Society (APS), 2022  DOI: 10.1103/physrevd.105.123023
  * [76] Tomohiro Harada, Chul-Moon Yoo, Tomohiro Nakama and Yasutaka Koga  “Cosmological long-wavelength solutions and primordial black hole formation”  In _Physical Review D_ 91.8 American Physical Society (APS), 2015  DOI: 10.1103/physrevd.91.084057
  * [77] Sam Young, Christian T. Byrnes and Misao Sasaki  “Calculating the mass fraction of primordial black holes”  In _JCAP_ 07, 2014, pp. 045  DOI: 10.1088/1475-7516/2014/07/045
  * [78] Qing Wang, Yi-Chen Liu, Bing-Yu Su and Nan Li  “Primordial black holes from the perturbations in the inflaton potential in peak theory”  In _Phys. Rev. D_ 104.8, 2021, pp. 083546  DOI: 10.1103/PhysRevD.104.083546
  * [79] Jens Chluba, Jan Hamann and Subodh P. Patil  “Features and New Physical Scales in Primordial Observables: Theory and Observation”  In _Int. J. Mod. Phys. D_ 24.10, 2015, pp. 1530023  DOI: 10.1142/S0218271815300232
  * [80] Tomohiro Nakama, Bernard Carr and Joseph Silk  “Limits on primordial black holes from $\mu$ distortions in cosmic microwave background”  In _Phys. Rev. D_ 97.4, 2018, pp. 043525  DOI: 10.1103/PhysRevD.97.043525
  * [81] Dan Hooper, Aurora Ireland, Gordan Krnjaic and Albert Stebbins  “Supermassive Primordial Black Holes From Inflation”, 2023  arXiv:2308.00756 [astro-ph.CO]
  * [82] Devanshu Sharma, Julien Lesgourgues and Christian T. Byrnes  “Spectral distortions from acoustic dissipation with non-Gaussian (or not) perturbations”, 2024  arXiv:2404.18474 [astro-ph.CO]
  * [83] Christian T. Byrnes, Julien Lesgourgues and Devanshu Sharma  “Robust $\mu$-distortion constraints on primordial supermassive black holes from non-Gaussian perturbations”, 2024  arXiv:2404.18475 [astro-ph.CO]
  * [84] A.. Iovino, G. Perna, A. Riotto and H. Veermäe  “Curbing PBHs with PTAs”, 2024  arXiv:2406.20089 [astro-ph.CO]
  * [85] Federico Bianchini and Giulio Fabbian  “CMB spectral distortions revisited: A new take on $\mu$ distortions and primordial non-Gaussianities from FIRAS data”  In _Phys. Rev. D_ 106.6, 2022, pp. 063527  DOI: 10.1103/PhysRevD.106.063527
  * [86] D.. Fixsen et al.  “The Cosmic Microwave Background Spectrum from the FullCOBEFIRAS Data Set”  In _The Astrophysical Journal_ 473.2 American Astronomical Society, 1996, pp. 576–587  DOI: 10.1086/178173
  * [87] Guillem Domènech  “Scalar Induced Gravitational Waves Review”  In _Universe_ 7.11, 2021, pp. 398  DOI: 10.3390/universe7110398
  * [88] Gabriella Agazie al  “The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background”  In _The Astrophysical Journal Letters_ 951.1 The American Astronomical Society, 2023, pp. L8  DOI: 10.3847/2041-8213/acdac6
  * [89] Gabriella Agazie al.  “The NANOGrav 15 yr Data Set: Observations and Timing of 68 Millisecond Pulsars”  In _The Astrophysical Journal Letters_ 951.1 The American Astronomical Society, 2023, pp. L9  DOI: 10.3847/2041-8213/acda9a
  * [90] J. al. Antoniadis  “The second data release from the European Pulsar Timing Array: III. Search for gravitational wave signals”  In _Astronomy & Astrophysics_ 678 EDP Sciences, 2023, pp. A50  DOI: 10.1051/0004-6361/202346844
  * [91] J.. Antoniadis  “The second data release from the European Pulsar Timing Array: I. The dataset and timing analysis”  In _Astronomy & Astrophysics_ 678 EDP Sciences, 2023, pp. A48  DOI: 10.1051/0004-6361/202346841
  * [92] J. al.  “The second data release from the European Pulsar Timing Array: V. Implications for massive black holes, dark matter and the early Universe”, 2023  arXiv:2306.16227 [astro-ph.CO]
  * [93] Daniel J. al. Reardon  “Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array”  In _The Astrophysical Journal Letters_ 951.1 American Astronomical Society, 2023, pp. L6  DOI: 10.3847/2041-8213/acdd02
  * [94] Andrew Zic al.  “The Parkes Pulsar Timing Array Third Data Release”, 2023  arXiv:2306.16230 [astro-ph.HE]
  * [95] Daniel J. al. Reardon  “The Gravitational-wave Background Null Hypothesis: Characterizing Noise in Millisecond Pulsar Arrival Times with the Parkes Pulsar Timing Array”  In _The Astrophysical Journal Letters_ 951.1 American Astronomical Society, 2023, pp. L7  DOI: 10.3847/2041-8213/acdd03
  * [96] Heng et al. Xu  “Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I”  In _Research in Astronomy and Astrophysics_ 23.7 IOP Publishing, 2023, pp. 075024  DOI: 10.1088/1674-4527/acdfa5
  * [97] Zhi-Qiang You, Zhu Yi and You Wu  “Constraints on primordial curvature power spectrum with pulsar timing arrays”  In _Journal of Cosmology and Astroparticle Physics_ 2023.11 IOP Publishing, 2023, pp. 065  DOI: 10.1088/1475-7516/2023/11/065
  * [98] Qin Fei  “Constraints on the Primordial Curvature Power Spectrum by Pulsar Timing Array Data: A Polynomial Parameterization Approach”, 2024  arXiv:2310.17199 [gr-qc]
  * [99] Yi-Fu Cai et al.  “Limits on scalar-induced gravitational waves from the stochastic background by pulsar timing array observations”  In _Science Bulletin_ 68.23 Elsevier BV, 2023, pp. 2929–2935  DOI: 10.1016/j.scib.2023.10.027
  * [100] Shyam Balaji, Guillem Domènech and Gabriele Franciolini  “Scalar-induced gravitational wave interpretation of PTA data: the role of scalar fluctuation propagation speed”, 2023  arXiv:2307.08552 [gr-qc]
  * [101] Gabriele Franciolini, Antonio Junior Iovino, Ville Vaskonen and Hardi Veermäe  “Recent Gravitational Wave Observation by Pulsar Timing Arrays and Primordial Black Holes: The Importance of Non-Gaussianities”  In _Physical Review Letters_ 131.20 American Physical Society (APS), 2023  DOI: 10.1103/physrevlett.131.201401
  * [102] Zhu Yi et al.  “Exploring the NANOGrav Signal and Planet-mass Primordial Black Holes through Higgs Inflation”, 2023  arXiv:2308.14688 [astro-ph.CO]
  * [103] Mohammad Ali Gorji, Misao Sasaki and Teruaki Suyama  “Extra-tensor-induced origin for the PTA signal: No primordial black hole production”, 2023  arXiv:2307.13109 [astro-ph.CO]
  * [104] Sayantan Choudhury, Kritartha Dey and Ahaskar Karde  “Untangling PBH overproduction in $w$-SIGWs generated by Pulsar Timing Arrays for MST-EFT of single field inflation”, 2023  arXiv:2311.15065 [astro-ph.CO]
  * [105] Gourab Bhattacharya et al.  “Evading no-go for PBH formation and production of SIGWs using Multiple Sharp Transitions in EFT of single field inflation”, 2023  arXiv:2309.00973 [astro-ph.CO]
  * [106] Keisuke Inomata, Kazunori Kohri and Takahiro Terada  “The Detected Stochastic Gravitational Waves and Subsolar-Mass Primordial Black Holes”, 2023  arXiv:2306.17834 [astro-ph.CO]
  * [107] Paul Frederik Depta, Kai Schmidt-Hoberg, Pedro Schwaller and Carlo Tasillo  “Do pulsar timing arrays observe merging primordial black holes?”, 2023  arXiv:2306.17836 [astro-ph.CO]
  * [108] Nilanjandev Bhaumik, Rajeev Kumar Jain and Marek Lewicki  “Ultralow mass primordial black holes in the early Universe can explain the pulsar timing array signal”  In _Physical Review D_ 108.12 American Physical Society (APS), 2023  DOI: 10.1103/physrevd.108.123532
  * [109] Adeela Afzal al  “The NANOGrav 15 yr Data Set: Search for Signals from New Physics”  In _The Astrophysical Journal Letters_ 951.1 The American Astronomical Society, 2023, pp. L11  DOI: 10.3847/2041-8213/acdc91
  * [110] Haipeng An et al.  “Phase transition during inflation and the gravitational wave signal at pulsar timing arrays”, 2023  arXiv:2308.00070 [astro-ph.CO]
  * [111] Gabriele Franciolini, Davide Racco and Fabrizio Rompineve  “Footprints of the QCD Crossover on Cosmological Gravitational Waves at Pulsar Timing Arrays”, 2023  arXiv:2306.17136 [astro-ph.CO]
  * [112] Andrea Addazi, Yi-Fu Cai, Antonino Marcianò and Luca Visinelli  “Have pulsar timing array methods detected a cosmological phase transition?”  In _Physical Review D_ 109.1 American Physical Society (APS), 2024  DOI: 10.1103/physrevd.109.015028
  * [113] Kohei Fujikura, Sudhakantha Girmohanta, Yuichiro Nakai and Motoo Suzuki  “NANOGrav Signal from a Dark Conformal Phase Transition”, 2023  arXiv:2306.17086 [hep-ph]
  * [114] Martin Wolfgang Winkler and Katherine Freese  “Origin of the Stochastic Gravitational Wave Background: First-Order Phase Transition vs. Black Hole Mergers”, 2024  arXiv:2401.13729 [astro-ph.CO]
  * [115] Lei Zu et al.  “Mirror QCD phase transition as the origin of the nanohertz Stochastic Gravitational-Wave Background”  In _Science Bulletin_ Elsevier BV, 2024  DOI: 10.1016/j.scib.2024.01.037
  * [116] Deng Wang  “Constraining Cosmological Phase Transitions with Chinese Pulsar Timing Array Data Release 1”, 2023  arXiv:2307.15970 [astro-ph.CO]
  * [117] Ricardo Z. Ferreira, Alessio Notari, Oriol Pujolàs and Fabrizio Rompineve  “Collapsing Domain Wall Networks: Impact on Pulsar Timing Arrays and Primordial Black Holes”, 2024  arXiv:2401.14331 [astro-ph.CO]
  * [118] Yann Gouttenoire and Edoardo Vitagliano  “Domain wall interpretation of the PTA signal confronting black hole overproduction”  In _arXiv e-prints_ , 2023, pp. arXiv:2306.17841  DOI: 10.48550/arXiv.2306.17841
  * [119] John Ellis, Marek Lewicki, Chunshan Lin and Ville Vaskonen  “Cosmic Superstrings Revisited in Light of NANOGrav 15-Year Data”, 2023  arXiv:2306.17147 [astro-ph.CO]
  * [120] George Lazarides, Rinku Maji and Qaisar Shafi  “Superheavy quasistable strings and walls bounded by strings in the light of NANOGrav 15 year data”  In _Physical Review D_ 108.9 American Physical Society (APS), 2023  DOI: 10.1103/physrevd.108.095041
  * [121] Ziwei Wang et al.  “The nanohertz stochastic gravitational wave background from cosmic string loops and the abundant high redshift massive galaxies”  In _Science China Physics, Mechanics & Astronomy_ 66.12 Springer ScienceBusiness Media LLC, 2023  DOI: 10.1007/s11433-023-2262-0
  * [122] Zhao Zhang et al.  “Nano-Hertz gravitational waves from collapsing domain walls associated with freeze-in dark matter in light of pulsar timing array observations”  In _Physical Review D_ 108.9 American Physical Society (APS), 2023  DOI: 10.1103/physrevd.108.095037
  * [123] John Ellis et al.  “Gravitational waves from supermassive black hole binaries in light of the NANOGrav 15-year data”  In _Physical Review D_ 109.2 American Physical Society (APS), 2024  DOI: 10.1103/physrevd.109.l021302
  * [124] Yan-Chen Bi, Yu-Mei Wu, Zu-Cheng Chen and Qing-Guo Huang  “Implications for the supermassive black hole binaries from the NANOGrav 15-year data set”  In _Science China Physics, Mechanics & Astronomy_ 66.12 Springer ScienceBusiness Media LLC, 2023  DOI: 10.1007/s11433-023-2252-4
  * [125] Emiko C. Gardiner, Luke Zoltan Kelley, Anna-Malin Lemke and Andrea Mitridate  “Beyond the Background: Gravitational Wave Anisotropy and Continuous Waves from Supermassive Black Hole Binaries”, 2023  arXiv:2309.07227 [astro-ph.HE]
  * [126] Hao Ma, Youjun Lu, Zhiwei Chen and Yunfeng Chen  “Diffractive lensing of nano-Hertz gravitational waves emitted from supermassive binary black holes by intervening galaxies”  In _Monthly Notices of the Royal Astronomical Society_ 524.2 Oxford University Press (OUP), 2023, pp. 2954–2965  DOI: 10.1093/mnras/stad2040
  * [127] Seppe Staelens and Gijs Nelemans  “The Astrophysical Gravitational Wave Background in the mHz band is likely dominated by White Dwarf binaries”, 2023  arXiv:2310.19448 [astro-ph.HE]
  * [128] Lukas T. Witkowski  “SIGWfast: a python package for the computation of scalar-induced gravitational wave spectra”  In _arXiv e-prints_ , 2022, pp. arXiv:2209.05296  DOI: 10.48550/arXiv.2209.05296
  * [129] Guillem Domenech  “Scalar Induced Gravitational Waves Review”  In _Universe_ 7.11 MDPI AG, 2021, pp. 398  DOI: 10.3390/universe7110398
  * [130] Gabriella et al. Agazie  “The NANOGrav 15 yr Data Set: Bayesian Limits on Gravitational Waves from Individual Supermassive Black Hole Binaries”  In _Astrophysical Journal, Letters_ 951.2, 2023, pp. L50  DOI: 10.3847/2041-8213/ace18a
  * [131] Brian C. Lacki and John F. Beacom  “Primordial Black Holes as Dark Matter: Almost All or Almost Nothing”  In _Astrophys. J. Lett._ 720, 2010, pp. L67–L71  DOI: 10.1088/2041-8205/720/1/L67
  * [132] Julian Adamek, Christian T. Byrnes, Mateja Gosenca and Shaun Hotchkiss  “WIMPs and stellar-mass primordial black holes are incompatible”  In _Phys. Rev. D_ 100.2, 2019, pp. 023506  DOI: 10.1103/PhysRevD.100.023506
  * [133] Bernard Carr, Florian Kuhnel and Luca Visinelli  “Black holes and WIMPs: all or nothing or something else”  In _Mon. Not. Roy. Astron. Soc._ 506.3, 2021, pp. 3648–3661  DOI: 10.1093/mnras/stab1930
  * [134] Mathieu Boudaud et al.  “In-depth analysis of the clustering of dark matter particles around primordial black holes. Part I. Density profiles”  In _JCAP_ 08, 2021, pp. 053  DOI: 10.1088/1475-7516/2021/08/053
  * [135] Kenji Kadota and Hiroyuki Tashiro  “Primordial black hole dark matter in the presence of p-wave WIMP annihilation”  In _JCAP_ 03.03, 2022, pp. 045  DOI: 10.1088/1475-7516/2022/03/045
  * [136] Estanis Utrilla Ginés, Olga Mena and Samuel J. Witte  “Revisiting constraints on WIMPs around primordial black holes”  In _Phys. Rev. D_ 106.6, 2022, pp. 063538  DOI: 10.1103/PhysRevD.106.063538
  * [137] Yury N. Eroshenko  “Evolution of the DM distribution function in the density spikes around PBHs”, 2024  arXiv:2407.10225 [astro-ph.CO]
  * [138] E. Mediavilla et al.  “Limits on the Mass and Abundance of Primordial Black Holes from Quasar Gravitational Microlensing”  In _The Astrophysical Journal Letters_ 836.2 American Astronomical Society, 2017, pp. L18  DOI: 10.3847/2041-8213/aa5dab
  * [139] Łukasz Wyrzykowski and Ilya Mandel  “Constraining the masses of microlensing black holes and the mass gap with Gaia DR2”  In _Astronomy & Astrophysics_ 636 EDP Sciences, 2020, pp. A20  DOI: 10.1051/0004-6361/201935842
  * [140] A. Kashlinsky, R.. Arendt, J. Mather and S.. Moseley  “Tracing the first stars with fluctuations of the cosmic infrared background”  In _Nature_ 438.7064 Springer ScienceBusiness Media LLC, 2005, pp. 45–50  DOI: 10.1038/nature04143
  * [141] Sam Young  “Computing the abundance of primordial black holes”, 2024  arXiv:2405.13259 [astro-ph.CO]
  * [142] Sam Young, Ilia Musco and Christian T. Byrnes  “Primordial black hole formation and abundance: contribution from the non-linear relation between the density and curvature perturbation”  In _JCAP_ 11, 2019, pp. 012  DOI: 10.1088/1475-7516/2019/11/012
  * [143] Sam Young and Marcello Musso  “Application of peaks theory to the abundance of primordial black holes”  In _Journal of Cosmology and Astroparticle Physics_ 2020.11 IOP Publishing, 2020, pp. 022–022  DOI: 10.1088/1475-7516/2020/11/022
  * [144] Christian T. Byrnes, Philippa S. Cole and Subodh P. Patil  “Steepest growth of the power spectrum and primordial black holes”  In _Journal of Cosmology and Astroparticle Physics_ 2019.06, 2019, pp. 028  DOI: 10.1088/1475-7516/2019/06/028
  * [145] Rongrong Zhai, Hongwei Yu and Puxun Wu  “Power spectrum with $k^{6}$ growth for primordial black holes”, 2024  arXiv:2308.09286 [gr-qc]
  * [146] Gianmassimo Tasinato  “An analytic approach to non-slow-roll inflation”  In _Phys. Rev. D_ 103.2, 2021, pp. 023535  DOI: 10.1103/PhysRevD.103.023535
  * [147] Pedro Carrilho, Karim A. Malik and David J. Mulryne  “Dissecting the growth of the power spectrum for primordial black holes”  In _Phys. Rev. D_ 100.10, 2019, pp. 103529  DOI: 10.1103/PhysRevD.100.103529
  * [148] H.. Ragavendra, Pankaj Saha, L. Sriramkumar and Joseph Silk  “Primordial black holes and secondary gravitational waves from ultraslow roll and punctuated inflation”  In _Phys. Rev. D_ 103.8, 2021, pp. 083510  DOI: 10.1103/PhysRevD.103.083510
  * [149] Sam Young  “Peaks and primordial black holes: the effect of non-Gaussianity”  In _Journal of Cosmology and Astroparticle Physics_ 2022.05 IOP Publishing, 2022, pp. 037  DOI: 10.1088/1475-7516/2022/05/037
  * [150] Ken’ichi Saikawa and Satoshi Shirai  “Primordial gravitational waves, precisely: the role of thermodynamics in the Standard Model”  In _Journal of Cosmology and Astroparticle Physics_ 2018.05 IOP Publishing, 2018, pp. 035–035  DOI: 10.1088/1475-7516/2018/05/035
  * [151] Katsuya T. Abe and Yuichiro Tada  “Translating nano-Hertz gravitational wave background into primordial perturbations taking account of the cosmological QCD phase transition”  In _Phys. Rev. D_ 108.10, 2023, pp. L101304  DOI: 10.1103/PhysRevD.108.L101304