Collapsing Domain Wall Networks: Impact on Pulsar Timing Arrays and Primordial Black Holes

Unstable domain wall (DW) networks in the early universe are cosmologically viable and can emit a large amount of gravitational waves (GW) before annihilating. As such, they provide an interpretation for the recent signal reported by Pulsar Timing Array (PTA) collaborations. A related important question is whether such a scenario also leads to significant production of Primordial Black Holes (PBH). We investigate both GW and PBH production using 3D numerical simulations in an expanding background, with box sizes up to $N=3240$, including the annihilation phase. We find that: i) the network decays exponentially, i.e. the false vacuum volume drops as $\sim \exp(-\eta^3)$, with $\eta$ the conformal time; ii) the GW spectrum is larger than traditional estimates by more than one order of magnitude, due to a delay between DW annihilation and the sourcing of GWs. We then present a novel semi-analytical method to estimate the PBH abundances: rare false vacuum pockets of super-Hubble size collapse to PBHs if their energy density becomes comparable to the background when they cross the Hubble scale. Smaller (but more abundant) pockets will instead collapse only if they are close to spherical. This introduces very large uncertainties in the final PBH abundance. The first phenomenological implication is that the DW interpretation of the PTA signal is compatible with observational constraints on PBHs, within the uncertainties. Second, in a different parameter region, the dark matter can be entirely in the form of asteroid-mass PBHs from the DW collapse. Remarkably, this would also lead to a GW background in the observable range of LIGO-Virgo-KAGRA and future interferometers, such as LISA and Einstein Telescope.


I. Introduction
Physical models that feature the formation of cosmic domain wall (DW) networks have typically been seen as problematic due to the so-called domain wall problem, the fact that the network tends to dominate the energy density of the universe [1].However, if domain walls are biased and annihilate, this problem turns into a virtue, as the network naturally tends to be an abundant component in the universe before its collapse, and is thus easier to probe.
Gravitational waves (GWs) are one of the potential signatures [2][3][4].The spectrum is stochastic and analyses with LIGO-Virgo O3 data already place constraints on the parameters of the network [5].Recently, the evidence for nano-Hz GWs at Pulsar Timing Arrays (PTAs) [6][7][8][9] brought renewed interest in this possibility: DW networks that annihilate around the QCD phase transition provide a good explanation of the signal and outperform several other models [10,11].In case of detection of GWs, it is also crucial to find additional signatures, that can help in selecting DWs over other early Universe sources.Dark radiation and collider signals may be some of such signatures [10], as well as the production of Primordial Black Holes (PBHs) from the collapsing network [12,13] (see also [14][15][16][17]), although the resulting PBH abundance is subject to large uncertainties.
Overall, these aspects motivate a detailed investigation of the evolution of a DW network during its collapse phase (see [18][19][20][21][22][23] for previous work) and of its gravitational relics, i.e.GWs (see [24][25][26][27][28][29]) and PBHs, which we aim to perform in this work.We consider the simple case of DWs with a double-well potential, with an energy per unit area (tension) σ, and where the potential is slightly tilted by a term of size ∆V , i.e. a bias such that one of the Z 2 -symmetric minima becomes a false vacuum [30] (see e.g.[31][32][33][34] for other mechanisms to have viable long-lived walls).
We simulate the corresponding DW network throughout the formation, scaling and annihilation regimes, using field theory numerical simulations in an expanding radiation-dominated background, with box size up to 3240 3 , while computing the GWs radiated throughout the evolution.
Until recently, most analyses assumed that GWs in this scenario are radiated until the pressure on the walls caused by the bias overcomes the self-acceleration due to the wall's tension, i.e. when ∆V = σH, where H is the Hubble rate. 1 Before this epoch, the network is in the socalled scaling regime, with most of its energy density in a fixed O(1) number of walls per Hubble patch.However, the collapsing network consists of large DWs of various shapes, which contain False Vacuum (FV) pockets, that shrink to small sizes after the scaling epoch.These last stages of evolution can certainly source GWs in addition to the ones from the so-called scaling epoch, and has been found in recent work [28].An order-one change in the estimate of the time scale for GW production can lead to an order-of-magnitude enhancement of the final GW signal.Our simulations improve on the determination of such time-scale, while also providing new insights into the properties of the network at the onset of annihilation and highlighting the role of the kinetic energy of the scalar field in the production of gravitational waves in the final stages of the network annihilation.
Our numerical results will also allow us to establish the time evolution of the decaying network, in particular of its FV pockets.Indeed, the collapse of the network takes some time: the abundance of Hubble-sized FV pockets at some point drops very quickly, but a small fraction of rare super-Hubble sized pockets survive for a longer time as they must cross the Hubble radius to annihilate.Their abundance dramatically decreases in time, but their likelihood to collapse into a BH grows instead simply because the Schwarzschild radius associated with the FV pocket grows faster in time than the Hubble length.This may result in a tiny population of BHs at formation, but potentially large at present if the network collapses in the early universe (this formation mechanism is similar to the one in [36,37] for isolated DWs, see also [38], with the crucial difference that for a network the collapse is in general far from spherical).
We provide an analytical understanding of this picture, which complements our numerical findings to provide an important step forward in the estimate of the PBH abundance.Nonetheless, large uncertainties in the final PBH abundance persist, due to an exponential sensitivity to parameters and the departure from spherical collapse, which at this point is still difficult to quantify.
With these new estimates of both GW and PBH relics, we then analyze the phenomenological consequences of a generic DW network that annihilates at different epochs in the early Universe.First, we assess the viability of the DW interpretation of the PTA signal in light of PBH production.Second, we discuss PBHs from DWs as a candidate for the observed dark matter, with a possible correlated GW signal at interferometers such as LIGO-Virgo-KAGRA (LVK) [39], LISA [40], Einstein Telescope (ET) [41] and Cosmic Explorer (CE) [42].
The paper is organized as follows.We summarize the evolution of annihilating DW networks in Sec.II, highlighting the important time scales in the problem.We present the results from our numerical simulations in Sec.III.We present an analytical model to account for the late FV pockets in Sec.IV.We conclude in Sec.V with a discussion on the phenomenological impact.

II. Domain Wall Networks: Scaling and Annihilation regimes
In this work, we focus on a simple model exhibiting DWs: a real scalar ϕ with a Z 2 symmetry ϕ → −ϕ and with a double well potential keeping in mind that several aspects of our discussion may also apply to other models (e.g. with more minima or from axion potentials).In this model the DW tensionthe energy per unit area of the walls -is σ = √ 2 λ 3 v 3 and the scalar mass squared in the minima is m 2 = 2λv 2 .In addition, we will assume a small bias term in the potential that breaks such a symmetry, of size ∆V , that will lead to annihilation of the walls.
We assume that a network of walls is formed by a phase transition in the early universe during radiation domination, i.e. we start with zero initial field plus small random fluctuations and the walls are formed via the so-called Kibble mechanism [1,43].The subsequent evolution follows a sequence of events, schematically represented in Fig. 1, and explained in detail below: (i) the network reaches the so-called scaling regime, (ii) the network starts annihilating when the bias term becomes relevant, (iii) a peak of gravitational waves is produced, (iv) some very rare domain walls that survive for a longer time may collapse and give rise to PBHs.

A. Scaling regime
Soon after formation, the DW network reaches a selfsimilar 'scaling' regime characterized by, on average, about one Hubble-sized DW per Hubble volume.If we consider a generic network of total comoving area A in a box of comoving volume V , its total energy is σAa 2 and so its total energy density is given by where a is the scale factor, H is the Hubble parameter and we introduced the so-called area parameter which is a dimensionless number related to the area density of the network.During scaling A is expected to be of order unity [26].
One of the remarkable signatures of the DW network is the stochastic spectrum of GWs that it creates.In the scaling regime, the spectrum Ω gw (k, t) ≡ ρ −1 c dρ gw /d log k, where ρ gw is the energy density in GWs and ρ c is the critical background density, as a function of the wave number k and cosmic time t, peaks at the Hubble scale and previous studies have found that in scaling its amplitude at the peak is given by [44] Ω (scaling)   gw where ϵ ≃ O( 1) is an efficiency factor extracted from the numerical simulations and α(t) ≡ ρ dw /ρ c is the fraction of the total energy density stored in the walls.The fraction α increases over time and thus one expects the GW spectrum to have a peak around the annihilation time of the network.

B. Annihilation Phase
The annihilation phase occurs roughly once the Hubble rate becomes smaller than the pressure acceleration, that is, for conformal time η ≳ η ∆V , defined by one then expects the network to start collapsing.So far in most literature it has been assumed that the peak of GW production occurs exactly at the time given by eq. 5, identified with the annihilation time of the network.As we will see, it is important to determine precisely both the annihilation time of the network and the time of emission of the peak of GWs, since the final GW spectrum grows as η −4 , as one can see by extrapolating the scaling properties in eqs. 2 and 4.
A useful quantity to monitor the degree of annihilation of the network is the fraction of volume occupied by the False Vacuum, F fv (η).In a Z 2 model, F fv = 1/2 in scaling.There is some literature on how F fv decays during the annihilation phase [18][19][20][21]32] (see also [22]) with no current consensus.In the following sections, we will see that the network decays exponentially in the annihilation phase according to with η ann possibly differing from η ∆V by O(1) factors, and p ≃ 3.

III. Numerical results
We now present the results of our lattice simulations of DW networks, obtained by means of the Cosmolattice code [45,46].We set initial conditions at the initial conformal time η i in radiation domination such that m = H(η i = 0), and a(η i = 0) = 1.We assign a white noise spectrum of small Gaussian fluctuations in Fourier space to the scalar field, while setting its homogeneous component at ϕ = 0.The number of gridpoints N 3 and the comoving box size L determine the comoving lattice spacing ∆x = L/N .We set ∆x and the duration η f of our simulations such that the physical lattice spacing a(η f )∆x at the end of the simulation is smaller than (or equal to) the domain wall width δ w ∼ m −1 , and impose that the simulation box contains at least one Hubble patch at the final time.From here on we set v = m = 1, which corresponds to the choice λ = 1/2 for the quartic coupling.With these choices, a(η) = 1 + η, H(η) = (1 + η) −2 and the domain wall tension is simply σ = 2/3.The maximal dynamical range that can be probed with the simulation is then η f,max = √ N − 1, obtained by choosing L = √ N , such that at η f there is only one Hubble patch in the box and a(η f )∆x = δ w .The field evolution is performed using the leapfrog algorithm with a time step ∆η ≲ ∆x/2, such that the Courant criterion is well satisfied.More details about the numerical setup can be found in Appendix A.
Here we shall consider a cubic bias, i.e.V bias = qϕ 3 , rather than the linear term.The reason for this choice is only technical: a linear bias shifts the location of the maximum of the potential, thereby introducing a bias in the population of the two minima already at the early times of the simulation.Since the limited dynamical range that can be simulated requires a sizable ∆V , such a population bias would prevent the formation of a network and/or alter its evolution.The cubic bias allows to overcome this problem, as the maximum is not displaced, and for sufficiently small initial field fluctuations the system does not notice the asymmetry introduced by V bias in the initial steps. 2The size of 2 An alternative strategy is to use a time-dependent linear bias [28], which is initially negligible and becomes important only after a certain activation time.However, we have found that this technique introduces additional uncertainties in the final results, due to different possible choices of the activation time of the bias.
The scenario relevant to our analysis is that of a domain wall network that achieves the scaling regime sufficiently before it starts collapsing as a consequence of a pressure bias between the two vacua.For the potential adopted in this work, this occurs for ∆V /V (0) ≲ 0.007, corresponding to η ∆V ≳ 18 and thus we focus on this range of bias sizes.For larger bias, the network does not fully achieve scaling for a sufficient time before collapse.While this scenario can certainly occur, it is characterized by some residual dependence on initial conditions, which makes it difficult to extract general conclusions.
We show the evolution of the average field value together with its standard deviation as a function of conformal time in Fig. 2, as obtained in simulations with a box size (N = 3000) 3 and maximal time η f = 55.
Here we fixed initial conditions and varied the size of the bias potential, such that η ∆V = {19, 22, 25}.The following features can be clearly distinguished: first, the field is initially localized very close to the maximum of the potential, until around η ≃ 6 it relaxes to the two minima.Field oscillations are sizeable until η ≃ 10, when they are significantly diluted by Hubble friction.The standard deviation then begins to shrink after η ∼ 20 and for η ≳ 30 only the leftmost minimum is populated, signaling that the network is rapidly dissolving under the action of the bias.For comparison, the behavior for vanishing bias is also shown (dot-dashed curve).In this case, a tendency towards the right-most minimum occurs at around η ∼ 15, which is to be attributed to the relatively small (and decreasing) number of Hubble patches at those times, i.e. ∼ (L/η) 3 ≲ 50 at η ≳ 15.
The deviations in the network evolution in the presence of a bias are best understood by focusing on the following quantities.

A. False Vacuum Fraction
First, we look at the fraction of volume in the false vacuum F fv , which is numerically obtained as the fraction of the simulation grid points where ϕ > 0, shown in Fig. 3 (blue curves).As expected, initially F fv = 0.5 in all simulations.This remains approximately true for the simulation without a bias, although a slight deviation to larger values is observed at late times, corresponding to the aforementioned small number of Hubble patches near the end of the simulation.In simulations where a bias is included, F fv decreases rapidly after a time which depends on the size of the bias, obviously the later the smaller ∆V .In this work, we are mostly interested in false vacuum regions that are at least Hubble-sized, since their radius becomes equal to their Schwarzschild radius if their energy density becomes comparable to the background at Hubble crossing, as we will explain in more detail in Section IV.When the false vacuum fraction drops below the inverse number of Hubble volumes in the box, given by n H ≃ (L/(1 + η)) 3 and shown by the black dot-dashed curve in Fig. 3, such regions no longer exist in our simulation.For the bias sizes of interest, this occurs around η ≃ 30.At later times, the remaining false vacuum regions are in sub-Hubble structures.Correspondingly, a much steeper decrease of F fv is observed at late times than at early times, when super-Hubble false vacuum regions can still be present.Thus we attempt fitting the numerical results until the time at which F fv ≃ n −1 H with eq. ( 6) where η ann and p are fitting parameters.The result of this procedure is shown by the orange curves in Fig. 3.The fitting function above provides an excellent fit to the early time data, and we find p ≃ 3.3 − 3.5.
To investigate the dependence of these results on the simulation box, we increase the number of Hubble patches in our box by increasing L (and N to a smaller extent) at the price of slightly worsening the spatial resolution, and thus limiting the dynamical range of our simulations.Guided by the previous discussion, we choose L and η f such that δ w /(a∆x) ≃ 1 at time η f ≃ 35, corresponding to L ≃ 90, with N = 3240.These choices increase the number of Hubble patches by a factor of ≈ 4.5 with respect to the results in Fig. 3.The new results are reported in Fig. 4, together with the fitting curves.It can be appreciated that F fv now remains almost constant for the entire simulation time in the absence of bias, thereby confirming that the previously observed deviation is due to the limited number of Hubble volumes.The most relevant result of this improved analysis is a decrease of the inferred value of p (which is closer to analytical expectations discussed in Sec.IV).
We then perform several realizations with increased number of Hubble patches, for several values of bias size and changing the random seed that sets the initial conditions, to estimate numerical uncertainties in our results.We report results in Table I.Averaging over all realizations, we infer When comparing η ann to the rough expectation η ∆V : σH = ∆V , we find a slight delay η ann ≈ 1.5η ∆V . 3 3 In this regard, we disagree with [23] on the dependence of ηann with ∆V .3), as a function of conformal time.Note that at early times, η ≲ 10, there is a transient which does not carry physically relevant information.

B. Area parameter
In the absence of a bias, a domain wall network achieves a scaling regime where the area parameter A remains approximately constant with time.The energy density in the domain walls during the scaling regime, is determined by eq. ( 2), which is generally smaller than the total energy density stored by the scalar field.We show in Fig. 5 the area parameter extracted from our high resolution simulations, in the absence of a bias (solid magenta curve).As expected, it remains approximately constant for η ≳ 25, saturating for this particular realization to a value A ≃ 0.9.These findings agree roughly with those of [26].
The corresponding evolution in the presence of the bias is shown by the blue dashed, dotted and dot-dashed curves, for several bias sizes.One can appreciate that the biased network follows the unbiased one until η ≃ 14−17, depending on the size of the bias. 4

C. Energy density
The evolution of the energy density of the scalar field is shown in Fig. 6 (left) for unbiased (solid) and biased (dashed, dotted, dot-dashed) potential.These results are obtained in our longest simulations, corresponding to the choice L = √ N .As expected, in the unbiased case the total energy density in the scalar field remains approximately constant after η ≳ 30, with a final value ρ tot ≈ 3.5 σH.On the other hand, in the biased case the total energy density decreases sharply, together with the decrease of the vacuum contribution from the bias potential (orange curves), due to the annihilation of the network.The evolution of the three components of the energy density, i.e. kinetic, gradient and potential energy, is shown in Fig. 6 (right) for the unbiased (solid) potential as well as for an example case of biased potential (dashed) with the choice η ∆V = 22.
The following observations can be made: in the unbiased case, the gradient energy rapidly saturates to a constant value ρ grad ≃ 1.4 σH, while the potential energy density ρ pot decreases slowly until the end of our simulations (here by potential we denote only the Z 2 symmetric term, whose behavior is shown by the purple curve).The latter decrease may however partially be a numerical effect since it occurs almost at the end of the simulation where the domain wall width becomes comparable to the lattice spacing of the simulation.Overall, these two components make most of the energy density of the scalar field, in agreement with the intuition that most of the energy density is in domain walls, and in particular ρ grad + ρ pot ≃ 2.6 σH at the end of the simulation.On the other hand, the kinetic energy decreases rapidly until it saturates to an approximately constant value ρ kin ≃ 0.9 σH, thereby making a subleading contribution to the energy density.
The situation in the biased case is strikingly different, where deviations from the unbiased scenario occur around η ≳ 20: most importantly, the kinetic energy stops decreasing and quickly begins to increase, while the potential energy decreases sharply.The former overcomes the latter around η ≳ 34.This is very close to the value of the annihilation time η ann inferred from our fit of the false vacuum fraction, nicely confirming that this time scale indeed characterizes the annihilation of the network when the vacuum pressure from the bias term (shown by the dashed brown curve) accelerates the walls, thereby increasing the kinetic energy.The gradient energy initially remains constant, before sharply decreasing at η ≳ 40.This is easy to interpret, as the existence of the network is the source of gradient energy, which is thus quickly dissipated away once the walls annihilate.Eventually, at η ≳ 45 also the kinetic energy starts decreasing.We have checked that for η ≳ 45 both the potential and kinetic components decrease approximately like non-relativistic matter, as expected since we are working with a massive scalar field with m ≫ H at the end of the simulation.Overall, the kinetic energy dominates the energy density of the biased domain wall network at the end of our simulations.Notice also that the bias term very rapidly vanishes after η ≳ 40, corresponding to the exponential decay of false vacuum regions.
Our findings on the behavior of the energy density in the unbiased case may seem different from the common lore in the literature, which mostly adopts ρ dw ≃ 2AσH.In our simulations, we find that the total energy density in the scalar field, including all the simulation box, is roughly twice as large as the estimate above, ρ tot ≈ 3.5σH.One should nonetheless notice that: 1) the commonly adopted estimate is expected to apply only to static domain walls, i.e. it neglects the contribution of the kinetic energy, which in our simulations accounts for ρ kin ≈ 0.9σH; 2) the simulation box includes regions where |ϕ|> 1, which cannot be attributed to domain walls, but rather to scalar waves.The energy density in this region of field space is reported in Fig 7 , for a simulation with a large number of Hubble patches.Its size at the end of the simulation is ρ scal ≈ 0.8σH, mostly coming from kinetic and potential energy.Ignoring the kinetic part, it contributes ≈ 0.5σH.Therefore our larger total energy density in the scalar field is easily explained in terms of the two contributions above (plus a small unavoidable contribution from scalar waves in the region |ϕ|≤ 1). 5 We conclude that, according to our simulations, domain walls carry an energy density ρ dw ≃ 2.4σH.If interpreted in terms of a relativistic correction to the standard formula, i.e. ρ dw ≃ 2Aγ 2 σH, it implies γ ≈ 1.2.

D. Gravitational Waves
It has long been appreciated that a domain wall network acts as an efficient source of GWs during its evolution.Previous numerical calculations have mostly focused on the contribution from the scaling regime.In our simulations, we are able to extract the energy density radiated in GWs throughout the annihilating phase as well.Due to higher memory consumption, our simulations including GWs are limited to N = 2040, with a maximal simulation time η f ≲ 45.Additionally, to speed up the calculation, we only start the numerical computation of the GWs at the latest times (η > 35).This is justified since we find that, unlike previous estimates, most of the GWs are emitted during the annihilation epoch rather than in the scaling regime (see also [28]).
A simple estimate of the maximal energy density fraction (i.e. the energy density in GWs compared to that of the radiation background) is provided by the quadrupole formula, which gives Ω gw (η) ∼ 3/(32π)α 2 tot (η), where α tot ≡ ρ ϕ /(3H 2 M 2 p ).This would correspond to the case in which all the energy density in the scalar field sources GWs.We compare our results with this simple estimate in Fig. 8, for three different choices of bias size for which we are able to capture the full GW production (with same initial conditions as in all previous figures).For all our choices, we find that the GW energy density fraction reaches a peak at η ≳ 40.The efficiency factor with respect to the simple quadrupole estimate at peak production is ϵ ≡ Ω gw /(3/(32π)α 2 tot ) ≃ 0.5 − 0.6.Fig. 9 shows that GWs are compatible with being mostly sourced by the kinetic component of the energy density in the scalar field, with a subdominant contribution from the gradient component, for one example value of bias size.The GW energy density fraction is computed by integrating the numerically obtained GW spectrum, although the final result is dominated by the region around the peak of the spectrum.
Our results importantly point to a mild delay between the characteristic time of the annihilation epoch η ann and the time at which most GW production occurs, η gw .In our simulations, this is estimated to be Physically, our findings suggest that efficient GW production continues throughout the annihilation epoch, and so beyond the production during scaling studied in previous works, and that most of the relic abundance of GWs is determined by the late stages of the DW network collapse.Compared to the previous literature, the delay amounts to a factor of 1.3 × 1.4 ∼ 1.8 in the time of GW emission, which was previously estimated to be η ∆V , and so, following eq.( 4), an enhancement of 1.8 4 ∼ 10 in the total GW abundance.

IV. Semi-Analytical approach
Let us now turn to a different way to analyze the problem.The network annihilation process can be viewed as a transition from DWs in the scaling regime to, eventually, a collection of rare FV pockets.During scaling, in each Hubble patch, there is typically one Hubble-sized DW, some scalar radiation and, to a good approximation, no closed DWs.Most of the energy is stored in the form of large DWs and any DW is typically separated from neighbour DWs by the correlation length, which in scaling is set by the Hubble length.
Annihilation starts when the force per unit area from ∆V is larger than σH, pushing the DWs to reduce the volume in the FV.Let us assume that this effect turns on instantaneously at η ∆V .Since the typical separation between walls is also of order η ∆V , a Hubble-sized FV pocket takes a time of about η ∆V to shrink to zero, and so one expects that the fraction of the volume in FV becomes tiny after a delay of order η ∆V , that is, at FIG. 6: Left: Energy density in the scalar field as a function of conformal time (blue), normalized to the scaling behavior σH, without (solid) and with several choices of bias sizes (dashed, dotted, dot-dashed).In orange, the evolution of the average of V bias in the simulation box is shown.Right: Components of the energy densities of the scalar field as a function of conformal time, for unbiased (solid) and biased (dashed) potential, the latter for η ∆V = 22.The vertical gray line corresponds to the value of η ann inferred from fitting the false vacuum fraction according to eq. ( 6).In both figures, N = 3000 and L = √ N .η ≃ 2η ∆V , since after that time only very rare structures survive, i.e, the ones that started with super-Hubble size at η ∆V .As an additional confirmation the numerical results shown in section III A, this delay can be qualitatively reproduced by solving the equations of motion for a DW enclosing the FV pocket, in the socalled Nambu-Goto approximation, for some simple DW shapes (see Appendix B for details), as shown in Fig. 10.
This simple observation has two interesting consequences.First, we identify that η ∼ 2η ∆V is a 'maximal shrinking time', when most of the FV has shrunk to zero, and so it is natural to expect an additional contribution to GW production on top of the GW that the DWs have sourced during scaling.Qualitatively, this confirms the results presented in the previous section of η gw ∼ 2η ∆V .
Second, this also sets the time when we can start to picture the 'remainders' of the network in a simple and useful way: as an ensemble of FV pockets of different sizes, which are placed far apart so that we can treat them independently.One can call this a dilute gas of FV pockets.

A. False Vacuum fraction
This approximation allows to compute the FV volume fraction F fv and extrapolate it to times that are inaccessible with numerical methods.The basic idea is simple: FV pockets shrink in time (see Fig. 10), so that the lifespan of each pocket is determined by its initial size R 0 at η = η ∆V .Once the network is sufficiently fragmented we can approximate F fv (η) as a sum over an ensemble of pockets of different initial radii R 0 .Moreover, the relative weights in the sum R 0 are known because they are inherited from the scaling regime.
Indeed, since there are only 2 vacua (for a Z 2 model) that are essentially distributed randomly during scaling, the probability of finding a super-Hubble region of radius R 0 , where the field is in one vacuum, is expected to be with the correlation length at onset of annihilation, which is identified as The distribution eq.satisfies the wanted normalization: initially, a Hubble-sized region with R 0 = L ann , has 50% chance to be in either vacuum.
The FV volume fraction then can be obtained by adding up the (shrinking) volumes of all pockets with weights given by eq. ( 8) over R 0 > L ann , We assume here a flat integration measure in log R 0 for simplicity, but the results do not depend very dramatically on the measure choice.Even if this formula FIG.10: Evolution of FV pockets with spherical (solid), cylindrical (dashed) and planar (dotted) shapes as extracted from the Nambu Goto approximation (see Appendix B).Blue lines show the comoving radius R(η) ≡ R(η; R 0 ) as a function of conformal time for various initial radii R 0 , starting at η ∆V from rest.Lower R 0 pockets are more numerous.The dashed vertical gray line denotes when most of the pockets shrink to small sizes.This is expected to coincide with the peak of GW production η gw .The black solid line is the comoving Hubble length at each time.The inset shows the initial values of R 0 of (spherical or cylindrical) pockets that reach zero (purple) or the Hubble length (black) at a given time.
holds only for η > 2η ∆V , the overall normalization constant N can be fixed to have F fv = 1/2 when extrapolated to initial time η ∆V .
Generically, pockets shrink and vanish after some time.The 'trajectories' R(η; R 0 ), that we obtain by solving the dynamics in the Nambu-Goto approximation (see Appendix B) are basically triangular, they decrease to zero and remain zero.As shown in Fig. 10, the shrinking time is quite shape independent too.At any time η, one can track back the initial radius of the pocket that reaches R = 0 at that moment.We call this and show it in the inset of Fig. 10 (purple curves) 6 .By construction, then, the lower integration limit in eq. ( 9) can be replaced by R min 0 (η).Since the size distribution in eq. ( 8) is exponentially biased towards small R 0 , it is clear that F fv must be suppressed by 2 −(R min 0 /η ∆V ) 3 .The 6 Since there is a slight shape dependence, we take two values of R min 0 (η), from the spherical and the cylindrical pockets, to give a sense of 'theoretical' errors.
Nambu-Goto approximation also provides some useful information to narrow down the asymptotic behaviour at large η.Figures 10 and 13 suggest that the mock trajectory R(η; R 0 ) → R 0 − w(η − η ∆V ) with w an O(1) constant gives a reasonable approximation.With this, the asymptotic form follows.Ignoring the power law term, this allows to recognize the FV decay time introduced in eq. ( 6) as The value of 1/w can be read off from Fig. 13, which in the end results in η ann /η ∆V being around 1.3.This is in quite remarkable agreement with the numerical simulations, see Table I, representing an important nontrivial validation of the FV pocket picture.This is also manifest in Fig. 11, where we compare the fits from the numerical simulations (orange curves) to the analytic expression eq. ( 9) (blue curves).
In this picture, it is also possible to isolate the part of this FV fraction F hor fv given by Hubble sized pockets (or larger) at any time.It reduces to integrate from a higher value of R 0 , the one corresponding to pockets that enter the Hubble radius (rather than vanishing) at that time.We then introduce which is also shown in the inset of Fig. 10 (black line).Notice how this radius grows significantly than R min 0 (η), hence of the FV fraction Hubble sized (or bigger) is much smaller than the total F fv .
In this case we obtain Since R hor 0 grows linearly in η (see Fig. 10), we conclude that the FV fraction contained in super-Hubble pockets behaves asymptotically like eq. ( 6) with p = 3.
Note also that a very simple relation holds if we use the approximate mock trajectory linear in η given above, i.e.R hor 0 (η) = R min 0 (η) + η, which immediately implies that the fraction in Hubble-sized patches is: F hor fv ∼ P 0 (R hor 0 ) ∼ 2 −((1+w)η/η ∆V ) 3 with w ≈ 0.8, much more suppressed than the total fraction in FV pockets as can be appreciated in Fig. 11. while the total fraction in FV pockets is roughly F fv ∼ P 0 (R min 0 ) ∼ 2 −(w η/η ∆V ) 3 .These results are reminiscent of the scaling exp −η d , in d + 1 dimensions, suggested by [18].However, let us emphasize that the similarity is accidental.The analysis of [18] refers actually to an annihilation mechanism based on population bias, with ∆V = 0.In 2 + 1 dimensions, it was found in [21] that the analog of [18] actually did not 0.5 1.0 1.5 2.0 2.5 3.0 3.5 FIG. 11: FV fractions from numerical simulations and from analytic estimates as a function of conformal time.
Orange lines correspond to eq. ( 6) with η ann and p extracted from the fits to the simulations shown in Table I (only the first and third row, giving the lowest and highest fractions respectively).The remaining curves arise from the analytic FV pocket 'gas' approximation.Blue lines are the estimate of the full FV fraction from eq. ( 9) with w = 0.9 (solid) and 0.8 (dashed), which represents a justified margin of uncertainty (see Fig. 13).
The red curves are two estimates of the FV fraction contained in Hubble-sized (or larger) pockets, for a spherical (dashed) and a cylindrical (dotted) pocket.For α c = 1 these correspond to the intersection of the red curves with the dotted vertical line -a tiny number well outside the plot range.The FV fractions for other choices of the effective collapse criterion α c are given by the intersections of the vertical lines and the red curves.
apply for the population bias case, whereas it did work for pressure bias, with ∆V ̸ = 0. Our work extends the analysis of [21] to 3 + 1 dimensions (see also [23]) and clarifies the physical reasons behind this decay law for the pressure bias (∆V ̸ = 0) mechanism.

B. PBH formation
We are now ready to give some rough estimates of the abundance of PBHs produced during the network decay.As in [13], the natural strategy is to follow the various FV pockets as they shrink.Part of this evolution takes place while the pocket is super-Hubble sized.It is useful then to look at the 'figure of merit' defined by the ratio of the Schwarzschild radius of the pocket to the Hubble radius when the pocket size crosses the Hubble radius.This quantity actually coincides with the local overdensity of a FV with energy density ρ pocket , For α loc ≪ 1 the pocket needs to contract significantly after entering the Hubble radius, in order to form a PBH.This is less likely to happen if it is non-spherical and, since FV pockets descend from a DW network, asphericities can actually be large.For larger α loc , this is not so challenging and so one can expect PBHs to form in this way.Moreover, α loc grows in time, so some PBHs are certainly produced.Notice that the limit α loc → 1 is special: the pocket collapses to a BH as soon as it enters the cosmological Hubble radius.These BHs are actually expected to carry a baby-universe.For spherical symmetry, they have been considered in [36,37,47,48].We will consider both types of BHs, but we anticipate that baby-universe BHs should be much rarer than the ordinary ones.
The rest of the argument to estimate the PBH abundances is as follows.First, the overdensity produced by FV pockets that enter the Hubble radius at η scales like with α gw the average fraction of energy density in DWs at η gw , and the prefactor fixed by numerical simulations, see Appendix B. Second, as argued before, the collapsing network can only be acceptably approximated as an ensemble of pockets after η ≳ 2η ∆V .
Third, in principle one should do an analysis of how many of the different pockets actually manage to shrink enough to form BHs. This will depend on their degree of asphericity and angular momentum and it can be model dependent.Its estimate deserves a dedicated analysis of data from numerical simulations that is outside the scope of this work.The result of such an analysis should effectively result in a threshold of collapse, α c , such that (on average) pockets which reach α loc ≥ α c collapse into a BH.
At present we are unable to obtain a reliable estimate of α c .We thus consider a range of benchmark values that could be reasonable for this quantity.Since α c = 1 is the threshold for baby-universe BHs, α c < 1 corresponds to sub-Hubble BHs.
Fourth, we will identify the PBH abundance (of either type) with the FV fraction F hor fv evaluated at the time η PBH when α loc = α c is met.
We show in Fig. 11 the FV fraction expected to collapse into BHs according to this simplified criterion, with α c = 0.15 and 0.2.Clearly, since the collapse criterion is satisfied first with smaller α c , the sub-Hubble PBHs are much more abundant than the ones carrying baby universes.
As expected, the abundance is exponentially sensitive to α c .The abundance of baby-universe BHs (α c ≃ 1) is extremely suppressed (and given by the extrapolation of the red curves to the vertical dotted line) even for quite large α gw .On the other hand, the sub-Hubble PBHs that can be formed 'soon' could be much more abundant.The basic reason why sub-Hubble PBHs are more abundant is simply that the FV pockets they descend from are (extremely) more common.Unfortunately, with the current analysis, it is difficult to make any further quantitative statements due to the large uncertainty in α c .A more detailed study is left for the future.

V. Phenomenological Implications
We now proceed to discuss the phenomenological impact of our estimates of the spectrum of stochastic GWs and of the abundance of PBHs from the DW network.We base our results on the numerical output described in Sec.III and on the analytical understanding of the evolution of the false vacuum fraction described in Sec.IV.We start by estimating the minimal PBH abundance, formed from Hubble sized PBHs pockets that reach α loc = 1, using eq.( 11).The total abundance may be expressed as a function of the energy fraction of the network at the time of GW emission, α gw , and the background temperature at that time, T gw .The mass of these black holes is set roughly by the total energy in the Hubble volume and more precisely by eq.(B3).The abundance at a given mass is experimentally constrained by a wide variety of probes.In Fig. 12 we translate the bounds on the PBH abundance from [49] and [50] into bounds on the α gw − T gw plane (thick blue curve).We relate α loc with α gw using eq.(B4), which is a minor refinement of eq. ( 16).Constraint on α gw are obviously stronger for larger T gw , as PBHs redshift as matter for a longer time.As an interesting example of the typical mass of these PBHs, we show in green the boundaries of the asteroid mass range 10 −16 M ⊙ ≲ M PBH ≲ 10 −11 M ⊙ where the PBHs can account for the whole of the dark matter [51].
The remaining blue lines in Fig. 12 refer to tentative estimates of the abundance of PBHs, by assuming some benchmark values, α c = 0.3 and α c = 0.1, at Hubble crossing, which indicates how much the structure has to further shrink (without dissipation) to enter its Schwarzschild radius.
In the same plot, we also show the range of parameters where the GW signal from the DW network could be observed by different GW detectors, from SKA [52] at the lowest frequencies to LIGO-Virgo-KAGRA (LVK) [39], ET and CE at the highest ones.We have fixed the frequency of the GW spectrum at the peak to be dictated by the Hubble radius at T gw , and amplitude given by eq. ( 4), with efficiency ϵ = 0.6, as obtained from the numerical results, and assumed the spectrum to decrease as 1/ω for frequencies ω larger than the peak and ω 3 for smaller frequencies.This behavior corresponds to that observed in the scaling regime [26], and has been roughly confirmed recently also during the annihilation phase by [28].The GW spectra computed in our simulations roughly agree with those works, although a dedicated study is necessary to firmly establish the high frequency slope.
We then plotted the regions of parameters where the spectrum overlaps with the power law sensitivity curves derived in [53].In the same figure, we also show the current bounds obtained in [5] with LIGO-Virgo O3 data (LV), which already indirectly constraints the maximal PBH abundance for α c = 1, and the region of parameter space which could provide an interpretation of the PTA signal (red contours) [11].
Let us focus on the PTA region first, i.e. at T gw ≈ 1 − 10 GeV.We find the DW interpretation of the PTA signal to be overall compatible with constraints on PBHs for collapse thresholds α c ≳ 0.1.Interestingly, if α c ≲ 1, a fraction of PBH dark matter from DW collapse is expected, which is compatible with what is inferred from the BH merger rate measured by LIGO/Virgo.
Significant tension between PTA observations and astrophysical bounds on PBHs would instead arise only if the collapse threshold were even smaller, i.e. α c ≲ 0.1.The likelihood of such low thresholds remains to be assessed by future work.Our results differ drastically both from [17] (for example regarding the time dependence of the DW decay) and from the earlier contradictory in the first version of [54].
As mentioned above, a particularly interesting region of parameter space is the asteroid mass range, for 10 6 GeV ≲ T gw ≲ 10 9 GeV, where the totality of dark matter could be explained by PBHs from the network.This mass range is typically very hard to probe, given the particle-like size of their Schwarzschild radius.Crucially however, if PBHs originate from the DW network, a complementary GW signature of their existence is expected, that can be probed partially by LVK at design sensitivity and fully by ET and CE.
Beside these two interesting regions, the plot shows that in a large range of parameters, the annihilation of the network can be "heard" by different GW observatories, and that a non-negligible abundance of PBHs might also be expected if α gw > O(10 −2 ).
The interplay of GW and PBHs signatures described in this work is similar to the more studied scenario of PBH formation from the collapse of large adiabatic perturbations from inflation.
However, GWs from density perturbations as an interpretation of PTA data have been shown to be in tension with constraints from PBH overproduction [55,56].In contrast, the mechanism presented here remains viable according to our current understanding.Additionally, asteroid-mass PBH DM from inflationary perturbations has been argued to give a GW signal peaked in the LISA frequency band [57], whereas ground-based interferometers (LVK, ET, CE) are best-suited to indirectly probe PBH DM from DW collapse.The shaded regions are the constraints on the GW spectrum from LIGO-VIRGO O3 data (LV) and the prospects for detection with LIGO-Virgo-KAGRA design sensitivity (LVK), Einstein Telescope, Cosmic Explorer, LISA and SKA using the power-law integrated sensitivity curves from [53].Finally, the green band corresponds the asteroid mass range (10 −16 − 10 −11 M ⊙ ) [58,59] and the black dashed lines correspond to PBHs between 1 and 100 solar masses, for α c = 1.For the other values of α c the band moves slightly to the left.
The thin wall approximation provides additional insight on the DW network behaviour, especially in the annihilation phase where the remains of the network is a collection of separate FV pockets.This is a good approximation when the DW worldsheet curvature radius, R, is large compared to its width, which is set by the inverse scalar mass.(This is satisfied in most of the network during scaling, except for a small fraction of the total volume where collisions, interconnections or pinch-off events occur).It is possible to include the gravitational effect from the DWs themselves (see e.g.[36,37,60,61]) but we shall ignore this here.
In this approximation, the evolution of a FV pocket follows from the 'equation of motion' for the DW at its boundary, which reduces to the Nambu-Goto (NG) equation Here K is the extrinsic curvature of the DW worldsheet and we included a pressure term given by the bias ∆V , see e.g.[62].It is not easy to solve this equation in general.However, the equation simplifies for walls with higher symmetry.We are interested in the motion of a DW network where inevitably the DW shapes are random.However, once the network annihilation starts, the DW motion, in a way, simplifies.Indeed, the DWs are simply the boundaries of FV pockets, which shrink quite quickly and quite independently of their initial shape.
This can be illustrated by comparing 3 extreme cases that can be easily computed, where the shape of the DW is: i) spherical, ii) cylindrical and iii) planar.The NG equation then reduces to with n = 2, 1 for spherical or cylindrical DW of comoving radius R respectively.The case n = 0 corresponds to a planar wall placed at, say, z = R(η).Primes denote derivatives w.r.t.conformal time, and It is straightforward to integrate this equation and thus follow the evolution of a structure of certain initial comoving radius R 0 .Some representative examples are shown in Fig. 10.It is clear from the figure that the structures reach arbitrarily small size after a finite (conformal) time ∆η, that is the time lapse until R approaches 0. Of course eventually a small enough structure transfers its energy into scalar waves (which are not captured in the NG approximation).Note that if one prefers to define the collapse time as when R(η) reaches a small radius r c , the result would be the same so long as r c ≪ R 0 .
The trajectories shown in Figs. 10 are readily understood: closed DWs shrink under both the effect of the tension and the pressure ∆V , reaching relativistic speeds quite quickly.
Interestingly, the collapse time ∆η depends mostly on the initial size.As shown in Fig. 13, for large enough initial FV pockets, the collapse time approaches C R 0 , with a constant C in the range 1.15 − 1.2, quite independently of the pocket shape.
Of course, in a network the DW shapes are not symmetric.However, Fig. 13 signals a quite clear time scale in the network decay: the first stage of annihilation, (that is still far from the FV pocket gas picture) takes about one Hubble time.
One expects a 'burst' of GWs from the collapsing FV regions as they shrink because they are significantly nonspherical.This GW production time η gw is then expected to be near the collapse time because it is dominated by the most numerous pockets, with sizes of order η ∆V .In summary, the fact that Hubble-sized structures have to reach small sizes naturally leads to the expectation It is also easy to keep track of the energy of a given FV pocket, and how it evolves in time.Focusing on spherical symmetry for simplicity, the energy is ∆V R 3 (η)a 3 (η) + 4πσR 2 (η)a 2 (η)γ(η) , (B3) and because of the expansion it is not conserved.For super-Hubble pockets, first E grows (because both the FV region and the DW gain volume/area by the expansion).Only when they enter the Hubble radius the energy stabilizes.
We can keep track of the energy carried by the FV pocket that enters the Hubble radius at each time, by evaluating eq.(B3).The γ factor grows in time, but not enough to compare to the volume contribution.Thus, after about one Hubble time, E scales like the physical Hubble volume, ∼ η 6 .Solving numerically the NG equation (B1) we find that the expression E/E 0 ≃ (τ 6 + τ 5 + 2τ 4 )/4, with τ = η/η ∆V and E 0 the initial energy of the pocket, provides a good fit (within 5%) of the actual time dependence (the τ 5 term can be identified as the DW gamma factor).This leads to the following improvement of eq. ( 16) where the factor of 1.5 comes from rewriting α loc (η gw ) = ∆V /(3H(η gw ) 2 M 2 p ) in terms of the fraction of the energy density in the DW network at the same time α dw = ρ dw /ρ c , with ρ dw (η gw ) = 2.6 σH(η gw ) and τ gw ≃ 2. This is smaller than eq.( 16) for τ > τ gw , so eq.( 16) provides a conservative bound.

FIG. 1 :
FIG.1: Timeline of DW network evolution.After spending some time in the scaling regime, the bias (vacuum energy difference ∆V ) becomes effective at (conformal) time η ∆V .GW and PBH production occur somewhat later, at η gw and η pbh respectively.The decay of the False Vacuum volume fraction, F fv , is parameterized by a decay time η ann , slightly larger than η ∆V , and an exponent p.Both our numerical simulations and analytical model point to p = 3.

2 FIG. 2 :
FIG.2: Evolution of the average of the field in the simulation box, φ, plus (minus) its standard deviation σ ϕ , for several sizes of the bias ∆V , as a function of conformal time.

FIG. 3 :FIG. 4 :
FIG. 3: (Log of) Volume fraction in the simulation box in the false vacuum, as a function of conformal time.The blue curves show the numerical results from our simulations.We fitted such results for small η, before the curves cross the black dot-dashed curve, that shows the inverse of the number of Hubble volumes in the simulation box.The orange curves are the resulting fits.

FIG. 5 :
FIG.5: Area fraction, eq.(3), as a function of conformal time.Note that at early times, η ≲ 10, there is a transient which does not carry physically relevant information.

7 FIG. 8 :
FIG.8: Energy density fraction in GWs compared to quadrupole formula estimate from the total energy density available in the scalar field.Note that the initial Ω gw is very small because we start the computation of GWs only at η = 35 for computational reasons.

FIG. 9 :
FIG.9: Energy density fraction in GWs compared to quadrupole formula estimate from the total energy density available in the scalar field.

FIG. 12 :
FIG. 12: Constraints on the PBH abundance from the collapse of the DW network for different values of the collapse threshold α c (blue lines) in terms of the fraction of the energy density in the DW network at the time of GW emission (α gw ) and the background temperature at that point (T gw ).The region above the solid blue curve is conservatively excluded.The red ellipsis shows the region of parameters where the DW network interpretation explains the recent GW signals detected at PTA[11] at one and two standard deviations.The shaded regions are the constraints on the GW spectrum from LIGO-VIRGO O3 data (LV) and the prospects for detection with LIGO-Virgo-KAGRA design sensitivity (LVK), Einstein Telescope, Cosmic Explorer, LISA and SKA using the power-law integrated sensitivity curves from[53].Finally, the green band corresponds the asteroid mass range (10 −16 − 10 −11 M ⊙ )[58,59] and the black dashed lines correspond to PBHs between 1 and 100 solar masses, for α c = 1.For the other values of α c the band moves slightly to the left.

2 FIG. 13 :
FIG.13: Ratio of the initial radius R 0 and the (conformal) time to reach R(η) = 0, ∆η, of super-Hubble FV pockets.It depends mildly on the shape and it asymptotes to a constant value (dotted black curve).

TABLE I :
Mean and standard deviation of the parameters of the fitting function eq.(6), over several realizations of our lattice simulations (about 3 − 4 per each row), with different random seeds.