Primordial Black Holes and Wormholes from Domain Wall Networks

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 nearlydegenerate 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 ∼ t [13], and the fraction of the total universe energy density stored in DWs increases linearly with time ρ DW /ρ tot ∝ 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 bias between the different vacua, as we review in Sec.II.After a time t ann , the vacuum energy difference V bias counterbalances the pressure due to the wall tension σ, driving DWs towards annihilation before they can dominate the universe at a time t 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 ann ≳ t 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 ann ≲ t 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].
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 con- ann at the onset of the annihilation phase, which are dubbed "late-annihilators" in the Orange region, enter inside their Schwarzschild radius R ≲ 2GM in the Blue region, leading them to collapse into PBHs.DWs with radius larger than R baby ann become larger than the cosmological horizon H −1 b of an inflating universe with same vacuum energy V bias as their interior, in the Purple region, with H 2 b = V bias /3M 2 Pl .They become eternally-inflating baby-universes detaching from the parent universe through wormholes and leaving PBHs behind.
nected to an eternally-inflating baby-universe.Previous studies [4][5][6] incorrectly assumed that the radius of DWs grows as R ≃ t.Using energy arguments, Ref. [22] commented that the relation R ≃ t is not applicable during the annihilation phase.In this study, we provide a definitive resolution to the issue, stating that R ≃ 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 equa-arXiv:2311.07670v1[hep-ph] 13 Nov 2023 tions in the thin-shell limit.The main results are sketched in Fig. 1.The radius grows as R ∝ a or slower if superhorizon R ≳ t, but quickly shrink if sub-horizon R ≲ t.We find that PBHs enter their Schwarschild radius if they have a radius larger than R ≳ R PBH ann > 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 lateannihilators, continue to grow during the annihilation phase until they finally collapse when entering inside the cosmic horizon.We also derive a critical radius R baby ann larger than R PBH ann by about 30% above which DWs collapse into wormholes connected to eternally-inflating baby-universes.We find that the PBH formation threshold R PBH ann is highly sensitive to the surface energy of DWs, which can have Lorentz factors up to γ ∼ 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 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 babyuniverse production.

II. BIASED NETWORK EVOLUTION
Initially, the work of the DW surface tension σ, with equivalent pressure P T = C d σ/L toward straightening, is dampened by the friction pressure P V ≃ β T 4 .The quantity L is the correlation length of the DW network, C d and β 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 ≃ vt, one obtains that DWs start moving with relativistic velocity v ≃ O(0.1) after the time where GeV is the reduced Planck mass, g ⋆ the number of relativistic degrees of freedom and where we used Friedmann's equation ⋆ .The magnitude of the friction coefficient, denoted by β, varies depending on the particle physics model [23,24].In our study, we assume β ≪ 1 and briefly touch upon the implications of friction on PBH formation towards the end of the paper in Sec.V C. We assume that the energy bias V 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, where A is the area parameter.Based on numerical simulations, the area parameter for a Z 2 symmetric model is approximately A ≃ 0.8 ± 0.1 [25].In contrast, Z NDW symmetric models with N DW > 2, where cosmic strings attached to N DW DWs are present, tend to yield larger values A = (0.37 ± 0.04)N DW [26].DWs annihilate when the vacuum pressure P V = V bias dominates over the pressure P T = C d σ/L arising from their surface tension σ, after the time where the factor C d ≃ O(a few) can be inferred from numerical simulations [26] and which we set to C d ≃ 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 2 Pl H 2 ≃ V bias around the time which upon comparing to Eq. ( 3), leads to the condition In the absence of a bias, the DW network would dominate when 3M 2 Pl H 2 ≃ σ/L after the time

A. 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 C. We denote by R(t) = a(t)χ(t) the physical radius of a spherical DW and by χ(t) its comoving radius.We introduce the wall Lorentz factor γ = 1/ 1 − (a χ) 2 .The EoM of a spherical DW can be determined through the Israel junction conditions [7][8][9][27][28][29][30][31], 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] χ + (4  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 baby in Eq. ( 21) of formation of a baby-universe.RIGHT: Composition of the DW mass at the time tPBH when it enters its Schwarszchild radius.The bulk component only accounts for 25 % of the DW mass.Another 25% is given by gravitational binding energies M bind and 50% of the DW mass is given by its surface term M bdy .
where a(t) is the scale factor of the universe inside the DW, the EoM in Eq. ( 7) can be rewritten as where the dot now denotes the derivative with respect to τ , H is defined below, and α ann is the DW network energy fraction at t ann where we used ρ tot ≃ 3M 2 Pl /4t 2 and Eqs.(2,3,4,6).Upon approximating ρ rad (t ann ) ≃ ρ tot (t ann ), which is exact if t dom ≫ t ann , the dimensionless Hubble factor reads

B. PBH collapse time
Closed DWs collapse into PBHs at a time t PBH after that they shrink below their Schwarzschild radius, 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] where is a bulk term, is a boundary term, and with H b = V bias /3M 2 Pl , includes the repulsive surfacesurface 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 PBH ) = t PBH (PBHs cannot form faster than what causality allows) we find that PBHs form after the time, see the left panel of Fig. 2: which is around the time t dom in Eq. ( 4) when DW interiors become dominated by the vacuum energy V bias .If we had approximated the DW mass by its volume term M ≃ M bulk , we would have found that PBHs are twice longer to form t PBH ≃ 2t 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.

C. Late-annihilators
Typical DWs annihilate at t 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 ≳ 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 ≃ 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 ≃ t (black line).The radius R PBH ann at t ann of the smallest DW which collapse into PBH is plotted in the right panel of Fig. 4 as a function of α ann , for different approximations for the Misner-Sharp mass in Eq. (14).We also show in gray the approximation R(t) ∝ t1/2 only accounting for the Hubble flow and neglecting effects from the bias V bias and surface tension σ which become important when R approaches R ∼ t.The most precise approximation shown with the black line is well reproduced with the fitting function:

D. Baby-universe
If the radius R of a DW becomes larger than the cosmological horizon H −1 b of an inflating universe with vacuum energy density V bias , which would occur around the time t baby 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 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.[7,19].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 spacetime 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 babyuniverse form roughly around twice the time of PBH formation, as given in Eq. ( 18): Since 2GM bulk = H 2 b R 3 in Eq. ( 14), the time t 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 The gray line shows the would-be DW trajectory from completely neglecting effects from the bias V bias and surface tension σ, i.e. solely accounting for Hubble flow R(t) ∝ a(t).RIGHT: We show the minimal radius that spherical DW must have at tann in order to collapse into PBH after tPBH.The Schwarzschild radius is calculated with volume term only (purple), plus surface term ("bdy") (blue) and binding terms (black).Due to the relationship 2GM bulk = H −1 b , 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.approximated by the fitting function: which for the same reason as before corresponds to the purple line in the right panel of Fig. 4. The abundance of lateannihilators, 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].

A. 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 ∼ t.The correlation length is usually estimated from the network energy density [13,26,46,47]: where A is the area parameter introduced in Eq. ( 2).Using results from field-theory lattice simulations, for Z 2 -symmetric models one finds L/t ≃ 1.25 +0.18 −0.14 [26] or L/t ≃ 1.14 ± 0.04 [47], while Z NDW -symmetric models with N DW > 2 are associated to a shorter correlation length, L/t ≃ (2.70 ± 0.30)/N DW [26].The values of the field ϕ(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 [45,[48][49][50].
Assuming a 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 bias , the probability of false vacuum occupation reads [18] where T is the temperature, ∆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 bias ≪ 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 × N 3 .The probability that any selected lattice site is part of a s-cluster is given by P s = s × n s .We refer the reader to [43][44][45] for reviews on percolation theory.

B. 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 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 min .In doing so, any deviation from sphericity would add more mass than required for the collapse to occur.We define s ball (r min ) as the minimal number of lattice sites of size L to fully cover a ball of radius R min , with r min ≡ R min /L.The probability that all those sites are occupied is p s ball (rmin) 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 ball where c 0 = π/3 √ 2 ≃ 0.74, we obtain the number of sclusters large enough to contain a ball of radius r min (in unit of L): where we normalized by the total number of sites N 3 .Those cluster carry an energy: where ϵ is the energy of one occupied sites.This must be compared with the energy of the DW network: where the two terms accounts for the bulk and surface energy of the DW network which at t ann are exactly equal.We deduce that the energy fraction contained in closed DW large enough to contain a spherical ball of radius r min is given by: where we recall r min ≡ R min /L, c 0 ≃ 0.74 and p = 0.5.In App.B, we calculate the minimal number s latt ball (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: We observe that for small radii r ∈ [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 ∈ [10, 10 3 ].Conversely, for larger radii r ≳ 10, it approaches its lower limit, representing perfect packing.The stair-case behaviour of s latt ball (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 smooth ball (r) obtained under the double condition to interpolate the analytical sum s latt ball (r) in Eq. (B3) as best as possible while keeping it always smaller than its upper limit s smooth ball (r) ≤ ⌈2r⌉ 3 .
The resulting smooth interpolation shown with black line in the left panel of Fig. 5 reads: where Lr ≡ Log 10 (r) and r ≡ R/L.The energy fraction F(r min ) of closed DWs larger than a ball of radius r min is obtained from plugging s smooth ball (r 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 latt ball (r min ) in Eq. (B3) of App.B as well as its upper and lower limits in Eq. ( 29), resulting in lower and upper limits on F shown with red and blue lines respectively.Additionally, we depict with a gray line the energy fraction F of late-annihilators which would have been obtained if results from MCMC simulations presented in Eq. (A3) of App.A were plugged in Eq. ( 28).

V. PBH AND WORMHOLE ABUNDANCE
Now that we are in possession of the critical size R PBH min and R baby min 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.

A. PBHs
The PBH contribution to the DM abundance is where ρ DM (t) is the DM energy density, ρ coll (t) is the energy density stored in collapsing DWs and F coll is the collapsing Left: Minimum number of lattice sites, denoted as s 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 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. (A3) in Eq. ( 28), tends to overestimate the fraction 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.fraction at the onset of annihilation, 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 ⋆ (T )T 4 respectively and matching them at t unbias dom in Eq. ( 6) and matter-radiation equality when T eq ≃ 0.80 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: with and The different temperatures can be calculated from the characteristic times in Eqs.(3,4,6,13) using T ≃ 1.23(M Pl /t) 1/2 /g ⋆ (T ) 1/4 , with g ⋆ (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 X s, * with X the associated epoch in Eq. ( 34).In Eqs. ( 34) and ( 35), behind the sign ∼ 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 ann ∼ t dom ∼ t unbias dom , which together with Eq. ( 18), implies G ∼ 1 and R ∼ 1.The collapsing fraction F coll is given by the fraction of closed DWs in Eq. ( 28) with a radius larger than the critical threshold: where p = 0.5.The function s 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 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 PBH ann at t 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):  [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 [64][65][66].The brown area shows where DM would over-close the universe.Above the purple dashed line, more than one baby-universe f baby ≳ 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 tann ≳ t 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 O3-run of LIGO-Virgo-Kagra (LVK) [75].We use that no SGWB with a signal-to-noise ratio larger than SNR ≤ 5 has been detected in NG15 and SNR ≤ 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 eff ≲ 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 Z2-symmetric network, see Eq. ( 23).
The minimal PBH mass is given by the mass inside the Schwarszchild radius at horizon crossing R sch (t PBH ) = t PBH : where t PBH is given by Eq. ( 18) and M L ≃ 3.00 × 10 −6 M s where M L and M s are Earth and Sun masses.The PBH mass can also be calculated as M PBH ≃ 4 × 4πt 3 PBH V 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 ≫ R PBH ann 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.

B. Wormholes
The critical DW radius R baby min 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 PBH min 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 baby of baby universes formed in our past lightcone is simply: where the DW fraction F(R min ) is defined in Eqs. ( 28) and (30), and the baby threshold R baby ann is given in Eq. ( 22).The function N patches (T dom ) is the number of Hubble patches, at the time of collapse in Eqs. ( 21) and ( 4) when the temperature is approximately T dom , in our past light-cone.It reads: where h ≡ H 0 /100 km/s/Mpc with H 0 the Hubble constant today.We approximated g * (T ) ≃ g * ,s (T ) and used g * s (T 0 ) ≃ 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.

C. 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.
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 sch /t ∼ (t/2t dom ) 2 , possibly much smaller than t, where we assumed R sch = 2GM with M ∼ 4πt 3 /3.PBH formation well within the horizon can only occur if they are spherical enough to shrink by a factor (2t dom /t) 2 ≫ 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 PBH ann 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 min , their abundance F(R min ) being derived in Sec.IV.If R min = R PBH ann , those DWs will all unavoidably collapse into PBHs whatever their shape.This is because nonsphericity will participate to a mass larger than M (R 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 F(R 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 ≈ σ/ρ DW , as presented in Eq. ( 23).This approach benefits from previous lattice calculations, where the relationship L/t ≈ (2.70 ± 0.30)/N DW was established for the Z NDWsymmetric model [26].This result implies that the PBH abundance dramatically decreases as the number N DW of nearlydegenerate vacua increases.A more precise determination of the correlation length, along with its dependency on N DW , is left for future studies.Also the precise dependence of the lateannihilator 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 DW > 2 due to the substitution of p = 0.5 with p = 1/N 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 ∼ 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 bias is absent at the time of network formation.The presence of a large initial bias V bias already active at the time of network formation [37,50,80], i.e. such that t ann < t form , would impact the occupation probability p in Eq. ( 24) and reduce the fraction 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.

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 Ω 0 GW can be related to the GW power spectrum at the annihilation epoch Ω ann GW by: where D is the redshift factor: where Pl H 2 0 with H 0 ≃ 100 km/s/Mpc, g s, * (T 0 ) ≃ 3.38 and g ⋆ (T 0 ) ≃ 3.94, assuming N eff ≃ 3.045 [100,101].The GW spectrum produced by long-lived DWs annihilating at t ann follows from the quadrupole formula [25,[102][103][104]: The dimensionless quantity ϵ GW encodes the deviation from the quadrupole formula and is fitted on lattice simulations ϵ GW ≃ 0.7 ± 0.4 [25].We model the spectral function as where the IR slope is Ω GW ∝ f 3 to respect causality [105][106][107][108] and the UV slope is Ω GW ∝ 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]: 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 Signalto-Noise Ratios at SNR = 5 and SNR = 2 for NG15 and LVK O3, respectively, and (SNR = 1, T = 10 years) for SKA, as well as (SNR = 1, T = 20 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 Z N symmetry with N > 1.In presence of a vacuum energy difference V 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 sch ∝ 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 ≃ t.This seems to imply that there must be a time t PBH when most of the DWs enter their Schwarschild radius and collapse into PBHs.This time is closely related the time t dom ∼ t 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 ann ≲ t PBH to prevent DW domination.In this work, we find that only super-horizon DWs continue to grow like R ∝ a(t)  [74] and run O3 from LIGO-Virgo-Kagra (LVK) [75], assuming the Signal-to-Noise Ratio detection thresholds SNR = 5 and 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 SNR = 1 (equivalent to 1σ deviation from noise) after T = 20 years of observation (T = 10 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 ∈ [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.
or slower while sub-horizon DWs shrink.This implies that only super-horizon DWs which are larger than a given threshold, R ≥ R PBH ann , can enter the horizon after t PBH and collapse into PBHs.We call this population late-annihilators.We calculate the PBH formation threshold R PBH ann 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 α 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 M s ∼ 10 9 g.This occurs if the bias energy density is larger than V bias ≳ (10 12 GeV) 4 , see Figs. 6 and 7.
Finally, the µ-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 µ-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).

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 = 0.0501 s τ exp {−0.6299Ps σ (Ps σ + 1.6679)} , (A1) where τ ≃ 2.17, σ ≃ 0.48, P ≡ (p − p c )/p c .For p > p c , the average cluster radius R s is given by [45] For an occupation probability p = 0.5 and the percolation threshold of a cubic lattice p c = 0.311, Eqs.(A1) and (A2) simplify to n s ≃ 0.0501s −2.17 exp −0.2326s 0.96 − 0.6385s 0.48 , (A3) and R s = 0.863s 1/3 .Since p > p c , most of the lattice points belong to the infinite cluster.We note P ∞ 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 ∞ , or occupied in a finite cluster with a probability s sn s .Summing all the possibilities leads to [43] Plugging Eq. (A3) into Eq.(A4) for p = 0.5, we calculate P ∞ ≃ 0.92.For the cosmological system (N → +∞), it means that 92% of the false vacuum region (with vacuum energy V 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. (A3).However, the formula given in Eq. (A3) 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 ∝ R 3 s in Eq. (A2) 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. (A3), 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 B. 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.

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 ⌈r⌉ layers denoted by j ∈ 0..⌊r⌋ , 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 = r 2 − j 2 .Denoting the column number by i ∈ 0..⌊r⌋ , the number of cubes filling a row, shown in red color, is r 2 j − i 2 .We deduce the total number of cubes filling the ball for the configuration shown in the left panel of Fig. 8: (B1) 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 (Y Z) 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⌊r + 0.5⌋ + 1.The remaining number of cubes in those three orthogonal layer outside the 3axis cross is 4 ⌊r⌋ i=0 r 2 − (i + 0.5) 2 + 0.5 .Finally, accounting for the remaining cubes in the 8 octants outside the three internal layers give:   The function s latt ball (r) with r = R/L is shown with a green line in Fig. 5.We find that at small radius ≳ 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 latt ball (r) ≳ ⌈2r⌉ 3 .

(B4)
At large radius r ≫ 1, the number of cubes asymptotes the perfect-packing (lower) limit shown in blue for which the tiny cubes smoothly fit inside the ball: 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. (C4) is the BBN constraints, from the scenario in which DWs reheat to SM, in which case Eq. (C4) applies only if DWs annihilate below the neutrino decoupling temperature T ann ≲ 1 MeV.

Figure 1 .
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 PBH ann at the onset of the annihilation phase, which are dubbed "late-annihilators" in the Orange region, enter inside their Schwarzschild radius R ≲ 2GM in the Blue region, leading them to collapse into PBHs.DWs with radius larger than R baby ann become larger than the cosmological horizon H −1

Figure 2 .
Figure 2. LEFT: We show the time tPBH 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 baby in Eq. (21) of formation of a baby-universe.RIGHT: Composition of the DW mass at the time tPBH when it enters its Schwarszchild radius.The bulk component only accounts for 25 % of the DW mass.Another 25% is given by gravitational binding energies M bind and 50% of the DW mass is given by its surface term M bdy .

Figure 4 .
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 bias and surface tension σ, i.e. solely accounting for Hubble flow R(t) ∝ a(t).RIGHT: We show the minimal radius that spherical DW must have at tann in order to collapse into PBH after tPBH.The Schwarzschild radius is calculated with volume term only (purple), plus surface term ("bdy") (blue) and binding terms (black).Due to the relationship 2GM bulk = H −1 b , 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.

1 Figure 5 .
Figure 5. Left: Minimum number of lattice sites, denoted as s 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 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. (A3) in Eq. (28), tends to overestimate the fraction 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.

Figure 6 .
Figure 6.Exclusion constraints on DW network occupying an energy fraction αann defined in Eq.(11) at the onset of the annihilation phase driven by a bias energy density difference V bias between distinct vacua.PBHs are produced around the temperature V 1/4 bias 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 µ-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[64][65][66].The brown area shows where DM would over-close the universe.Above the purple dashed line, more than one baby-universe f baby ≳ 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 tann ≳ t 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 O3-run of LIGO-Virgo-Kagra (LVK)[75].We use that no SGWB with a signal-to-noise ratio larger than SNR ≤ 5 has been detected in NG15 and SNR ≤ 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 eff ≲ 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 Z2-symmetric network, see Eq.(23).

Figure 7 .
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 O3 from LIGO-Virgo-Kagra (LVK)[75], assuming the Signal-to-Noise Ratio detection thresholds SNR = 5 and 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 SNR = 1 (equivalent to 1σ deviation from noise) after T = 20 years of observation (T = 10 years for SKA)[91] and the perfect subtraction of all the stochastic astrophysical GW foregrounds, see e.g.Fig.6in[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 ∈ [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.

Figure 8 .
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 ball,1 (r) and s ball,2 (r).

2 30 g ⋆ (T )T 4 .
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 where ρ γ is the photon number density.We introduce the DW energy fraction in unit of radiation energy density at annihilation temperatureα DW (T ) = ρ DW (T ) π (C2)where T is the SM photon temperature.From Eq. (C1) and Eq.(C2), the maximal DW contribution to N eff occurs at an-BBN bound ∆N eff ≲ 0.3[77,78], the effective number of extra relativistic degrees of freedom must be evaluated below the neutrino decoupling temperature where g ⋆ (T < T dec ) ≡ 2 + (7/8) • 6 • (4/11) 4/3 ≃ 3.36.Hence, we obtain ∆N eff = 7.4 α GW ≲ 0.4, (C4)