# Primordial black holes as dark matter: Interferometric tests of phase transition origin

Iason Baldes, María Olalla Olea-Romacho

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

To cite this version:

Iason Baldes, María Olalla Olea-Romacho. Primordial black holes as dark matter: Interferometric tests of phase transition origin. Journal of High Energy Physics, 2024, 2024 (01), pp.133. 10.1007/JHEP01(2024)133. hal-04179145

HAL Id: hal-04179145

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

Submitted on 23 Aug 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
RECEIVED: September 21, 2023

REVISED: December 22, 2023

ACCEPTED: January 10, 2024

PUBLISHED: January 23, 2024

# Primordial black holes as dark matter: interferometric tests of phase transition origin

Iason Baldes and María Olalla Olea-Romacho

Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, 24 rue Lhomond, F-75005 Paris, France

E-mail: iasonbaldes@gmail.com, mariaolalla.olearomacho@phys.ens.fr

**ABSTRACT:** We show that primordial black holes — in the observationally allowed mass window with $f_{PBH} = 1$ — formed from late nucleating patches in a first order phase transition imply upcoming gravitational wave interferometers will see a large stochastic background arising from the bubble collisions. As an example, we use a classically scale invariant $B-L$ model, in which the right handed neutrinos explain the neutrino masses and leptogenesis, and the dark matter consists of primordial black holes. The conclusion regarding the gravitational waves is, however, expected to hold model independently for black holes coming from such late nucleating patches.

**KEYWORDS:** Cosmology of Theories BSM, Early Universe Particle Physics, Phase Transitions in the Early Universe

ARXIV ePRINT: 2307.11639
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>A simple argument</td><td>3</td></tr><tr><td>3</td><td>The model</td><td>5</td></tr><tr><td>3.1</td><td>Particles and charges</td><td>5</td></tr><tr><td>3.2</td><td>Radiative symmetry breaking at finite temperature</td><td>5</td></tr><tr><td>3.3</td><td>Bubble nucleation and percolation</td><td>7</td></tr><tr><td>4</td><td>Primordial black hole production</td><td>10</td></tr><tr><td>4.1</td><td>Background evolution</td><td>10</td></tr><tr><td>4.2</td><td>Late patch evolution</td><td>12</td></tr><tr><td>4.3</td><td>Collapse and fractional abundance</td><td>13</td></tr><tr><td>4.4</td><td>Application to example model</td><td>16</td></tr><tr><td>5</td><td>Gravitational waves</td><td>17</td></tr><tr><td>5.1</td><td>The spectra</td><td>17</td></tr><tr><td>5.2</td><td>Signal-to-noise ratio</td><td>21</td></tr><tr><td>6</td><td>Validity of assumptions</td><td>23</td></tr><tr><td>6.1</td><td>Rapid decay of the condensate excitations</td><td>23</td></tr><tr><td>6.2</td><td>Pressure on the bubble wall</td><td>24</td></tr><tr><td>6.3</td><td>Reliability of the effective potential</td><td>25</td></tr><tr><td>7</td><td>Comments on leptogenesis</td><td>25</td></tr><tr><td>8</td><td>Conclusion</td><td>26</td></tr><tr><td>A</td><td>Towards a non-monochromatic mass distribution</td><td>27</td></tr></table>

## 1 Introduction

Putting aside the possibility being observationally misled by a theory of modified Newtonian dynamics [1–6], the microscopic nature of the cosmic fluid labelled dark matter (DM) is an open question. One possibility is for all of the DM to be made up of primordial black holes (PBHs) [7–9]. In this paper, we adopt the standard notation of denoting the PBH-to-DM density ratio, and assume throughout

$$f_{\text{PBH}} \equiv \frac{\rho_{\text{PBH}}}{\rho_{\text{DM}}} = 1. \qquad (1.1)$$

Observationally, this is allowed for approximately monochromatic PBH distributions with mass in the range [10],

$$10^{-16} M_{\odot} \lesssim M_{\text{PBH}} \lesssim 10^{-10} M_{\odot}. \qquad (1.2)$$
Here the lower limit comes from the would-be-flux of particles coming from Hawking evap-
oration — which would have altered the CMB [11, 12], produced too much extragalactic
background light [13–15], been detected by Voyager [16] or SPI/INTEGRAL [17–19] — and
the upper limit comes from microlensing observations [20].

In order for PBHs to be produced in the early universe, physics beyond the standard model (SM) of particles and cosmology is required. One possibility is for a feature in the inflaton potential to give an enhancement in the amplitude of overdensities at small scales [21, 22]. When the overdensities re-enter the Hubble horizon, these will collapse into PBHs, provided the density contrast is sufficiently large for the given equation-of-state and shape of the curvature power spectrum [23–27, 27–33]. As the evidence for inflation is rather compelling — albeit circumstantial — this possibility has garnered a good deal of attention. In such scenarios, gravitational waves (GWs) are produced at second order in perturbation theory, when the enhanced curvature perturbations at small scales re-enter the Hubble horizon. Such GWs may be detected at future interferometers [34–39].

An alternative is for the overdensities which collapse into PBHs to be produced at a later stage in the evolution of the universe. This could occur if there was a sufficiently strong first-order phase transition (PT) between the end of slow roll inflation and big bang nucleosynthesis. This possibility has also been the topic of a number of papers which have appeared in the last four decades, albeit more sporadically than those regarding the inflation scenario, and a number of mechanisms allowing for PBH formation during a PT have been identified [40–57].

Among the PT mechanisms that trigger PBH formation, the formation of large false
vacuum remnants from stochastically late-nucleated patches in supercooled phase transitions
stands as an intriguing possibility [42, 43, 47, 48, 52–57]. This is an unavoidable process if
sufficiently strong PTs take place. The remnants lead to Hubble sized overdensities being
present in the radiation following post-PT reheating, most closely emulating the inflation
mechanism. The concurrent production of smaller PBHs has also recently been considered
in [55], using a Hoop conjecture criterion, but the spectrum remains relatively peaked.<fnref n="1" />

The reheating following bubble percolation in the PT leads to an initially inhomogenous
state, not modelled in General Relativistic numerical simulations of PBH formation [24–
27, 27–32], but the large number of bubbles present and the fact that an $O(1)$ overdensity
is achieved, leads us to assume that collapse is rather plausible. Therefore, for the scope
of this paper, we operate under the assumption that this mechanism functions as expected,
along the lines of [47, 48, 52–54, 56, 57]. Nevertheless, we warn the reader that the collapse
criterion in this scenario, which will necessarily feature a certain level of anisotropy, turbulent
inhomogeneity, non-sphericity, and initial velocity vectors for the radiation fluid differing
from the inflationary overdensity scenario, is far from certain.

Strong first order PTs also produce a stochastic background of GWs [64, 65]. The production of the PBHs in the PT will lead to GWs from the bubble collisions, and these GWs are observable in upcoming interferometers, as we explain next using a simple argument. In the rest of the subsequent paper we provide a more detailed calculation.

<footnote n="1">¹For PBHs produced from PTs during or at the end of inflation see [58–63].</footnote>
## 2 A simple argument

In the delayed patch mechanism, the mass of the PBH is set by the energy contained inside a Hubble volume during the PT,

$$M_{\text{PBH}} \sim \frac{M_{\text{Pl}}^3}{T_{\text{RH}}^2}, \qquad (2.1)$$

where $M_{Pl} \sim 10^{19}$ GeV is the Planck mass, and $T_{RH}$ is the reheating temperature following the completion of the PT (we assume rapid reheating). The observationally allowed mass window for which $f_{\text{PBH}} = 1$ is possible, eq. (1.2), constrains the reheating temperature to lie between [56]

$$10 \text{ TeV} \lesssim T_{\text{RH}} \lesssim 10^4 \text{ TeV}. \qquad (2.2)$$

The abundance is controlled by the time the PT takes to complete compared to the Hubble rate. To achieve an abundance of $f_{\text{PBH}} = 1$, we require the inverse timescale of the transition to be [56],

$$\beta \equiv \frac{1}{\Gamma_{\text{bub}}} \frac{d\Gamma_{\text{bub}}}{dt} \approx 8H, \qquad (2.3)$$

where $\Gamma_{\text{bub}}$ is the bubble nucleation rate per unit volume and $H$ is the Hubble rate during the PT.

The amplitude of GWs from bubble collisions during the PT scales as $\propto R^2 \sim v_w^2/\beta^2$, where $R$ is the bubble radius and $v_w$ is the wall velocity. For the PTs of interest here, it is a good approximation to set $v_w \simeq 1$, and we do so throughout the rest of this paper. The precise amplitude and spectral shape of the GWs has been addressed in a number of works [66–77]. In our calculations we shall use the recent determinations, relevant for our type of PTs, in which the peak amplitude redshifted to today is [73, 75, 77]

$$\Omega_{\text{GW}} \sim 10^{-6} \left(\frac{H}{\beta}\right)^2. \qquad (2.4)$$

Said GW estimates still feature uncertainties for the very strong PTs we will be interested in, due to the expansion of the universe during the PT itself, see [78, 79]. So in this regard our results should be considered preliminary, albeit promising, until the GW predictions can be further refined. The peak frequency of the GW signal is similarly set by the inverse bubble radius at the collision time, which is then redshifted to today,

$$f_{\text{peak}} \sim 1 \text{ mHz} \left(\frac{\beta}{H}\right) \left(\frac{T_{\text{RH}}}{100 \text{ TeV}}\right). \qquad (2.5)$$

From the above, we see that in such a $f_{\text{PBH}} = 1$ from a PT scenario, we will have a peak amplitude $\Omega_{\text{GW}} \sim 10^{-8}$, i.e. well above stochastic astrophysical foregrounds, somewhere in the frequency range $10^{-3}$ Hz $\lesssim f_{\text{peak}} \lesssim 1$ Hz. The GWs at lower peak frequencies are detectable by LISA [80]. Those at higher peak frequencies are similarly detectable by the Einstein Telescope (ET) [81]. A mid-frequency interferometer such as BDECIGO [82] would
be sensitive to the entire allowed PBH mass range and could in any case help confirm and characterise the signal.

Similar conclusions hold if we replace a detector with a variant of broadly similar equivalence, such as ET with Cosmic Explorer [83] or BDECIGO with AEDGE [84]. Note we consider only the GWs from the bubble wall collisions, the calculation and addition of GWs sourced from the enhanced curvature perturbations at small scales in the context of the PT mechanism [85], is left for future work.

The rest of the paper is dedicated to calculating the PBH production and GW signal within a simple beyond the SM scenario, motivated by neutrino masses and leptogenesis, for which the PBHs will act as a DM candidate. We will scan over the parameters giving $M_{\text{PBH}}$ in the allowed window. The model we choose is a classically scale invariant, gauged $B-L$, extension of the SM. Such close-to-conformal models can provide the necessary supercooled PTs [86–88]. Even though fine-tuning is required to get the cosmological constant close-to-zero today, this seems to be the case almost generically when dealing with the vacuum energy density in a cosmological context. We therefore put questions of fine-tuning aside, which may anyway similarly be present in other models triggering PBH formation, and simply explore the scenario of achieving the observed $\rho_{\text{DM}}$.

We find it beneficial to conduct an explicit calculation within the aforementioned beyond the SM physics scenario. The calculation in such a model serves to illuminate several facets of the subject matter. For instance, many previous studies (e.g. [47, 56]) assume a model independent approximation for the bubble nucleation rate

$$ \Gamma_{\text{bub}}(t) = H^4(t_n)e^{\beta(t-t_n)}, \qquad (2.6) $$

where $t_n$ is the bubble nucleation time, by definition when $\Gamma_{\text{bub}} = H^4$, and $\beta$ acts as the coefficient of the first order Taylor approximation of the full nucleation rate. One may wonder, however, whether the second and possibly higher derivatives will effectively play a role in setting $f_{\text{PBH}}$, because the probability of a Hubble patch collapse is very sensitive to the precise behaviour of $\Gamma_{\text{bub}}$. Nevertheless, our results below show the approximate form yields accurate predictions for the required $\beta$, at least in close-to-conformal potentials. Thus confirming the applicability of the GW signal based on the requirements on $\beta$ for achieving $f_{\text{PBH}} = 1$, as outlined in our discussion around eqs. (2.1)–(2.5).

Let us now also mention some work which previously headed along this path, but with different overall scope, e.g. [47, 48, 53, 89, 90].<fnref n="2" /> Indeed, already in 1984, Witten was discussing the concomitant aspects of a strong QCD phase transition, GWs, production of PBHs, and future pulsar timing array measurements [64]. Compared to the aforementioned papers, we provide a thorough scan over the expected PBH mass range coming from the PT, giving a more detailed picture of the expected GW signal. We also further develop the late nucleating patch formalism [47, 56], showing the calculation can be recast from using time as a variable to a conveniently chosen temperature, here that of the false vacuum plasma.

<footnote n="2">²As this paper was in the final stages of preparation, an analytic approach to estimating the PT parameters, GW and PBH production, partly in the context of the $B-L$ model, was presented in [57].</footnote>
Our focus here is on $f_{\text{PBH}} = 1$ which results in a scale of new physics beyond collider reach. The interplay of micro-lensing hints for $0 < f_{\text{PBH}} < 1$ at larger PBH masses, GW signals from supercooled PTs, and collider physics has instead recently been studied in [91].

# 3 The model

## 3.1 Particles and charges

As our example, we take a classically conformal $B-L$ extension of the SM, in which we add three right handed neutrinos, $N_i$ with $Q_{B-L} = -1$, to cancel off anomalies [86, 92, 93]. To break the $B-L$ symmetry, we add a complex scalar, $\rho$, with $Q_{B-L} = -2$. Ultimately, the radiative breaking of the $B-L$ symmetry induces a negative mass term for the electroweak Higgs field via a portal coupling $\lambda_{\rho h}$, thus eventually breaking electroweak symmetry. The kinetic term for the new scalar reads

$$ \mathcal{L} \supset (D_{\mu} \rho)^* (D^{\mu} \rho), \qquad (3.1) $$

where $D_{\mu} \equiv \partial_{\mu} - 2ig_{B-L}Z'_{\mu}$ is the covariant derivative with gauge coupling $g_{B-L}$. This will generate a mass for the gauge boson once $\rho$ gains a vacuum expectation value (vev) radiatively, $\langle \rho \rangle = v_{\rho}/\sqrt{2}$. At the classical level, the scalar potential reads

$$ V(H, \rho) = \lambda_{\rho} |\rho|^4 + \lambda_{\rho h} |\rho|^2 |H|^2 + \lambda_h |H|^4, \qquad (3.2) $$

where $H$ is the SM Higgs doublet. In addition to the SM Yukawas we also have

$$ \mathcal{L} \supset y_{\nu ij} \overline{l}_{Li} \tilde{H} N_j + \frac{1}{2} y_{Ni} \rho \overline{N} N^c + \text{H.c.}, \qquad (3.3) $$

where $l_{Li}$ are the SM lepton doublets and $\tilde{H} \equiv i\sigma_2 H^*$. The second term in eq. (3.3) will give Majorana masses to the $N_i$. These can then decay in a CP violating manner via the first Yukawa coupling realising leptogenesis, on which we will briefly comment later. The neutrino oscillation data [94–101] is explained through a type-I seesaw [102–105] (here at relatively low scales $\lesssim 10^8$ GeV). In the regime we are interested in, there is a significant hierarchy between $v_{\rho}$ and the electroweak vev $v_{\text{ew}} \simeq 246$ GeV, which implies a small cross quartic coupling,

$$ \lambda_{\rho h} \simeq - \left( \frac{m_h}{v_{\rho}} \right)^2, \qquad (3.4) $$

where $m_h \simeq 125$ GeV is the SM-like Higgs mass. This hierarchy situates the scenario beyond the reach of current collider constraints.

## 3.2 Radiative symmetry breaking at finite temperature

Our treatment of the field theory in this section, leading to the approximation of the effective potential, $V_{\text{eff}}$, follows [87]. At one loop the beta function of the $B-L$ Higgs quartic is

$$ \beta_{\lambda_{\rho}} \equiv \frac{d\lambda_{\rho}}{d\log \mu} = \frac{1}{(4\pi)^2} (96g_{B-L}^4 - y_{Ni}^2 + 2\lambda_{\rho h}^2 + 20\lambda_{\rho}^2 + \lambda_{\rho}[2y_{Ni}^2 - 48g_{B-L}^2]), \qquad (3.5) $$
where the sum over the $y_{Ni}$ is understood ($i = 1, 2, 3$). Consistency requires $\beta_{\lambda_{\rho}} > 0$, so that $\lambda_{\rho}$ turns negative when running down from high scales, triggering radiative symmetry breaking. Making the here justifiable approximation of the symmetry breaking as dominantly in a single field direction, following the approach of Gildener and Weinberg [106], we write the one loop zero temperature potential as

$$V_0(\rho) = \beta_{\lambda_{\rho}} \frac{\rho^4}{4} \left( \log \left[ \frac{\rho}{v_{\rho}} \right] - \frac{1}{4} \right), \qquad (3.6)$$

where from now we denote $v_{\rho}$ as the vev at the zero-temperature minimum and $\rho$ denotes the classical field value. At the $B-L$ symmetry breaking scale, $\lambda_{\rho} \simeq 0$, because the coupling is switching signs. We also require a small $\lambda_{\rho h}$, as mentioned above. Therefore we can approximate the beta function as

$$\beta_{\lambda_{\rho}} \approx \frac{1}{(4\pi)^2} (96g_{B-L}^4 - y_{Ni}^2). \qquad (3.7)$$

The mass of the physical scalar after symmetry breaking is $m_{\rho} \simeq \sqrt{\beta_{\lambda}} v_{\rho}$. Its field dependent mass gives only a small contribution to the thermal corrections in the effective potential. The dominant thermal corrections instead come from the field dependent mass of the $B-L$ gauge boson,

$$M_{Z'} = 2g_{B-L}\rho. \qquad (3.8)$$

and the heavy Majorana neutrinos,

$$M_{Ni} = \frac{y_{Ni}\rho}{\sqrt{2}}. \qquad (3.9)$$

Their thermal contributions are

$$V_T(\rho, T) = \frac{T^4}{2\pi^2} \left( 3J_B \left[ \frac{M_{Z'}^2}{T^2} \right] + 2J_F \left[ \frac{M_{Ni}^2}{T^2} \right] \right), \qquad (3.10)$$

where the thermal functions are

$$J_{B/F}(x) = \pm \int_{0}^{\infty} dk \, k^2 \log[1 \mp e^{-\sqrt{k^2+x}}], \qquad (3.11)$$

for bosons and fermions respectively. To improve the perturbative analysis, we include the corrections due to the resummation of daisy diagrams [107],

$$V_{\text{daisy}}(\rho, T) = \frac{T}{12\pi} (M_{Z'}^3 - [M_{Z'}^2 + \Pi_{Z'}]^{3/2}), \qquad (3.12)$$

where the thermal mass of the longitudinal component of the $B-L$ gauge boson is

$$\Pi_{Z'} = 4g_{B-L}^2 T^2. \qquad (3.13)$$

The full effective potential we consider is therefore

$$V_{\text{eff}}(\rho, T) = V_0(\rho) + V_T(\rho, T) + V_{\text{daisy}}(\rho, T), \qquad (3.14)$$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The evolution with temperature of the effective potential for a benchmark point with $g_{B-L} = 1$, $y_{Ni} = 1$ and $v_{\rho} = 100$ TeV. The barrier separating the two minima means the evolution from the false to the true vacuum corresponds to a first order phase transition. The barrier is always present when the finite temperature corrections are added to the radiatively induced symmetry breaking potential.</figcaption></figure>

after the aforementioned approximations have been made. To be consistent with observations, showing a small vacuum energy in the present day, a cosmological constant term,

$$ \Lambda_{\text{vac}} = \frac{\beta_{\lambda_{\rho}} v_{\rho}^4}{16} \qquad (3.15) $$

needs to be added to the potential, so that $V_0(v_\rho) = 0$. The gravitationally important terms, $\Lambda_{\text{vac}}$, together with the Planck mass, $M_{\text{Pl}}$, must be considered external to our classically scale invariant sector. An example showing the temperature evolution of the potential is shown in figure 1. Nucleation is possible once the false and true minima become degenerate at the critical temperature (time), $T_c$ ($t_c$).

### 3.3 Bubble nucleation and percolation

The bubble nucleation rate per unit volume can be approximated as [108–110]

$$ \Gamma_{\text{bub}} = A^4 e^{-S}, \qquad (3.16) $$

where $T$ is the temperature, $A$ is a mass dimension one pre-factor we elucidate later, and $S$ is the Euclidean bounce action. At zero temperature the configuration minimizing the action, $S \equiv S_4$, is $O(4)$ symmetric and

$$ S_4 = 2\pi^2 \int_0^\infty dr \, r^3 \left( \frac{1}{2} \left[ \frac{d\rho}{dr} \right]^2 + V_{\text{eff}}(\rho) \right) \qquad (3.17) $$

In this work we can limit ourselves PT dominated by single field directions. The action $S_4$ corresponds to quantum tunneling through the potential barrier. The equation of motion is

$$ \frac{d^2 \rho}{dr^2} + \frac{3}{r} \frac{d \rho}{dr} = \frac{d V_{\text{eff}}}{d \rho}, \qquad (3.18) $$

with boundary conditions

$$ \left. \frac{d\rho}{dr} \right|_{r=0} = 0, \quad \text{and} \quad \lim_{r \to \infty} \rho(r) = 0. \qquad (3.19) $$
At finite temperature the field instead becomes periodic in $1/T$ in the time coordinate. The configuration minimizing the action is $O(3)$ symmetric. Furthermore, at sufficiently high temperatures, the minimum action configuration becomes constant in the time direction and

$$S \equiv \frac{S_3}{T} = \frac{4\pi}{T} \int_0^\infty dr \, r^2 \left( \frac{1}{2} \left[ \frac{d\rho}{dr} \right]^2 + V_{\text{eff}}(\rho, T) \right). \qquad (3.20)$$

This represents bubble formation through classical field excitation over the barrier. The corresponding equation of motion is then

$$\frac{d^2\rho}{dr^2} + \frac{2}{r}\frac{d\rho}{dr} = \frac{dV_{\text{eff}}}{d\rho}, \qquad (3.21)$$

with the same boundary conditions as above. The solution with a non-trivial periodic bounce in the time coordinate, corresponding to quantum tunneling at finite temperature, is more difficult to evaluate in practical calculations. It is eventually, however, well approximated at low $T$ by $S_4$ but calculated with $V_{\text{eff}}$ including any finite temperature corrections. For simplicity, we therefore take a rather standard estimate for our nucleation rate

$$\Gamma_{\text{bub}} \approx \text{Max} \left[ \frac{1}{R_c^4} \left( \frac{S_4}{2\pi} \right)^2 e^{-S_4}, \quad T^4 \left( \frac{S_3}{2\pi T} \right)^{3/2} e^{-S_3/T} \right], \qquad (3.22)$$

where $R_c \sim 1/T$ is the bubble radius in the low $T$ limit of radiative symmetry breaking, and we have used standard estimates of the pre-factor (for discussion regarding uncertainties see [111]). For the PTs we study, we have verified that $\Gamma_{\text{bub}}$ is given by the $S_3$ action.

The nucleation temperature in an average Hubble patch, $T_n$, or the equivalent nucleation time, $t_n$, is defined as when the nucleation rate equals the Hubble rate,

$$\Gamma_{\text{bub}} = H^4. \qquad (3.23)$$

In practice, we compute the nucleation temperature using the Hubble rate in the false vacuum,

$$H \equiv H_{\text{false}} = \sqrt{\frac{8\pi}{3M_{\text{Pl}}^2} \left( \frac{g_* \pi^2}{30} T^4 + \Lambda_{\text{vac}} \right)}, \qquad (3.24)$$

where $T$ is the false vacuum temperature, and $g_*$ are the radiation degrees-of-freedom. In our scenario, $g_* \simeq 116$ prior to the PT, and $g_* \simeq 113$ just after the PT. Henceforth, any instances of the Hubble parameter and the scale factor presented without an index should be understood as these quantities evaluated in the false vacuum phase. We will be interested in strong supercooled phase transitions for which the universe becomes vacuum dominated. This occurs at a temperature

$$T_{\text{infl}} = \left( \frac{30\Lambda_{\text{vac}}}{g_*(T_{\text{infl}})\pi^2} \right)^{1/4}, \qquad (3.25)$$

For the case we will be interested in, namely efficient reheating and negligible changes in degrees-of-freedom between the two phases, we have $T_{\text{RH}} \simeq \text{Max}[T_{\text{infl}}, T_p]$, where $T_{\text{RH}}$ is the reheating temperature following the completion of the phase transition, and $T_p$ is the bubble percolation temperature, which we now review how to calculate.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Left: temperatures characterising the phase transition as a function of the gauge coupling. We display the following temperatures: critical, $T_c$, nucleation, $T_n$, percolation, $T_p$, and the temperature at which thermal inflation starts, $T_{\text{infl}}$. Right: the dimensionless phase transition latent heat, $\alpha$, and inverse timescale normalised to the Hubble rate, $\beta_H$.</figcaption></figure>

In the vacuum dominated regime, the temperature (time) at which the bubbles percolate, $T_p$ ($t_p$), may be appreciably different from $T_n$ ($t_n$). It is useful to introduce the comoving radius of a bubble at false vacuum temperature $T$, nucleated at some higher temperature $T'$,

$$r(t, t') = \int_{t'}^{t} \frac{d\tilde{t}}{a(\tilde{t})} = r(T, T') = \int_{T}^{T'} \frac{d\tilde{T}}{\tilde{T}H(\tilde{T})a(\tilde{T})}, \quad (3.26)$$

where $a$ is the scale factor, we have assumed $v_w \simeq 1$ shortly after nucleation, and have used $dT = -THdt$. The physical radius is then

$$R(T, T') = a(T)r(T, T') = a(T) \int_{T}^{T'} \frac{d\tilde{T}}{\tilde{T}H(\tilde{T})a(\tilde{T})} = \frac{1}{T} \int_{T}^{T'} \frac{d\tilde{T}}{H(\tilde{T})}, \quad (3.27)$$

where we have assumed a constant $g_*$ and so used $a(T)T = a(\tilde{T})\tilde{T}$. The expected volume of true-vacuum bubbles per comoving volume (double counting overlapping regions and including fictitious nucleations in true vacuum) is given by [112–116]

$$I = \frac{4\pi}{3} \int_{t_c}^{t} dt' \Gamma_{\text{bub}}(t') a^3(t') r^3(t, t') = \frac{4\pi}{3} \int_{T}^{T_c} dT' \frac{\Gamma_{\text{bub}}(T')}{T'^4 H(T')} \left( \int_{T}^{T'} \frac{dT}{H(\tilde{T})} \right)^3. \quad (3.28)$$

To find the probability of a given point in the comoving volume to be in the false vacuum, we need to exclude fictitious nucleations and avoid double counting overlapped regions, and the appropriate expression is given by [112–116]

$$P = e^{-I}. \quad (3.29)$$

The change in the physical false vacuum volume, $V_{\text{false}} = a^3 P(T)$, normalized to the Hubble rate, is captured by the equation [116]

$$\frac{1}{H V_{\text{false}}} \frac{dV_{\text{false}}}{dt} = 3 + T \frac{dI}{dT}. \quad (3.30)$$

We define the percolation temperature $T_p$ as the highest temperature for which both $I(T) > 1$ and $d \log V_{\text{false}}/dt < -H$ hold (sometimes slightly less stringent conditions are employed,
e.g. see discussion in [117], the precise choice does not affect our qualitative results below). Numerically, we find the first condition, $I(T) > 1$, occurs later and hence determines the percolation condition in our parameter space.

To characterise the strength of the phase transition, we find the free energy difference between the true and false vacua normalized to the radiation density,

$$ \alpha = \frac{1}{\rho_{\text{rad}}} \left( 1 - T \frac{\partial}{\partial T} \right) \left. \left( V_{\text{eff}}(0, T) - V_{\text{eff}}(\phi_n, T) \right) \right|_{T=T_*} , \quad (3.31) $$

where $T_*$ is set to either the nucleation or percolation temperature.<fnref n="3" /> The inverse timescale of the phase transition normalized to the Hubble rate reads

$$ \beta_H \equiv \left. \frac{\beta}{H} \right|_{t=t_*} = \frac{1}{H \Gamma_{\text{bub}}} \left. \frac{d\Gamma_{\text{bub}}}{dt} \right|_{t=t_*} = - \left. \frac{T}{\Gamma_{\text{bub}}} \frac{d\Gamma_{\text{bub}}}{dT} \right|_{T=T_*} . \quad (3.32) $$

By numerically calculating the bubble action $S$ for a sufficiently large number of points over a suitable temperature range, we obtain a good approximation of $\Gamma_{\text{bub}}(T)$, which allows us to find all the above quantities. Phase transition parameters as a function of $g_{B-L}$ are shown in figure 2. Clearly in the limit of a small $\beta_{\lambda\rho}$ we have strong supercooling, with $\beta_H \lesssim 8$ and $\alpha \gg 1$, which allows for cosmologically significant PBH production, e.g. as found in [56]. We now turn to calculating the PBH abundance.

# 4 Primordial black hole production

## 4.1 Background evolution

We first consider an average Hubble patch. The total energy density, $\rho$, consist of the vacuum energy density, $\rho_{\text{vac}}$, and the radiation, $\rho_{\text{rad}}$. The latter is made up of the plasma of relativistic particles together with the bubble walls moving at $v_w \simeq 1$ [56]. The Friedmann equations are given by

$$ H_{\text{bkg}}^2 = \frac{8\pi}{3} \frac{\rho_{\text{vac}} + \rho_{\text{rad}}}{M_{\text{Pl}}^2}, \quad (4.1a) $$

$$ \frac{d\rho_{\text{rad}}}{dt} = -4H_{\text{bkg}}\rho_{\text{rad}} - \frac{d\rho_{\text{vac}}}{dt}, \quad (4.1b) $$

where $H_{\text{bkg}}$ denotes the Hubble rate of the average background patch. We can relate $\rho_{\text{vac}}$ to the bubble nucleation rate through

$$ \rho_{\text{vac}} = \Lambda_{\text{vac}}e^{-I}. \quad (4.2) $$

To solve the Friedmann equations, we find it useful to switch coordinates, from time to false vacuum temperature. We then have

$$ H_{\text{bkg}}^2 = \frac{8\pi}{3} \frac{\rho_{\text{vac}} + \rho_{\text{rad}}}{M_{\text{Pl}}^2}, \quad (4.3a) $$

$$ \frac{d\rho_{\text{rad}}}{dT} = \frac{4H_{\text{bkg}}\rho_{\text{rad}}}{HT} - \frac{d\rho_{\text{vac}}}{dT}, \quad (4.3b) $$

<footnote n="3">³When eventually calculating the entropy dilution factor below, evaluating $\alpha$ at $T_* = T_p$ gives a better estimate, due to the additional expansion between $T_n$ and $T_p$. The choice is irrelevant when determining the GW signal for the PTs of interest as $\alpha \gg 1$. On the other hand, when approximating the mean bubble radius using $\beta_H$, we numerically find evaluating $\beta_H$ at $T_* = T_n$ gives a more accurate result, after comparing to the mean radius extracted from the bubble distribution itself.</footnote>
where $T$ is the temperature in the false vacuum, and which we solve to find $\rho_{\text{rad}}(T)$. Note $\rho_{\text{vac}}$, through the quantity $I$ which contains the bubble radius, effectively also contains the scale factor. Hence the equations should be solved self-consistently by evaluating $I$ by taking into account the scale factor at each time step, which in turn can be found through $H_{\text{bkg}}$. In order to solve the equations, however, a simplification is possible. First note that before percolation, $H_{\text{bkg}}$ is for practical purposes given by its corresponding false vacuum value. After percolation, the $\rho_{\text{vac}}$ source term for the radiation is, by definition, approximately negligible. We therefore first evaluate $I$ numerically using the scale factors and Hubble rate in the false vacuum, $H$, and use this to solve the Friedmann equations including $H_{\text{bkg}}$ of the background patch in the first equation. We then check whether the approximation is a good one by using the resulting scale factor in an updated determination of $I$. The resulting shift in $T_p$ is tiny — we will quantify it below — and so we confirm the approximation is justified.

For our calculations, we will also need the scale factor as a function of the false vacuum temperature. Ignoring changes in degrees-of-freedom, for the scale factor in a false vacuum patch we have the standard relation,

$$ \frac{a(T)}{a(T')} = \frac{T'}{T}. \qquad (4.4) $$

For the scale factor in the background patch that has nucleated, we instead solve

$$ \frac{da_{\text{bkg}}}{dT} = -\frac{a_{\text{bkg}}}{T} \frac{H_{\text{bkg}}}{H}. \qquad (4.5) $$

The equation is solved starting from a sufficiently high temperature where $H_{\text{bkg}} = H$ and $a_{\text{bkg}}(T) = a(T)$. Note, for bookkeeping purposes, the temperatures which appear here should be understood as the temperature in the false vacuum patch, not the temperature in the nucleating background patch, which contains nucleated bubbles and partially reheated plasma. Using our solution, it is convenient to define

$$ A(T, T') \equiv \frac{a_{\text{bkg}}(T)T}{a_{\text{bkg}}(T')T'}, \qquad (4.6) $$

which changes from unity to take into account deviations from the false vacuum relation, eq. (4.4), for the background scale factor.

In terms of the false vacuum temperatures and false vacuum $H$, the bubble radius in the background patch is then more precisely given by

$$ R_{\text{bkg}}(T, T') = \frac{1}{T} \int_{T}^{T'} \frac{A(T, \tilde{T}) d\tilde{T}}{H(\tilde{T})}. \qquad (4.7) $$

It can also be of interest to consider the bubble density as a function of radius [116, 117], which rewritten as a function of the false vacuum temperature is given by

$$ \frac{dn_{\text{bkg}}}{dR_{\text{bkg}}}(T, R_{\text{bkg}}) = \frac{dT'}{dR_{\text{bkg}}} \frac{T^3 \Gamma_{\text{bub}}(T') P(T')}{A(T, T')^3 T'^4 H(T')}. \qquad (4.8) $$

where here $T' \equiv T'(R)$ is understood as being the false vacuum temperature at which a bubble of physical size $R_{\text{bkg}}$ was nucleated, found by numerically inverting eq. (4.7), from which one also finds the corresponding derivative.
## 4.2 Late patch evolution

The above describes the evolution of an average background nucleating patch. We are now interested in the evolution of a late nucleating patch. We assume no bubble nucleates in the patch until some temperature $T_i < T_n$. By modifying the high temperature terminal in eq. (3.28), we find the vacuum energy in the late patch,

$$I_{\text{late}} = \frac{4\pi}{3} \int_{T}^{T_{i}} dT' \frac{\Gamma_{\text{bub}}(T')}{T'^{4} H(T')} \left( \int_{T}^{T'} \frac{dT}{H(\tilde{T})} \right)^{3}. \quad (4.9)$$

Again we use the same approximation as before, in assuming $H = H_{\text{false}}$, when calculating $I_{\text{late}}$. The probability of finding a point in the false vacuum in the late patch is then

$$P_{\text{late}} = e^{-I_{\text{late}}}. \quad (4.10)$$

Thus the vacuum energy density in the late patch is simply

$$\rho_{\text{vaclate}} = \Lambda_{\text{vac}} e^{-I_{\text{late}}}. \quad (4.11)$$

We then also solve the Friedmann equations for the late patch

$$H_{\text{late}}^{2} = \frac{8\pi}{3} \frac{\rho_{\text{vaclate}} + \rho_{\text{radlate}}}{M_{\text{Pl}}^{2}}, \quad (4.12a)$$

$$\frac{d\rho_{\text{radlate}}}{dT} = \frac{4H_{\text{late}}\rho_{\text{radlate}}}{HT} - \frac{d\rho_{\text{vaclate}}}{dT}. \quad (4.12b)$$

Eventually we will also need the scale factor in the late patch, found by solving

$$\frac{da_{\text{late}}}{dT} = -\frac{a_{\text{late}}}{T} \frac{H_{\text{late}}}{H}. \quad (4.13)$$

Similarly to what we did for the background patch, it is convenient to define

$$B(T, T') \equiv \frac{a_{\text{late}}(T)T}{a_{\text{late}}(T')T'}, \quad (4.14)$$

to parametrize deviations from the false vacuum scale factor/temperature relation. The bubble size in the late patch is then given by

$$R_{\text{late}}(T, T') = \frac{1}{T} \int_{T}^{T'} \frac{B(T, \tilde{T}) d\tilde{T}}{H(\tilde{T})}. \quad (4.15)$$

The bubble spectrum of the late patch is given by

$$\frac{dn_{\text{late}}}{dR_{\text{late}}} (T, R_{\text{late}}) = \frac{dT'}{dR_{\text{late}}} \frac{T^3 \Gamma_{\text{bub}}(T') P_{\text{late}}(T')}{B(T, T')^3 T'^4 H(T')} \Theta(T_i - T'). \quad (4.16)$$

where the $\Theta$ function takes into account that no bubbles formed at $T' > T_i$ in the late patch.
## 4.3 Collapse and fractional abundance

For each choice of $T_i$ we can also calculate the contrast density in radiation

$$ \delta(T) = \frac{\rho_{\text{radlate}} - \rho_{\text{rad}}}{\rho_{\text{rad}}}. \qquad (4.17) $$

This reaches a maximum shortly after late patch percolation, as the energy density in the background patch has began to become diluted a little earlier, while the late patch energy density is still constant due to the vacuum. After late patch percolation, $\delta$ decreases again as $H_{\text{late}} > H_{\text{bkg}}$ leading to a faster redshifting. We define $T_{\delta\text{max}}$ as the temperature at which $\delta$ is maximized, $\delta_{\text{max}} \equiv \delta(T_{\delta\text{max}})$. The smaller $T_i$ the larger $\delta_{\text{max}}$.

We assume a patch collapses if $\delta_{\text{max}}$ reaches a threshold value, $\delta_c$, for critical collapse. Calculations, some based on full general relativistic simulations, indicate $0.4 \le \delta_c \le 0.66$ in the context of overdensities from inflation re-entering the Hubble horizon [23–27, 27–33]. The precise value of $\delta_c$ depends on the shape of the overdensity, but its effects can be captured by considering the maximum of the so-called compaction function and its curvature, and thus be related to the power spectrum of curvature perturbations. Note the calculations have been performed assuming spherical symmetry and typically also isotropic pressure (although see [118]). We expect departures from these assumptions in the PT scenario, due to the nature of the Hubble patch just after bubble percolation. Note non-sphericity is expected to increase $\delta_c$ [119–123]. The value of $\delta_c$ will have an effect on which input parameters return $f_{\text{PBH}} = 1$, but the strong GW signal will not be sensitive to our precise choice. In this work, we follow previous PT literature [56], and take $\delta_c = 0.45$ to aide comparison. The late patch nucleation temperature, $T_i$, is thus from now on fixed by requiring the maximal density contrast, $\delta_{\text{max}}$, be equal to the critical collapse threshold, $\delta_c$.

We assume the collapse occurs at $T_{\delta\text{max}}$ with corresponding time $t_{\delta\text{max}}$. The probability of the horizon size patch, which eventually collapses into a PBH, having no bubbles at $T_i$ is given by [47]

$$ P_{\text{no bub}} = \text{Exp} \left[ - \int_{t_c}^{t_i} dt' \Gamma_{\text{bub}}(t') a_{\text{late}}(t')^3 V_{\text{no bub}} \right] \qquad (4.18a) $$

$$ = \text{Exp} \left[ - \int_{T_i}^{T_c} \frac{dT' \Gamma_{\text{bub}}(T')}{T' H(T')} a_{\text{late}}(T')^3 V_{\text{no bub}} \right], \qquad (4.18b) $$

where the temperature and Hubble rate are those of the false vacuum, and the volume factor,

$$ V_{\text{no bub}} = \frac{4\pi}{3} \left[ \frac{1}{a_{\text{late}}(t_{\delta\text{max}})H_{\text{late}}(t_{\delta\text{max}})} \right]^3, \qquad (4.19) $$

represents the comoving volume of the Hubble sized patch at the start of the collapse. However, the calculation of $\delta$ for the critical collapse also assumes no larger bubbles — from background patches surrounding the late nucleating patch — enter into the collapsing volume before $\delta_{\text{max}}$ is attained. Thus, it has been advocated [56], that the collapse probability is better estimated as

$$ P_{\text{coll}} = \text{Exp} \left[ - \int_{t_c}^{t_i} dt' \Gamma_{\text{bub}}(t') a_{\text{late}}(t')^3 V_{\text{coll}} \right] \qquad (4.20a) $$

$$ = \text{Exp} \left[ - \int_{T_i}^{T_c} \frac{dT' \Gamma_{\text{bub}}(T')}{T' H(T')} a_{\text{late}}(T')^3 V_{\text{coll}} \right], \qquad (4.20b) $$
where the volume factor is

$$V_{\text{coll}} = \frac{4\pi}{3} \left[ \frac{1}{a_{\text{late}}(t_{\text{delta}})H_{\text{late}}(t_{\text{delta}})} + r(t_{\text{delta}}, t')} \right]^3 \quad (4.21a)$$

$$= \frac{4\pi}{3} \left[ \frac{1}{a_{\text{late}}(T_{\text{delta}})H_{\text{late}}(T_{\text{delta}})} + \int_{T_{\text{delta}}}^{T'} \frac{d\tilde{T}}{\tilde{T}H(\tilde{T})a_{\text{bkg}}(\tilde{T})} \right]^3 . \quad (4.21b)$$

Here, in the evaluation of $r(t_{\delta\text{max}}, t')$, eq. (3.26), we use the background value for the scale factor. The eventual PBH mass is estimated as the energy inside the sound horizon of the collapsing patch at $t_{\delta\text{max}}$ (e.g. see [124])

$$M_{\text{PBH}} = \frac{4\pi(\rho_{\text{radlate}} + \rho_{\text{vaclate}})c_s^3}{3H_{\text{late}}^3} = \frac{c_s^3 M_{\text{Pl}}^2}{2H_{\text{late}}}, \quad (4.22)$$

where $c_s = 1/\sqrt{3}$ is the sound speed (although included, $\rho_{\text{vaclate}}$ is negligible at $t_{\delta\text{max}}$). The current study is limited to the above monochromatic estimate, we leave for future work the full derivation of a detailed mass spectrum taking into account the full critical collapse phenomenon [125–129]. Some preliminary estimates, showing the monochromatic estimate is a good one, are given in appendix A.

We now wish to find, $f_{\text{PBH}}$, the PBH-to-DM fraction at late times. To do this we note that PBH formation occurs a little after bubble percolation, around time $t_{\delta\text{max}}$ when the false vacuum temperature is $T_{\delta\text{max}}$, and we denote the corresponding background temperature as $T_{\text{bkg form}}$ (note $T_{\text{bkg form}} \approx T_{\text{RH}} \approx T_{\text{infl}}$ with small differences found numerically and taken into account in our results). The ratio of densities at this time is

$$\frac{\rho_{\text{PBH}}(t_{\delta\text{max}})}{\rho_{\text{rad}}(t_{\delta\text{max}})} \approx c_s^3 P_{\text{coll}} \frac{H_{\text{bkg}}(t_{\delta\text{max}})}{H_{\text{late}}(t_{\delta\text{max}})}, \quad (4.23)$$

where the ratio of Hubble rates is a correction which takes into account that $M_{\text{PBH}}$ is set by the late rather than background patch density. Later, at matter radiation equality, $T_{\text{eq}} \simeq 0.8$ eV, we have $\rho_{\text{rad}} = \rho_{\text{m}} = \rho_{\text{DM}}(\Omega_{\text{DM}} + \Omega_{\text{B}})/\Omega_{\text{DM}} \simeq 1.2\rho_{\text{DM}}$. Between PBH formation, at $T \simeq T_{\text{bkg form}}$, and $T_{\text{eq}}$, one has $\propto 1/a^3$ dilution so

$$\frac{\rho_{\text{PBH}}(T_{\text{eq}})}{\rho_{\text{PBH}}(T_{\text{bkg form}})} = \left( \frac{a_{\text{bkg}}(T_{\text{bkg form}})}{a_{\text{bkg}}(T_{\text{eq}})} \right)^3 . \quad (4.24)$$

From entropy conservation we have

$$\frac{a_{\text{bkg}}(T_{\text{bkg form}})}{a_{\text{bkg}}(T_{\text{eq}})} = \left( \frac{g_{*s}(T_{\text{eq}})}{g_{*s}(T_{\text{bkg form}})} \right)^{1/3} \frac{T_{\text{eq}}}{T_{\text{bkg form}}}, \quad (4.25)$$

where $g_{*s}$ are the entropic degrees-of-freedom. The radiation instead redshifts as $\propto 1/a^4$, but it also receives reheating contributions, so

$$\frac{\rho_{\text{rad}}(T_{\text{eq}})}{\rho_{\text{rad}}(T_{\text{bkg form}})} = \frac{g_{*}(T_{\text{eq}})T_{\text{eq}}^4}{g_{*}(T_{\text{bkg form}})T_{\text{bkg form}}^4} \quad (4.26)$$

Combing all the above we have

$$f_{\text{PBH}} = \frac{\rho_{\text{PBH}}}{\rho_{\text{DM}}} \approx c_s^3 P_{\text{coll}} \frac{\Omega_{\text{DM}} + \Omega_{\text{B}}}{\Omega_{\text{DM}}} \frac{g_{*s}(T_{\text{eq}})}{g_{*}(T_{\text{eq}})} \frac{T_{\text{bkg form}}}{T_{\text{eq}}} \frac{H_{\text{bkg}}(t_{\delta\text{max}})}{H_{\text{late}}(t_{\delta\text{max}})}, \quad (4.27)$$

where we have used $g_{*}(T_{\text{bkg form}}) = g_{*s}(T_{\text{bkg form}})$, and we remind the reader that $g_{*s}(T_{\text{eq}}) = 3.91$ and $g_{*}(T_{\text{eq}}) = 3.36$.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Characterisation of several important quantities for PBH formation and the phase transition dynamics for a benchmark scenario with $v_ρ = 10^3$ TeV, $y_{Ni} = 0.01$, and $g_{B-L} = 0.2958$ that predict $f_{PBH} = 1$, with $M_{PBH} \simeq 4.5 \times 10^{-12} M_⊙$. We depict the behaviour of the nucleation rate comparing the full calculation of the nucleation rate with the approximation often used in the literature [47, 56] (top left), energy densities (top right) and density contrast (middle left, with $T_i \simeq 53.6$ GeV) as a function of the false vacuum temperature $T$. We also show the density contrast as a function of $T_i$ (middle right). Finally we display the bubble distribution at percolation (bottom left), where RH is the bubble size as a fraction of the Hubble length. Weighting the distribution by $R^3$ we see the differences in the large bubbles between the background and late patches (bottom right).</figcaption></figure>
## 4.4 Application to example model

Having developed the above machinery, we now calculate the PBH fraction for our example $B-L$ model, and search for parameter space in which $f_{\text{PBH}} = 1$. In figure 3, we characterise PBH formation and the strongly first-order phase transition for a benchmark point for which $v_\rho = 10^3$ TeV, $y_{Ni} = 0.01$, and $g_{B-L} = 0.2958$. These values were judiciously chosen so as to return $f_{\text{PBH}} = 1$. The PBH mass is $M_{\text{PBH}} \simeq 4.5 \times 10^{-12} M_\odot$. In the top left plot we show our full determination of $\Gamma_{\text{bub}}(T)$ and compare it with the approximate form, eq. (2.6), rewritten in terms of temperature

$$ \Gamma_{\text{bub}}(T) \approx H(T_n)^4 \left(\frac{T_n}{T}\right)^{\beta_H(T_n)}. \qquad (4.28) $$

Deviations between the full and approximate expression are visibly present at $T = T_i \approx 54$ GeV which returns the required contrast for PBH formation. We will eventually show the resulting modest differences in the $\beta_H$ required to achieve $f_{\text{PBH}} = 1$ below.

In figure 3 top right we first show the evolution of energy densities with the false vacuum temperature $T$. The dashed lines in the plot represent the vacuum, radiation, and total energy densities of an average Hubble patch, while the solid lines illustrate these same quantities for a late-nucleating patch. At sufficiently high false vacuum temperatures, the radiation energy densities for both types of Hubble patches align, and the same is observed with the vacuum energy densities. As nucleation starts within an average Hubble patch, its radiation energy density begins to increase relative to that of the late-nucleating patch. Once the percolation threshold of the true vacuum phase is achieved, there is a drastic drop in the vacuum energy density, leading the average background patch to become dominated by radiation. Concurrently, the energy density within this patch begins to dilute due to expansion, whilst the late-nucleating patch remains vacuum dominated. Once the percolation threshold is reached within the late-nucleating patch, an excess of radiation energy density relative to the surrounding background Hubble patches follows, triggering the collapse into a PBH.

In the middle row of figure 3, we show the dependence of $\delta_{\text{max}}$ on $T_i$ (right) and the dependence of the density contrast with the false vacuum temperature $T$ (left). We observe how we start off from a homogeneous state at sufficiently high false vacuum temperatures. The density contrast drops to its minimum value approximately when bubbles start nucleating in the background Hubble patches and increases again when nucleation starts in the delayed patch. Ultimately, it reaches a maximum shortly after late patch percolation, as the energy density in the background patch has began to become diluted a little earlier, while the late patch energy density is still constant due to the vacuum.

Finally, in the bottom row we show the bubble distribution of the average background patch and the late patch with $\delta_{\text{max}} = 0.45$ at their respective percolation temperatures. The late patch features approximately twice the number of bubbles as the background patch. Due to the first bubbles being delayed, the distribution of bubbles in the late patch features a cut in the distribution, here at roughly half the Hubble length. In the right plot, we weight the distribution by $R^3$ to show the volume occupied by the bubble. Therein we see a significant volume is occupied by rare large bubbles, especially for the background patch, which explains why the late patch is filled instead by roughly twice as many small bubbles.
Now we turn to figure 4. For a given choice of $v_ρ$, we can scan over $g_{B-L}$ to find where $f_{PBH} = 1$ (we always set $y_i = 0.01$ as an example). The required $g_{B-L}$ to achieve $f_{PBH} = 1$ is shown in figure 4 (top left) together with the corresponding value of the $β_{λρ}$ function (top right).

In the middle panel of figure 4 we show various temperatures characterizing the PT. The percolation temperature and also $T_{δmax}$ remains above QCD confinement, so we do not need to take into account QCD enhancement of $Γ_{bub}$, present at lower temperatures [130–133]. A smaller $β_{λρ}$ corresponds to a smaller vacuum energy difference between the two phases, which strengthens the phase transition [134]. Furthermore, smaller values of $g_{B-L}$ imply a weaker running of couplings which, in turn, delays the transition time and makes the derivative of the bubble action $S$ smaller [86].

As mentioned above, before solving the Friedmann equations, we first fix $H = H_{false}$ in our calculation of $I(T)$. To check this is a good approximation, we then iterate, using our $H_{bkg}$ determination to calculate a new $I(T)$, and use this to find an updated percolation temperature $T_{pnew}$. The ratio of the two is displayed in figure 4 middle right, showing only very small changes in $T_p$, thus justifying the approximation. In figure 4 bottom we show the resulting $M_{PBH}$ and various measures of the bubble radius at collision.

In the bottom row on the left, we show how the PBH mass decreases as we go to larger vevs, which corresponds to smaller physical Hubble volumes during the PT. In the right panel, we show the measures of the bubble radius at percolation, normalized to the Hubble length (mean, approximated mean, and radius corresponding to the mode of the bubble volume distribution). As the inverse timescale of the transition does not need to vary much from $β_H ≈ 8$ in order to give $f_{PBH} = 1$ over the entire allowed window, the typical bubble sizes when normalized to Hubble are also always similar.

# 5 Gravitational waves

## 5.1 The spectra

We next turn to the spectra we employ for the gravitational wave signal defined as

$$h^2\Omega_{GW}(f) = \frac{h^2}{\rho_c} \frac{d\rho_{GW}}{d\log f}, \qquad (5.1)$$

where $ρ_{GW}$ is the energy density in GWs and $ρ_c$ is the critical density.

As spectra from dedicated studies are often given at production rather than today, we briefly review how to take into account the redshifting. We note that GWs redshift like radiation but are not reheated, and $Ω_{rad}(T_{RH}) = ρ_{rad}/ρ_c ≈ 1$ just after the PT. The above implies an amplitude today [66, 73]

$$h^2\Omega_{GW} = h^2\Omega_{rad}(T_0) \left( \frac{g_{*s}(T_0)}{g_{*s}(T_{RH})} \right)^{4/3} \frac{g_{*}(T_{RH})}{g_{*}(T_0)} \Omega_{GW*} \qquad (5.2a)$$

$$= \frac{g_{*s}(T_0)^{4/3} \pi^2 T_0^4}{30} \frac{8\pi}{3H_{100}^2 M_{Pl}^2} \frac{\Omega_{GW*}}{g_{*s}(T_{RH})^{1/3}} \qquad (5.2b)$$

$$= 7.64 \times 10^{-5} g_{*s}(T_{RH})^{-1/3} \Omega_{GW*}, \qquad (5.2c)$$
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Top: the required gauge coupling (left) and associated beta function (right) giving $f_{\text{PBH}} = 1$ as a function of the $B-L$ vev, $v_\rho$. Shaded regions are excluded by PBH constraints. Middle: various temperatures characterizing the PT as a function of $v_\rho$. Note the new percolation temperature $T_{\text{pnew}}$ following our iterative check shows only a per mille correction with respect to $T_p$ (right). Bottom: the resulting PBH mass (left) and the bubble radii at collision (right), namely the mean radius $R_{\text{mean}}$, the approximate radius, $R_{\text{approx}} = \pi^{1/3}/\beta(T_n)$ [115], and the radius at which the bubble energy (and volume) distribution is maximized, $R_{\text{max}}$ [117].</figcaption></figure>
where $\Omega_{GW*}$ is the GW signal just after the PT, $T_0 \simeq 0.23$ meV is the CMB temperature today, and $H_{100} \equiv 100$ km/s/Mpc. The frequency of the signal is redshifted as

$$f_0 = \left( \frac{g_{*s}(T_0)}{g_{*s}(T_{RH})} \right)^{1/3} \frac{T_0}{T_{RH}} f_* \qquad (5.3)$$

where $f_*$ is the frequency just after the PT. Concretely, for the reheating temperature, we take

$$T_{\text{RH}} = \left( \frac{90 M_{\text{Pl}}^2 H_{\text{bkg}} (t_p)^2}{g_*(T_{\text{RH}}) 8 \pi^3} \right)^{1/4}, \qquad (5.4)$$

which follows directly from the first Friedmann equation. We are interested in supercooled PTs with highly relativistic bubble walls. We will use three determinations of the spectrum, showing our conclusions are hardly affected by the precise choice.

* The first estimate comes from $(3+1)$ dimensional lattice simulations of thick wall bubble collisions by Cutting, Escartin, Hindmarsh, and Weir [75]. The spectrum is found to be

$$h^2 \Omega_{\text{thw}}(f) = 4.38 \times 10^{-9} \left( \frac{100}{g_*(T_{\text{RH}})} \right)^{1/3} \left( \frac{10}{\beta_H} \right)^2 S_{\text{thw}}(f), \qquad (5.5)$$

where the spectral shape is given by

$$S_{\text{thw}}(f) = \frac{2.902 \tilde{f}_{\text{thw}}^{2.16} f^{0.742}}{2.16 \tilde{f}_{\text{thw}}^{2.902} + 0.742 f^{2.902}}, \qquad (5.6)$$

and the peak frequency is

$$\tilde{f}_{\text{thw}} = 32 \text{ mHz} \left( \frac{g_*(T_{\text{RH}})}{100} \right)^{1/6} \left( \frac{\beta_H}{10} \right) \left( \frac{T_{\text{RH}}}{10^2 \text{ TeV}} \right). \qquad (5.7)$$

In the above, we have used the central values of the fit provided in [75] for the thickest walls of their simulations, and taken the bubble diameter to be $D_{\text{bub}}H \simeq (8\pi)^{1/3}\beta_H^{-1}$ [115] in converting the spectrum in [75] to be in terms of $\beta_H$. From figure 4 we see this diameter is in good agreement with the mean diameter extracted from the bubble distribution when using $\beta_H(T_n)$, which we subsequently use in calculating the GW spectrum. Note if we were to instead use $2R_{\text{max}}$, the diameter at which the bubble energy distribution is maximized, as advocated in [117], our GW signals would be stronger.

* The second estimate comes from the hybrid simulations by Vaskonen and Lewicki [77]. Therein, the anisotropic stress induced in a bubble collision is first determined in a $(1+1)$ dimensional simulation. This is then used as a source at points at which walls collide in a $(3+1)$ dimensional lattice simulation in the thin walled limit. The advantage here is that the lower dimensional simulation allows one to study the effect of non-trivial scalar gradients and associated gauge field production during the bubble collision. The simulations find differences in the spectra between the non-gauged [76] and gauged cases [77], of U(1) symmetry breaking. Of course we use the gauged example here, where the spectrum is

$$h^2 \Omega_{\text{hyb}}(f) = 5.93 \times 10^{-9} \left( \frac{100}{g_*(T_{\text{RH}})} \right)^{1/3} \left( \frac{10}{\beta_H} \right)^2 S_{\text{hyb}}(f), \qquad (5.8)$$
with the shape,

$$S_{\text{hyb}}(f) = \frac{695}{\left[ 2.41 \left( \frac{f}{\tilde{f}_{\text{hyb}}} \right)^{-0.557} + \left( \frac{f}{\tilde{f}_{\text{hyb}}} \right)^{0.574} \right]^{4.20}} \quad (5.9)$$

and peak frequency

$$\tilde{f}_{\text{hyb}} = 22 \text{ mHz} \left( \frac{g_*}{100} \right)^{1/6} \left( \frac{\beta_H}{10} \right) \left( \frac{T_{\text{RH}}}{10^2 \text{ TeV}} \right). \quad (5.10)$$

* Reassuringly also semi-analytic methods which avoid the need for running lattice simulations have been developed [72, 73]. We use the results of the bulk flow from Konstandin [73] with amplitude

$$h^2 \Omega_{\text{bulk}}(f) = 1.06 \times 10^{-8} \left( \frac{100}{g_*(T_{\text{RH}})} \right)^{1/3} \left( \frac{10}{\beta_H} \right)^2 S_{\text{bulk}}(f), \quad (5.11)$$

and spectral shape

$$S_{\text{bulk}}(f) = \frac{3 \tilde{f}_{\text{bulk}}^{2.1} f^{0.9}}{2.1 \tilde{f}_{\text{bulk}}^{3} + 0.9 f^{3}}, \quad (5.12)$$

with peak frequency at

$$\tilde{f}_{\text{bulk}} = 21 \text{ mHz} \left( \frac{g_*(T_{\text{RH}})}{100} \right)^{1/6} \left( \frac{\beta_H}{10} \right) \left( \frac{T_{\text{RH}}}{10^2 \text{ TeV}} \right). \quad (5.13)$$

Note all the above results are given in the limit $\alpha \gg 1$, $v_w \simeq 1$, and the spectral shape functions have been normalized to return unity at their respective peak frequencies. As the finite cosmological horizon is not taken into account in the above simulations, we also impose the correct scaling $\Omega_{\text{GW}} \propto f^3$, for super-horizon modes at the time of the PT [135–139]. These correspond to the frequencies at IR tail of the spectrum, below

$$f_{\text{horizon}}^{\text{PT}} = 2.6 \text{ mHz} \left( \frac{g_*(T_{\text{RH}})}{100} \right)^{1/6} \left( \frac{T_{\text{RH}}}{10^2 \text{ TeV}} \right), \quad (5.14)$$

as measured today. In figure 5 we show the resulting GW spectra for two example parameter points with $f_{\text{PBH}} = 1$. For practical purposes here, we see the GW spectra from the three estimates are all rather similar. They are all above foregrounds and detectable at some upcoming interferometer.

One word of caution is in order. The above estimates do not take into account the Hubble expansion during the transition itself. The issue has been studied in ref. [78], which extended the analysis to an expanding background, but only for PTs in a radiation dominated background and using the older envelope approximation. The results indicated a suppression of an order-of-magnitude in the signal for PTs with $\beta_H \approx 10$. Note even with such a suppression, our GW signals would still be easily detectable, although this effect adds considerable theory uncertainty. To better gain a handle on the GW signal, the bulk flow analysis or the simulations will also have to be modified to take into account expansion, in particular for a vacuum dominated background. Our analysis further motivates such efforts.
<figure><img src="image_7.png" /><figcaption>Figure 5. GW spectra for two PTs returning $f_{\text{PBH}} = 1$. For each of the PTs, we show three estimates of the spectra — given the macroscopic PT parameters — available in the literature. Also shown are power-law-integrated sensitivity curves (SNR = 10) for current and future interferometers together with astrophysical foregrounds. Other key parameters associated with these example PTs can be read off figure 4.</figcaption></figure>

## 5.2 Signal-to-noise ratio

We now wish to estimate the signal-to-noise ratio and also show foregrounds will not mask our signal. A similar method to what we describe in this section, up to some small updates, has also been used in [140]. The signal-to-noise ratio is given by [141–145]

$$ \text{SNR} = \sqrt{t_{\text{obs}} \int \left( \frac{\Omega_{\text{GW}}^2}{\Omega_{\text{sens}}^2 + 2\Omega_{\text{GW}}\Omega_{\text{sens}} + 2\Omega_{\text{GW}}^2} \right) df}, \quad (5.15) $$

where $t_{\text{obs}}$ is the observation time and $\Omega_{\text{sens}}$ the sensitivity of the interferometer. For these quantities we take without deviation the choices made in [140] and use this to also calculate the power-law-integrated sensitivity curves shown in figure 5. Note the second and third terms in the denominator of the integral in eq. (5.15) means the SNR is saturated for large GW signals. We do not want to count signals below astrophysical foregrounds as detectable, so as a simple and naive estimate, we define a foreground limited signal-to-noise ratio in which we only count GW signals above the astrophysical contamination. Namely, we compute,

$$ \text{SNR}_{\text{FGL}} = \sqrt{t_{\text{obs}} \int \left( \frac{\text{Max}[0, \Omega_{\text{GW}}(\nu) - \Omega_{\text{FG}}(\nu)]^2}{\Omega_{\text{sens}}^2 + 2\Omega_{\text{GW}}\Omega_{\text{sens}} + 2\Omega_{\text{GW}}^2} \right) df}, \quad (5.16) $$

As estimates of the foregrounds we take:

* The component of the GW galactic white dwarf binaries which are not subtractable after four years of LISA observations, with approximate form given in [146, 147]. The same quantity has also been estimated in [148], with a higher peak value, but at lower frequencies. For this latter estimate we use the analytic fit from [149]. To be conservative we simply include both estimates.

* The upper value of the extra-galactic white-dwarf binary estimate from [150], a broadly similar estimate is derived in [151].

* The central value binary neutron star and binary black hole signal inferred from the observations in [152], extrapolated to lower frequencies assuming the $f^{2/3}$ scaling, as is
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Top: the latent heat $\alpha$ and inverse timescale of the transitions $\beta/H$ returning $f_{\text{PBH}} = 1$ as a function of $M_{\text{PBH}}$. Also shown is a comparison to the results of [56], which uses the approximate $\Gamma_{\text{bub}}$ ansatz, so some differences are expected. Middle: the peak amplitude and peak frequency of the three GW spectrum models as a function of $M_{\text{PBH}}$. Bottom: the signal-to-noise ratio (SNR), and signal-to-noise ratio above foregrounds (SNR$_{\text{FGL}}$), as a function of $M_{\text{PBH}}$. The signal can be detected by BDECIGO or a combination of LISA and ET over the entire allowed range of $M_{\text{PBH}}$ with $f_{\text{PBH}} = 1$. The shaded band encompasses the three different determinations of the GW signal. Additional theory and experimental forecasting uncertainty means the true error band is larger.</figcaption></figure>
appropriate for a stationary ensemble of binaries with circular orbits losing energy
solely through GW emission [153].

* The continuous super-massive binary black hole signal estimate of [151] will not mask our signal, but we display it in our plots.

To summarize, in figure 6 we show the latent heat $\alpha$, $\beta_H$, peak amplitude $\Omega_{GW}(f_{peak})$, the peak frequency $f_{peak}$, SNR, and SNR$_{FGL}$ as a function of $M_{PBH}$ (the vev $v_{B-L}$ and coupling are fixed by $M_{PBH}$ and the requirement $f_{PBH} = 1$).

In the top row of figure 6 we see, as expected, PBH formation occurs for $\alpha \gg 1$ and $\beta_H(T_n) \lesssim 8$, with slightly slower transitions required for more massive PBHs, in order to compensate the smaller enhancement in $\rho_{PBH}/\rho_{rad}$ between PBH formation time and matter-radiation equality. We also compare to the required $\beta_H$ found for $f_{PBH} = 1$ in ref. [56], which used the approximate ansatz eq. (2.6) for the nucleation rate. We obtain the required abundance with a slightly higher $\beta_H$, because the bubble nucleation rate is somewhat suppressed in the full calculation compared to the approximate method. Note the collapse probability, $P_{coll}$, is very sensitive to $\beta_H$, thus a small shift in the latter can compensate for changes in $\Gamma_{bub}$ between the full and approximate expressions.

We thus confirm the approximate method is valid, at least in close-to-conformal potentials as studied here, which do not have a zero-temperature barrier. (If a zero temperature barrier exists, then there is the parameter dependent risk of having the field become permanently trapped in the false vacuum [154, 155]. Then PBH formation would not occur simply because the background patches themselves do not percolate. Such models may require a more careful examination. Nevertheless, bubble sizes and hence GW signals are still necessarily large in successfully percolating models which produce PBHs and feature more dramatic decreases in $\Gamma_{bub}$ below $T_n$, see [55] for work along these lines.)

In the middle row, we see the amplitude in $\Omega_{GW}$ is large and almost constant over the range of $M_{PBH}$, as $\beta_H$ is almost constant. The peak frequency covers three orders of magnitude, ranging over the sensitivities for LISA, BDECIGO, and ET, as anticipated through eqs. (2.1) and (2.5).

In the bottom row of figure 6, we observe that the typical supercooled phase transitions in the classically conformal $B-L$ model explaining all the dark matter in our universe in the form of PBHs will give rise to extremely strong GWs signals detectable by LISA, BDECIGO, and ET. These experiments will be able to probe the entire parameter space of the model producing $f_{PBH} = 1$. Current limits do not constrain this PBH formation mechanism in the given range of $M_{PBH}$ [156, 157].

# 6 Validity of assumptions

## 6.1 Rapid decay of the condensate excitations

In the above, we have assumed rapid reheating, so that $T_{RH} \simeq T_{infl}$. To check whether this is a valid assumption in our model, we can compare the decay rate of the scalar $\rho$ with the Hubble rate at the end of the PT. The first partial width we consider is of $\rho$ decaying
<figure><img src="image_9.png" /><figcaption>Figure 7. Left: the decay rates compared to the Hubble rate. Right: the LO and the γ enhanced NLO retarding pressure on the wall compared to the driving pressure, i.e. the vacuum energy difference.</figcaption></figure>

to the EW Higgs doublet

$$ \Gamma_{\rho \to H^{\dagger}H} = \frac{\lambda_{\rho h}^2 v_{\rho}^2}{8\pi m_{\rho}} \text{Re} \left[ \sqrt{1 - \frac{4m_H^2}{m_{\rho}^2}} \right]. \qquad (6.1) $$

Here we take into account the thermal mass of the Higgs, which begins to become present as the bath is reheated,

$$ m_H^2 \approx \left( \frac{\lambda_h^2}{2} + 3\frac{g_2^2}{16} + \frac{g_Y^2}{16} + \frac{\lambda_{h\rho}}{12} \right) T^2 \approx 0.4T^2. \qquad (6.2) $$

The above width is tiny, due to the suppression of $\lambda_{\rho h}^2$, and possible kinematic factors as $m_H(T)$ becomes large compared to $m_\rho$. But we can also consider the partial width to the $N_i$,

$$ \Gamma_{\rho \to N_i N_i} = \frac{y_{N_i}^2 m_{\rho}}{32\pi} \left(1 - \frac{2M_{N_i}^2}{m_{\rho}^2}\right) \text{Re} \left[ \sqrt{1 - \frac{4M_{N_i}^2}{m_{\rho}^2}} \right]. \qquad (6.3) $$

As $M_{Z'} \approx 10T_{\text{RH}}$, the thermal masses of the $N_i$ are negligible following the PT. For standard assumptions regarding the typical Yukawa coupling size — as inferred in from the type-I seesaw and the atmospheric mass squared difference of the active neutrinos — we can assume the $N_i$ decay rapidly compared to the Hubble rate, i.e. the strong washout regime [158]. We can thus simply compare eq. (6.3) to the Hubble rate to determine whether we can avoid a long lived $\rho$ or other matter domination following the PT. In figure 7 left we show both the partial widths compared to the Hubble rate for our parameter choice, showing our assumptions of rapid decay, giving $T_{\text{RH}} \simeq T_{\text{infl}}$ is justified. Note the $Z'$ will also rapidly decay to SM fermions.

## 6.2 Pressure on the bubble wall

For supercooled transitions, such as those studied here, leading order pressure on the wall as particles gain their mass [159],

$$ P_{\text{LO}} \approx \frac{T_{\rho}^2}{24} \left( 3M_{Z'}^2 + m_{\rho}^2 + M_{N_i}^2 \right), \qquad (6.4) $$
is generically insufficient to stop the wall from accelerating. The presence of gauge bosons which change their mass across the wall, however, leads to a NLO retarding pressure which increases linearly with $\gamma$, the bubble wall Lorentz factor [160]. We therefore check whether the vacuum energy density is converted into the wall energy, or whether sufficient soft-gauge bosons are produced at the wall, for the energy to instead be stored in the latter before wall collision. Such a situation has been modelled in [161]. Whether this would affect the GW signal appreciably is a separate question, macroscopically both cases seem to be captured by the bulk flow model (e.g. see [161] and discussion in [162]). Soft gauge boson production leads to a retarding pressure on the wall $P_{\text{NLO}} \propto \gamma T_p^3 M_{Z'}$ [160, 163]. Here we use the leading-log determination for this term given in [164],

$$P_{\text{NLO}} \approx \frac{\kappa \zeta(3) (Q_{B-L}^{\text{eff}})^2 \alpha_{B-L} \gamma M_{Z'} T_p^3}{\pi^3} \log \left( \frac{v_{\rho}}{T_p} \right), \quad (6.5)$$

where $\kappa \approx 4$, $(Q_{B-L}^{\text{eff}})^2 = 40$ is an effective charge factor summing over the bath degrees-of-freedom with $3/4$ (1) weighting for fermions (scalsars), $\alpha_{B-L} \equiv g_{B-L}^2/4\pi$, and we have used the thermal mass of the $Z'$ as an IR cutoff in the logarithm. The maximum Lorentz factor is obtained at collision for an effectively run-away wall, i.e. approximately zero retarding pressure, and is given by [164]

$$\gamma \approx \frac{1}{3} \frac{R_{\text{coll}}}{R_{\text{nuc}}} \approx \frac{1}{3} \frac{\pi^{1/3}}{\beta_H H} \left( \frac{T_n}{10} \right), \quad (6.6)$$

where $R_{\text{nuc}} \approx 10/T_n$ is the bubble radius at nucleation and $R_{\text{coll}} \approx \pi^{1/3}/\beta$ is the bubble radius at collision (we have checked both analytic estimates through numerical determination). For consistency of the above $\gamma$ estimate we require $P_{\text{NLO}} < \Lambda_{\text{vac}}$. We calculate these quantities and display the results in figure 7 right, showing the wall is in the effective run-away regime at collision.

## 6.3 Reliability of the effective potential

We now briefly discuss the reliability of the approximations we have made in deriving $V_{\text{eff}}$. To answer this question requires the application of more refined techniques such as the renormalization group improved effective potential [111, 165]. Indeed, in the context of the particle-DM dilution mechanisms, there are known examples in which going to the RGE improved $V_{\text{eff}}$ changes the qualitative outcome of the mechanism itself [166]. This is due to changes to the details of the nucleation, percolation, and reheating temperatures, although the PTs remain very strong. Such an investigation is left for future work.

## 7 Comments on leptogenesis

In figure 8 we show the dilution factor due to the entropy production following the PT. As the dilution factor is quite large, $\sim 10^7 - 10^9$, it seems unlikely that we could use a mechanism which relies on the dilute false vacuum plasma as a source of the baryon asymmetry. For example, the mass gain mechanism when the $N_i$ suddenly gain their mass during the PT and then decay to generate the BAU [167–171]. Similar conclusions most likely also hold for
<figure><img src="image_10.png" /><figcaption>Figure 8. The entropy production (dilution) factor following the PT.</figcaption></figure>

the related mechanism in [172]. Electroweak style baryogenesis mechanisms suffer from the same problem in the current context, together with a completely negligible yield for $v_w \simeq 1$, due to lack of diffusion back into the symmetric phase [173].

If instead, heavy particles are produced in the reheating process itself, following along the lines of [174, 175], then these may have a larger abundance compared to the entropy density, and act as a source of the asymmetry. But note our bubble collisions are expected to be inelastic, which suppresses heavy particle production [176].

Anyway, for our choice of example Yukawa coupling $y_i = 0.01$, the asymmetry will be generated by the usual thermal leptogenesis around $T \approx M_N$ after the PT (for our parameter choices $T_{RH} \simeq 10 M_{Ni}$). In either case, resonant leptogenesis is generically required [177–180], as the heavy neutrinos are much below the Davidson-Ibarra bound [181] (also see [182]). Note the heavy neutrinos are efficiently produced in the early universe due to their Yukawa coupling to $\rho$ and their gauge interactions.

## 8 Conclusion

In this work we have studied the GW signal from PBH formation during a PT. We have assumed $f_{PBH} = 1$, which limits the range for the peak frequency of the GWs, and calculated the required underlying parameters of the model to obtain the PBH abundance. We used the properties of the PT to find the expected GW signal from bubble collisions over the allowed parameter space, showing it can be detected by BDECIGO or a combination of LISA and ET, independently of the precise $M_{PBH}$, over the entire range allowed for $f_{PBH} = 1$. Substituting one of the experiments by a broadly equivalent one will, of course, not affect the conclusions. An important caveat is that the GW predictions for strong phase transitions may overestimate the signal because expansion during the transition is not taken into account. This motivates further refinement of such estimates, along the lines of [78], but extended to vacuum dominated transitions and beyond the envelope approximation. We note a suppression in $\Omega_{GW}$ by an order of magnitude for $\beta_H \approx 10$, would still return SNRs of $10^2$ and above over the $f_{PBH} = 1$ parameter space, so prospects remain very promising.
The classically scale invariant $B-L$ model we considered does not contain an automatically stable DM candidate, although given a small enough Yukawa coupling, $N_1$ could play the role of the possibility of a ~keV scale sterile neutrino (provided an appropriate mechanism was also in play in order to set its relic abundance correctly, as it carries gauge charge). Such sterile neutrino scenarios are by now by now heavily constrained [183]. Additional field content can of course be included to act as DM [184]. We have instead proposed using the PBH generation mechanism from the strong PT realised with the close-to-conformal potential to obtain the observed $\Omega_{DM}$. Leptogenesis in this model should occur sometime after the PT, through resonant enhancement of the CP violation, as in the usual low-scale thermal leptogenesis [177–180]. Furthermore, it would also be of interest to re-examine models of particle DM, in which the relic abundance relies on supercooled PTs [87, 140, 166, 185–187], to check compatibility with PBH production.

Our calculations in a specific model, taking into account the full temperature dependence of the bubble nucleation rate, have shown only modestly small changes to the inverse timescale of the transition $\beta_H$ required to achieve $f_{PBH} = 1$, compared to model independent calculation which used an approximate ansatz for the nucleation rate. As $\beta_H$ sets the typical bubble size at collision and therefore also controls the properties of the GW signal, our results regarding the sensitivities of future detectors to such a late patch mechanism are independent of the underlying particle physics model.

Models of slow roll inflation producing PBHs also result in detectable GWs, in the same frequency range as the PT, due to the enhanced power spectrum on small scales leading to significant anisotropic stress [34–39]. The additional power also present on small scales from the PT, would also induce anisotropic stress and a further GW source. Accordingly, this paper justifies the further development of computational techniques in order to pin down the expected PBH spectrum in the late patch mechanism (beyond the monochromatic approximation primarily used here), the bubble collision signal, together with any additional GWs at second order in perturbation theory, in light of testing $f_{PBH} = 1$ in upcoming interferometers. Along with the GW signal, the PBHs can be searched for using improved lensing studies, through future MeV telescopes [188, 189], and the 21-cm absorption signal [190, 191].

## Acknowledgments

We thank Yann Gouttenoire, Thomas Konstandin, Florian Kühnel, and Filippo Sala for useful discussions.

**Funding information.** This work was supported by the European Union’s Horizon 2020 research and innovation programme under grant agreement No 101002846, ERC CoG “CosmoChart”.

## A Towards a non-monochromatic mass distribution

In the main text, we used a monochromatic approximation for the PBH mass,

$$M_{\text{PBH}} = \frac{4\pi(\rho_{\text{radlate}} + \rho_{\text{vaclate}})c_s^3}{3H_{\text{late}}^3} = \frac{c_s^3 M_{\text{Pl}}^2}{2H_{\text{late}}} \quad (\text{A.1})$$
We also found

$$f_{\text{PBH}} = \frac{\Omega_{\text{DM}} + \Omega_{\text{B}} g_{*s}(T_{\text{eq}}) T_{\text{bkg form}}}{\Omega_{\text{DM}} g_{*}(T_{\text{eq}}) T_{\text{eq}}} \frac{\rho_{\text{PBH}}(T_{\text{bkg form}})}{\rho_{\text{rad}}(T_{\text{bkg form}})}, \quad (A.2)$$

where $(\Omega_{\text{DM}} + \Omega_{\text{B}})/\Omega_{\text{DM}} \simeq 1.2$. We then used

$$\frac{\rho_{\text{PBH}}(T_{\text{bkg form}})}{\rho_{\text{rad}}(T_{\text{bkg form}})} \approx P_{\text{coll}} \frac{M_{\text{PBH}}}{M_{\text{Hor}}} \approx c_s^3 P_{\text{coll}} \frac{H_{\text{bkg}}(t_{\delta\text{max}})}{H_{\text{late}}(t_{\delta\text{max}})}, \quad (A.3)$$

where $M_{\text{Hor}} = M_{\text{Pl}}^2/(2H_{\text{bkg}})$ is the energy inside the Hubble horizon of a background patch. In the monochromatic approximation $P_{\text{coll}}(T_i)$, given in eq. (4.20b), is evaluated at $T_i$ chosen so that $\delta_{\text{max}} = \delta_c$.

We now assume instead critical collapse, as is known from numerical codes relevant for the inflation scenario, in which for $\delta_m > \delta_c$ the PBH mass is found to be [125, 192]

$$M_{\text{PBH}} = \mathcal{K} M_{\text{Hor}} (\delta_m - \delta_c)^{\gamma_c}, \quad (A.4)$$

where $1 \lesssim \mathcal{K} \lesssim 10$ is an efficiency factor, and $\gamma_c \simeq 0.36$ is the critical exponent assuming radiation domination (both found from simulations). In the inflationary scenario, $\delta_m$ is the overdensity as measured at the maximum of the compactification function, at cosmological horizon crossing time, e.g. see [27]. We now make assume something similar holds in the PT mechanism and replace $\delta_m \to \delta_{\text{max}}$ in the above equation. We warn we are making a big leap making such an assumption for critical collapse in the PT mechanism, and do so simply out of curiousity regarding the resulting PBH spectrum, mainly as an academic exercise. Note in using eq. (A.4) we also do not include the small correction taking into account the different Horizon masses of the collapsing and background patches as we have done for our monochromatic approximation in eq. (A.1).

From eq. (A.4) we see that patches with $\delta_{\text{max}} = \delta_c$ fail to form a black hole, but instead a whole spectrum of PBH masses is generated from patches with larger $\delta_{\text{max}}$, coming from progressively smaller $T_i$. The probability of keeping large regions bubble free also falls as $T_i$ is reduced. (Roughly speaking, $P_{\text{coll}}$ is the probability of keeping a sufficiently large region bubble free sufficiently long, it therefore includes the probability of staying bubble free even longer. In other words it acts as a cumulative probability function.) A plot of $P_{\text{coll}}(T_i)$ is shown in figure 9.

We now wish to calculate the resulting PBH mass spectra assuming eq. (A.4) holds. The first spectrum typically considered in the literature is normalized to the PBH density and is conventionally defined as [192]

$$\Psi(M_{\text{PBH}}) \equiv \frac{1}{\rho_{\text{PBH}}} \frac{d\rho_{\text{PBH}}}{dM_{\text{PBH}}}. \quad (A.5)$$

In our scenario we have

$$\Psi(M_{\text{PBH}}) = \frac{M_{\text{PBH}}}{\mathcal{N}} \frac{dP_{\text{coll}}}{dM_{\text{PBH}}} = \frac{M_{\text{PBH}}}{\mathcal{N}} \frac{dP_{\text{coll}}}{dT_i} \frac{dT_i}{dM_{\text{PBH}}} = \frac{1}{\mathcal{N}\gamma_c} \left( \frac{M_{\text{PBH}}}{\mathcal{K}M_{\text{Hor}}} \right)^{1/\gamma_c} \frac{dP_{\text{coll}}}{dT_i} \frac{dT_i}{d\delta_{\text{max}}}, \quad (A.6)$$

where the normalization factor is given by

$$\mathcal{N} = \int_{0}^{\infty} dM'_{\text{PBH}} M'_{\text{PBH}} \frac{dP_{\text{coll}}}{dM'_{\text{PBH}}}. \quad (A.7)$$
<figure><img src="image_11.png" /><figcaption>Figure 9. Left: the cumulative probability $P(T_i)$ for the example parameter point, $v_ρ = 10^3$ TeV, $y_{Ni} = 0.01$, and $g_{B-L} = 0.2958$. Right: for convenience we again show $δ_{max}(T_i)$ for the same parameter point, where $δ_{max} = 0.45$ for $T_i = 53.6$ GeV (for a larger range of $T_i$ see figure 3, middle row, right).</figcaption></figure>

Note $dP_{coll}/dM_{PBH} < 0$ given our definition of $P_{coll}$ in eq. (4.20b), but the factor appears both in the numerator and the denominator, $N$, so $Ψ$ is positive. Integrating $Ψ(M_{PBH})$ over the entire PBH mass range yields unity. In the above, we have related the spectrum to the collapse probability calculated in eq. (4.20b), and we have used that $M_{Hor}$ and $T_{bkg}$ form are approximately independent of $M_{PBH}$ for PBHs originating from overdensities with a preferred length scale (in the current context of the late patch mechanism this is simply the Hubble length during the PT). This is further justified by the strongly falling $P_{coll}$ with decreasing $T_i$, as shown in figure 9.

The second commonly used spectrum is instead normalized to the observed DM density and is defined as [192]

$$f(M_{\text{PBH}}) = \frac{M_{\text{PBH}}}{\rho_{\text{DM}}} \frac{d\rho_{\text{PBH}}}{dM_{\text{PBH}}}, \qquad (A.8)$$

so that

$$f_{\text{PBH}} = \int_{0}^{\infty} dM_{\text{PBH}} \frac{f(M_{\text{PBH}})}{M_{\text{PBH}}} \quad \text{and} \quad f(M_{\text{PBH}}) = f_{\text{PBH}} M_{\text{PBH}} \Psi(M_{\text{PBH}}). \qquad (A.9)$$

In our scenario we find

$$f(M_{\text{PBH}}) = - \frac{\Omega_{\text{DM}} + \Omega_B g_{*s}(T_{\text{eq}}) \frac{T_{\text{bkg form}}}{T_{\text{eq}}} \frac{M_{\text{PBH}}^2}{M_{\text{Hor}}} \frac{dP_{\text{coll}}}{dM_{\text{PBH}}}}{\Omega_{\text{DM}} g_{*}(T_{\text{eq}}) T_{\text{eq}}} \qquad (A.10a)$$

$$= - \frac{\Omega_{\text{DM}} + \Omega_B g_{*s}(T_{\text{eq}}) \frac{T_{\text{bkg form}}}{T_{\text{eq}}} \frac{M_{\text{PBH}}^2}{M_{\text{Hor}}} \frac{dP_{\text{coll}}}{dT_i} \frac{dT_i}{dM_{\text{PBH}}}}{\Omega_{\text{DM}} g_{*}(T_{\text{eq}}) T_{\text{eq}}} \qquad (A.10b)$$

$$= - \frac{\Omega_{\text{DM}} + \Omega_B g_{*s}(T_{\text{eq}}) \frac{T_{\text{bkg form}}}{T_{\text{eq}}} \frac{1}{\gamma_c K^{1/\gamma_c}} \left(\frac{M_{\text{PBH}}}{M_{\text{Hor}}}\right)^{1+1/\gamma_c} \frac{dP_{\text{coll}}}{dT_i} \frac{dT_i}{d\delta_{\text{max}}}}{\Omega_{\text{DM}} g_{*}(T_{\text{eq}}) T_{\text{eq}}} \qquad (A.10c)$$

where we have again used that $M_{Hor}$ and $T_{bkg}$ form are approximately independent of $M_{PBH}$ for PBHs originating from overdensities with a preferred length scale. We now include an overall negative sign because $dP_{coll}/dM_{PBH} < 0$. For either of the spectra of interest, we can evaluate the last two differential factors numerically as functions of $M_{PBH}$, and thus find the desired quantity.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** The PBH spectra for the example parameter point, $v_ρ = 10^3$ TeV, $y_{Ni} = 0.01$, and $g_{B-L} = 0.2958$, now taking into account that $δ_{max}$ is a function of $T_i$, and assuming the critical collapse scaling with $K = 1$. Upon integration we find $\langle M_{PBH} \rangle \simeq 5.4 \times 10^{-12} M_⊙$ and an integrated total $f_{PBH} \simeq 1.05$, compared with $M_{PBH} = 4.5 \times 10^{-12} M_⊙$ and $f_{PBH} \simeq 1.02$ using the monochromatic approximation. The grey dashed lines show the expected IR scalings, $\Psi(M_{PBH}) \propto M_{PBH}^{1/\gamma_c}$ and $f_{PBH}(M_{PBH}) \propto M_{PBH}^{1+1/\gamma_c}$ respectively [192], coming from values of $T_i$ close to which $δ_{max} = 0.45$. These IR scalings can be understood from eqs. (A.6) and (A.10c).</figcaption></figure>

We now return to our example parameter point used in figure 3 and calculate the PBH mass spectra assuming the critical collapse phenomenon holds. The resulting PBH spectra are shown in figure 10. Clearly the spectra are strongly peaked at $M_{PBH} \approx c_{sound}^3 M_{Hor} \approx 0.2 M_{Hor}$, the value assumed for the monochromatic approximation. Furthermore, the total integrated PBH density normalized to the observed DM density remains largely unchanged. Thus the monochromatic approximation is justified. We can also compare the resulting spectra to existing limits as done in figure 11, but only to gain a rough idea, as the limits are strictly only valid in the monochromatic approximation. Still it is reassuring that the spectra are strongly enough peaked that the IR tails do not intersect with the low mass constraints, as long as $\langle M_{PBH} \rangle \gtrsim \mathcal{O}(10^{-16}) M_⊙$. The UV cut-off is much steeper and so even less of a concern.

Even assuming the critical collapse relation holds for the late patch mechanism, the above is still not a complete estimate of the full PBH spectrum. This is because: (i) we have ignored the possibility of PBHs being formed at slightly different times with different $M_{Hor}$, and (ii) the collapsing volume is not limited to eq. (4.21b), but may take a whole range of values (for us, horizon and superhorizon at the time of the PT, but one may also imagine that sub-horizon patches with larger densities can also collapse, as in ref. [55]). Given the theoretical uncertainties involved in applying the critical collapse to the PT scenario, however, together with the computational time required for the numerical evaluation, we have not yet evaluated this full expression. Nevertheless, the peaked nature of the PBH spectrum is expected to hold also in this more general case, because the probability of obtaining larger $V_{coll}$ will also fall rapidly.

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.
<figure><img src="image_13.png" /></figure>

**Figure 11.** The PBH spectra for two example parameter points, both with $y_{Ni} = 0.01$ and $\mathcal{K} = 1$, compared to limits on monochromatic spectra from the CMB [11, 12], extra galactic background light [13], galactic gamma ray background [15], Voyager $e^{\pm}$ [16], and Subaru Hyper Suprime-Cam (HSC) lensing [20]. The GWs from the same two parameter points are shown in figure 5.

References

[1] M. Milgrom, *A modification of the Newtonian dynamics as a possible alternative to the hidden mass hypothesis*, Astrophys. J. 270 (1983) 365 [INSPIRE].

[2] J. Bekenstein and M. Milgrom, *Does the missing mass problem signal the breakdown of Newtonian gravity?*, Astrophys. J. 286 (1984) 7 [INSPIRE].

[3] J.D. Bekenstein, *Relativistic gravitation theory for the MOND paradigm*, Phys. Rev. D 70 (2004) 083509 [Erratum ibid. 71 (2005) 069901] [astro-ph/0403694] [INSPIRE].

[4] F. Lelli, S.S. McGaugh, J.M. Schombert and M.S. Pawlowski, *One law to rule them all: the radial acceleration relation of galaxies*, Astrophys. J. 836 (2017) 152 [arXiv:1610.08981] [INSPIRE].

[5] C. Skordis and T. Zlosnik, *New relativistic theory for modified Newtonian dynamics*, Phys. Rev. Lett. 127 (2021) 161302 [arXiv:2007.00082] [INSPIRE].

[6] D.B. Thomas, A. Mozaffari and T. Zlosnik, *Consistent cosmological structure formation on all scales in relativistic extensions of MOND*, JCAP 06 (2023) 006 [arXiv:2303.00038] [INSPIRE].

[7] Y.B. Zel'dovich and I.D. Novikov, *The hypothesis of cores retarded during expansion and the hot cosmological model*, Sov. Astron. 10 (1967) 602 [INSPIRE].

[8] S. Hawking, *Gravitationally collapsed objects of very low mass*, Mon. Not. Roy. Astron. Soc. 152 (1971) 75 [INSPIRE].

[9] B.J. Carr and S.W. Hawking, *Black holes in the early universe*, Mon. Not. Roy. Astron. Soc. 168 (1974) 399 [INSPIRE].

[10] B. Carr and F. Kuhnel, *Primordial black holes as dark matter candidates*, SciPost Phys. Lect. Notes 48 (2022) 1 [arXiv:2110.02821] [INSPIRE].

[11] S.K. Acharya and R. Khatri, *CMB and BBN constraints on evaporating primordial black holes revisited*, JCAP 06 (2020) 018 [arXiv:2002.00898] [INSPIRE].
[12] J. Chluba, A. Ravenni and S.K. Acharya, Thermalization of large energy release in the early universe, Mon. Not. Roy. Astron. Soc. 498 (2020) 959 [arXiv:2005.11325] [INSPIRE].

[13] B.J. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, New cosmological constraints on primordial black holes, Phys. Rev. D 81 (2010) 104019 [arXiv:0912.5297] [INSPIRE].

[14] G. Ballesteros, J. Coronado-Blázquez and D. Gaggero, X-ray and gamma-ray limits on the primordial black hole abundance from Hawking radiation, Phys. Lett. B 808 (2020) 135624 [arXiv:1906.10113] [INSPIRE].

[15] B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, Constraints on primordial black holes, Rept. Prog. Phys. 84 (2021) 116902 [arXiv:2002.12778] [INSPIRE].

[16] M. Boudaud and M. Cirelli, Voyager 1e$^{+}$ and 1e$^{-}$ further constrain primordial black holes as dark matter, Phys. Rev. Lett. 122 (2019) 041104 [arXiv:1807.03075] [INSPIRE].

[17] W. DeRocco and P.W. Graham, Constraining primordial black hole abundance with the galactic 511 keV line, Phys. Rev. Lett. 123 (2019) 251102 [arXiv:1906.07740] [INSPIRE].

[18] R. Laha, Primordial black holes as a dark matter candidate are severely constrained by the galactic center 511 keV γ-ray line, Phys. Rev. Lett. 123 (2019) 251101 [arXiv:1906.09994] [INSPIRE].

[19] R. Laha, J.B. Muñoz and T.R. Slatyer, INTEGRAL constraints on primordial black holes and particle dark matter, Phys. Rev. D 101 (2020) 123514 [arXiv:2004.00627] [INSPIRE].

[20] H. Niikura et al., Microlensing constraints on primordial black holes with Subaru/HSC Andromeda observations, Nature Astron. 3 (2019) 524 [arXiv:1701.02151] [INSPIRE].

[21] B.J. Carr and J.E. Lidsey, Primordial black holes and generalized constraints on chaotic inflation, Phys. Rev. D 48 (1993) 543 [INSPIRE].

[22] P. Ivanov, P. Naselsky and I. Novikov, Inflation and primordial black holes as dark matter, Phys. Rev. D 50 (1994) 7173 [INSPIRE].

[23] B.J. Carr, The primordial black hole mass spectrum, Astrophys. J. 201 (1975) 1 [INSPIRE].

[24] M. Shibata and M. Sasaki, Black hole formation in the Friedmann universe: formulation and computation in numerical relativity, Phys. Rev. D 60 (1999) 084002 [gr-qc/9905064] [INSPIRE].

[25] I. Musco, J.C. Miller and L. Rezzolla, Computations of primordial black hole formation, Class. Quant. Grav. 22 (2005) 1405 [gr-qc/0412063] [INSPIRE].

[26] T. Harada, C.-M. Yoo and K. Kohri, Threshold of primordial black hole formation, Phys. Rev. D 88 (2013) 084051 [Erratum ibid. 89 (2014) 029903] [arXiv:1309.4201] [INSPIRE].

[27] I. Musco, Threshold for primordial black holes: dependence on the shape of the cosmological perturbations, Phys. Rev. D 100 (2019) 123524 [arXiv:1809.02127] [INSPIRE].

[28] I. Musco, V. De Luca, G. Franciolini and A. Riotto, Threshold for primordial black holes. II. A simple analytic prescription, Phys. Rev. D 103 (2021) 063538 [arXiv:2011.03014] [INSPIRE].

[29] C. Germani and I. Musco, Abundance of primordial black holes depends on the shape of the inflationary power spectrum, Phys. Rev. Lett. 122 (2019) 141302 [arXiv:1805.04087] [INSPIRE].

[30] A. Escrivà, C. Germani and R.K. Sheth, Universal threshold for primordial black hole formation, Phys. Rev. D 101 (2020) 044022 [arXiv:1907.13311] [INSPIRE].

[31] A. Escrivà and A.E. Romano, Effects of the shape of curvature peaks on the size of primordial black holes, JCAP 05 (2021) 066 [arXiv:2103.03867] [INSPIRE].
[32] A. Escrivà, E. Bagui and S. Clesse, *Simulations of PBH formation at the QCD epoch and comparison with the GWTC-3 catalog*, JCAP 05 (2023) 004 [arXiv:2209.06196] [INSPIRE].

[33] I.D. Stamou, *Exploring critical overdensity thresholds in inflationary models of primordial black holes formation*, Phys. Rev. D 108 (2023) 063515 [arXiv:2306.02758] [INSPIRE].

[34] T. Nakama, J. Silk and M. Kamionkowski, *Stochastic gravitational waves associated with the formation of primordial black holes*, Phys. Rev. D 95 (2017) 043511 [arXiv:1612.06264] [INSPIRE].

[35] J. Garcia-Bellido, M. Peloso and C. Unal, *Gravitational wave signatures of inflationary models from primordial black hole dark matter*, JCAP 09 (2017) 013 [arXiv:1707.02441] [INSPIRE].

[36] R.-G. Cai, S. Pi and M. Sasaki, *Gravitational waves induced by non-Gaussian scalar perturbations*, Phys. Rev. Lett. 122 (2019) 201101 [arXiv:1810.11000] [INSPIRE].

[37] N. Bartolo et al., *Primordial black hole dark matter: LISA serendipity*, Phys. Rev. Lett. 122 (2019) 211301 [arXiv:1810.12218] [INSPIRE].

[38] N. Bartolo et al., *Testing primordial black holes as dark matter with LISA*, Phys. Rev. D 99 (2019) 103521 [arXiv:1810.12224] [INSPIRE].

[39] W. Qin et al., *Planck constraints and gravitational wave forecasts for primordial black hole dark matter seeded by multifield inflation*, Phys. Rev. D 108 (2023) 043508 [arXiv:2303.02168] [INSPIRE].

[40] S.W. Hawking, I.G. Moss and J.M. Stewart, *Bubble collisions in the very early universe*, Phys. Rev. D 26 (1982) 2681 [INSPIRE].

[41] M. Crawford and D.N. Schramm, *Spontaneous generation of density perturbations in the early universe*, Nature 298 (1982) 538 [INSPIRE].

[42] H. Kodama, M. Sasaki and K. Sato, *Abundance of primordial holes produced by cosmological first order phase transition*, Prog. Theor. Phys. 68 (1982) 1979 [INSPIRE].

[43] S.D.H. Hsu, *Black holes from extended inflation*, Phys. Lett. B 251 (1990) 343 [INSPIRE].

[44] I.G. Moss, *Singularity formation from colliding bubbles*, Phys. Rev. D 50 (1994) 676 [INSPIRE].

[45] M.Y. Khlopov, R.V. Konoplich, S.G. Rubin and A.S. Sakharov, *Formation of black holes in first order phase transitions*, hep-ph/9807343 [INSPIRE].

[46] M. Lewicki and V. Vaskonen, *On bubble collisions in strongly supercooled phase transitions*, Phys. Dark Univ. 30 (2020) 100672 [arXiv:1912.00997] [INSPIRE].

[47] J. Liu et al., *Primordial black hole production during first-order phase transitions*, Phys. Rev. D 105 (2022) L021303 [arXiv:2106.05637] [INSPIRE].

[48] K. Hashino, S. Kanemura and T. Takahashi, *Primordial black holes as a probe of strongly first-order electroweak phase transition*, Phys. Lett. B 833 (2022) 137261 [arXiv:2111.13099] [INSPIRE].

[49] C. Gross, G. Landini, A. Strumia and D. Teresi, *Dark matter as dark dwarfs and other macroscopic objects: multiverse relics?*, JHEP 09 (2021) 033 [arXiv:2105.02840] [INSPIRE].

[50] M.J. Baker, M. Breitbach, J. Kopp and L. Mittnacht, *Detailed calculation of primordial black hole formation during first-order cosmological phase transitions*, arXiv:2110.00005 [INSPIRE].

[51] K. Kawana and K.-P. Xie, *Primordial black holes from a cosmic phase transition: the collapse of Fermi-balls*, Phys. Lett. B 824 (2022) 136791 [arXiv:2106.00111] [INSPIRE].
[52] S. He, L. Li, Z. Li and S.-J. Wang, *Gravitational waves and primordial black hole productions from gluodynamics by holography*, arXiv:2210.14094 [insPIRE].

[53] K. Hashino, S. Kanemura, T. Takahashi and M. Tanaka, *Probing first-order electroweak phase transition via primordial black holes in the effective field theory*, Phys. Lett. B 838 (2023) 137688 [arXiv:2211.16225] [insPIRE].

[54] K. Kawana, T.H. Kim and P. Lu, *PBH formation from overdensities in delayed vacuum transitions*, Phys. Rev. D 108 (2023) 103531 [arXiv:2212.14037] [insPIRE].

[55] M. Lewicki, P. Toczek and V. Vaskonen, *Primordial black holes from strong first-order phase transitions*, JHEP 09 (2023) 092 [arXiv:2305.04924] [insPIRE].

[56] Y. Gouttenoire and T. Volansky, *Primordial black holes from supercooled phase transitions*, arXiv:2305.04942 [insPIRE].

[57] A. Salvio, *Supercooling in radiative symmetry breaking: theory extensions, gravitational wave detection and primordial black holes*, JCAP 12 (2023) 046 [arXiv:2307.04694] [insPIRE].

[58] J. Garriga, A. Vilenkin and J. Zhang, *Black holes and the multiverse*, JCAP 02 (2016) 064 [arXiv:1512.01819] [insPIRE].

[59] H. Deng, J. Garriga and A. Vilenkin, *Primordial black hole and wormhole formation by domain walls*, JCAP 04 (2017) 050 [arXiv:1612.03753] [insPIRE].

[60] H. Deng and A. Vilenkin, *Primordial black hole formation by vacuum bubbles*, JCAP 12 (2017) 044 [arXiv:1710.02865] [insPIRE].

[61] A. Kusenko et al., *Exploring primordial black holes from the multiverse with optical telescopes*, Phys. Rev. Lett. 125 (2020) 181304 [arXiv:2001.09160] [insPIRE].

[62] A. Ashoorioon, A. Rostami and J.T. Firouzjaee, *Examining the end of inflation with primordial black holes mass distribution and gravitational waves*, Phys. Rev. D 103 (2021) 123512 [arXiv:2012.02817] [insPIRE].

[63] C. Animali and V. Vennin, *Primordial black holes from stochastic tunnelling*, JCAP 02 (2023) 043 [arXiv:2210.03812] [insPIRE].

[64] E. Witten, *Cosmic separation of phases*, Phys. Rev. D 30 (1984) 272 [insPIRE].

[65] C.J. Hogan, *Gravitational radiation from cosmological phase transitions*, Mon. Not. Roy. Astron. Soc. 218 (1986) 629 [insPIRE].

[66] A. Kosowsky, M.S. Turner and R. Watkins, *Gravitational waves from first order cosmological phase transitions*, Phys. Rev. Lett. 69 (1992) 2026 [insPIRE].

[67] A. Kosowsky, M.S. Turner and R. Watkins, *Gravitational radiation from colliding vacuum bubbles*, Phys. Rev. D 45 (1992) 4514 [insPIRE].

[68] A. Kosowsky and M.S. Turner, *Gravitational radiation from colliding vacuum bubbles: envelope approximation to many bubble collisions*, Phys. Rev. D 47 (1993) 4372 [astro-ph/9211004] [insPIRE].

[69] M. Kamionkowski, A. Kosowsky and M.S. Turner, *Gravitational radiation from first order phase transitions*, Phys. Rev. D 49 (1994) 2837 [astro-ph/9310044] [insPIRE].

[70] C. Caprini, R. Durrer and G. Servant, *Gravitational wave generation from bubble collisions in first-order phase transitions: an analytic approach*, Phys. Rev. D 77 (2008) 124015 [arXiv:0711.2593] [insPIRE].
[71] S.J. Huber and T. Konstandin, *Gravitational wave production by collisions: more bubbles*, JCAP 09 (2008) 022 [arXiv:0806.1828] [INSPIRE].

[72] R. Jinno and M. Takimoto, *Gravitational waves from bubble dynamics: beyond the envelope*, JCAP 01 (2019) 060 [arXiv:1707.03111] [INSPIRE].

[73] T. Konstandin, *Gravitational radiation from a bulk flow model*, JCAP 03 (2018) 047 [arXiv:1712.06869] [INSPIRE].

[74] D. Cutting, M. Hindmarsh and D.J. Weir, *Gravitational waves from vacuum first-order phase transitions: from the envelope to the lattice*, Phys. Rev. D 97 (2018) 123513 [arXiv:1802.05712] [INSPIRE].

[75] D. Cutting, E.G. Escartin, M. Hindmarsh and D.J. Weir, *Gravitational waves from vacuum first order phase transitions II: from thin to thick walls*, Phys. Rev. D 103 (2021) 023531 [arXiv:2005.13537] [INSPIRE].

[76] M. Lewicki and V. Vaskonen, *Gravitational wave spectra from strongly supercooled phase transitions*, Eur. Phys. J. C 80 (2020) 1003 [arXiv:2007.04967] [INSPIRE].

[77] M. Lewicki and V. Vaskonen, *Gravitational waves from colliding vacuum bubbles in gauge theories*, Eur. Phys. J. C 81 (2021) 437 [Erratum ibid. 81 (2021) 1077] [arXiv:2012.07826] [INSPIRE].

[78] H. Zhong, B. Gong and T. Qiu, *Gravitational waves from bubble collisions in FLRW spacetime*, JHEP 02 (2022) 077 [arXiv:2107.01845] [INSPIRE].

[79] T. Bringmann et al., *Does NANOGrav observe a dark sector phase transition?*, JCAP 11 (2023) 053 [arXiv:2306.09411] [INSPIRE].

[80] LISA collaboration, *Laser Interferometer Space Antenna*, arXiv:1702.00786 [INSPIRE].

[81] M. Maggiore et al., *Science case for the Einstein Telescope*, JCAP 03 (2020) 050 [arXiv:1912.02622] [INSPIRE].

[82] S. Kawamura et al., *Current status of space gravitational wave antenna DECIGO and B-DECIGO*, PTEP 2021 (2021) 05A105 [arXiv:2006.13545] [INSPIRE].

[83] M. Evans et al., *A horizon study for cosmic explorer: science, observatories, and community*, arXiv:2109.09882 [INSPIRE].

[84] AEDGE collaboration, *AEDGE: Atomic Experiment for Dark matter and Gravity Exploration in space*, EPJ Quant. Technol. 7 (2020) 6 [arXiv:1908.00802] [INSPIRE].

[85] J. Liu et al., *Constraining first-order phase transitions with curvature perturbations*, Phys. Rev. Lett. 130 (2023) 051001 [arXiv:2208.14086] [INSPIRE].

[86] R. Jinno and M. Takimoto, *Probing a classically conformal B-L model with gravitational waves*, Phys. Rev. D 95 (2017) 015020 [arXiv:1604.05035] [INSPIRE].

[87] T. Hambye, A. Strumia and D. Teresi, *Super-cool dark matter*, JHEP 08 (2018) 188 [arXiv:1805.01473] [INSPIRE].

[88] C. Marzo, L. Marzola and V. Vaskonen, *Phase transition and vacuum stability in the classically conformal B-L model*, Eur. Phys. J. C 79 (2019) 601 [arXiv:1811.11169] [INSPIRE].

[89] K.-P. Xie, *Pinning down the primordial black hole formation mechanism with gamma-rays and gravitational waves*, JCAP 06 (2023) 008 [arXiv:2301.02352] [INSPIRE].

[90] I.K. Banerjee and U.K. Dey, *Probing the origin of primordial black holes through novel gravitational wave spectrum*, JCAP 07 (2023) 024 [arXiv:2305.07569] [INSPIRE].
[91] Y. Gouttenoire, *Primordial black holes from conformal Higgs*, arXiv:2311.13640 [insPIRE].

[92] S. Iso, N. Okada and Y. Orikasa, *Classically conformal B-L extended Standard Model*, Phys. Lett. B 676 (2009) 81 [arXiv:0902.4050] [insPIRE].

[93] S. Iso, N. Okada and Y. Orikasa, *The minimal B-L model naturally realized at TeV scale*, Phys. Rev. D 80 (2009) 115007 [arXiv:0909.0128] [insPIRE].

[94] R. Davis Jr., D.S. Harmer and K.C. Hoffman, *Search for neutrinos from the sun*, Phys. Rev. Lett. 20 (1968) 1205 [insPIRE].

[95] SNO collaboration, *Measurement of the rate of $\nu_e + d \to p + p + e^-$ interactions produced by $^8$B solar neutrinos at the Sudbury Neutrino Observatory*, Phys. Rev. Lett. 87 (2001) 071301 [nucl-ex/0106015] [insPIRE].

[96] SUPER-KAMIOKANDE collaboration, *Evidence for oscillation of atmospheric neutrinos*, Phys. Rev. Lett. 81 (1998) 1562 [hep-ex/9807003] [insPIRE].

[97] K2K collaboration, *Measurement of neutrino oscillation by the K2K experiment*, Phys. Rev. D 74 (2006) 072003 [hep-ex/0606032] [insPIRE].

[98] MINOS collaboration, *First direct observation of muon antineutrino disappearance*, Phys. Rev. Lett. 107 (2011) 021801 [arXiv:1104.0344] [insPIRE].

[99] DAYA BAY collaboration, *Observation of electron-antineutrino disappearance at Daya Bay*, Phys. Rev. Lett. 108 (2012) 171803 [arXiv:1203.1669] [insPIRE].

[100] RENO collaboration, *Observation of reactor electron antineutrino disappearance in the RENO experiment*, Phys. Rev. Lett. 108 (2012) 191802 [arXiv:1204.0626] [insPIRE].

[101] I. Esteban et al., *The fate of hints: updated global analysis of three-flavor neutrino oscillations*, JHEP 09 (2020) 178 [arXiv:2007.14792] [insPIRE].

[102] P. Minkowski, $\mu \to e\gamma$ at a rate of one out of $10^9$ muon decays?, Phys. Lett. B 67 (1977) 421 [insPIRE].

[103] T. Yanagida, *Horizontal gauge symmetry and masses of neutrinos*, Conf. Proc. C 7902131 (1979) 95 [insPIRE].

[104] M. Gell-Mann, P. Ramond and R. Slansky, *Complex spinors and unified theories*, Conf. Proc. C 790927 (1979) 315 [arXiv:1306.4669] [insPIRE].

[105] R.N. Mohapatra and G. Senjanovic, *Neutrino mass and spontaneous parity nonconservation*, Phys. Rev. Lett. 44 (1980) 912 [insPIRE].

[106] E. Gildener and S. Weinberg, *Symmetry breaking and scalar bosons*, Phys. Rev. D 13 (1976) 3333 [insPIRE].

[107] P.B. Arnold and O. Espinosa, *The effective potential and first order phase transitions: beyond leading-order*, Phys. Rev. D 47 (1993) 3546 [Erratum ibid. 50 (1994) 6662] [hep-ph/9212235] [insPIRE].

[108] S.R. Coleman, *The fate of the false vacuum. 1. Semiclassical theory*, Phys. Rev. D 15 (1977) 2929 [Erratum ibid. 16 (1977) 1248] [insPIRE].

[109] C.G. Callan Jr. and S.R. Coleman, *The fate of the false vacuum. 2. First quantum corrections*, Phys. Rev. D 16 (1977) 1762 [insPIRE].

[110] A.D. Linde, *Decay of the false vacuum at finite temperature*, Nucl. Phys. B 216 (1983) 421 [Erratum ibid. 223 (1983) 544] [insPIRE].
[111] D. Croon et al., Theoretical uncertainties for cosmological first-order phase transitions, JHEP 04 (2021) 055 [arXiv:2009.10080] [INSPIRE].

[112] A.H. Guth and S.H.H. Tye, Phase transitions and magnetic monopole production in the very early universe, Phys. Rev. Lett. 44 (1980) 631 [Erratum ibid. 44 (1980) 963] [INSPIRE].

[113] A.H. Guth and E.J. Weinberg, Cosmological consequences of a first order phase transition in the SU(5) grand unified model, Phys. Rev. D 23 (1981) 876 [INSPIRE].

[114] A.H. Guth and E.J. Weinberg, Could the universe have recovered from a slow first order phase transition?, Nucl. Phys. B 212 (1983) 321 [INSPIRE].

[115] K. Enqvist, J. Ignatius, K. Kajantie and K. Rummukainen, Nucleation and bubble growth in a first order cosmological electroweak phase transition, Phys. Rev. D 45 (1992) 3415 [INSPIRE].

[116] M.S. Turner, E.J. Weinberg and L.M. Widrow, Bubble nucleation in first order inflation and other cosmological phase transitions, Phys. Rev. D 46 (1992) 2384 [INSPIRE].

[117] J. Ellis, M. Lewicki and J.M. No, On the maximal strength of a first-order electroweak phase transition and its gravitational wave signal, JCAP 04 (2019) 003 [arXiv:1809.08242] [INSPIRE].

[118] I. Musco and T. Papanikolaou, Primordial black hole formation for an anisotropic perfect fluid: initial conditions and estimation of the threshold, Phys. Rev. D 106 (2022) 083017 [arXiv:2110.05982] [INSPIRE].

[119] R.K. Sheth, H.J. Mo and G. Tormen, Ellipsoidal collapse and an improved model for the number and spatial distribution of dark matter haloes, Mon. Not. Roy. Astron. Soc. 323 (2001) 1 [astro-ph/9907024] [INSPIRE].

[120] T.W. Baumgarte and P.J. Montero, Critical phenomena in the aspherical gravitational collapse of radiation fluids, Phys. Rev. D 92 (2015) 124065 [arXiv:1509.08730] [INSPIRE].

[121] K. Clough and E.A. Lim, Critical phenomena in non-spherically symmetric scalar bubble collapse, arXiv:1602.02568 [INSPIRE].

[122] F. Kühnel and M. Sandstad, Ellipsoidal collapse and primordial black hole formation, Phys. Rev. D 94 (2016) 063514 [arXiv:1602.04815] [INSPIRE].

[123] M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, Primordial black holes — perspectives in gravitational wave astronomy, Class. Quant. Grav. 35 (2018) 063001 [arXiv:1801.05235] [INSPIRE].

[124] T. Fujita, M. Kawasaki, K. Harigaya and R. Matsuda, Baryon asymmetry, dark matter, and density perturbation from primordial black holes, Phys. Rev. D 89 (2014) 103501 [arXiv:1401.1909] [INSPIRE].

[125] M.W. Choptuik, Universality and scaling in gravitational collapse of a massless scalar field, Phys. Rev. Lett. 70 (1993) 9 [INSPIRE].

[126] J.C. Niemeyer and K. Jedamzik, Near-critical gravitational collapse and the initial mass function of primordial black holes, Phys. Rev. Lett. 80 (1998) 5481 [astro-ph/9709072] [INSPIRE].

[127] A.M. Green and A.R. Liddle, Critical collapse and the primordial black hole initial mass function, Phys. Rev. D 60 (1999) 063509 [astro-ph/9901268] [INSPIRE].

[128] I. Musco, J.C. Miller and A.G. Polnarev, Primordial black hole formation in the radiative era: investigation of the critical nature of the collapse, Class. Quant. Grav. 26 (2009) 235001 [arXiv:0811.1452] [INSPIRE].
[129] F. Kühnel, C. Rampf and M. Sandstad, *Effects of critical collapse on primordial black-hole mass spectra*, Eur. Phys. J. C **76** (2016) 93 [arXiv:1512.00488] [INSPIRE].

[130] E. Witten, *Cosmological consequences of a light Higgs boson*, Nucl. Phys. B **177** (1981) 477 [INSPIRE].

[131] S. Iso, P.D. Serpico and K. Shimada, *QCD-electroweak first-order phase transition in a supercooled universe*, Phys. Rev. Lett. **119** (2017) 141301 [arXiv:1704.04955] [INSPIRE].

[132] B. von Harling and G. Servant, *QCD-induced electroweak phase transition*, JHEP **01** (2018) 159 [arXiv:1711.11554] [INSPIRE].

[133] L. Sagunski, P. Schicho and D. Schmitt, *Supercool exit: gravitational waves from QCD-triggered conformal symmetry breaking*, Phys. Rev. D **107** (2023) 123512 [arXiv:2303.02450] [INSPIRE].

[134] G.C. Dorsch, S.J. Huber, K. Mimasu and J.M. No, *The Higgs vacuum uplifted: revisiting the electroweak phase transition with a second Higgs doublet*, JHEP **12** (2017) 086 [arXiv:1705.09186] [INSPIRE].

[135] R. Durrer and C. Caprini, *Primordial magnetic fields and causality*, JCAP **11** (2003) 010 [astro-ph/0305059] [INSPIRE].

[136] C. Caprini, R. Durrer, T. Konstandin and G. Servant, *General properties of the gravitational wave spectrum from phase transitions*, Phys. Rev. D **79** (2009) 083519 [arXiv:0901.1661] [INSPIRE].

[137] G. Barenboim and W.-I. Park, *Gravitational waves from first order phase transitions as a probe of an early matter domination era and its inverse problem*, Phys. Lett. B **759** (2016) 430 [arXiv:1605.03781] [INSPIRE].

[138] R.-G. Cai, S. Pi and M. Sasaki, *Universal infrared scaling of gravitational wave background spectra*, Phys. Rev. D **102** (2020) 083528 [arXiv:1909.13728] [INSPIRE].

[139] A. Hook, G. Marques-Tavares and D. Racco, *Causal gravitational waves as a probe of free streaming particles and the expansion of the universe*, JHEP **02** (2021) 117 [arXiv:2010.03568] [INSPIRE].

[140] I. Baldes, Y. Gouttenoire, F. Sala and G. Servant, *Supercool composite dark matter beyond 100 TeV*, JHEP **07** (2022) 084 [arXiv:2110.13926] [INSPIRE].

[141] B. Allen and J.D. Romano, *Detecting a stochastic background of gravitational radiation: signal processing strategies and sensitivities*, Phys. Rev. D **59** (1999) 102001 [gr-qc/9710117] [INSPIRE].

[142] H. Kudoh, A. Taruya, T. Hiramatsu and Y. Himemoto, *Detecting a gravitational-wave background with next-generation space interferometers*, Phys. Rev. D **73** (2006) 064006 [gr-qc/0511145] [INSPIRE].

[143] E. Thrane and J.D. Romano, *Sensitivity curves for searches for gravitational-wave backgrounds*, Phys. Rev. D **88** (2013) 124032 [arXiv:1310.5300] [INSPIRE].

[144] C. Caprini et al., *Reconstructing the spectral shape of a stochastic gravitational wave background with LISA*, JCAP **11** (2019) 017 [arXiv:1906.09244] [INSPIRE].

[145] D. Brzeminski, A. Hook and G. Marques-Tavares, *Precision early universe cosmology from stochastic gravitational waves*, JHEP **11** (2022) 061 [arXiv:2203.13842] [INSPIRE].

[146] T. Robson, N.J. Cornish and C. Liu, *The construction and use of LISA sensitivity curves*, Class. Quant. Grav. **36** (2019) 105011 [arXiv:1803.01944] [INSPIRE].
[147] K. Schmitz, *LISA sensitivity to gravitational waves from sound waves*, *Symmetry* **12** (2020) 1477 [arXiv:2005.10789] [INSPIRE].

[148] A. Lamberts et al., *Predicting the LISA white dwarf binary population in the Milky Way with cosmological simulations*, *Mon. Not. Roy. Astron. Soc.* **490** (2019) 5888 [arXiv:1907.00014] [INSPIRE].

[149] G. Boileau et al., *Prospects for LISA to detect a gravitational-wave background from first order phase transitions*, *JCAP* **02** (2023) 056 [arXiv:2209.13277] [INSPIRE].

[150] A.J. Farmer and E.S. Phinney, *The gravitational wave background from cosmological compact binaries*, *Mon. Not. Roy. Astron. Soc.* **346** (2003) 1197 [astro-ph/0304393] [INSPIRE].

[151] P.A. Rosado, *Gravitational wave background from binary systems*, *Phys. Rev. D* **84** (2011) 084004 [arXiv:1106.5795] [INSPIRE].

[152] KAGRA et al. collaborations, *Upper limits on the isotropic gravitational-wave background from advanced LIGO and advanced Virgo's third observing run*, *Phys. Rev. D* **104** (2021) 022004 [arXiv:2101.12130] [INSPIRE].

[153] E.S. Phinney, *A practical theorem on gravitational wave backgrounds*, astro-ph/0108028 [INSPIRE].

[154] T. Biekötter et al., *Fate of electroweak symmetry in the early universe: non-restoration and trapped vacua in the N2HDM*, *JCAP* **06** (2021) 018 [arXiv:2103.12707] [INSPIRE].

[155] T. Biekötter et al., *The trap in the early universe: impact on the interplay between gravitational waves and LHC physics in the 2HDM*, *JCAP* **03** (2023) 031 [arXiv:2208.14466] [INSPIRE].

[156] KAGRA et al. collaborations, *Upper limits on the isotropic gravitational-wave background from advanced LIGO and advanced Virgo's third observing run*, *Phys. Rev. D* **104** (2021) 022004 [arXiv:2101.12130] [INSPIRE].

[157] C. Badger et al., *Probing early universe supercooled phase transitions with gravitational wave data*, *Phys. Rev. D* **107** (2023) 023511 [arXiv:2209.14707] [INSPIRE].

[158] W. Buchmuller, P. Di Bari and M. Plumacher, *Leptogenesis for pedestrians*, *Annals Phys.* **315** (2005) 305 [hep-ph/0401240] [INSPIRE].

[159] D. Bodeker and G.D. Moore, *Can electroweak bubble walls run away?*, *JCAP* **05** (2009) 009 [arXiv:0903.4099] [INSPIRE].

[160] D. Bodeker and G.D. Moore, *Electroweak bubble wall speed limit*, *JCAP* **05** (2017) 025 [arXiv:1703.08215] [INSPIRE].

[161] R. Jinno, H. Seong, M. Takimoto and C.M. Um, *Gravitational waves from first-order phase transitions: ultra-supercooled transitions and the fate of relativistic shocks*, *JCAP* **10** (2019) 033 [arXiv:1905.00899] [INSPIRE].

[162] I. Baldes, M. Dichtl, Y. Gouttenoire and F. Sala, *Bubbletrons*, arXiv:2306.15555 [INSPIRE].

[163] A. Azatov and M. Vanvlasselaer, *Bubble wall velocity: heavy physics effects*, *JCAP* **01** (2021) 058 [arXiv:2010.02590] [INSPIRE].

[164] Y. Gouttenoire, R. Jinno and F. Sala, *Friction pressure on relativistic bubble walls*, JHEP **05** (2022) 004 [arXiv:2112.07686] [INSPIRE].

[165] T. Prokopec, J. Rezacek and B. Świeżewska, *Gravitational waves from conformal symmetry breaking*, *JCAP* **02** (2019) 009 [arXiv:1809.11129] [INSPIRE].
[166] M. Kierkla, A. Karam and B. Swiezewska, Conformal model for gravitational waves and dark matter: a status update, JHEP 03 (2023) 007 [arXiv:2210.07075] [INSPIRE].

[167] I. Baldes et al., Baryogenesis via relativistic bubble expansion, Phys. Rev. D 104 (2021) 115029 [arXiv:2106.15602] [INSPIRE].

[168] P. Huang and K.-P. Xie, Leptogenesis triggered by a first-order phase transition, JHEP 09 (2022) 052 [arXiv:2206.04691] [INSPIRE].

[169] A. Dasgupta, P.S.B. Dev, A. Ghoshal and A. Mazumdar, Gravitational wave pathway to testable leptogenesis, Phys. Rev. D 106 (2022) 075027 [arXiv:2206.07032] [INSPIRE].

[170] D. Borah, A. Dasgupta and I. Saha, Leptogenesis and dark matter through relativistic bubble walls with observable gravitational waves, JHEP 11 (2022) 136 [arXiv:2207.14226] [INSPIRE].

[171] E.J. Chun et al., Bubble-assisted leptogenesis, JHEP 09 (2023) 164 [arXiv:2305.10759] [INSPIRE].

[172] A. Azatov, M. Vanylasselaer and W. Yin, Baryogenesis via relativistic bubble walls, JHEP 10 (2021) 043 [arXiv:2106.14913] [INSPIRE].

[173] J.M. Cline and K. Kainulainen, Electroweak baryogenesis at high bubble wall velocities, Phys. Rev. D 101 (2020) 063525 [arXiv:2001.00568] [INSPIRE].

[174] A. Masiero and A. Riotto, Cosmic $\Delta B$ from lepton violating interactions at the electroweak phase transition, Phys. Lett. B 289 (1992) 73 [hep-ph/9206212] [INSPIRE].

[175] A. Katz and A. Riotto, Baryogenesis and gravitational waves from runaway bubble collisions, JCAP 11 (2016) 011 [arXiv:1608.00583] [INSPIRE].

[176] A. Falkowski and J.M. No, Non-thermal dark matter production from the electroweak phase transition: multi-TeV WIMPs and ‘baby-zillas’, JHEP 02 (2013) 034 [arXiv:1211.5615] [INSPIRE].

[177] M. Flanz, E.A. Paschos and U. Sarkar, Baryogenesis from a lepton asymmetric universe, Phys. Lett. B 345 (1995) 248 [hep-ph/9411366] [INSPIRE].

[178] M. Flanz, E.A. Paschos, U. Sarkar and J. Weiss, Baryogenesis through mixing of heavy Majorana neutrinos, Phys. Lett. B 389 (1996) 693 [hep-ph/9607310] [INSPIRE].

[179] A. Pilaftsis, CP violation and baryogenesis due to heavy Majorana neutrinos, Phys. Rev. D 56 (1997) 5431 [hep-ph/9707235] [INSPIRE].

[180] A. Pilaftsis and T.E.J. Underwood, Resonant leptogenesis, Nucl. Phys. B 692 (2004) 303 [hep-ph/0309342] [INSPIRE].

[181] S. Davidson and A. Ibarra, A lower bound on the right-handed neutrino mass from leptogenesis, Phys. Lett. B 535 (2002) 25 [hep-ph/0202239] [INSPIRE].

[182] T. Hambye, Leptogenesis at the TeV scale, Nucl. Phys. B 633 (2002) 171 [hep-ph/0111089] [INSPIRE].

[183] A. Boyarsky et al., Sterile neutrino dark matter, Prog. Part. Nucl. Phys. 104 (2019) 1 [arXiv:1807.07938] [INSPIRE].

[184] G. Arcadi et al., Probing a dark sector with collider physics, direct detection, and gravitational waves, Phys. Lett. B 848 (2024) 138382 [arXiv:2307.06376] [INSPIRE].

[185] I. Baldes and C. Garcia-Cely, Strong gravitational radiation from a simple dark matter model, JHEP 05 (2019) 190 [arXiv:1809.01198] [INSPIRE].
[186] I. Baldes, Y. Gouttenoire and F. Sala, String fragmentation in supercooled confinement and implications for dark matter, JHEP 04 (2021) 278 [arXiv:2007.08440] [INSPIRE].

[187] X.-R. Wong and K.-P. Xie, Freeze-in of WIMP dark matter, Phys. Rev. D 108 (2023) 055035 [arXiv:2304.00908] [INSPIRE].

[188] A. Ray, R. Laha, J.B. Muñoz and R. Caputo, Near future MeV telescopes can discover asteroid-mass primordial black hole dark matter, Phys. Rev. D 104 (2021) 023516 [arXiv:2102.06714] [INSPIRE].

[189] K. Agashe et al., Correlating gravitational wave and gamma-ray signals from primordial black holes, Phys. Rev. D 105 (2022) 123009 [arXiv:2202.04653] [INSPIRE].

[190] S. Clark et al., 21 cm limits on decaying dark matter and primordial black holes, Phys. Rev. D 98 (2018) 043006 [arXiv:1803.09390] [INSPIRE].

[191] A.K. Saha and R. Laha, Sensitivities on nonspinning and spinning primordial black hole dark matter with global 21 cm troughs, Phys. Rev. D 105 (2022) 103026 [arXiv:2112.10794] [INSPIRE].

[192] G. Franciolini, Primordial black holes: from theory to gravitational wave observations, Ph.D. thesis, Geneva U., Geneva, Switzerland (2021) [arXiv:2110.06815] [INSPIRE].