Primordial black hole dark matter from catastrogenesis with unstable pseudo-Goldstone bosons

We propose a new scenario for the formation of asteroid-mass primordial black holes (PBHs). Our mechanism is based on the annihilation of the string-wall network associated with the breaking of a $U(1)$ global symmetry into a discrete $Z_N$ symmetry. If the potential has multiple local minima ($N>1$) the network is stable, and the annihilation is guaranteed by a bias among the different vacua. The collapse of the string-wall network is accompanied by catastrogenesis, a large production of pseudo-Goldstone bosons (pGBs) -- e.g. axions, ALPs, or majorons -- gravitational waves, and PBHs. If pGBs rapidly decay into products that thermalize, as predicted e.g. in the high-quality QCD axion and heavy majoron models, they do not contribute to the dark matter population, but we show that PBHs can constitute 100\% of the dark matter. The gravitational wave background produced by catastrogenesis with heavy unstable axions, ALPs, or majorons could be visible in future interferometers.


Introduction
The only known interaction that dark matter (DM) has with ordinary matter is gravitational [1,2]. This is why black holes can be good DM candidates if they are formed in the early Universe, i.e. are primordial black holes (PBHs) [3][4][5][6][7][8]. Moreover, they need to be stable, so that they do not evaporate through Hawking's radiation before the present time [9][10][11][12][13][14], and they need to avoid current constraints from microlensing [15][16][17], as well as many other bounds collected in e.g. Ref. [18]. There is a window of PBH masses, between 10 −16 and 10 −10 M (roughly equivalent to asteroid masses), for which PBHs can constitute 100% of DM. A necessarily incomplete list of proposed formation scenarios for PBHs in the asteroid-mass range includes density perturbations in the early Universe (see e.g. [6,[19][20][21]), bubble collisions [22,23], the collapse of cosmic strings [24], scalar field dynamics [25][26][27][28], long-range interactions [29], and collapse of domain walls [30,31] or vacuum bubbles [32,33] in multi-field inflationary scenarios. We propose here a novel production scenario for asteroid-mass PBHs based on the mechanism we have dubbed catastrogenesis in our previous work [34]. Our mechanism does not depend on inflationary physics, nor on primordial density perturbations, and crucially it can be linked to other unsolved problems in particle physics, e.g. the strong CP problem.
We assume the existence of a global U (1) symmetry whose spontaneous breaking is associated with the existence of Goldstone bosons. A generic feature of these models is that the global symmetry is also broken explicitly, so that the bosons are massive, i.e. pseudo-Goldstone bosons (pGBs). This happens in a plethora of extensions of the Standard Model of particle physics, including the original axion model [35][36][37], invisible axion (also called QCD axion) models [38][39][40], generic axion-like particle (ALP) models [41], and (singlet) majoron models [42][43][44][45][46][47][48]. Intriguingly, some models predict a larger-than-expected mass for the QCD axion [49][50][51][52], including the so-called "high-quality QCD axion" [53] and previous models. (For a review on earlier attempts of model building predicting heavy QCD axions, see Section 6.7 of Ref. [54].) Heavy majorons, even of mass in the TeV range, have been considered as well (see e.g. [44,47]). If the pGBs are unstable and decay in the early Universe, they could not constitute the DM, and thus such models might lack a DM candidate. This can be the case for axions and ALPs, as well as for majorons which could get a mass from soft breaking terms or from gravitational effects (see e.g. [43][44][45][46][47]).
We deal with bosons that are heavier that a GeV and decay very early into Standard Model products which rapidly thermalize. Axions would decay into photons and charged fermions. Majorons could decay mostly into two relatively heavy right-handed neutrinos, which in turn decay very fast into pions, charged fermion pairs, and active neutrinos. Therefore, such pGBs cannot be the DM in the models we consider, but catastrogenesis [34] may lead to a different DM candidate, namely PBHs. We consider explicit breaking potentials which admit a number N > 1 of minima along the orbit of vacua. In this case, if inflation happens before the spontaneous symmetry breaking, as we assume, a cosmological problem may arise: a stable string-wall network forms and could eventually come to dominate the Universe leading to an unacceptable cosmology [79,80]. Both axion [81] and majoron [45] models can result in N > 1.
The cosmological problem is solved by adding an additional explicit breaking in the Lagrangian [81][82][83], possibly due to the effect of Planck suppressed operators [34,84,85], which produces an energy difference (a bias) between the minima that leads to a unique vacuum. As first proposed in 1974 in Ref. [79], this small bias results in the annihilation of the string-wall network, and only one vacuum is eventually left. The annihilation is accompanied by catastrogenesis [34], i.e. the production of gravitational waves (GWs), pGBs, and PBHs. Our main result is that for a range of the string-wall network annihilation temperatures, one can form asteroid-mass PBHs and potentially a GW signal in future observatories. As we will see, both the abundance of PBHs and the amplitude of the stochastic GW background produced depend crucially on the details of the annihilation.
Although we focus on models featuring pGBs (that we will call, as customary in the recent literature, ALPs), our mechanism also applies in principle to the breaking of any discrete symmetry with multiple vacua, when the wall system annihilation produces unstable particles, such as flavons [86], whose decay products thermalize.
The paper is structured as follows. In Sec. 2 we introduce the generic particle model we assume and describe its cosmology. In Sec. 3 we present our estimates for the PBH relic abundance and find the range of parameters in which they could account for the whole of the DM. In Sec. 4 we show that for the same range of parameters a signal in GWs could be observable in future detectors. Sec. 5 explains in detail restrictions imposed by consistency conditions of our model. Finally the last section contains some concluding remarks.

ALP models and their cosmology
Here we briefly review well-known aspects of the models we study. A more detailed account can be found e.g. in Ref. [34] (see also [87]). We consider a complex scalar field φ = |φ|e iθ whose phase after spontaneous symmetry breaking at an energy scale V is related to the ALP field a as a = θ V . The generic potential V (φ) for the scalar field includes three terms, one with a spontaneously broken global U (1) symmetry and two others which break the symmetry explicitly, The first term is U (1) invariant and leads to the spontaneous breaking of this symmetry during a phase transition in the early Universe, at which the magnitude of the field becomes |φ| V and cosmic strings are produced due to the random distribution of the phase θ in different correlation volumes. Assuming that the bosons have the same average energy of the visible sector (which is possible in many inflation models), this phase transition happens at a temperature T V . Upper bounds on the energy scale of inflation impose V 10 16 GeV [88,89]. A short time after its formation, the string system enters into a "scaling" regime, in which the population of strings in a Hubble volume tends to remain of O(1) (see e.g. Ref. [80] and references therein).
The second term in Eq. (2.1) explicitly breaks the U (1) symmetry of the first term into a Z N discrete subgroup, producing N degenerate vacua along the previous U (1) orbit of minima |φ| V . Close to each minimum the ALP field a has a mass We assume for simplicity that temperature corrections to m a are negligible. We do not expect that these corrections will be important in our scenario since, as we explain below, catastrogenesis is independent of the temperature dependence as it happens when the ALP mass has reached its asymptotic present value. The equation of motion of the field a in the expanding Universe is that of a harmonic oscillator with damping term 3Hȧ, where H = (2t) −1 is the Hubble expansion rate during the radiation dominated epoch. While 3H m a the field does not change with time. As H decreases with time, at a temperature T w when H(T w ) m a /3, regions of the Universe with different values of θ evolve to different minima separated by domain walls of mass per unit surface (or surface tension) Here f σ is a model dependent dimensionless parameter. For N = 2 the model is solvable analytically and f σ 6. Whenever choosing specific values of N and f σ will be required, we will assume N = 6 and f σ /N 1. The temperature T w at which walls appear, i.e. H(T w ) m a /3, is where we have used that for a radiation dominated Universe the Hubble expansion rate is Here, M P = 1.22 × 10 19 GeV is the Planck mass and g is the energy density number of degrees of freedom [90]. A short time after walls form, when friction of the walls with the surrounding medium is negligible, the string-wall system enters into another scaling regime in which the linear size of the walls is the cosmic horizon size t. Therefore, its energy density at time t is Friction on walls, due to scattering with particles in the surrounding medium, is negligible when the energy of these particles is much larger than the ALP mass m a [91,92], which is the case for most if not all the models we consider. This has been confirmed e.g. in Ref. [93].
There the impact of friction on the parameter space of the high-quality QCD axion was studied (see their appendix D), and it was concluded that even when friction was important their results were only mildly different from those obtained neglecting friction. The third term in Eq. (2.1), proposed in Ref. [81], is assumed to be much smaller than the second one, i.e. b 1. It breaks explicitly the Z N symmetry, raising all vacua with respect to the lowest energy one by a bias of order The motion of the walls after they are formed is determined by two quantities. The surface tension σ tends to straighten out curved walls to the horizon scale as it produces a pressure p T ρ wall , while the bias produces a volume pressure p V V bias [83]. This latter pressure tends to accelerate the walls towards their higher energy adjacent vacuum, converting the higher energy vacuum into the lower energy one, which releases energy that fuels the wall motion.
Assuming that p V p T when walls form (i.e. b 1), as p T decreases with time it becomes similar to and then smaller than p V . At this point the bias drives the walls (and strings bounding them) to annihilation. Taking p T p V , i.e. ρ wall σ/t ann V bias b v 4 , as the condition for wall annihilation, which defines the temperature T ann at which the string-wall system annihilates, At this point the energy stored in the string-wall system goes almost entirely into nonrelativistic or mildly relativistic ALPs (since the wall thickness is m −1 a ) [94], into GWs, and into PBHs.
The GWs produced at annihilation (see Appendix A and e.g. Ref. [34] for more details) have a characteristic spectrum peaked at and peak energy density (see e.g. [95][96][97][98]) Following Ref. [98], we include in Eq. (2.11) a dimensionless factor GW 10 -20 for N = 6 (see Fig. 8 of Ref. [98]) found in numerical simulations that parameterizes the efficiency of GW production. In the following we conservatively assume GW = 10. Eq. (2.11) is also (see Appendix A) the usually quoted maximum of the GW energy density spectrum at present as a function of the wave-number at present k, which for the present scale factor R 0 = 1 coincides with the comoving wave-number, or of the frequency f = k/(2π). This is defined as 86,99]. GWs are also emitted by the string network before walls appear. The dominant source of these waves are loops continuously formed by string fragmentation. From the spectra computed in Refs. [100][101][102], also based on the quadrupole formula, a simple good fit can be derived for the energy spectrum of GWs emitted by strings during the radiation dominated era, namely (2.12) This spectrum has a low frequency cutoff f st cut at the frequency of GWs emitted the latest possible by strings, when the horizon size is largest, namely when walls are formed, which is given by Eq. (2.10) using T w instead of T ann , i.e.
which can be written in terms of the ALP mass as (2.14) Given that we consider only ALPs with mass m a 1 GeV, the emission of GWs by strings before walls appear only contributes at high frequencies, f 82 Hz. Comparing Eqs. (2.12) and (2.11) we are going to show that only in a restricted parameter space GWs from strings could be observable by the future GW Einstein Telescope [103].
In the next section we will argue that when ALPs are unstable and decay fast in the early Universe, the string-wall system annihilation could produce PBHs that account for the whole (or part) of the dark matter, while also producing GWs with peak frequency from 10 −5 to 10 2 Hz and density Ω GW h 2 > 10 −15 , which corresponds to an observable signal in future GW detectors.

Primordial black holes
In N = 1 models the walls appear in the form of ribbons flanked by two strings which become narrower due to surface tension and disappear very fast. By contrast, in the latest stages of wall annihilation in N > 1 models, closed walls are expected to arise and collapse in an approximately spherically symmetric way. In this case, some fraction of the closed walls could shrink to their Schwarzschild radius R Sch (t) = 2GM (t) = 2M (t)/M 2 P and collapse into PBHs [104]. Here M (t) is the mass within the collapsing closed wall at time t. Considering that during the scaling regime the typical linear size of the walls is the horizon size, which we take to be t, PBH formation could happen if the ratio is close to one, i.e. p(t) 1. Simulations of the annihilation process [105] can then be used to estimate at which temperature T T ann this could happen. As we mentioned above, the annihilation process starts at T ann when the contribution of the volume energy density to the mass within a closed wall of radius t becomes as important as the contribution of the wall energy density. Shortly after, the volume density term dominates over the surface term, and the volume pressure accelerates the walls towards each other. Close to the start of the annihilation process, the mass within a closed wall as a function of the lifetime t of the Universe is Therefore, the ratio p(t) increases with time (as t 2 once the volume term becomes dominant).
At T > T ann the surface term dominates and at T < T ann the volume term dominates. If at the moment t ann when annihilation starts p(t ann ) is close to 1, PBHs would form immediately.
If p(t ann ) 1 instead, PBHs could only form later, at a time t > t ann corresponding to a temperature T < T ann , for which p(t ) = 1. Temperature and time are related by H = 1/2t (see Eq. (2.5)), assuming radiation domination. If p(t ann ) 1, at T only a small fraction of the original string-wall system still remains.
The wall surface energy density σ/t decreases with time with respect to the volume energy density due to the bias V bias which is constant in time. Although the volume energy density is negligible initially, both contributions become similar at T ann , i.e. V bias σ/t ann , thus Eq. (3.2) implies As expected of all quantities depending only on the annihilation process, M (t ann ) and p(T ann ) only depend on the parameters V bias and σ. After t ann , the volume contribution to the density rapidly becomes dominant over the surface contribution, and When t 3 t ann we can we neglect the second term in Eqs. (3.5) and (3.6), consequently We are here assuming that the characteristic linear dimension of the walls continues to be close to t even after annihilation is underway. At some point along this process one expects larger deviations from this scaling behavior, but detailed simulations that are not available at the moment would be needed to assess when and how this departure happens. Using Eq. (3.7) we find the PBH formation temperature T as the temperature at which p(T ) = 1, which in terms of T ann is This relation defines T and its corresponding time t = 1/2H(T ). The PBH mass is then M (t ), (3.9) Eqs. (3.8) and (3.9) show that the temperature of PBH formation T depends only on V bias (or, equivalently, M PBH ) (3.10) From Eq. (3.9) and Eq. (3.1) one finds (3.11) Small deviations from a spherically symmetric collapse, as well as effects of some angular momentum, could make the probability of forming a PBH at temperature T somewhat smaller than p(T ), which we could attempt to parameterize as p(T ) β with a real positive exponent β. A modification of the probability of this type does not have any effect on our estimates of PBH formation, since when p = 1 also p β = 1. A large enough deviation from the spherical shape could prevent the formation of a PBH, since the degree of asymmetry may decrease initially during the contraction but increase in the late stages of the collapse [106]. The details of the collapse process become more important if p(T ann ) 1, since a longer evolution is needed to potentially reach p = 1 (Ref. [104] suggests that in this case asphericities, energy losses or angular momentum might make more difficult the formation of a black hole). However, a large deviation from sphericity is unlikely, as shown in the context of the collapse of vacuum bubbles produced during inflation [32]. Thus, at least for some portion of the walls, the departure from spherical symmetry during collapse is likely to be small. In this case, the PBH density at formation is given by the energy density in the wall system when PBHs form times the probability of PBH formation, namely ρ PBH (T ) p β (T )ρ wall (T ), (3.12) and the fraction f PBH of the total DM density ρ DM in PBHs is where we used that at the moment of PBH formation temperature p(T ) = 1. Notice that f PBH is the same at present as it was at T (since the PBH and DM densities redshift in the same manner).
Simulations of the process of the string-wall system annihilation for N > 1 models have been performed [105], which provided measurements of the times at which the area density of walls (area per unit volume A/V ) is 10% and 1% of what it would have been without a bias. Notice that this ratio is the same in comoving (as given in Ref. [105]) or physical coordinates. We call these times t(10%) and t(1%), and the corresponding temperatures T (10%) and T (1%). In Ref. [105], the pressure due to the walls is parameterized as where A(t) is a function which accounts for deviations from scaling, and find that A(t) is close to 1. Assuming that the wall contribution to the energy density is dominant, a simple approximation as a power law for the evolution of the wall energy density with temperature after annihilation starts, seems to work well to extract the real positive exponent α from the mentioned simulations [105]. Table VI and Fig. 4 of Ref. [105] show that the t(1%)/t(10%) = [T (10%)/T (1%)] 2 ratio takes up central values close to 2, actually from 1.7 to 1.5, under different assumptions. If the energy of the string-wall system is still dominated by the contribution of the walls until t(1%), i.e. the energy density is proportional to A/V , ρ wall = Aσ/V , then (3.15) and the central values measured of the t(1%)/t(10%) ratio translate into values of the exponent α roughly between 9 and 14. In our figures we include also the value 7 used in Ref. [104]. Taking into account the systematic errors quoted in Table VI of Ref. [105], α could be as large as 21. However, the volume contribution to the energy density of the string-wall system may not be negligible, which introduces a further uncertainty in the determination of α. To get an estimate of this uncertainty, we can proceed assuming that the volume energy is dominant in the string-wall system, i.e. that its density is proportional to A 3/2 . In fact, since the simulation volume in Ref. [105] is the same in both the cases with and without a bias (in the latter case the evolution follows the scaling solution), the ratio of the area densities can be identified with the ratio of the area A in the evolution with bias and the characteristic area t 2 in the scaling case within a Hubble volume. Thus [A(t(10%))] 1/2 √ 0.10 t(10%), and similarly [A(t(1%))] 1/2 √ 0.01 t(1%). Moreover, in a Hubble volume, if the V bias contribution dominates over the surface energy, the wall-system density is ρ wall V bias [A(t)] 3/2 /t 3 , and T (10%) T (1%) α ρ wall (T (10%)) ρ wall (T (1%)) 10 3/2 . (3.16) This leads to values of α larger by a factor 3/2 with respect to our previous estimates minus two (i.e. from (3/2)7 to (3/2)19). Thus, α could be as large as 28. Using Eq. (3.14), Eq. (3.13) becomes In this equation ρ DM (T ) is easily related by redshift to the present dark matter density, and ρ wall (T ann ) can be related to the present radiation density, in the following manner.
Had the string-wall system not annihilated, its energy density would have continued evolving in the scaling regime, with ρ wall σ/t, until the time t wd at which ρ wall (t wd ) = ρ rad (t wd ), thus ρ wall (t ann ) = t wd t ann ρ rad (t wd ) = H(T ann ) H(T wd ) ρ rad (t wd ). (3.20) Considering the redshift of the radiation density to the present, Eq. (3.17) becomes or, using Eq. (3.19), In the following we neglect the possible change of energy degrees of freedom between T ann and T , since these temperatures are always close to each other. Fig. 1 shows upper limits on the fraction of DM in PBH f PBH as a function of the PBH mass M PBH . The constraints correspond to Hawking's evaporation bounds [9,12,13,[107][108][109][110], microlensing bounds [15,16,33,111,112], and CMB bounds [113]. All bounds on the PBH abundance are taken from Ref. [18], except for the ones from the Subaru Hyper Suprime-Cam (HSC) data [16,33]. Additional constraints, independent of cosmology, relies on dwarf galaxy heating [114,115]. The parameter space will be further probed by upcoming experiments [116,117]. The lines of constant annihilation temperature T ann implied by Eq. The range of T ann for which PBHs can constitute 100% of the DM is better seen as the orange band in Fig. 2. The yellow band in the same figure shows the range of T ann in which, according to Eq. (3.23) f PBH = 10 −2 , a DM fraction that is just below all the observational bounds for 10 −10 M < M PBH < 1 M . Such PBHs could potentially account for some putative events, as found by HSC [16] and other [15] microlensing observation, as well as LIGO observations [118].
As we explain in Sec. 5, the consistency requirement of having walls formed well after strings appear, i.e. T w V , implies that a lower limit on m a translates into a region of this plot being excluded. Assuming T w < 10 −2 V , we obtain the limits shown in Fig. 2 as red lines. Each lower limit on the ALP mass (from 1 GeV to 1 TeV) excludes the region above and to Two combined consistency conditions, T ann < 10 −2 T w and T w < 10 −2 V , reject the regions above the thick blue line in Fig. 2. Finally, requiring the string-wall system not to dominate the energy density of the Universe, i.e. T wd < T ann , rejects the region below the dashed green line, thus it does not affect any of our regions of interest.
Notice that the some of the above mentioned consistency conditions depend on T w , and thus on m a . Therefore, they might be affected by a temperature dependence of m a which we have neglected. However, the results depending only on the string-wall system annihilation, i.e. the present density of PBHs and GWs due to catastrogenesis, are independent of this assumption.

Potentially observable GWs
Since the GWs due to the string-wall system and PBHs are produced at annihilation, their present abundance depends only on two parameters, V bias and σ, which determine when annihilation happens. It is convenient here to choose the two independent parameters to be instead the GW peak frequency f peak (which depends only on T ann through Eq. (2.10)) and  Fig. 1). Each lower limit on m a excludes the region above and to the right of the corresponding red line, imposing for consistency of the model that T w < 10 −2 V (i.e. m a can only be larger than the quoted limit to the left of the lines, and can be both smaller on both sides). Other consistency conditions, T ann < 10 −2 T w < 10 −4 V and T wd < T ann , reject the regions above the thick blue line and below the dashed green line, respectively. the fraction of DM in PBHs, f PBH , and write Ω GW as . (4.1) Fig. 3 shows the expected present GW density produced by the string-wall system as a function of f peak (lower abscissa axis) or T ann (upper abscissa axis). The region where PBHs can be all of the DM (f PBH = 1) and the one where they can constitute only a DM subcomponent (f PBH = 10 −2 ) are shown respectively in orange and yellow. We consider a range of α values, 7 to 28 (black solid lines). The predictions are compared to the current upper limits (solid contours) and the expected reach (dashed contours) of several GW detectors. We include the projected sensitivities of the space-based experiments TianQin [119], Taiji [120], and the Laser Interferometer Space Antenna (LISA) [121] in green, the reach of Figure 3: Present GW density Ω GW h 2 emitted by the string-wall system predominantly at annihilation, as a function of f peak (lower abscissa axis) or T ann (upper abscissa axis) for f PBH = 1 (orange region) or f PBH = 10 −2 (yellow region) expected for different values of the power index α, between 7 and 28 (solid black lines). We also show upper limits (solid contours) and the expected reach (dashed contours) of existing and future GW detectors, respectively. The quoted lower limits on m a allow only the regions above and to the right of the respective slanted red lines (meaning that m a can be larger than the quoted values only to the right side of the lines, and can be smaller on both sides) if T w /V = 10 −2 (for smaller values of this ratio the allowed regions shrink). Gray dotted lines indicate constant PBH mass values. The consistency condition T ann < 10 −2 T w < 10 −4 V rejects the region below and to the right of the thick blue line.
the Atom Interferometer Observatory and Network (AION) [122], the Atomic Experiment for Dark Matter and Gravity Exploration in Space (AEDGE) [123], the Deci-hertz Interferometer Gravitational wave Observatory (DECIGO) [124], and the Big Bang Observer (BBO) [125] in blue. Finally, we show in red the reach of the ground based experiments Einstein Telescope (ET) (projection) [103] and in grey limits and future reach of the Laser Interferometer Gravitational-Wave Observatory (LIGO) [126]. The cyan band corresponds to the 95% C.L. upper limit on the effective number of degrees of freedom during CMB emission from Planck, and other data [127], which imposes Ω GW h 2 < 10 −6 .
Also shown in Fig. 3  (4.2) Figure 4: Two examples of combined GW spectra from wall annihilation and from string emission: for a model with α = 28 the spectrum from wall annihilation (thick red solid lines) is observable in several future GW detectors, and the corresponding spectrum from strings (thick dark red lines) for the m a values indicated would be observable by ET only for m a 1 GeV; for a more typical example of a model with lower α value, α = 12, the GWs from wall annihilation (in violet) are observable, but the corresponding spectrum from strings (in blue) is not observable.
As explained in the next section, the lower limits quoted on the ALP mass (between 1 GeV and 1 TeV) allow only the regions above and to the right of the respective slanted red lines (meaning that m a cannot be larger than the quoted limit to the left) if T w /V = 10 −2 .
For smaller values of the T w /V ratio the allowed regions shrink, i.e. the slanted red lines move to the right. Fig. 3 clearly demonstrates that a GW signal due to the annihilation of the cosmic walls is within the reach of several experiments for values of α above about 12, both if PBHs constitute the whole of the DM or just a small fraction of it.
Contrary to the GW signal from wall annihilation, the GW signal emitted by strings before walls appear depends on the ALP mass value. Therefore, we do not show the signal from strings in Fig. 3. Eq. (2.14) shows that for m a > 1 GeV strings contribute to the GW predicted spectrum only above f st cut = 82 Hz, and the cutoff frequency moves rapidly to higher frequencies as m 1/2 a , moving from LIGO frequency range to the range which will be probed by ET. However, Eq. (2.12) shows that the GW energy density due to cosmic strings is always below 10 −9 , since we only consider V < 10 16 GeV, thus putting the signal out of the reach of LIGO. The signal is potentially within the reach of ET only for models with very large values of α. Fig. 4 shows two examples of GW spectra from wall annihilation and from string emission. One of them with α = 28 has a large peak amplitude from wall annihilation, with a spectrum ∼ f 3 below the peak and ∼ 1/f above it (thick red lines), that would be largely observable in several future GW detectors. The corresponding spectrum from strings, shown in dark red for several values of m a , would be observable by ET only for m a 1 GeV. The other spectra correspond to more typical examples with lower value of α, α = 12. In this case the GWs from wall annihilation (shown in violet) could still be observed by DECIGO and BBO, but the GWs from strings (in blue, each for the indicated m a value) could not.
Similar predictions for GWs from wall annihilation were found in Ref. [93] for the highquality QCD axion model. Because our consistency conditions are tighter and we include the production of asteroid-mass PBHs, we find a less prominent GW signal. Assuming for example α = 7, we find that above 1 Hz the GW amplitude cannot be larger than Ω GW h 2 10 −19 .

Self consistency conditions
Consistency of our models requires that the annihilation of the string-wall system happens before walls would become dominant in the Universe and well after the formation of domain walls, which in turn appear sufficiently after the formation of the cosmic string network, i.e. that 10 16 GeV > V T w T ann > T wd . Conservatively, we will require T w /V 10 −2 and T ann /T w 10 −2 as well. Our model has three independent parameters, which at the level of the Lagrangian are V , v V , and b 1. Instead of v, we could alternatively use m a through Eq. (2.2), with m a V . But for our problem it is convenient to choose as the independent parameters T ann and M PBH (on which anything dependent only on the annihilation depends), as well as m a . In terms of these parameters, (5.7) Thus having b small enough is the same as having T ann T w . Extracting m a from the ratio T ann /T w one gets, Thus, requiring T w /V < 10 −2 and T ann /T w < 10 −2 , this last expression imposes an M PBH dependent upper limit on T ann shown as the blue thick line in Fig. 2 and Fig. 3.
Finally, let us ensure that wall-domination never happens before annihilation. Had the walls not annihilated previously, the walls would dominate over radiation when σ/t wd ρ(T wd ) or (5.10) Thus, H(T wd ) < H(T ann ) 0.33 (T ann /T w ) 2 m a means that This condition is automatically fulfilled for T ann /T w = 10 −2 since we require V < 10 16 GeV.
For smaller values of this ratio we check now that Eq.

Concluding remarks
Many extensions of the Standard Model of particle physics feature U (1) symmetries whose breaking is associated with the existence of heavy pGBs. There are QCD axion, ALP, and majoron models of this type. If the pGBs are unstable and decay in the early Universe, they cannot constitute the DM, and thus such models might lack a particle DM candidate.
Potentials which admit a number N > 1 of minima along the orbit of vacua imply the formation of a string-wall network whose annihilation is accompanied by catastrogenesis [34], i.e. the production of GWs, pGBs, and PBHs. We find that for annihilation temperatures 10 4 GeV T ann 10 9 GeV, one can form asteroid-mass PBHs (10 −16 M M PBH 10 −10 M ) that can constitute the entirety of the dark matter. For slightly smaller annihilation temperatures, 1 GeV T ann 10 5 GeV, the PBHs produced have mass 10 −10 M M PBH 10 M , and could potentially account for HSC [16] and other [15] microlensing observations, and they can account for part of the mergers observed by LIGO [118]. The resulting PBH masses depend on α, an exponent that parameterizes the annihilation rate of the string-wall network and needs to be determined from simulations of this process.
The GW signal produced by catastrogenesis together with PBHs can have an amplitude as large as Ω GW h 2 10 −7 and can be observable in future GW detectors. For a particularly favorable set of parameters, the GWs emitted from strings alone, before walls appear, could be observable as well.
Notice that although in this paper we do not study the breaking of purely discrete symmetries, in principle our mechanism also applies to this case. The main ingredient is the existence of multiple vacua. Therefore, the breaking of a discrete symmetry can also result in a scenario similar to the one described here, as far as the wall system annihilates mostly into unstable particles, such as flavons [86], whose decay products thermalize.
There are several caveats concerning the interpretation of our results. The production of PBHs relies on the formation of closed domain walls which must retain a high degree of sphericity during the last stages of their collapse. The amplitude of the stochastic GW background produced by catastrogenesis is delicately sensitive to the precise evolution of the string-wall network at the very end of its existence. Clarifying both these issues requires dedicated simulations of the process of annihilation of the network, which we hope will be tackled in the future. The collapse of string-wall networks associated with heavy ALPs is an open problem worth examining in further detail. quadrupole moment of the walls is Q ij E w t 2 and ... Q ij σt. Therefore P Gσ 2 t 2 . The energy density ∆ρ GW emitted by the string-wall system in a time interval ∆t is then ∆ρ GW (t) P ∆t/t 3 Gσ 2 ∆t/t. Thus, in a time interval equal to the Hubble time ∆t t, ∆ρ GW (t) Gσ 2 independently of the emission time t. This energy density is redshifted so its contribution to the present-day GW energy density is ∆ρ GW (t)(R(t)/R 0 ) 4 Gσ 2 (R(t)/R 0 ) 4 , where R(t) and R 0 = 1 are the scale factors of the Universe at time t and at present.
This argument clearly shows that the largest contribution to the present GW energy density spectrum emitted by the string-wall network, i.e. its peak amplitude, corresponds to the latest emission time, t t ann , i.e. ρ GW | peak Gσ 2 (R(t ann )/R 0 ) 4 . Defining as usual Ω GW h 2 | peak = ρ GW | peak (h 2 /ρ c ), this leads to Eq. (2.11). Here ρ c is the present critical density, h is the reduced Hubble constant, and entropy conservation implies that g s (t)[R(t)T (t)] 3 is constant, and at present the number of entropy degrees of freedom is g s (t 0 ) 3.93 [90].
What is usually quoted is the maximum of the GW energy density spectrum at time t as a function of the wave-number at present k, which for R 0 = 1 coincides with the comoving wave-number (or the frequency f = k/(2π)), defined as (see e.g. Refs. [86,99]). This is also given by Eq. Using from above that dρ GW (t) Gσ 2 (dt/t) = Gσ 2 dln(t), we find Eq. (A.2). Thus the peak amplitude of this GW spectrum at present, for t = t 0 , coincides with the result in Eq. (2.11). Since the peak GW density of walls is emitted at annihilation, its present frequency is f peak R(t ann )H(t ann ), as given in Eq. (2.10).
The spectrum of the GWs emitted by cosmic walls cannot be determined from the previous considerations alone, but has been computed numerically for N > 1 in Ref. [98]. Figure 6 of Ref. [98] shows that there is a peak at k| peak R(t f )m a , a bump at the scale k R(t f )H(t f ) where t f is the latest time in their simulation, and the spectral slope changes at these two scales. Causality requires that for frequencies below the peak, k < k peak corresponding to super-horizon wavelengths at t ann , the spectrum goes as k 3 . This is a white noise spectrum characteristic of the absence of causal correlations [128]. At frequencies above the peak, k > k peak , the spectrum depends instead on the particular GW production model assumed. The spectrum 1/k was found analytically for a source that is not correlated at different times, i.e. that consists of a series of short events [128] (see also Ref. [129]). Ref. [98] obtained roughly a 1/k spectrum for k > k peak , with approximate slope and height of the secondary bump that depend on N .
Similarly to what is done for the string-wall system, the estimates of the GW density emitted by the string network before walls appear are based on the quadrupole formula (see e.g. Refs. [100][101][102] and references therein). The energy of the string network is E st µH −1 , where µ is the string mass per unit length. Thus ... Q ij µ, and the power emitted in GWs is P Gµ 2 . Using the same assumptions as for the walls, the energy density ∆ρ GW for strings is ∆ρ st GW (t) Gµ 2 (∆t). The resulting spectrum has been computed in Refs. [100][101][102], and the simple expression in Eq. (2.12) provides a good fit to all of them, for frequencies f > f st cut . The cutoff frequency f st cut in Eq. (2.14) corresponds to the frequency emitted at T w when walls appear and the string network ceases to exist.