# Primordial Black Holes and Wormholes from Domain Wall Networks

Yann Gouttenoire  0000-0003-2225-6704<EMAIL_ADDRESS>School of
Physics and Astronomy, Tel-Aviv University, Tel-Aviv 69978, Israel  Edoardo
Vitagliano  0000-0001-7847-1281<EMAIL_ADDRESS>Racah
Institute of Physics, Hebrew University of Jerusalem, Jerusalem 91904, Israel

###### Abstract

Domain walls (DWs) are topological defects originating from phase transitions
in the early universe. In the presence of an energy imbalance between distinct
vacua, enclosed DW cavities shrink until the entire network disappears. By
studying the dynamics of thin-shell bubbles in General Relativity, we
demonstrate that closed DWs with sizes exceeding the cosmic horizon tend to
annihilate later than the average. This delayed annihilation allows for the
formation of large overdensities, which, upon entering the Hubble horizon,
eventually collapse to form Primordial Black Holes (PBHs). We rely on 3D
percolation theory to calculate the number density of these late-annihilating
DWs, enabling us to infer the abundance of PBHs. A key insight from our study
is that DW networks with the potential to emit observable Gravitational Waves
are also likely to yield detectable PBHs. Additionally, we study the
production of wormholes connected to baby-universes and conclude on the
possibility to generate a multiverse.

## I. INTRODUCTION

The study of primordial black holes (PBHs) has become a vibrant field of
exploration since the 2015 discovery of gravitational waves from solar-mass
black holes mergers [1]. Any observation of black holes with masses below that
of the Sun would be a smoking-gun for the gravitational collapse of large
density fluctuation, present in the initial cosmic plasma [2]. An intriguing
possibility involves the inhomogeneities generated by domain walls (DWs)
either under the form of networks [3, 4, 5, 6] or nucleated during inflation
[7, 8, 9, 10, 11, 12].

DW networks are topological defects that form during cosmological phase
transitions in the early universe when nearly-degenerate vacua are present.
Shortly after its formation, the DW network evolves in a scaling regime during
which the correlation length approximately equals the horizon size $L\sim t$
[13], and the fraction of the total universe energy density stored in DWs
increases linearly with time $\rho_{\rm DW}/\rho_{\rm tot}\propto t$ [13]. A
universe dominated by DWs would lead to a universe either primarily filled
with black holes or to eternally-inflating universes, depending on whether the
observer location lies within a true vacuum or false vacuum region [14]. A
graceful exit from this daunting scenario could occur if there is an energy
bias $V_{\rm bias}$ between the different vacua, as we review in Sec. II.
After a time $t_{\rm ann}$, the vacuum energy difference $V_{\rm bias}$
counterbalances the pressure due to the wall tension $\sigma$, driving DWs
towards annihilation before they can dominate the universe at a time $t_{\rm
dom}$ [15, 16, 17, 18]. During this annihilation phase, closed DWs shrink and
under specific conditions, could enter within their Schwarzschild radius and
form PBHs [4, 7, 8, 9], a process dubbed “catastrogenesis” in [5, 6]. As
explained above, if $t_{\rm ann}\gtrsim t_{\rm dom}$, the universe is filled
with PBHs and wormholes connected to baby-universes [14]. Per continuity, we
conclude that there must exist a region in the ballpark $t_{\rm ann}\lesssim
t_{\rm dom}$, where DWs networks produce just enough PBHs to explain dark
matter (DM), or just enough wormholes to generate a multiverse in our past
lightcone [19, 20, 21].

Figure 1:  Domain wall late-annihilation mechanism. During annihilation (Green
region), closed DWs with sub-horizon size shrink exponentially fast while
super-horizon DWs continue to grow due to Hubble flow. Closed DW
configurations with a radius larger than $R_{\rm ann}^{\rm PBH}$ at the onset
of the annihilation phase, which are dubbed “late-annihilators” in the Orange
region, enter inside their Schwarzschild radius $R\lesssim 2GM$ in the Blue
region, leading them to collapse into PBHs. DWs with radius larger than
$R_{\rm ann}^{\rm baby}$ become larger than the cosmological horizon
$H_{b}^{-1}$ of an inflating universe with same vacuum energy $V_{\rm bias}$
as their interior, in the Purple region, with $H_{b}^{2}=V_{\rm bias}/3M_{\rm
Pl}^{2}$. They become eternally-inflating baby-universes detaching from the
parent universe through wormholes and leaving PBHs behind.

Correcting for mistakes of previous work [4, 5, 6], we calculate the abundance
of PBHs, and for the first time the number of those which hid behind their
horizon a wormhole connected to an eternally-inflating baby-universe. Previous
studies [4, 5, 6] incorrectly assumed that the radius of DWs grows as $R\simeq
t$. Using energy arguments, Ref. [22] commented that the relation $R\simeq t$
is not applicable during the annihilation phase. In this study, we provide a
definitive resolution to the issue, stating that $R\simeq t$ is never valid,
neither during the annihilation phase nor at any other period. In Sec. III, we
solve the equation of motion (EoM) of a single spherical DW in an expanding
universe, derived from Einstein equations in the thin-shell limit. The main
results are sketched in Fig. 1. The radius grows as $R\propto a$ or slower if
superhorizon $R\gtrsim t$, but quickly shrink if sub-horizon $R\lesssim t$. We
find that PBHs enter their Schwarschild radius if they have a radius larger
than $R\gtrsim R_{\rm ann}^{\rm PBH}>t$ just before the annihilation phase.
While most of DWs shrink and decay away into scalar waves, those DWs with
super-horizon size, which we call late-annihilators, continue to grow during
the annihilation phase until they finally collapse when entering inside the
cosmic horizon. We also derive a critical radius $R_{\rm ann}^{\rm baby}$
larger than $R_{\rm ann}^{\rm PBH}$ by about $30\%$ above which DWs collapse
into wormholes connected to eternally-inflating baby-universes. We find that
the PBH formation threshold $R_{\rm ann}^{\rm PBH}$ is highly sensitive to the
surface energy of DWs, which can have Lorentz factors up to $\gamma\sim 10$,
and to gravitational binding energies. These contributions, which constitute
$75\%$ of the PBH mass, were overlooked in prior studies [4, 5, 6]. In Sec.
IV, we employ for the first time percolation theory to estimate the fraction
of closed DWs which are large enough to host a false vacuum ball of radius
$R_{\rm min}$. In Sec. V, we deduce the abundance of PBHs and wormholes. In
Sec. VI, we review the Stochastic Gravitational Wave Background (SGWB)
produced by annihilating DW networks. We point the complementarity between GW
detectability by future observatories, PBHs, and baby-universe production.

## II. BIASED NETWORK EVOLUTION

Initially, the work of the DW surface tension $\sigma$, with equivalent
pressure $\mathcal{P}_{T}=C_{d}\,\sigma/L$ toward straightening, is dampened
by the friction pressure $\mathcal{P}_{V}\simeq\,\beta\,T^{4}$. The quantity
$L$ is the correlation length of the DW network, $C_{d}$ and $\beta$ are
dimensionless parameters setting respectively the strength of the
straightening pressure and of the DW-plasma interactions. Estimating DWs to
have typical curvature radius $L\simeq vt$, one obtains that DWs start moving
with relativistic velocity $v\simeq\mathcal{O}(0.1)$ after the time

$t_{\rm rel}\simeq 2\times
10^{-3}\left(\frac{v}{0.1}\right)\left(\frac{100}{g_{\star}}\right)\frac{\beta}{C_{d}}\frac{M_{\rm
Pl}^{2}}{\sigma}.$ (1)

where $M_{\rm Pl}=1/\sqrt{8\pi G}\simeq 2.435\times 10^{18}\,\rm GeV$ is the
reduced Planck mass, $g_{\star}$ the number of relativistic degrees of freedom
and where we used Friedmann’s equation $T=1.2\sqrt{M_{\rm
Pl}/t}/g_{\star}^{1/4}$. The magnitude of the friction coefficient, denoted by
$\beta$, varies depending on the particle physics model [23, 24]. In our
study, we assume $\beta\ll 1$ and briefly touch upon the implications of
friction on PBH formation towards the end of the paper in Sec. V.3. We assume
that the energy bias $V_{\text{bias}}$ is initially insignificant at the time
of DW formation and remains negligible when they begin to move without
friction. Under these conditions, numerical simulations have shown that the
energy density of DWs quickly reaches the scaling regime,

$\rho_{\rm DW}=\frac{\sigma}{L},\qquad L=t/\mathcal{A},$ (2)

where $\mathcal{A}$ is the area parameter. Based on numerical simulations, the
area parameter for a $\mathcal{Z}_{2}$ symmetric model is approximately
$\mathcal{A}\simeq 0.8\pm 0.1$ [25]. In contrast, $\mathcal{Z}_{N_{\rm DW}}$
symmetric models with $N_{\rm DW}>2$, where cosmic strings attached to $N_{\rm
DW}$ DWs are present, tend to yield larger values $\mathcal{A}=(0.37\pm
0.04)N_{\rm DW}$ [26]. DWs annihilate when the vacuum pressure
$\mathcal{P}_{V}=V_{\rm bias}$ dominates over the pressure
$\mathcal{P}_{T}=C_{d}\,\sigma/L$ arising from their surface tension $\sigma$,
after the time

$t_{\rm ann}\leavevmode\nobreak\ \simeq\leavevmode\nobreak\
C_{d}\,\mathcal{A}\,\frac{\sigma}{V_{\rm bias}},$ (3)

where the factor $C_{d}\simeq\mathcal{O}(\rm a\leavevmode\nobreak\ few)$ can
be inferred from numerical simulations [26] and which we set to $C_{d}\simeq
3$. To prevent entering the catastrophic scenario of DW domination, DW
annihilation must proceed before the network dominates the energy budget of
the universe, occurring when $3M_{\rm Pl}^{2}H^{2}\simeq V_{\rm bias}$ around
the time

$t_{\rm dom}\leavevmode\nobreak\ \simeq\leavevmode\nobreak\
\frac{\sqrt{3}M_{\rm Pl}}{2\sqrt{V_{\rm bias}}}\leavevmode\nobreak\
=\leavevmode\nobreak\ 1.5\times 10^{-11}{\leavevmode\nobreak\ \rm
s}\left(\frac{300\leavevmode\nobreak\ \rm GeV}{V^{1/4}_{\rm
bias}}\right)^{\\!2},$ (4)

which upon comparing to Eq. (3), leads to the condition

$V_{\rm bias}\leavevmode\nobreak\ \gtrsim\leavevmode\nobreak\
\frac{4}{3}\left(\mathcal{A\,}C_{d}\frac{\sigma}{M_{\rm
Pl}}\right)^{\\!2}\simeq(300\leavevmode\nobreak\ \rm
GeV)^{4}\mathcal{A\,}^{\\!2}C_{d}^{2}\left(\frac{\sigma^{1/3}}{58\leavevmode\nobreak\
\rm PeV}\right)^{\\!6}.$ (5)

In the absence of a bias, the DW network would dominate when $3M_{\rm
Pl}^{2}H^{2}\simeq\sigma/L$ after the time

$t_{\rm dom}^{{\rm unbias}}\leavevmode\nobreak\ =\leavevmode\nobreak\
\frac{3M_{\rm Pl}^{2}}{4\mathcal{A}\sigma}\leavevmode\nobreak\
=\leavevmode\nobreak\ \frac{1.5\times 10^{-11}{\leavevmode\nobreak\ \rm
s}}{\mathcal{A}}\left(\frac{58\leavevmode\nobreak\ \rm
PeV}{\sigma^{1/3}}\right)^{\\!3}.$ (6)

Figure 2:  LEFT: We show the time $t_{\rm PBH}$ in Eq. (18) at which a
spherical DW collapses into PBH, for different levels of approximation of the
DW mass, either accounting for volume (“bulk”), surface (“bdy”) and
gravitational binding (“bind”) terms in black, or neglecting binding term in
blue or neglecting both binding and surface terms in purple. The purple line,
where only the vacuum energy is included, also corresponds to the time $t_{\rm
baby}$ in Eq. (21) of formation of a baby-universe. RIGHT: Composition of the
DW mass at the time $t_{\rm PBH}$ when it enters its Schwarszchild radius. The
bulk component only accounts for $25\leavevmode\nobreak\ \%$ of the DW mass.
Another $25\%$ is given by gravitational binding energies $M_{\rm bind}$ and
$50\%$ of the DW mass is given by its surface term $M_{\rm bdy}$.

## III. TIME EVOLUTION OF A SPHERICAL DW

### III.1. Equation of motion

The DW network contains closed configurations of size $R(t)$ which we treat as
spherical for simplicity. We discuss the implication of such approximation
later in Sec. V.3. We denote by $R(t)=a(t)\chi(t)$ the physical radius of a
spherical DW and by $\chi(t)$ its comoving radius. We introduce the wall
Lorentz factor $\gamma=1/\sqrt{1-(a\dot{\chi})^{2}}$. The EoM of a spherical
DW can be determined through the Israel junction conditions [27, 28, 29, 30,
31, 7, 8, 9], named after Werner Israel [32]. These conditions equate the
discontinuity in the extrinsic curvature across the wall to the wall surface
energy-momentum tensor [33]. For a shrinking DW, one obtains [9]

$\ddot{\chi}+(4-3a^{2}\dot{\chi}^{2})H\dot{\chi}+\frac{2}{a^{2}\chi}(1-a^{2}\dot{\chi}^{2})\\\
=-\left(\frac{V_{\rm bias}}{\sigma}-6\pi
G\sigma\right)\frac{(1-a^{2}\dot{\chi}^{2})^{3/2}}{a},$ (7)

where $a(t)$ is the scale factor of the universe inside the DW,

$\dot{a}(t)=H(t)a(t),\qquad H^{2}(t)=\frac{\rho_{\rm rad}(t)+V_{\rm
bias}}{3M_{\rm Pl}^{2}},$ (8)

with $\rho_{\rm rad}(t)=\rho_{\rm rad}(t_{\rm ann})\left[a(t_{\rm
ann})/a(t)\right]^{4}$. Introducing the dimensionless quantities

$\tilde{\chi}(\tau)\equiv\chi(t/t_{\rm ann})/t_{\rm ann},\qquad\tau\equiv
t/t_{\rm ann},$ (9)

the EoM in Eq. (7) can be rewritten as

$\ddot{\tilde{\chi}}+(4-3a^{2}\dot{\tilde{\chi}}^{2})\tilde{H}\dot{\tilde{\chi}}+\frac{2}{a^{2}\tilde{\chi}}(1-a^{2}\dot{\tilde{\chi}}^{2})\\\
=-\left(\mathcal{A}C_{d}-\frac{9\alpha_{\rm
ann}}{16\mathcal{A}}\right)\frac{(1-a^{2}\dot{\tilde{\chi}}^{2})^{3/2}}{a},$
(10)

where the dot now denotes the derivative with respect to $\tau$, $\tilde{H}$
is defined below, and $\alpha_{\rm ann}$ is the DW network energy fraction at
$t_{\rm ann}$

$\alpha_{\rm ann}\equiv\frac{\rho_{\rm DW}}{\rho_{\rm tot}}\Big{|}_{t=t_{\rm
ann}}\simeq\left(\frac{t_{\rm ann}}{t_{\rm dom}^{{\rm
unbias}}}\right)^{2}\simeq C_{d}^{-1}\left(\frac{t_{\rm ann}}{t_{\rm
dom}}\right)^{2},$ (11)

where we used $\rho_{\rm tot}\simeq 3M_{\rm Pl}^{2}/4t^{2}$ and Eqs. (2, 3, 4,
6). Upon approximating $\rho_{\rm rad}(t_{\rm ann})\simeq\rho_{\rm tot}(t_{\rm
ann})$, which is exact if $t_{\rm dom}\gg t_{\rm ann}$, the dimensionless
Hubble factor reads

$\tilde{H}(\tau)=\sqrt{\left(a(t_{\rm ann})/a(t)\right)^{4}+C_{d}\alpha_{\rm
ann}^{2}}/2.$ (12)

### III.2. PBH collapse time

Closed DWs collapse into PBHs at a time $t_{\rm PBH}$ after that they shrink
below their Schwarzschild radius,

$R(t_{\rm PBH})=2GM(t_{\rm PBH}).$ (13)

The mass-energy $M$ contained within the spherical DW is given by the Misner-
Sharp mass [34, 35]. The later can be decomposed as [9]

$M=M_{\rm bulk}+M_{\rm bdy}+M_{\rm bind},$ (14)

where

$M_{\rm bulk}=\frac{4\pi}{3}V_{\rm bias}R^{3},$ (15)

is a bulk term,

$M_{\rm bdy}=4\pi\sigma\gamma R^{2}$ (16)

is a boundary term, and

$M_{\rm bind}=-8\pi^{2}G\sigma^{2}R^{3}+4\pi\sigma
H_{b}\sqrt{\gamma^{2}-1}R^{3},$ (17)

with $H_{b}=\sqrt{V_{\rm bias}/3M_{\rm Pl}^{2}}$, includes the repulsive
surface-surface gravitational binding energy [36] and the attractive surface-
volume gravitational binding energy, with the former typically being dominant.
From numerically integrating the EoM in Eq. (10) and solving for roots of Eq.
(13) assuming $R(t_{\rm PBH})=t_{\rm PBH}$ (PBHs cannot form faster than what
causality allows) we find that PBHs form after the time, see the left panel of
Fig. 2:

$t_{\rm PBH}\simeq\frac{1}{2}\sqrt{\frac{3M_{\rm Pl}^{2}}{V_{\rm
bias}}}=\frac{t_{\rm ann}}{\sqrt{C_{d}\alpha_{\rm ann}}}\simeq t_{\rm dom},$
(18)

which is around the time $t_{\rm dom}$ in Eq. (4) when DW interiors become
dominated by the vacuum energy $V_{\rm bias}$. If we had approximated the DW
mass by its volume term $M\simeq M_{\rm bulk}$, we would have found that PBHs
are twice longer to form $t_{\rm PBH}\simeq 2t_{\rm PBH}$, see the left panel
of Fig. 2. Such difference can impact the PBH abundance by many orders of
magnitude. In fact the bulk term composes only $25\%$ of the DW total mass at
the time of collapse, see the right panel of Fig. 2. The large surface term,
composing $50\%$ of the mass budget, can be attributed to the large Lorentz
factor shown in Fig. 3.

Figure 3:  Lorentz factor $\gamma=1/\sqrt{1-(a\dot{\chi})^{2}}$ of shrinking
DW at the time $t_{\rm PBH}$ of collapse into PBH.

Figure 4:  LEFT: Time evolution of the radius $R(t)$ of a spherical DW from
solving the full DW EoM for different initial sizes. DWs collapses into PBHs
if they cross their Schwarzschild radius in the blue region. The dark orange
line shows the trajectory of the smallest spherical DW collapsing into PBH.
The gray line shows the would-be DW trajectory from completely neglecting
effects from the bias $V_{\rm bias}$ and surface tension $\sigma$, i.e. solely
accounting for Hubble flow $R(t)\propto a(t)$. RIGHT: We show the minimal
radius that spherical DW must have at $t_{\rm ann}$ in order to collapse into
PBH after $t_{\rm PBH}$. The Schwarzschild radius is calculated with volume
term only (purple), plus surface term (“bdy”) (blue) and binding terms
(black). Due to the relationship $2GM_{\rm bulk}=H_{b}^{-1}$, the purple line
is also the critical radius to form a baby universe. The fitting formula for
black and purple lines are given in Eqs. (19) and (22) respectively.

### III.3. Late-annihilators

Typical DWs annihilate at $t_{\rm ann}$ in Eq. (3) when the volume term in Eq.
(15) dominates over the surface term in Eq. (16). This argument however
neglect Hubble expansion which becomes important for $R\gtrsim t$. In the
absence of analytical solution, we solve Eqs. (10) numerically and display the
DW trajectory in the left panel of Fig. 4. We find that DWs grow when their
radius is super-horizon $R>t$, and quickly shrink under the work of their own
surface tension and vacuum pressure when entering the horizon $R<t$. The
trajectory of the smallest spherical DWs which can collapse into PBHs is shown
in dark orange. It enters its Schwarzschild radius at Hubble crossing $R\simeq
t$. Larger DWs enter their Schwarzschild radius while $R>t$ (blue region) but
they only start to collapse after they become causally-connected when $R\simeq
t$ (black line). The radius $R_{\rm ann}^{\rm PBH}$ at $t_{\rm ann}$ of the
smallest DW which collapse into PBH is plotted in the right panel of Fig. 4 as
a function of $\alpha_{\rm ann}$, for different approximations for the Misner-
Sharp mass in Eq. (14). We also show in gray the approximation $R(t)\propto
t^{1/2}$ only accounting for the Hubble flow and neglecting effects from the
bias $V_{\rm bias}$ and surface tension $\sigma$ which become important when
$R$ approaches $R\sim t$. The most precise approximation shown with the black
line is well reproduced with the fitting function:

$\frac{R_{\rm ann}^{\rm PBH}}{t_{\rm ann}}\simeq 0.780\,{\rm
Log_{10}^{2}}\left(\alpha_{\rm ann}\right)-0.618\,{\rm
Log_{10}}\left(\alpha_{\rm ann}\right)+0.407,$ (19)

Closed DW configurations with a radius larger than $R\gtrsim R_{\rm ann}^{\rm
PBH}$ at the onset of the annihilation phase at $t_{\rm ann}$ are the
population of late-annihilators which collapse into PBHs after $t_{\rm
PBH}$.111Their existence was postulated in [4, 37] under the name of late-
birds. They are the DW network analogs of late-bloomers which are formed in
the context of first-order phase transitions [38, 39].

### III.4. Baby-universe

If the radius $R$ of a DW becomes larger than the cosmological horizon
$H_{b}^{-1}$ of an inflating universe with vacuum energy density $V_{\rm
bias}$, which would occur around the time $t_{\rm baby}$

$R(t_{\rm baby})\simeq H_{b}^{-1},\qquad H_{b}\equiv\sqrt{V_{\rm bias}/3M_{\rm
Pl}^{2}},$ (20)

then the DW radius grows exponentially with time, see the left panel of Fig.
4.This leads to the formation of an eternally inflating universe with a vacuum
energy denoted by $V_{\rm bias}$ called child-universe in [40] and baby-
universe in [41, 42].

The baby-universe detaches from the parent universe through a wormhole, which
closes off within a timescale comparable to its light crossing time. Observers
outside the DW will continue to experience the power-law Hubble expansion
characteristic of a radiation-dominated universe. They will perceive the DW as
shrinking until it collapses into a PBH, see e.g. [19, 7]. The observation of
the DW collapsing into a PBH, while at the same time the coordinate $R(t)$ is
exponentially growing, appears paradoxical. This seeming contradiction finds
resolution upon recognizing that the growing DW radius $R(t)$ is in fact
defined in the parallel exterior region of the maximally-extended
Schwarszchild space-time diagram [29]. Consequently, a wormhole forms between
the growing DW and the external observer, as illustrated by the space-time
diagrams C, D and E in [29]. The surface area along the wormhole throat
interpolates between the growing DW surface observed from the inside, and the
surface area perceived from the outside, as depicted in Fig. 11 of [7]. The
same figure shows that the throat quickly pinches off, forming a BH horizon on
both sides, one in the inflating baby-universe and one in the parent
radiation-dominated universe. Another consequence of the wormhole geometry,
demonstrated in App. C of [29], is that in spite of having a growing radius
$R(t)$, the DW is subject to a proper acceleration which is oriented inward,
toward the expanding baby-universe.

Numerically, we find that wormholes leading to the baby-universe form roughly
around twice the time of PBH formation, as given in Eq. (18):

$t_{\rm baby}\simeq 2t_{\rm PBH}\simeq 2t_{\rm dom}.$ (21)

Since $2GM_{\rm bulk}=H_{b}^{2}R^{3}$ in Eq. (14), the time $t_{\rm baby}$
actually corresponds to the purple line in the left panel of Fig. 2. We find
that the minimal DW radius at the onset of the annihilation phase for the
latter to collapse into a wormhole is well approximated by the fitting
function:

$\frac{R_{\rm ann}^{\rm baby}}{t_{\rm ann}}\simeq 0.951\,{\rm
Log_{10}^{2}}\left(\alpha_{\rm ann}\right)-0.860\,{\rm
Log_{10}}\left(\alpha_{\rm ann}\right)+1.15,$ (22)

which for the same reason as before corresponds to the purple line in the
right panel of Fig. 4. The abundance of late-annihilators, whose size are
larger than either Eq. (19) or Eq. (22), depending on whether they form a
simple PBH or a PBH + wormhole, is determined in the next section using
percolation theory.

## IV. PERCOLATION THEORY IN 3D

The goal of this section is to estimate the size distribution of closed DWs.
In anticipation of future, appropriate numerical simulations to accurately
compute this distribution, we here propose to discretize the DW network on a
lattice and to use principles and results from percolation theory [43, 44,
45].

### IV.1. Discretization

A DW network can be viewed as a collection of domains of constant field
configuration whose typical size is set by the correlation length $L$, which
in the scaling regime becomes of the order of the horizon size $L\sim t$. The
correlation length is usually estimated from the network energy density [13,
26, 46, 47]:

$\frac{L}{t}\simeq\frac{\sigma}{\rho_{\rm DW}t}=\mathcal{A}^{-1},$ (23)

where $\mathcal{A}$ is the area parameter introduced in Eq. (2). Using results
from field-theory lattice simulations, for $\mathcal{Z}_{2}$-symmetric models
one finds $L/t\simeq 1.25^{+0.18}_{-0.14}$ [26] or $L/t\simeq 1.14\pm 0.04$
[47], while $\mathcal{Z}_{N_{\rm DW}}$-symmetric models with $N_{\rm DW}>2$
are associated to a shorter correlation length, $L/t\simeq(2.70\pm
0.30)/N_{\rm DW}$ [26]. The values of the field $\phi(r,t)$ giving rise to the
DW networks are correlated within region of size $L$ and becomes quickly
uncorrelated at larger distances. This motivates the modeling of the DW
network as a 3D periodic lattice containing $N^{3}$ sites with lattice spacing
set equal to the correlation length $L$ [48, 45, 49, 50].

Assuming a $\mathcal{Z}_{2}$-symmetric model, every individual lattice sites
can be in one of two states, either occupied or empty, corresponding to the
field sitting in the false or true vacuum, respectively. Each site is occupied
or empty entirely randomly, independently of the state of its neighbors, with
a probability $p$ or $1-p$, respectively. In presence of an initial bias
$V_{\rm bias}$, the probability of false vacuum occupation reads [18]

$\frac{p}{1-p}=\exp\left(-\Delta F/T\right)\simeq\exp\left(-V_{\rm
bias}/V_{0}\right),$ (24)

where $T$ is the temperature, $\Delta F$ and $V_{0}$ are the free energy
difference and energy density barrier between the two minima. We assume that
the potential bias is negligible at the time of DW formation $V_{\rm bias}\ll
V_{0}$ such that we have $p=0.5$. The occupied sites are either isolated from
one another or they form small groups of neighbors. These groups are called
clusters [43, 44, 45]. An $s$-cluster is defined as a group of $s$ occupied
lattice sites connected by nearest neighbor distances. Above some critical
probability $p_{c}$ ($=0.311$ for a cubic lattice), the system undergoes a
phase transition from having no infinite cluster when $p<p_{c}$ to having one
infinite cluster when $p>p_{c}$; in other words, the system becomes percolated
by one continuous path of occupied sites. This infinite cluster is unique.
Apart from the infinite cluster, there will also be many finite clusters
dispersed throughout the lattice. We define $n_{s}$ as the total number of
finite clusters of size $s$ divided by the total number of lattice sites
$N^{3}$. Hence, the total number of such clusters is $N_{s}=n_{s}\times
N^{3}$. The probability that any selected lattice site is part of a
$s$-cluster is given by $P_{s}=s\times n_{s}$. We refer the reader to [43, 44,
45] for reviews on percolation theory.

Figure 5:  Left: Minimum number of lattice sites, denoted as $s_{\rm ball}$,
required to completely cover a spherical region of radius $R$. The green line
illustrates the lattice result, which is given by the analytical summations
derived in App. B. To eliminate the unphysical stair-step behavior, a smooth
interpolation is applied, as depicted by the black line, with the additional
condition that it must not surpass the upper limit represented by the red
line, see Eq. (30). As the radius $R$ increases, the number of lattice sites
approaches the perfect-packing limit, indicated by the blue line. Right:
Energy fraction of the network retained in closed DWs, which are sufficiently
large to contain a false vacuum sphere with radius $R$. It is obtained by
injecting the function $s_{\rm ball}(R)$ shown in the left panel in Eq. (28).
The black line represents the most accurate result, achieving a smooth
interpolation of the lattice data depicted in green, while consistently
exceeding the lower limit, marked by the red line. In contrast, the MCMC
result, presented in App. A and shown in gray, obtained from injecting Eq.
(48) in Eq. (28), tends to overestimate the fraction $\mathcal{F}$ of DWs on
the verge to collapse into PBHs as it also includes DWs that are significantly
non-spherical or that possess internal holes. It is important to note the
anti-correlation between the legends of the two panels: the lower limit in the
left panel corresponds to the upper limit in the right panel, and vice versa
for the upper limit. Instead, the colors match across both panels. $L$ is the
network correlation length.

### IV.2. Number of late-annihilators

In principle, the number $n_{s}$ of s-cluster has been calculated with Markov
Chain Monte Carlo (MCMC) simulations [45, 51, 52, 53, 54] together with their
typical radius $R_{s}$. We address this topic in App. A. However, we conclude
that since these clusters tend to be irregular in shape and contain internal
holes, we cannot rely on existing results from MCMC simulations to determine
the fraction of DWs that collapse into PBHs. Therefore, we follow a different
route which we now outline in detail. In Sec. III, we concluded that DWs
collapsing into PBHs are the ones which are larger than a given radius $R_{\rm
min}$ in Eqs. (19) or (22). We suppose that we can relax the assumption on
their spherical shape as long as they are large enough to contain a fully-
occupied ball of radius $R_{\rm min}$. In doing so, any deviation from
sphericity would add more mass than required for the collapse to occur. We
define $s_{\rm ball}(r_{\rm min})$ as the minimal number of lattice sites of
size $L$ to fully cover a ball of radius $R_{\rm min}$, with $r_{\rm
min}\equiv R_{\rm min}/L$. The probability that all those sites are occupied
is $p^{s_{\rm ball}(r_{\rm min})}$ where $p=0.5$ is the occupation
probability. Multiplying it by the highest number of balls which can be packed
on the full lattice $c_{0}N^{3}/s_{\rm ball}$ where $c_{0}=\pi/3\sqrt{2}\simeq
0.74$, we obtain the number of s-clusters large enough to contain a ball of
radius $r_{\rm min}$ (in unit of $L$):

$n_{\rm sph}(r_{\rm min})\leavevmode\nobreak\ =\leavevmode\nobreak\
\frac{c_{0}}{s_{\rm ball}(r_{\rm min})}p^{s_{\rm ball}(r_{\rm min})},$ (25)

where we normalized by the total number of sites $N^{3}$. Those cluster carry
an energy:

$E_{\rm sph}=\epsilon\times s_{\rm ball}(r_{\rm min})\times n_{\rm sph}(r_{\rm
min}),$ (26)

where $\epsilon$ is the energy of one occupied sites. This must be compared
with the energy of the DW network:

$E_{\rm DW}=\epsilon(p+p)=\epsilon,$ (27)

where the two terms accounts for the bulk and surface energy of the DW network
which at $t_{\rm ann}$ are exactly equal. We deduce that the energy fraction
contained in closed DW large enough to contain a spherical ball of radius
$r_{\rm min}$ is given by:

$\mathcal{F}(r_{\rm min})\leavevmode\nobreak\ =\leavevmode\nobreak\
\frac{E_{\rm sph}}{E_{\rm DW}}\leavevmode\nobreak\ =\leavevmode\nobreak\
s_{\rm ball}n_{\rm sph}\leavevmode\nobreak\ =\leavevmode\nobreak\
c_{0}\,p^{s_{\rm ball}(r_{\rm min})},$ (28)

where we recall $r_{\rm min}\equiv R_{\rm min}/L$, $c_{0}\simeq 0.74$ and
$p=0.5$. In App. B, we calculate the minimal number $s_{\rm ball}^{\rm
latt}(r)$ of cubes of unit length required to fully cover a ball of radius
$r$, under the form of analytical sums. The result is plotted with a green
line in the left panel of Fig. 5 against its upper and lower limit shown with
red and blue lines:

$\left\lceil\frac{4\pi r^{3}}{3}\right\rceil\leavevmode\nobreak\
<\leavevmode\nobreak\ s_{\rm ball}^{\rm latt}(r)\leavevmode\nobreak\
\leq\leavevmode\nobreak\ \left\lceil 2r\right\rceil^{3}.$ (29)

We observe that for small radii $r\in[1,5]$, the value adheres to its upper
limit, representing a cube of length $2R$. This is due to the low number of
lattice sites $\in[10,10^{3}]$. Conversely, for larger radii $r\gtrsim 10$, it
approaches its lower limit, representing perfect packing. The stair-case
behaviour of $s_{\rm ball}^{\rm latt}(r)$ arises because of the modeling of
the field configuration on a discrete lattice. It should not reflect any
physical property of the real system that should instead show a smooth
behavior. This motivates the introduction of a smooth version $s_{\rm
ball}^{\rm smooth}(r)$ obtained under the double condition to interpolate the
analytical sum $s_{\rm ball}^{\rm latt}(r)$ in Eq. (52) as best as possible
while keeping it always smaller than its upper limit $s_{\rm ball}^{\rm
smooth}(r)\leavevmode\nobreak\ \leq\leavevmode\nobreak\ \left\lceil
2r\right\rceil^{3}$.

The resulting smooth interpolation shown with black line in the left panel of
Fig. 5 reads:

$\textrm{Log}_{10}(s_{\rm ball}^{\rm
smooth})=-0.129\textrm{Lr}^{2}+2.88\textrm{Lr}+0.949,$ (30)

where $\textrm{Lr}\equiv\textrm{Log}_{10}{(r)}$ and $r\equiv R/L$. The energy
fraction $\mathcal{F}(r_{\rm min})$ of closed DWs larger than a ball of radius
$r_{\rm min}$ is obtained from plugging $s_{\rm ball}^{\rm smooth}(r_{\rm
min})$ in Eq. (30) into Eq. (28), and is shown with a black line in the right
panel of Fig. 5. It is compared to the expression obtained from injecting the
analytical sum $s_{\rm ball}^{\rm latt}(r_{\rm min})$ in Eq. (52) of App. B as
well as its upper and lower limits in Eq. (29), resulting in lower and upper
limits on $\mathcal{F}$ shown with red and blue lines respectively.
Additionally, we depict with a gray line the energy fraction $\mathcal{F}$ of
late-annihilators which would have been obtained if results from MCMC
simulations presented in Eq. (48) of App. A were plugged in Eq. (28).

## V. PBH AND WORMHOLE ABUNDANCE

Now that we are in possession of the critical size $R_{\rm min}^{\rm PBH}$ and
$R_{\rm min}^{\rm baby}$ beyond which closed DW collapse into PBH and
wormholes, derived in Sec. III, and of the DW size distribution derived in
Sec. IV, we can now proceed to the derivation of the PBH and wormhole
abundance.

### V.1. PBHs

The PBH contribution to the DM abundance is

$f_{\rm PBH}=\frac{\rho_{\rm coll}(t_{\rm PBH})}{\rho_{\rm DM}(t_{\rm
PBH})}=\frac{\rho_{\rm DW}(t_{\rm ann})}{\rho_{\rm DM}(t_{\rm
PBH})}\frac{\rho_{\rm coll}(t_{\rm PBH})}{\rho_{\rm coll}(t_{\rm
ann})}\mathcal{F}_{\rm coll},$ (31)

where $\rho_{\rm DM}(t)$ is the DM energy density, $\rho_{\rm coll}(t)$ is the
energy density stored in collapsing DWs and $\mathcal{F}_{\rm coll}$ is the
collapsing fraction at the onset of annihilation,

$\mathcal{F}_{\rm coll}\equiv\frac{\rho_{\rm coll}(t_{\rm ann})}{\rho_{\rm
DW}(t_{\rm ann})}.$ (32)

The first factor in right side of Eq. (31) can be evaluated from evolving DWs,
DM and radiation energy densities like $t^{-1}$, $a^{-3}$ and
$g_{\star}(T)T^{4}$ respectively and matching them at $t_{\rm dom}^{{\rm
unbias}}$ in Eq. (6) and matter-radiation equality when $T_{\rm eq}\simeq
0.80\leavevmode\nobreak\ \rm eV$. The second factor in Eq. (31) can be
evaluated from using that DWs collapsing into PBHs have super-horizon size and
therefore their energy evolves approximately as $a^{-1}$. One obtains:

$f_{\rm PBH}\leavevmode\nobreak\ =\leavevmode\nobreak\
\mathcal{G}\times\mathcal{R}\times\frac{T_{\rm dom}}{T_{\rm
eq}}\times\mathcal{F}_{\rm coll},$ (33)

with

$\mathcal{G}\equiv\left(\frac{g_{s\star}^{\rm eq}}{g_{\star}^{\rm
eq}}\right)\left(\frac{g_{s\star}^{\rm PBH}}{g_{s\star}^{\rm
ann}}\right)^{\\!\\!1/3}\left(\frac{(g_{\star}^{\rm ann}g^{\rm
dom}_{\star})^{\rm 1/4}}{g_{s\star}^{\rm PBH}/(g^{\rm
PBH}_{\star})^{\\!1/2}}\right)\sim 1,$ (34)

and

$\mathcal{R}\simeq\leavevmode\nobreak\ \left(\frac{t_{\rm PBH}}{t_{\rm
dom}^{{\rm unbias}}}\right)\left(\frac{t_{\rm dom}}{t_{\rm
ann}}\right)^{\\!1/2}\sim 1.$ (35)

The different temperatures can be calculated from the characteristic times in
Eqs. (3,4,6,13) using $T\simeq 1.23(M_{\rm Pl}/t)^{1/2}/g_{\star}(T)^{1/4}$,
with $g_{\star}(T)$ and $g_{*s}(T)$ the number of relativistic degrees of
freedom appearing in energy and entropy density respectively, written as
$g_{*}^{X}$ and $g_{s,*}^{X}$ with $X$ the associated epoch in Eq. (34). In
Eqs. (34) and (35), behind the sign $\sim$ we anticipated that PBH formation
reaches maximum efficiency when DWs have the longest lifespan. This happens
when they annihilate just before dominating the universe $t_{\rm ann}\sim
t_{\rm dom}\sim t_{\rm dom}^{{\rm unbias}}$, which together with Eq. (18),
implies $\mathcal{G}\sim 1$ and $\mathcal{R}\sim 1$. The collapsing fraction
$\mathcal{F}_{\rm coll}$ is given by the fraction of closed DWs in Eq. (28)
with a radius larger than the critical threshold:

$\mathcal{F}_{\rm coll}=\mathcal{F}(R_{\rm ann}^{\rm
PBH}/L)\leavevmode\nobreak\ \simeq\leavevmode\nobreak\ p^{s_{\rm ball}(R_{\rm
ann}^{\rm PBH}/L)},$ (36)

where $p=0.5$. The function $s_{\rm ball}(r=R/L)$, given in Eq. (30), is the
number of correlated regions within a ball of radius $R$. The network
correlation length $L$ must be evaluated at $t_{\rm ann}$ just before the
annihilation stage starts, thus in the scaling regime. Possible values for $L$
are discussed below Eq. (23). The critical radius $R_{\rm ann}^{\rm PBH}$ at
$t_{\rm ann}$ beyond which DWs are expected to collapse into PBHs is given by
the fitting function in Eq. (19). For rapid use, we propose the following
fitting function for Eq. (36):

$\mathcal{F}_{\rm
coll}\simeq\exp\left[-\frac{a}{\ell^{b}}\left(\frac{1}{\alpha_{\rm
ann}}\right)^{\dfrac{c}{\ell^{d}}}\right],$ (37)

with $a\simeq 0.979$, $b\simeq 2.30$, $c\simeq 1.50$, and $d=0.195$,
successfully tested against the numerical calculation for $\ell\equiv
L/t\in[0.2,2]$. The PBH abundance in shown with brown lines in Fig. 6 and Fig.
7.

Figure 6:  Exclusion constraints on DW network occupying an energy fraction
$\alpha_{\rm ann}$ defined in Eq. (11) at the onset of the annihilation phase
driven by a bias energy density difference $V_{\rm bias}$ between distinct
vacua. PBHs are produced around the temperature $V_{\rm bias}^{1/4}$ shown on
the bottom x-axis with the PBH mass shown on the top x-axis. The colorful
regions are excluded due to the presence of PBHs. The two yellow areas on the
left side indicate regions excluded by the Cosmic Microwave Background (CMB)
either due to $\mu$-distortion caused by inhomogeneities accompanying the
production of very large PBH shown in dashed [55, 56, 57], or due to the
accretion of large PBHs shown in solid [58, 59, 60]. The region in cyan is
ruled out due to limits on the black hole merging rate from LIGO-Virgo-Kagra
(LVK) [61]. The purple region is excluded by the constraints from MACHO [62],
Eros [63], OGLE [64], and HSC [65] microlensing experiments. The red and
purple dotted horizontal ellipses show the best-fit regions that can explain
the anomalies observed in OGLE and HSC microlensing data [65, 64, 66]. The
brown area shows where DM would over-close the universe. Above the purple
dashed line, more than one baby-universe $f_{\rm baby}\gtrsim 1$ is produced
in our past lightcone. We only show it in the region not excluded by PBH
overproduction. At smaller masses, PBH evaporate within a universe lifetime
[67, 68]. Since we do not observe the presence of Hawking radiation, neither
in terms of cosmic-ray fluxes [69, 70], nor in terms of modification of the
ionisation fraction in CMB [71, 72, 73], nor in terms of modificaton of the
abundance of light elements produced during Big-Bang Nucleosynthesis (BBN)
[69], we can exclude the red, yellow and green regions on the right side
respectively. The gray region labelled “DW domination” defined by $t_{\rm
ann}\gtrsim t_{\rm dom}$, see Eq. (11), indicates where the bias vacuum energy
becomes larger than the radiation density of the universe. Such region is
expected to lead to a PBH-dominated universe as mentioned in the introduction.
Additionally, the blue-gray regions indicate where stochastic GW background
(SGWB) produced from DW annihilation are excluded by NANOGrav 15-year (NG15)
data [74] and $O_{3}$-run of LIGO-Virgo-Kagra (LVK) [75]. We use that no SGWB
with a signal-to-noise ratio larger than $\rm SNR\leq 5$ has been detected in
NG15 and $\rm SNR\leq 2$ for LVK. The blue-filled black dotted vertical
ellipse shows the $90\%$ favored region which can explain the SGWB recently
detected in NG15, using the Bayesian analysis of [76]. The regions labeled
“BBN” in gray, are excluded by the constraints on number of effective degrees
of freedom $N_{\rm eff}\lesssim 0.4$ [77, 78]. Two scenarios must be
distinguished according to whether the DW network annihilates into degrees of
freedom from the Standard Model (SM) or from a dark sector, see App. C for the
details. Lastly, the DW network correlation length has been set to $L=0.8t$,
which is a rather conservative assumption for $\mathcal{Z}_{2}$-symmetric
network, see Eq. (23).

The minimal PBH mass is given by the mass inside the Schwarszchild radius at
horizon crossing $R_{\rm sch}(t_{\rm PBH})=t_{\rm PBH}$:

$M_{\rm PBH}\simeq\frac{t_{\rm PBH}}{2G}\simeq
M_{\mathTerra}\left(\frac{217\leavevmode\nobreak\ \rm GeV}{V_{\rm
bias}^{1/4}}\right)^{2},$ (38)

where $t_{\rm PBH}$ is given by Eq. (18) and $M_{\mathTerra}\simeq 3.00\times
10^{-6}\leavevmode\nobreak\ M_{\mathSun}$ where $M_{\mathTerra}$ and
$M_{\mathSun}$ are Earth and Sun masses. The PBH mass can also be calculated
as $M_{\rm PBH}\simeq 4\times 4\pi t_{\rm PBH}^{3}V_{\rm bias}/3$, which is
the bulk mass inside the horizon multiplied by a factor $4$ accounting for the
surface and gravitational binding energies, see the right panel of Fig. 2. In
principle, PBH produced above the threshold $R\gg R_{\rm ann}^{\rm PBH}$ have
larger mass than the minimal value in Eq. (38). The production of heavier PBH
is exponentially suppressed by the late-annihilator fraction in Eq. (28), thus
we expect a nearly-monochromatic mass distribution for PBHs produced by DW
networks. We leave its precise calculation for future works. In Fig. 6, we
recast the usual PBHs astrophysical constraints in the parameter space of DW
networks.

### V.2. Wormholes

The critical DW radius $R_{\rm min}^{\rm baby}$ in Eq. (22) beyond which the
collapse into PBHs is also associated to the creation of a wormhole connected
to an inflating baby-universe is about $30\%$ times larger than the critical
radius $R_{\rm min}^{\rm PBH}$ in Eq. (19) for forming simple PBHs. Due to the
exponential suppression in Eq. (28), this implies that the wormhole production
rate is significantly lower than the PBH production rate. In spite of this
suppression, is it possible that at least one baby universe has been created
within our past light cone? The number $f_{\rm baby}$ of baby universes formed
in our past lightcone is simply:

$f_{\rm baby}=\mathcal{N}_{\rm patches}(T_{\rm dom})\times\mathcal{F}(R_{\rm
ann}^{\rm baby}),$ (39)

where the DW fraction $\mathcal{F}(R_{\rm min})$ is defined in Eqs. (28) and
(30), and the baby threshold $R_{\rm ann}^{\rm baby}$ is given in Eq. (22).
The function $\mathcal{N}_{\rm patches}(T_{\rm dom})$ is the number of Hubble
patches, at the time of collapse in Eqs. (21) and (4) when the temperature is
approximately $T_{\rm dom}$, in our past light-cone. It reads:

$\displaystyle\mathcal{N}_{\rm patches}$ $\displaystyle=\left(\frac{a_{\rm
dom}H_{\rm dom}}{a_{0}H_{0}}\right)^{3}$ $\displaystyle\simeq\frac{3.3\times
10^{38}}{h^{3}}\left(\frac{g_{\star}(T_{\rm
dom})}{100}\right)^{1/2}\left(\frac{T_{\rm dom}}{100\leavevmode\nobreak\ \rm
GeV}\right)^{3},$ (40)

where $h\equiv H_{0}/100\leavevmode\nobreak\ \rm km/s/Mpc$ with $H_{0}$ the
Hubble constant today. We approximated $g_{*}(T)\simeq g_{*,s}(T)$ and used
$g_{*s}(T_{0})\simeq 3.38$. The DW network parameter space producing at least
one wormhole connected to baby universe is indicated with the dashed purple
line in Fig. 6 and Fig. 7.

### V.3. Further considerations

Finally, we discuss additional effects which could impact the PBHs abundance
predicted in this work, as well as aspects reserved for future investigation.

Collapse of sub-horizon PBHs. In this study, we only account for PBHs formed
from DWs entering their Schwarzchild radius at the moment of their horizon
entry or before. We do not account for the possibility for DWs to enter their
Schwarschild radius much after their horizon entry, at a radius $R_{\rm
sch}/t\sim(t/2t_{\rm dom})^{2}$, possibly much smaller than $t$, where we
assumed $R_{\rm sch}=2GM$ with $M\sim 4\pi t^{3}/3$. PBH formation well within
the horizon can only occur if they are spherical enough to shrink by a factor
$(2t_{\rm dom}/t)^{2}\gg 1$ to become fully contained within their
Schwarschild sphere. The PBH abundance from sub-horizon collapse is highly
sensitive to the distribution of DW shapes. For this reason, in this analysis
we made the conservative choice to not account for this additional
contribution, which we leave for future research.

Deviation from spherical symmetry. In Sec. III, we calculate the critical
radius $R_{\rm ann}^{\rm PBH}$ beyond which DWs collapse into PBH, see e.g.
Eq. (14). The fact that DWs are in general not spherical raises the question
whether this treatment is appropriate to describe the real physical system.
However, we apply this criteria on closed DW which are large enough to contain
a ball of radius $R_{\rm min}$, their abundance $\mathcal{F}(R_{\rm min})$
being derived in Sec. IV. If $R_{\rm min}=R_{\rm ann}^{\rm PBH}$, those DWs
will all unavoidably collapse into PBHs whatever their shape. This is because
non-sphericity will participate to a mass larger than $M(R_{\rm min})$ in Eq.
(14) and therefore make it easier to collapse into PBHs. Instead, our analysis
is in fact conservative with respect to non-sphericity since it misses DW
configurations which are not accounted by $\mathcal{F}(R_{\rm min})$ due to
their non-spherical shapes in spite of being successful at forming PBHs.

Correlation length. The PBH abundance shows an exponential sensitivity to the
correlation length $L$ of the DW network, as indicated in Eq. (37). In this
study, we suggest estimating $L$ through the DW network energy density, with
$L\approx\sigma/\rho_{\rm DW}$, as presented in Eq. (23). This approach
benefits from previous lattice calculations, where the relationship
$L/t\approx(2.70\pm 0.30)/N_{\rm DW}$ was established for the
$\mathcal{Z}_{N_{\rm DW}}$-symmetric model [26]. This result implies that the
PBH abundance dramatically decreases as the number $N_{\rm DW}$ of nearly-
degenerate vacua increases. A more precise determination of the correlation
length, along with its dependency on $N_{\rm DW}$, is left for future studies.
Also the precise dependence of the late-annihilator fraction as a function of
the number of vacua is left for future studies. The question at hand is
whether the fraction of late annihilators experiences a significant reduction
for $N_{\text{DW}}>2$ due to the substitution of $p=0.5$ with
$p=1/N_{\text{DW}}$ in Eq. (28).

Friction. The assumption of a significant source of friction – which is
considered negligible in this study but may become relevant in certain
scenarios [23, 24] – would slow down the motion of DW, thereby reducing the
correlation length $L$, and preventing it from keeping pace with the cosmic
horizon $L\sim t$ [79]. This could subsequently lower the resulting PBH
abundance, see e.g. the right panel of Fig. 5. Nevertheless, DWs that are
larger than the cosmic horizon should remain unaffected by friction,
suggesting that PBH formation should still occur in the presence of friction.
In any case, whether or not friction is active, as mentioned in the
introduction if the DW network were to dominate the energy density of the
universe, one should anticipate a substantial production of PBHs, including
those that host a baby universe.

Presence of an initial bias. Our analysis assumes that $V_{\rm bias}$ is
absent at the time of network formation. The presence of a large initial bias
$V_{\rm bias}$ already active at the time of network formation [50, 80, 37],
i.e. such that $t_{\rm ann}<t_{\rm form}$, would impact the occupation
probability $p$ in Eq. (24) and reduce the fraction $\mathcal{F}$ of late-
annihilators in Eq. (36) and the resulting PBH abundance. Also models of DW
with time-decreasing surface tension [81, 82, 83] would give different
results.

Figure 7:  Observability of stochastic GW background (SGWB) produced by
annihilating DW networks (blue) compared to PBH abundance (brown). Solid blue
lines indicate existing GW constraints from NANOGrav 15-year data release
(NG15) [74] and run $O_{3}$ from LIGO-Virgo-Kagra (LVK) [75], assuming the
Signal-to-Noise Ratio detection thresholds $\rm SNR=5$ and $\rm SNR=2$,
respectively. The dashed blue lines indicate the future prospects from SKA
[84], LISA [85, 86, 87] and ET/CE [88, 89, 90]. The low opacity regions show
the most optimistic prospects, assuming a low detection threshold $\rm SNR=1$
(equivalent to $1\sigma$ deviation from noise) after $T=20\leavevmode\nobreak\
\rm years$ of observation ($T=10\leavevmode\nobreak\ \rm years$ for SKA) [91]
and the perfect subtraction of all the stochastic astrophysical GW
foregrounds, see e.g. Fig. 6 in [92]. Instead the higher opacity regions show
the more conservative prospects assuming no astrophysical foreground
subtraction. The thickness of the brown bands represents the uncertainty on
the correlation length $L$ of the DW network $L/t\in[0.6,1]$, the solid line
showing the central value $L=0.8t$. Above the purple line, at least one
eternally inflating baby-universe is produced in our past light cone, assuming
$L=0.8t$. The regions labelled “BBN” and “DW domination” are the same as in
Fig. 6.

## VI. GRAVITATIONAL WAVES

During the annihilation process, DWs are driven to relativistic speed and
radiate GWs (see e.g. [16, 93, 94, 95, 96, 97, 98]). The GW power spectrum
today $\Omega_{\rm GW}^{0}$ can be related to the GW power spectrum at the
annihilation epoch $\Omega_{\rm GW}^{\rm ann}$ by:

$\displaystyle\Omega_{\rm GW}^{0}h^{2}=\mathcal{D}\,\Omega_{\rm GW}^{\rm
ann},$ (41)

where $\mathcal{D}$ is the redshift factor:

$\displaystyle\mathcal{D}$ $\displaystyle=\Omega_{\rm
rad}^{0}h^{2}\left(\frac{g_{*}(T_{\rm
ann})}{g_{\star}(T_{0})}\right)\left(\frac{g_{s,*}(T_{0})}{g_{s,*}(T_{\rm
ann})}\right)^{4/3}$ $\displaystyle\simeq 1.62\times
10^{-5}\left(\frac{g_{*}(T_{\rm
ann})}{106.75}\right)\left(\frac{106.75}{g_{s,*}(T_{\rm ann})}\right)^{4/3},$
(42)

where $\Omega_{\rm
rad}^{0}h^{2}=\pi^{2}g_{\star}(T_{0})T_{0}^{4}/30/\rho_{0}\simeq 4.18\times
10^{-5}$ with $T_{0}\simeq 2.73\leavevmode\nobreak\ \rm K$ [99],
$\rho_{0}=3M_{\rm Pl}^{2}H_{0}^{2}$ with $H_{0}\simeq 100\leavevmode\nobreak\
\rm km/s/Mpc$, $g_{s,*}(T_{0})\simeq 3.38$ and $g_{\star}(T_{0})\simeq 3.94$,
assuming $N_{\rm eff}\simeq 3.045$ [100, 101]. The GW spectrum produced by
long-lived DWs annihilating at $t_{\rm ann}$ follows from the quadrupole
formula [102, 103, 25, 104]:

$\Omega_{\rm GW}^{\rm ann}=\epsilon_{{\rm
GW}}\frac{G\mathcal{A}^{2}\sigma^{2}}{\rho_{\rm rad}(t_{\rm
ann})}S(f)=\frac{3}{32\pi}\epsilon_{{\rm GW}}\alpha_{\rm ann}^{2}S(f).$ (43)

The dimensionless quantity $\epsilon_{{\rm GW}}$ encodes the deviation from
the quadrupole formula and is fitted on lattice simulations $\epsilon_{{\rm
GW}}\simeq 0.7\pm 0.4$ [25]. We model the spectral function as

$S(f)=\frac{2}{\left(f/f_{\rm peak}\right)+\left(f_{\rm peak}/f\right)^{3}},$
(44)

where the IR slope is $\Omega_{\rm GW}\propto f^{3}$ to respect causality
[105, 106, 107, 108] and the UV slope is $\Omega_{\rm GW}\propto f^{-1}$ as
suggested by lattice simulations results [25] (though more complicated spectra
are possible, see e.g. [109]). The peak frequency is given the Hubble factor
at the time of annihilation redshifted until today [25]:

$f_{\rm peak}=\frac{a(t_{\rm ann})}{a(t_{0})}H(t_{\rm ann})\simeq
1.6\leavevmode\nobreak\ {\rm mHz}\left(\frac{g_{\star}(T_{\rm
ann})}{106.75}\right)^{1/2}\\\ \times\left(\frac{106.75}{g_{*s}(T_{\rm
ann})}\right)^{1/3}\left(\frac{T_{\rm ann}}{10\leavevmode\nobreak\ \rm
TeV}\right).$ (45)

In Fig. 6, we show the exclusion regions due to the limit on Stochastic GW
Background (SGWB) by Pulsar Timing Arrays NANOGrav [74], and earth-based
interometers LIGO-Virgo-Kagra (run O3) [75]. In Fig. 7, we also present the
potential for detection by future pulsar timing array SKA [84], as well as
forthcoming space-based interferometers LISA [85, 86], ET [88, 89], and CE
[90]. We have set the Signal-to-Noise Ratios at $\rm SNR=5$ and $\rm SNR=2$
for NG15 and LVK O3, respectively, and $\rm(SNR=1,\leavevmode\nobreak\
T=10\leavevmode\nobreak\ years)$ for SKA, as well as
$\rm(SNR=1,\leavevmode\nobreak\ T=20\leavevmode\nobreak\ years)$ for LISA, ET,
and CE, where $T$ is the observation time. The power-law integrated curves
(PLIC) are adopted from [91]. For CE and ET, we choose the minimum of the two
PLIC. We explore two different scenarios, depending on whether the expected
astrophysical foreground can be subtracted or not. We take into account the
galactic white dwarf binaries, referencing either the model from [110, 111] or
the one from [86], as well as extragalactic supermassive black hole binaries
from [112], and extragalactic compact binaries (comprising of neutron stars
and black holes) fitted on LIGO O3 data [75]. We refer to Fig. 6 in [92] for a
visualisation of the different GW astrophysical foregrounds across
frequencies.

## VII. CONCLUSION

Domain Wall (DW) networks could have formed in the early universe after the
spontaneous breaking of a $\mathcal{Z}_{N}$ symmetry with $N>1$. In presence
of a vacuum energy difference $V_{\rm bias}$ lifting the degeneracy between
the $N$ vacua, DWs are driven toward annihilating each other. In order to be
viable, the DW network must annihilate before occupying a significant energy
fraction of the universe. The Schwarzschild radius of DWs grows with radius as
$R_{\rm sch}\propto R^{2}$ or $R^{3}$ according to whether surface or bulk
energy dominates. The average size of DWs, which is set by the correlation
length $L$ of the network, grows with cosmic time as $L\simeq t$. This seems
to imply that there must be a time $t_{\rm PBH}$ when most of the DWs enter
their Schwarschild radius and collapse into PBHs. This time is closely related
the time $t_{\rm dom}\sim t_{\rm PBH}$ when the universe becomes dominated by
the DW network energy density. This is the recipe followed by previous papers
[4, 5, 6] to calculate the PBH abundance from DW networks. However, as first
pointed in [22], the DW growth cannot be applied during the annihilation phase
which necessarily precedes PBH formation $t_{\rm ann}\lesssim t_{\rm PBH}$ to
prevent DW domination. In this work, we find that only super-horizon DWs
continue to grow like $R\propto a(t)$ or slower while sub-horizon DWs shrink.
This implies that only super-horizon DWs which are larger than a given
threshold, $R\geq R_{\rm ann}^{\rm PBH}$, can enter the horizon after $t_{\rm
PBH}$ and collapse into PBHs. We call this population late-annihilators. We
calculate the PBH formation threshold $R_{\rm ann}^{\rm PBH}$ accounting for
volume, surface, and gravitational binding energies arising from Einstein
equations. We introduce a new formalism borrowed from percolation theory in 3D
to calculate the abundance of late-annihilators which we find to be
exponentially suppressed for small DW energy fraction $\alpha_{\rm ann}$ or
small network correlation length $L$, see Eq. (37).

We are able to translate PBH exclusion constraints to the parameter space of
DW networks in Fig. 6. We find that the newly-constrained parameter space
overlaps with the regions which are either currently probed by GW
observatories or will be probed in the future, see Fig. 7. We conclude that
Bayesian analysis of DW network interpretation of stochastic GW background in
current and future observatories should account for constraints from PBH
production as initiated in our companion paper [76].

Another novelty of the present study is to have determined the conditions for
DWs to continue expanding forever as seen from an observer located inside.
Those eternally inflating baby-universes are connected to their parent
universe through a wormhole. The later pinches off in one light crossing time,
leaving a PBH behind, with no residual evidence for the baby-universe genesis.
Those eternally-inflating baby universes are in principle allowed to undergo
nucleation of bubbles of any other vacua permitted by the underlying particle
physics [19, 20, 21]. Hence, baby-universes will evolve into complex space-
time structures, comprising a multitude of eternally-inflating regions
interconnected by wormholes, giving rise to a Multiverse. We find that the
possibility to have a baby-universe formed in our past lightcone is excluded
by PBH over-production, except in the region where PBH evaporate before the
onset of BBN if they are lighter than $10^{-24}\leavevmode\nobreak\
M_{\mathSun}\sim 10^{9}\leavevmode\nobreak\ \rm g$. This occurs if the bias
energy density is larger than $V_{\rm bias}\gtrsim(10^{12}\leavevmode\nobreak\
\rm GeV)^{4}$, see Figs. 6 and 7.

Finally, the $\mu$-distortion constraints shown in dashed yellow region in
Fig. 6 assume that inhomogeneities sourced by the DW networks are not highly
non-Gaussian [57] in contrast to the claim in [4]. Those $\mu$-distortion
constraints could jeopardize the possibility for PBH to be produced from QCD
axion cosmological scenarios as initially proposed in [4]. We leave more
quantitative conclusions for future studies.

## Acknowledgments

YG thanks Aleksandr Azatov, Marek Lewicki, Ken’ichi Saikawa, Peera
Simakachorn, and Ville Vaskonen for helpful discussions, and is grateful to
the Azrieli Foundation for the award of an Azrieli Fellowship. EV acknowledges
support by the European Research Council (ERC) under the European Union’s
Horizon Europe research and innovation programme (grant agreement No.
101040019).

Figure 8:  We consider two possible way to pave a ball with cubes, according
to whether the ball center is at the center or a corner of a cube, whose 2D
projections are shown in the left and right panel respectively.. The
corresponding numbers of cubes fully covering the ball are denoted $s_{\rm
ball,1}(r)$ and $s_{\rm ball,2}(r)$.

##  Appendix A MCMC simulations

The number of $s$-cluster has been calculated with Markov Chain Monte Carlo
(MCMC) simulations in both condensed matter [51, 52, 53, 54] and cosmological
context [45]. One finds [45]

$n_{s}=\frac{0.0501}{s^{\tau}}\exp\left\\{-0.6299\mathcal{P}s^{\sigma}(\mathcal{P}s^{\sigma}+1.6679)\right\\},$
(46)

where $\tau\simeq 2.17$, $\sigma\simeq 0.48$,
$\mathcal{P}\equiv(p-p_{c})/p_{c}$. For $p>p_{c}$, the average cluster radius
$R_{s}$ is given by [45]

$R_{s}/L=0.702|p-p_{c}|^{-0.124}s^{1/3}.$ (47)

For an occupation probability $p=0.5$ and the percolation threshold of a cubic
lattice $p_{c}=0.311$, Eqs. (46) and (47) simplify to

$n_{s}\simeq 0.0501s^{-2.17}\exp{-0.2326s^{0.96}-0.6385s^{0.48}},$ (48)

and $R_{s}=0.863s^{1/3}$. Since $p>p_{c}$, most of the lattice points belong
to the infinite cluster. We note $P_{\infty}$ the probability that an occupied
site belongs to this percolating cluster. A given lattice point can either be
empty with a probability $1-p$, occupied in the percolating cluster with a
probability $pP_{\infty}$, or occupied in a finite cluster with a probability
$\sum_{s}sn_{s}$. Summing all the possibilities leads to [43]

$1-p+pP_{\infty}+\sum_{s}sn_{s}=1.$ (49)

Plugging Eq. (48) into Eq. (49) for $p=0.5$, we calculate $P_{\infty}\simeq
0.92$. For the cosmological system ($N\to+\infty$), it means that $92\%$ of
the false vacuum region (with vacuum energy $V_{\rm bias}$) is contained
inside an infinite DW spanning the whole universe and only $8\%$ of the false
vacuum region is contained inside finite, closed, DWs. The number of these
closed DWs is represented by $n_{s}$ in Eq. (48). However, the formula given
in Eq. (48) is not suitable for estimating the number of DWs collapsing into
PBHs due to two considerations.

Firstly, the analysis conducted in Sec. III presumes that the DWs are
spherical. The scaling relationship between cluster volume and radius
$s\propto R_{s}^{3}$ in Eq. (47) suggests that clusters exhibit a droplet-like
shape, as elaborated in [43]. Nonetheless, when observing the simulated shapes
of DWs on a lattice, e.g. as in Fig. 1 of [13], it becomes evident that
typical closed DWs can deviate significantly from a spherical shape at larger
sizes. This deviation can be characterized by their inner radius being
potentially much smaller than their average radius $R_{s}$. Secondly, at
larger sizes, clusters tend to encompass a considerable number of internal
holes, similar to Swiss cheese [43]. In light of these factors, we deduce that
the number distribution $n_{s}$, derived from MCMC simulations in Eq. (48),
which accounts for clusters that are significantly irregular in shape and
contain internal holes, would overestimate the number of DWs on the verge of
collapse. Therefore, in this work, we determine the fraction of DWs collapsing
into PBHs with our own methods as indicated in Sec. IV.2. The outcomes
displayed in the right panel of Fig. 5 indeed indicate that results from MCMC
simulations found in the literature would significantly overestimate the
fraction of DWs undergoing collapse.

##  Appendix B Number of cubes to cover a ball

In this appendix, we calculate the minimal number of cubes of length $L$
required to fully cover a ball of radius $R$, or equivalently the number of
cubes of unit length covering a ball of radius $r=R/L$. We propose the two
configurations shown in Fig. 8 where the center of the ball is either a face
or a center of a cube. We start with the later. We can divide the ball into
$8$ octants, which we subsequently divide into $\lceil r\rceil$ layers denoted
by $j\in\llbracket 0\mkern 1.5mu..\mkern 1.5mu\lfloor r\rfloor\rrbracket$, one
of them being shown in the left panel of Fig. 8. The radius of the circle
coinciding with the boundary of the ball at a layer $j$ is
$r_{j}=\sqrt{r^{2}-j^{2}}$. Denoting the column number by $i\in\llbracket
0\mkern 1.5mu..\mkern 1.5mu\lfloor r\rfloor\rrbracket$, the number of cubes
filling a row, shown in red color, is
$\left\lceil\sqrt{r_{j}^{2}-i^{2}}\right\rceil$. We deduce the total number of
cubes filling the ball for the configuration shown in the left panel of Fig.
8:

$s_{\rm ball,1}(r)=8\sum_{j=0}^{\lfloor
r\rfloor}\sum_{i=0}^{\lfloor\sqrt{r^{2}-j^{2}}\rfloor}\left\lceil\sqrt{r^{2}-i^{2}-j^{2}}\right\rceil.$
(50)

We now pass to the next configuration where the ball center coincides with the
center of a cube. We start by isolating the 3 orthogonal layers along the
planes $(XY)$, $(XZ)$ and $(YZ)$ passing through the center of the ball, one
of them being shown in the right panel of Fig. 8. The number of cubes in the
central tri-axis blue cross is $6\lfloor r+0.5\rfloor+1$. The remaining number
of cubes in those three orthogonal layer outside the 3-axis cross is
$4\sum_{i=0}^{\lfloor
r\rfloor}\left\lceil\sqrt{r^{2}-(i+0.5)^{2}}+0.5\right\rceil$. Finally,
accounting for the remaining cubes in the 8 octants outside the three internal
layers give:

$s_{\rm ball,2}(r)=6\lfloor r+0.5\rfloor+1+4\sum_{i=0}^{\lfloor
r\rfloor}\left\lceil\sqrt{r^{2}-(i+0.5)^{2}}+0.5\right\rceil+8\sum_{j=0}^{\lfloor
r\rfloor}\sum_{i=0}^{\lfloor\sqrt{r^{2}-j^{2}}\rfloor}\left\lceil\sqrt{r^{2}-(i+0.5)^{2}-(j+0.5)^{2}}\right\rceil.$
(51)

The minimal number of cubes to fully cover a ball obtained from the analytical
sums of this appendix is

$s_{\rm ball}^{\rm latt}(r)=\textrm{Min}\left[s_{\rm ball,1}(r),s_{\rm
ball,2}(r)\right].$ (52)

The function $s_{\rm ball}^{\rm latt}(r)$ with $r=R/L$ is shown with a green
line in Fig. 5. We find that at small radius $\gtrsim 1$, it is close to its
upper limit shown in red for which the ball is replaced by its outer-cube of
length $2r$

$s_{\rm ball}^{\rm latt}(r)\gtrsim\left\lceil 2r\right\rceil^{3}.$ (53)

At large radius $r\gg 1$, the number of cubes asymptotes the perfect-packing
(lower) limit shown in blue for which the tiny cubes smoothly fit inside the
ball:

$s_{\rm ball}^{\rm latt}(r)\lesssim\left\lceil\frac{4\pi
r^{3}}{3}\right\rceil.$ (54)

##  Appendix C BBN bound

Domain walls (DWs) form a component of the total energy density of the
universe. As such, they contribute to increase the expansion rate of the
universe which makes neutron freeze-out earlier, increase the $n/p$ ratio
which in turn increases the Helium abundance [113]. The presence of DWs can be
described in terms of an extra number of neutrino species

$N_{\rm eff}=\frac{8}{7}\left(\frac{\rho_{\rm
DW}}{\rho_{\gamma}}\right)\left(\frac{11}{4}\right)^{4/3},$ (55)

where $\rho_{\gamma}$ is the photon number density. We introduce the DW energy
fraction in unit of radiation energy density at annihilation temperature

$\alpha_{\rm DW}(T)=\frac{\rho_{\rm
DW}(T)}{\frac{\pi^{2}}{30}g_{\star}(T)T^{4}}.$ (56)

where $T$ is the SM photon temperature. From Eq. (55) and Eq. (56), the
maximal DW contribution to $N_{\rm eff}$ occurs at annihilation temperature

$\Delta N_{\rm
eff}(T)=\frac{8}{7}\left(\frac{{g_{*}(T)}}{2}\right)\left(\frac{11}{4}\right)^{4/3}\alpha_{\rm
DW}(T),$ (57)

To apply the BBN bound $\Delta N_{\rm eff}\lesssim 0.3$ [77, 78], the
effective number of extra relativistic degrees of freedom must be evaluated
below the neutrino decoupling temperature where $g_{\star}(T<T_{\rm
dec})\equiv 2+(7/8)\cdot 6\cdot(4/11)^{4/3}\simeq 3.36$. Hence, we obtain

$\Delta N_{\rm eff}=7.4\leavevmode\nobreak\ \alpha_{\rm
GW}\leavevmode\nobreak\ \lesssim\leavevmode\nobreak\ 0.4,$ (58)

which is slightly different from [114]. As discussed in the main text, we must
distinguish the scenario in which DWs reheat to dark radiation, in which case
Eq. (58) is the BBN constraints, from the scenario in which DWs reheat to SM,
in which case Eq. (58) applies only if DWs annihilate below the neutrino
decoupling temperature $T_{\rm ann}\lesssim 1\leavevmode\nobreak\ \rm MeV$.

## References

  * Abbott  _et al._ [2019] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Gwtc-1: a Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by Ligo and Virgo during the First and Second Observing Runs, Phys. Rev. X 9, 031040 (2019), arXiv:1811.12907 [astro-ph.HE] .
  * Carr and Hawking [1974] B. J. Carr and S. W. Hawking, Black holes in the early Universe, Mon. Not. Roy. Astron. Soc. 168, 399 (1974).
  * Vachaspati [2017] T. Vachaspati, Lunar Mass Black Holes from QCD Axion Cosmology, (2017), arXiv:1706.03868 [hep-th] .
  * Ferrer  _et al._ [2019] F. Ferrer, E. Masso, G. Panico, O. Pujolas, and F. Rompineve, Primordial Black Holes from the QCD axion, Phys. Rev. Lett. 122, 101301 (2019), arXiv:1807.01707 [hep-ph] .
  * Gelmini  _et al._ [2023a] G. B. Gelmini, A. Simpson, and E. Vitagliano, Catastrogenesis: DM, GWs, and PBHs from ALP string-wall networks, JCAP 02, 031, arXiv:2207.07126 [hep-ph] .
  * Gelmini  _et al._ [2023b] G. B. Gelmini, J. Hyman, A. Simpson, and E. Vitagliano, Primordial black hole dark matter from catastrogenesis with unstable pseudo-Goldstone bosons, (2023b), arXiv:2303.14107 [hep-ph] .
  * Deng  _et al._ [2017] H. Deng, J. Garriga, and A. Vilenkin, Primordial black hole and wormhole formation by domain walls, JCAP 04, 050, arXiv:1612.03753 [gr-qc] .
  * Deng and Vilenkin [2017] H. Deng and A. Vilenkin, Primordial black hole formation by vacuum bubbles, JCAP 12, 044, arXiv:1710.02865 [gr-qc] .
  * Deng [2020] H. Deng, Primordial black hole formation by vacuum bubbles. Part II, JCAP 09, 023, arXiv:2006.11907 [astro-ph.CO] .
  * Liu  _et al._ [2020] J. Liu, Z.-K. Guo, and R.-G. Cai, Primordial Black Holes from Cosmic Domain Walls, Phys. Rev. D 101, 023513 (2020), arXiv:1908.02662 [astro-ph.CO] .
  * Kusenko  _et al._ [2020] A. Kusenko, M. Sasaki, S. Sugiyama, M. Takada, V. Takhistov, and E. Vitagliano, Exploring Primordial Black Holes from the Multiverse with Optical Telescopes, Phys. Rev. Lett. 125, 181304 (2020), arXiv:2001.09160 [astro-ph.CO] .
  * Ge  _et al._ [2023] S. Ge, J. Guo, and J. Liu, A New Mechanism for Primordial Black Hole Formation from QCD Axion, (2023), arXiv:2309.01739 [hep-ph] .
  * Press  _et al._ [1989] W. H. Press, B. S. Ryden, and D. N. Spergel, Dynamical Evolution of Domain Walls in an Expanding Universe, Astrophys. J. 347, 590 (1989).
  * Vilenkin and Shellard [2000] A. Vilenkin and E. P. S. Shellard,  _Cosmic Strings and Other Topological Defects_ (Cambridge University Press, 2000).
  * Kibble [1976] T. W. B. Kibble, Topology of Cosmic Domains and Strings, J. Phys. A 9, 1387 (1976).
  * Vilenkin [1981] A. Vilenkin, Gravitational Field of Vacuum Domain Walls and Strings, Phys. Rev. D 23, 852 (1981).
  * Sikivie [1982] P. Sikivie, Of Axions, Domain Walls and the Early Universe, Phys. Rev. Lett. 48, 1156 (1982).
  * Gelmini  _et al._ [1989] G. B. Gelmini, M. Gleiser, and E. W. Kolb, Cosmology of Biased Discrete Symmetry Breaking, Phys. Rev. D 39, 1558 (1989).
  * Garriga  _et al._ [2016] J. Garriga, A. Vilenkin, and J. Zhang, Black holes and the multiverse, JCAP 02, 064, arXiv:1512.01819 [hep-th] .
  * Garriga and Vilenkin [1998] J. Garriga and A. Vilenkin, Recycling universe, Phys. Rev. D 57, 2230 (1998), arXiv:astro-ph/9707292 .
  * Linde [2017] A. Linde, A brief history of the multiverse, Rept. Prog. Phys. 80, 022001 (2017), arXiv:1512.01203 [hep-th] .
  * Gelmini and Hyman [2023] G. B. Gelmini and J. Hyman, Catastrogenesis with unstable ALPs as the origin of the NANOGrav 15 yr gravitational wave signal, (2023), arXiv:2307.07665 [hep-ph] .
  * Blasi  _et al._ [2023a] S. Blasi, A. Mariotti, A. Rase, A. Sevrin, and K. Turbang, Friction on ALP domain walls and gravitational waves, JCAP 04, 008, arXiv:2210.14246 [hep-ph] .
  * Blasi  _et al._ [2023b] S. Blasi, A. Mariotti, A. Rase, and A. Sevrin, Axionic domain walls at Pulsar Timing Arrays: QCD bias and particle friction, (2023b), arXiv:2306.17830 [hep-ph] .
  * Hiramatsu  _et al._ [2014] T. Hiramatsu, M. Kawasaki, and K. Saikawa, On the estimation of gravitational wave spectrum from cosmic domain walls, JCAP 02, 031, arXiv:1309.5001 [astro-ph.CO] .
  * Kawasaki  _et al._ [2015] M. Kawasaki, K. Saikawa, and T. Sekiguchi, Axion dark matter from topological defects, Phys. Rev. D 91, 065014 (2015), arXiv:1412.0789 [hep-ph] .
  * Berezin  _et al._ [1983] V. A. Berezin, V. A. Kuzmin, and I. I. Tkachev, THIN WALL VACUUM DOMAINS EVOLUTION, Phys. Lett. B 120, 91 (1983).
  * Berezin  _et al._ [1987] V. A. Berezin, V. A. Kuzmin, and I. I. Tkachev, Dynamics of Bubbles in General Relativity, Phys. Rev. D 36, 2919 (1987).
  * Blau  _et al._ [1987] S. K. Blau, E. I. Guendelman, and A. H. Guth, The Dynamics of False Vacuum Bubbles, Phys. Rev. D 35, 1747 (1987).
  * Maeda [1986] K.-i. Maeda, Bubble dynamics in the expanding universe, Gen. Rel. Grav. 18, 931 (1986).
  * Tanahashi and Yoo [2015] N. Tanahashi and C.-M. Yoo, Spherical Domain Wall Collapse in a Dust Universe, Class. Quant. Grav. 32, 155003 (2015), arXiv:1411.7479 [gr-qc] .
  * Israel [1966] W. Israel, Singular hypersurfaces and thin shells in general relativity, Nuovo Cimento B Serie 44, 1 (1966).
  * Israel [1966] W. Israel, Singular hypersurfaces and thin shells in general relativity, Nuovo Cim. B 44S10, 1 (1966), [Erratum: Nuovo Cim.B 48, 463 (1967)].
  * Misner and Sharp [1964] C. W. Misner and D. H. Sharp, Relativistic equations for adiabatic, spherically symmetric gravitational collapse, Phys. Rev. 136, B571 (1964).
  * Hayward [1996] S. A. Hayward, Gravitational energy in spherical symmetry, Phys. Rev. D 53, 1938 (1996), arXiv:gr-qc/9408002 .
  * Ipser and Sikivie [1984] J. Ipser and P. Sikivie, The Gravitationally Repulsive Domain Wall, Phys. Rev. D 30, 712 (1984).
  * Pujolas and Zahariade [2023] O. Pujolas and G. Zahariade, Domain wall annihilation: A QFT perspective, Phys. Rev. D 107, 123527 (2023), arXiv:2212.11204 [hep-th] .
  * Gouttenoire and Volansky [2023] Y. Gouttenoire and T. Volansky, Primordial Black Holes from Supercooled Phase Transitions, (2023), arXiv:2305.04942 [hep-ph] .
  * Gouttenoire [2023] Y. Gouttenoire, First-Order Phase Transition Interpretation of Pulsar Timing Array Signal Is Consistent with Solar-Mass Black Holes, Phys. Rev. Lett. 131, 171404 (2023), arXiv:2307.04239 [hep-ph] .
  * Sato  _et al._ [1982] K. Sato, H. Kodama, M. Sasaki, and K.-i. Maeda, Multiproduction of Universes by First Order Phase Transition of a Vacuum, Phys. Lett. B 108, 103 (1982).
  * Giddings and Strominger [1988] S. B. Giddings and A. Strominger, Axion Induced Topology Change in Quantum Gravity and String Theory, Nucl. Phys. B 306, 890 (1988).
  * Hawking [1990] S. W. Hawking, Baby universes, in  _1st Alexander Friedmann International Seminar on Gravitation and Cosmology_ (1990) pp. 81–92.
  * Stauffer [1979] D. Stauffer, Scaling theory of percolation clusters, Phys. Rept. 54, 1 (1979).
  * Essam [1980] J. W. Essam, Percolation theory, Reports on Progress in Physics 43, 833 (1980).
  * Lalak  _et al._ [1995a] Z. Lalak, B. A. Ovrut, and S. Thomas, Large scale structure as a critical phenomenon, Phys. Rev. D 51, 5456 (1995a).
  * Oliveira  _et al._ [2005] J. C. R. E. Oliveira, C. J. A. P. Martins, and P. P. Avelino, The Cosmological evolution of domain wall networks, Phys. Rev. D 71, 083509 (2005), arXiv:hep-ph/0410356 .
  * Martins  _et al._ [2016] C. J. A. P. Martins, I. Y. Rybak, A. Avgoustidis, and E. P. S. Shellard, Extending the Velocity-Dependent One-Scale Model for Domain Walls, Phys. Rev. D 93, 043534 (2016), arXiv:1602.01322 [hep-ph] .
  * Lalak and Ovrut [1993] Z. Lalak and B. A. Ovrut, Domain walls, percolation theory and Abell clusters, Phys. Rev. Lett. 71, 951 (1993).
  * Lalak  _et al._ [1995b] Z. Lalak, S. Lola, B. A. Ovrut, and G. G. Ross, Large scale structure from biased nonequilibrium phase transitions: Percolation theory picture, Nucl. Phys. B 434, 675 (1995b), arXiv:hep-ph/9404218 .
  * Coulson  _et al._ [1996] D. Coulson, Z. Lalak, and B. A. Ovrut, Biased domain walls, Phys. Rev. D 53, 4237 (1996).
  * Sykes  _et al._ [1976] M. F. Sykes, D. S. Gaunt, and M. Glen, Percolation processes in three dimensions, Journal of Physics A Mathematical General 9, 1705 (1976).
  * Flammang [1977] A. Flammang, Percolation cluster sizes and perimeters in three dimensions, Zeitschrift für Physik B Condensed Matter 28, 47 (1977).
  * Hoshen  _et al._ [1979] J. Hoshen, D. Stauffer, G. Bishop, R. Harrison, and G. Quinn, Monte carlo experiments on cluster size distribution in percolation, Journal of Physics A General Physics 12 (1979).
  * Nakanishi and Stanley [1980] H. Nakanishi and H. E. Stanley, Scaling studies of percolation phenomena in systems of dimensionality two to seven: Cluster numbers, Phys. Rev. B 22, 2466 (1980).
  * Kohri  _et al._ [2014] K. Kohri, T. Nakama, and T. Suyama, Testing scenarios of primordial black holes being the seeds of supermassive black holes by ultracompact minihalos and CMB $\mu$-distortions, Phys. Rev. D 90, 083514 (2014), arXiv:1405.5999 [astro-ph.CO] .
  * Nakama  _et al._ [2016] T. Nakama, T. Suyama, and J. Yokoyama, Supermassive black holes formed by direct collapse of inflationary perturbations, Phys. Rev. D 94, 103522 (2016), arXiv:1609.02245 [gr-qc] .
  * Nakama  _et al._ [2018] T. Nakama, B. Carr, and J. Silk, Limits on primordial black holes from $\mu$ distortions in cosmic microwave background, Phys. Rev. D 97, 043525 (2018), arXiv:1710.06945 [astro-ph.CO] .
  * Ali-Haïmoud and Kamionkowski [2017] Y. Ali-Haïmoud and M. Kamionkowski, Cosmic microwave background limits on accreting primordial black holes, Phys. Rev. D 95, 043534 (2017), arXiv:1612.05644 [astro-ph.CO] .
  * Poulin  _et al._ [2017a] V. Poulin, P. D. Serpico, F. Calore, S. Clesse, and K. Kohri, Cmb Bounds on Disk-Accreting Massive Primordial Black Holes, Phys. Rev. D 96, 083524 (2017a), arXiv:1707.04206 [astro-ph.CO] .
  * Serpico  _et al._ [2020] P. D. Serpico, V. Poulin, D. Inman, and K. Kohri, Cosmic microwave background bounds on primordial black holes including dark matter halo accretion, Phys. Rev. Res. 2, 023204 (2020), arXiv:2002.10771 [astro-ph.CO] .
  * De Luca  _et al._ [2020] V. De Luca, G. Franciolini, P. Pani, and A. Riotto, Primordial Black Holes Confront LIGO/Virgo data: Current situation, JCAP 06, 044, arXiv:2005.05641 [astro-ph.CO] .
  * Alcock  _et al._ [2000] C. Alcock _et al._ (MACHO), The MACHO project: Microlensing results from 5.7 years of LMC observations, Astrophys. J. 542, 281 (2000), arXiv:astro-ph/0001272 .
  * Tisserand  _et al._ [2007] P. Tisserand _et al._ (EROS-2), Limits on the Macho Content of the Galactic Halo from the EROS-2 Survey of the Magellanic Clouds, Astron. Astrophys. 469, 387 (2007), arXiv:astro-ph/0607207 .
  * Niikura  _et al._ [2019a] H. Niikura, M. Takada, S. Yokoyama, T. Sumi, and S. Masaki, Constraints on Earth-mass primordial black holes from OGLE 5-year microlensing events, Phys. Rev. D 99, 083503 (2019a), arXiv:1901.07120 [astro-ph.CO] .
  * Niikura  _et al._ [2019b] H. Niikura _et al._ , Microlensing constraints on primordial black holes with Subaru/HSC Andromeda observations, Nature Astron. 3, 524 (2019b), arXiv:1701.02151 [astro-ph.CO] .
  * Sugiyama  _et al._ [2023] S. Sugiyama, M. Takada, and A. Kusenko, Possible evidence of axion stars in HSC and OGLE microlensing events, Phys. Lett. B 840, 137891 (2023), arXiv:2108.03063 [hep-ph] .
  * Hawking [1974] S. W. Hawking, Black hole explosions, Nature 248, 30 (1974).
  * Hawking [1975] S. W. Hawking, Particle Creation by Black Holes, Commun. Math. Phys. 43, 199 (1975), [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * Carr  _et al._ [2010] B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, New cosmological constraints on primordial black holes, Phys. Rev. D 81, 104019 (2010), arXiv:0912.5297 [astro-ph.CO] .
  * Boudaud and Cirelli [2019] M. Boudaud and M. Cirelli, Voyager 1 $e^{\pm}$ Further Constrain Primordial Black Holes as Dark Matter, Phys. Rev. Lett. 122, 041104 (2019), arXiv:1807.03075 [astro-ph.HE] .
  * Poulin  _et al._ [2017b] V. Poulin, J. Lesgourgues, and P. D. Serpico, Cosmological Constraints on Exotic Injection of Electromagnetic Energy, JCAP 03, 043, arXiv:1610.10051 [astro-ph.CO] .
  * Stöcker  _et al._ [2018] P. Stöcker, M. Krämer, J. Lesgourgues, and V. Poulin, Exotic energy injection with ExoCLASS: Application to the Higgs portal model and evaporating black holes, JCAP 03, 018, arXiv:1801.01871 [astro-ph.CO] .
  * Poulter  _et al._ [2019] H. Poulter, Y. Ali-Haïmoud, J. Hamann, M. White, and A. G. Williams, CMB constraints on ultra-light primordial black holes with extended mass distributions, (2019), arXiv:1907.06485 [astro-ph.CO] .
  * Agazie  _et al._ [2023] G. Agazie _et al._ (NANOGrav), The NANOGrav 15 yr Data Set: Detector Characterization and Noise Budget, Astrophys. J. Lett. 951, L10 (2023), arXiv:2306.16218 [astro-ph.HE] .
  * Abbott  _et al._ [2021] R. Abbott _et al._ (KAGRA, Virgo, LIGO Scientific), Upper limits on the isotropic gravitational-wave background from Advanced LIGO and Advanced Virgo’s third observing run, Phys. Rev. D 104, 022004 (2021), arXiv:2101.12130 [gr-qc] .
  * Gouttenoire and Vitagliano [2023] Y. Gouttenoire and E. Vitagliano, Domain wall interpretation of the PTA signal confronting black hole overproduction, (2023), arXiv:2306.17841 [gr-qc] .
  * Pitrou  _et al._ [2018] C. Pitrou, A. Coc, J.-P. Uzan, and E. Vangioni, Precision big bang nucleosynthesis with improved Helium-4 predictions, Phys. Rept. 754, 1 (2018), arXiv:1801.08023 [astro-ph.CO] .
  * Dvorkin  _et al._ [2022] C. Dvorkin _et al._ , The Physics of Light Relics, in  _2022 Snowmass Summer Study_ (2022) arXiv:2203.07943 [hep-ph] .
  * Martins and Shellard [1996] C. J. A. P. Martins and E. P. S. Shellard, String evolution with friction, Phys. Rev. D 53, 575 (1996), arXiv:hep-ph/9507335 .
  * Hindmarsh [1996] M. Hindmarsh, Analytic scaling solutions for cosmic domain walls, Phys. Rev. Lett. 77, 4495 (1996), arXiv:hep-ph/9605332 .
  * Babichev  _et al._ [2022] E. Babichev, D. Gorbunov, S. Ramazanov, and A. Vikman, Gravitational shine of dark domain walls, JCAP 04 (04), 028, arXiv:2112.12608 [hep-ph] .
  * Ramazanov  _et al._ [2022] S. Ramazanov, E. Babichev, D. Gorbunov, and A. Vikman, Beyond freeze-in: Dark matter via inverse phase transition and gravitational wave signal, Phys. Rev. D 105, 063530 (2022), arXiv:2104.13722 [hep-ph] .
  * Babichev  _et al._ [2023] E. Babichev, D. Gorbunov, S. Ramazanov, R. Samanta, and A. Vikman, NANOGrav spectral index $\gamma=3$ from melting domain walls, (2023), arXiv:2307.04582 [hep-ph] .
  * Janssen  _et al._ [2015] G. Janssen _et al._ , Gravitational wave astronomy with the SKA, PoS AASKA14, 037 (2015), arXiv:1501.00127 [astro-ph.IM] .
  * Audley  _et al._ [2017] H. Audley _et al._ (LISA), Laser Interferometer Space Antenna, (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Robson  _et al._ [2019] T. Robson, N. J. Cornish, and C. Liu, The construction and use of LISA sensitivity curves, Class. Quant. Grav. 36, 105011 (2019), arXiv:1803.01944 [astro-ph.HE] .
  * Auclair  _et al._ [2023] P. Auclair _et al._ (LISA Cosmology Working Group), Cosmology with the Laser Interferometer Space Antenna, Living Rev. Rel. 26, 5 (2023), arXiv:2204.05434 [astro-ph.CO] .
  * Punturo  _et al._ [2010] M. Punturo _et al._ , The Einstein Telescope: A third-generation gravitational wave observatory, Class. Quant. Grav. 27, 194002 (2010).
  * Maggiore  _et al._ [2020] M. Maggiore _et al._ , Science Case for the Einstein Telescope, JCAP 03, 050, arXiv:1912.02622 [astro-ph.CO] .
  * Reitze  _et al._ [2019] D. Reitze _et al._ , Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO, Bull. Am. Astron. Soc. 51, 035 (2019), arXiv:1907.04833 [astro-ph.IM] .
  * Schmitz [2021] K. Schmitz, New Sensitivity Curves for Gravitational-Wave Signals from Cosmological Phase Transitions, JHEP 01, 097, arXiv:2002.04615 [hep-ph] .
  * Baldes  _et al._ [2023] I. Baldes, M. Dichtl, Y. Gouttenoire, and F. Sala, Bubbletrons, (2023), arXiv:2306.15555 [hep-ph] .
  * Preskill  _et al._ [1991] J. Preskill, S. P. Trivedi, F. Wilczek, and M. B. Wise, Cosmology and broken discrete symmetry, Nucl. Phys. B 363, 207 (1991).
  * Chang  _et al._ [1999] S. Chang, C. Hagmann, and P. Sikivie, Studies of the motion and decay of axion walls bounded by strings, Phys. Rev. D 59, 023505 (1999), arXiv:hep-ph/9807374 .
  * Gleiser and Roberts [1998] M. Gleiser and R. Roberts, Gravitational waves from collapsing vacuum domains, Phys. Rev. Lett. 81, 5497 (1998), arXiv:astro-ph/9807260 .
  * Hiramatsu  _et al._ [2013] T. Hiramatsu, M. Kawasaki, K. Saikawa, and T. Sekiguchi, Axion cosmology with long-lived domain walls, JCAP 01, 001, arXiv:1207.3166 [hep-ph] .
  * Gelmini  _et al._ [2021a] G. B. Gelmini, A. Simpson, and E. Vitagliano, Gravitational waves from axionlike particle cosmic string-wall networks, Phys. Rev. D 104, 061301 (2021a), arXiv:2103.07625 [hep-ph] .
  * Li  _et al._ [2023] Y. Li, L. Bian, R.-G. Cai, and J. Shu, Gravitational waves radiated from axion string-wall networks, (2023), arXiv:2311.02011 [astro-ph.CO] .
  * Workman  _et al._ [2022] R. L. Workman _et al._ (Particle Data Group), Review of Particle Physics, PTEP 2022, 083C01 (2022).
  * Mangano  _et al._ [2005] G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti, and P. D. Serpico, Relic neutrino decoupling including flavor oscillations, Nucl. Phys. B 729, 221 (2005), arXiv:hep-ph/0506164 .
  * de Salas and Pastor [2016] P. F. de Salas and S. Pastor, Relic neutrino decoupling with flavour oscillations revisited, JCAP 07, 051, arXiv:1606.06986 [hep-ph] .
  * Hiramatsu  _et al._ [2010] T. Hiramatsu, M. Kawasaki, and K. Saikawa, Gravitational Waves from Collapsing Domain Walls, JCAP 05, 032, arXiv:1002.1555 [astro-ph.CO] .
  * Kawasaki and Saikawa [2011] M. Kawasaki and K. Saikawa, Study of gravitational radiation from cosmic domain walls, JCAP 09, 008, arXiv:1102.5628 [astro-ph.CO] .
  * Saikawa [2017] K. Saikawa, A review of gravitational waves from cosmic domain walls, Universe 3, 40 (2017), arXiv:1703.02576 [hep-ph] .
  * Durrer and Caprini [2003] R. Durrer and C. Caprini, Primordial Magnetic Fields and Causality, JCAP 11, 010, arXiv:astro-ph/0305059 .
  * Caprini  _et al._ [2009] C. Caprini, R. Durrer, T. Konstandin, and G. Servant, General Properties of the Gravitational Wave Spectrum from Phase Transitions, Phys. Rev. D 79, 083519 (2009), arXiv:0901.1661 [astro-ph.CO] .
  * Cai  _et al._ [2020] R.-G. Cai, S. Pi, and M. Sasaki, Universal Infrared Scaling of Gravitational Wave Background Spectra, Phys. Rev. D 102, 083528 (2020), arXiv:1909.13728 [astro-ph.CO] .
  * Hook  _et al._ [2021] 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, 117, arXiv:2010.03568 [hep-ph] .
  * Gelmini  _et al._ [2021b] G. B. Gelmini, S. Pascoli, E. Vitagliano, and Y.-L. Zhou, Gravitational wave signatures from discrete flavor symmetries, JCAP 02, 032, arXiv:2009.01903 [hep-ph] .
  * Lamberts  _et al._ [2019] A. Lamberts, S. Blunt, T. B. Littenberg, S. Garrison-Kimmel, T. Kupfer, and R. E. Sanderson, Predicting the Lisa White Dwarf Binary Population in the Milky Way with Cosmological Simulations, Mon. Not. Roy. Astron. Soc. 490, 5888 (2019), arXiv:1907.00014 [astro-ph.HE] .
  * Boileau  _et al._ [2022] G. Boileau, A. C. Jenkins, M. Sakellariadou, R. Meyer, and N. Christensen, Ability of Lisa to Detect a Gravitational-Wave Background of Cosmological Origin: the Cosmic String Case, Phys. Rev. D 105, 023510 (2022), arXiv:2109.06552 [gr-qc] .
  * Rosado [2011] P. A. Rosado, Gravitational Wave Background from Binary Systems, Phys. Rev. D 84, 084004 (2011), arXiv:1106.5795 [gr-qc] .
  * Baumann [2022] D. Baumann,  _Cosmology_ (Cambridge University Press, 2022).
  * Ferreira  _et al._ [2023] R. Z. Ferreira, A. Notari, O. Pujolas, and F. Rompineve, Gravitational waves from domain walls in Pulsar Timing Array datasets, JCAP 02, 001, arXiv:2204.04228 [astro-ph.CO] .