Catastrogenesis with unstable ALPs as the origin of the NANOGrav 15 yr gravitational wave signal

In post-inflation axion-like particle (ALP) models, a stable domain wall network forms if the model's potential has multiple minima. This system must annihilate before dominating the Universe's energy density, producing ALPs and gravitational waves (a process we dub"catastrogenesis,"or"creation via annihilation"). We examine the possibility that the gravitational wave background recently reported by NANOGrav is due to catastrogenesis. For the case of ALP decay into two photons, we identify the region of ALP mass and coupling, just outside current limits, compatible with the NANOGrav signal.


INTRODUCTION
The fundamental importance of gravitational waves (GWs) as messengers of the pre-Big Bang Nucleosynthesis (BBN) era, a yet unknown epoch of the Universe from which we do not yet have any other remnants, cannot be overestimated.The NANOGrav pulsar timing array collaboration has recently reported the observation of a stochastic gravitational wave background [1] in 15 years of data, and has examined its possible origin in terms of new physics [2].They showed that the pre-BBN annihilation of cosmic walls provides a good fit to their signal, both as the sole source and in combination with a background from a population of inspiraling supermassive black hole binaries (SMBHBs), which is expected to be its primary conventional physics origin [2].The annihilation produces a peaked spectrum, whose peak frequency f peak is given by the inverse of the cosmic horizon ≃ t ann at annihilation redshifted to the present.In their fit to the wall annihilation model NANOGrav finds [2] a peak frequency and a peak energy density with coefficients c f and c Ω of order 1.In particular c Ω ≃ 1, while c f can have larger values.
Here we consider the annihilation of a U (1) pseudo Nambu-Goldstone boson stable string-wall system as the origin of the NANOGrav signal, based on our previous recent work [3][4][5], to which we refer often in the following.
We assume that the spontaneous symmetry breaking happens after inflation, in which case cosmic strings appear during the spontaneous breaking transition, and a system of cosmic walls bounded by strings is produced when the explicit breaking becomes dynamically relevant, when t ≃ m −1 a .The cosmic strings then enter into a "scaling" regime, in which the number of strings per Hubble volume remains of order 1 (see e.g.[35] and references therein).The subsequent evolution of the stringwall system depends crucially on the number of minima of the potential after the explicit breaking, which may be just one minimum, N = 1, or several, N > 1.With N = 1, "ribbons" of walls bounded by strings surrounded by true vacuum form, which shrink very fast due to the pull of the walls on the strings, leading to the immediate annihilation of the string-wall system (see e.g.[36]).
We concentrate on the N > 1 case, where the U (1) symmetry is broken into a discrete Z N symmetry, in which each string connects to N walls forming a stable string-wall system.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.Thus its energy density is ρ walls ≃ σ/t, arXiv:2307.07665v3[hep-ph] 28 Nov 2023 where σ is the energy density per unit area of the walls.The energy density in this system grows faster with time than the radiation density, and would come to dominate the energy density of the Universe, leading to an unacceptable cosmology [37], unless it annihilates earlier.
If the Z N is also an approximate symmetry, then there is a "bias," a small energy difference between the N minima, which chooses one of them to be that with minimum energy.The energy difference between two vacua at both sides of each wall accelerates each wall toward its adjacent higher-energy vacuum, which drives the domain walls to their annihilation [37] (see also e.g.[38]).As in our previous recent work [3][4][5], we adopt the Z N explicit breaking term in the scalar potential originally proposed for QCD axions [39,40], and parameterized as V bias ≃ ϵ b v 4 , with a dimensionless positive coefficient ϵ b ≪ 1.For small enough ϵ b values, ALPs are dominantly produced when the string-wall system annihilates, together with GWs, a process that we named "catastrogenesis" [4], after the Greek word καταστροφή, for "overturn" or "annihilation." The emission of GWs by the initial system of cosmic strings ends when walls appear.Thus, there is a lowfrequency cutoff of 82 (m a /GeV) 1/2 Hz [36,41,42], corresponding to the inverse of the cosmic horizon when walls appear, redshifted to the present.This is much higher than the relevant frequencies for m a ≃ GeV, so strings do not contribute to the NANOGrav signal in this model.
We assume radiation domination during the times of interest.In this case, the present peak GW density is related to the temperature at annihilation T ann by where g ⋆ and g s⋆ are the energy and entropy density numbers of degrees of freedom.Thus, Eq. ( 1) also gives T ann in terms of c f T ann ≃ 82.5 c f MeV 16.5 g ⋆ (T ann ) 1/2 g s⋆ (T ann ) 16.5 while in terms of the parameters of our model it is The peak energy density is where f σ is a parameter entering into the definition of the energy per unit area of the walls, σ ≃ f σ v 2 V /N , and f σ ≃ 6 for most assumed potentials.We include in Eq. ( 6) a dimensionless factor ϵ GW found in numerical simulations (e.g.[43]).When needing to fix its value we use ϵ GW = 0.7 as adopted in the NANOGrav fit [2] following [44] (in our earlier work we took instead ϵ GW = 10, using Fig. 8 of [43]).Since g ⋆ = g s⋆ for T > 1 MeV, we set them equal in the following.We address the reader to [3,4] for the derivation of these equations.
Our previous results [3,4] show that the requirement that the ALP density not exceed that of dark matter, Ω a h 2 ≲ 0.12, implies so the model cannot produce the NANOGrav signal with stable ALPs.Thus we concentrate on ALPs that are unstable and decay into SM products that thermalize early enough to leave no trace by the time of Big Bang Nucleosynthesis (BBN), such as we considered in [5].To escape existing laboratory, astrophysical, and cosmological limits on ALP decays into SM products, these ALPs must have a mass m a in the GeV range or higher, depending on the decay mode (see e.g.[5] and references therein).Similar or related models have been studied recently in relation to pulsar timing array data, e.g.[45][46][47][48][49][50][51].[50] considered the same type of models we study here, but with the purpose of excluding parameter regions disfavored by the NANOGrav 15 yr data, which they analyzed independently.Our purpose is instead to try to explain the signal, and thus we stay away from the disfavored region (shown in gray in the lower left panel of Fig. 12 of [2] and the right panel of Fig. 2 of [50]).

UNSTABLE ALP MODELS THAT CAN PRODUCE THE NANOGRAV SIGNAL
In [5] we assumed m a was sufficiently larger than 1 GeV for ALPs decaying into SM particles to comfortably escape existing experimental limits.However, we need to be more nuanced here and explore the viability of somewhat lighter ALPs.The reason is that requiring the walls to form at least one order of magnitude in temperature after strings appear, combined with upper limits on T ann determined by NANOGrav to explain the signal, impose m a ≲ 1.8 GeV, as we are going to show now.
Walls appear when the Hubble parameter is H(T w ) ≃ m a /3, i.e. when the temperature is Thus T w depends only on m a (g ⋆ (T w ) ≃ 105, since T w > 100 GeV).As in our previous papers, we consider m a to be temperature independent.A temperature dependence would not affect the annihilation process, which happens late enough for m a to have reached its present constant value in any case, but could affect T w .
Combining Eqs. ( 2) and ( 6) fixes the ratio V 2 /ϵ b , and Eqs. ( 4) and ( 5) fix the product ϵ b m a .Thus, given Eqs.( 1) and ( 2), we obtain V as a function of m a , ) and consequently m a in terms of the ratio .
(10) We require T w /V ≲ 0.1 so walls form at least one order of magnitude in temperature after strings appear.Larger values of N are favorable to allow larger m a .To our knowledge, upper limits on N have been studied only in QCD axion models [52] in which N = 20 is possible.For axions coupled to gluons, N is given by the color anomaly coefficient.Similarly, non-perturbative effects in a dark sector [25,28], generically lead to N > 1 for ALPs, thus possibly to similarly large N values.We will thus adopt N = 20.Replacing also f σ = 6 and ϵ GW = 0.7, An upper limit on c f thus provides an upper limit on m a and vice versa (since the NANOGrav fit prefers c Ω ≃ 1 [2]).Looking in Fig. 12 of [2], the range of annihilation temperatures (called T ⋆ in that paper) where the NANOGrav signal can be explained by the annihilation of domain walls into SM products (DW-SM, the model most similar to ours), we can see that T ann ≲ 1 GeV (close to the upper boundary of the red region in the figure).By Eq. ( 4), this corresponds to c f ≲ 15 (taking into account the rapid change of g ⋆ values for temperatures in the 100 MeV range, g ⋆ (1 GeV) ≃ 70).Through Eq. (11), this implies m a ≲ 1.8 GeV.A more conservative upper limit on the annihilation temperature is the upper boundary of the 95% credible interval including a SMBHB contribution quoted in the text of [2], 843 MeV.This implies through Eq. ( 4) (with g ⋆ (0.84 GeV) ≃ 68) c f ≲ 13, and through Eq. ( 11) m a ≲ 1.5 GeV.
Let us now consider the experimental limits on ALPs coupled to photons through a Lagrangian term where F µν is the electromagnetic field tensor and F µν its dual, |c γγ | is a dimensionless coupling constant, and f a = V /N is given by Eq. ( 9) divided by N , and is independent of N , thus . (13) Or using Eq.(11) in Eq. ( 13), So a larger c f (thus also a larger T ann ) makes the coupling smaller.
Assuming |c γγ | ≲ 1, these upper limits on 1/f a translate into upper limits on the ALP coupling to photons as a function of m a .These limits constitute the upper boundary of the gray region in Fig. 1.Fig. 1 also shows relevant regions rejected by the most up-to-date limits on ALPs.Notice that if |c γγ | > 1, the region extends upward (as indicated by the dashed lines) where experimental limits (not only astrophysical limits) become important.
The value of |c γγ | depends on the completion of the ALP model.It has been extensively studied only for the QCD axion, where |c γγ | ≃ α EM /8π ≃ 2.9 × 10 −4 in the simplest models.However, |c γγ | can be many orders of magnitude, even exponentially, larger in some models (see e.g. the "Axions and Other Similar Particles" review in [53] or [34,[54][55][56][57]).Notice that with the |c γγ | value in the simplest QCD axion models, the upper boundary of our region of compatibility (gray) in Fig. 1 would move to the dot-dashed line in the region excluded by BBN limits (yellow), i.e. the region of compatibility would not exist.
To compute the density of the string-wall system with respect to that of radiation at annihilation, we consider that, had the system not annihilated, its energy density ρ walls ≃ σ/t would have continued to grow until the moment we call wall-domination t wd , at which it becomes as large as the radiation energy, ρ walls (t wd ) ≃ ρ rad (t wd ).The temperature of wall-domination is (see [4,5]) Besides, ρ walls (t ann )/ρ walls (t wd ) ≃ t wd /t ann , and the ratio of radiation densities at wall-domination and annihilation is given by the ratio of g ⋆ T 4 at each temperature.Combining these equations we find Using Eq. ( 4) and Eq. ( 9) in Eq. ( 17), we find which shows that this ratio is always < 1 for the annihilation temperatures we consider.Since practically all the density in the string-wall system goes into nonrelativistic (or quasi-nonrelativistic) ALPs at annihilation, considering the redshift of the ALP and radiation densities until ALPs decay at temperature T decay , As we mentioned above, T decay can be very close to T ann , so ALPs do not get to matter dominate in our model and the decays happen early enough for the products to thermalize long before BBN.Otherwise, there would be a period of ALP matter domination before ALPs decay, which is in principle not problematic since the decays happen much before BBN, but would be a scenario deserving further study.The range of c f values, 2.9 to 15, that we have found above corresponds to a peak frequency range through Eq. ( 1).In Fig. 2 we indicate two approximate spectra, with the maximum and minimum f peak in the mentioned range.Frequencies f < f peak correspond to superhorizon wavelengths at annihilation, so causality requires a ∼ f 3 dependence [68] for wavelengths that enter into the horizon during radiation domination, see e.g.[69][70][71].For f > f peak the spectrum depends instead on the particular production model.[43] finds a roughly 1/f dependence (although the approximate slope slightly depends on N ), which we use for Fig. 2. In Fig. 2 the rough signal region of NANOGrav, as well as the limits and reach of other GW observatories, is shown.

POSSIBILITY OF PRIMORDIAL BLACK HOLE FORMATION
The formation of primordial black holes (PBHs) during the process of annihilation of the string-wall system is an exciting possible aspect of ALP models with N > 1.We recently dealt with the possibility of producing "asteroidmass" PBHs, in the range in which they could constitute all of the dark matter, in [5].If formed, the PBH mass in the models in the present paper would be in the range of 0.1 to a few solar masses, but PBH abundance would be too large to be allowed, and this would reject these models.
However, the formation of PBHs is uncertain.The argument for PBH formation, first presented in [84] for QCD axions, is that in the latest stages of wall annihilation in N > 1 models (t > t ann ) closed walls could arise and collapse in an approximately spherically symmetric way.In this case, if the characteristic linear size of the walls continues to grow with time after annihilation starts, some fraction of the closed walls could reach their Schwarzschild radius R Sch and collapse into PBHs.The figure of merit used is p(t) = R Sch /t = 2M (t)/(t M 2 P ), where M P is the Planck mass and M (t) is the mass within the collapsing closed wall at time t.Reaching p(t) = 1 would indicate the formation of PBHs.This definition is based on the fact that while walls are in the scaling regime, the linear size of the walls L is close to the horizon size (L ≃ t).
Annihilation starts when surface tension of the walls, which produces a pressure p T ≃ σ/t that decreases with time (which tends to rapidly straighten out curved walls and peak frequencies f peak = 2.9 × 10 −8 Hz and 1.5 × 10 −7 Hz, which are the minimum and maximum frequencies we found (see text).Allowed spectra would have peak frequencies in between these two.Also shown are the approximate NANOGrav 15 yr signal [2] (in purple) and limits (solid line boundaries) or reach (dashed line boundaries) of other GW detectors: the European Pulsar Timing Array (EPTA) [72] and the Square Kilometre Array (SKA) [73] in purple; the space-based experiments TianQin [74], Taiji [75], and the Laser Interferometer Space Antenna (LISA) [76] in green; the Atom Interferometer Observatory and Network (AION) [77], the Atomic Experiment for Dark Matter and Gravity Exploration in Space (AEDGE) [78], the Deci-hertz Interferometer Gravitational wave Observatory (DECIGO) [79], and the Big Bang Observer (BBO) [80] in blue; the ground-based Einstein Telescope (ET) in red [81]; and the Laser Interferometer Gravitational-wave Observatory (LIGO) in gray [82].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 [83], which imposes Ω GW h 2 < 10 −6 .
to the horizon scale t), is compensated by the volume pressure p V ≃ V bias (which tends instead to accelerate the walls toward their higher-energy adjacent vacuum).In our model, p V ≪ p T when walls form.At a later time, when p T ≃ p V , the bias drives the walls (and the strings bounding them) to annihilate within a Hubble time.This defines t ann ≃ σ/V bias , after which V bias dominates the energy density.After annihilation starts, L ≃ t is no longer guaranteed.
We have checked that for the models in this paper, we always have p(t ann ) < 1.If L continues being close to t for t > t ann , then p(t > t ann ) ≃ V bias L 3 /L grows with time as t 2 and eventually reaches 1.However, if L decreases with time at some point after annihilation starts, the figure of merit may never reach 1.Based on the simple power-law parameterization we used in our previous recent work [4,5] for the evolution of the energy density after annihilation starts, namely ρ walls (T )/ρ walls (T ann ) ≃ (T /T ann ) α (with a parameter α that needs to be extracted from simulations), we can make a naive estimate of how the characteristic linear wall size L within a Hubble volume t 3 evolves with time.In Appendix A we show how this naive estimate requires α < 6 for L to ever become larger than t ann .The only simulations of the annihilation process available [85] find α ≥ 7 [4,5].On the other hand, they also seem to indicate that the evolution of the string-wall system continues being close to that in the scaling regime for some time.Therefore, more detailed simulations of the annihilation process are required to elucidate the appearance of PBHs.
In addition, a large enough departure from spherical symmetry due to angular momentum or vacua with different energy on different sides of the collapsing closed wall could prevent the formation of PBHs.Since the formation of PBHs is such an uncertain consequence of ALP models with N > 1, we do not use this feature to reject any of these models.

CONCLUSIONS
We pointed out that the recently confirmed stochastic gravitational wave background could be due to pseudo Nambu-Goldstone bosons, whose existence could only be revealed through their decays and this background.In particular, we examined unstable ALP models which can produce the recent NANOGrav 15 yr signal.ALP models have a complex cosmology in which a stable system of walls bounded by strings develops (for N > 1), and nonrelativistic ALPs and gravitational waves are produced when the cosmic string-wall system annihilates (a process we dubbed "catastrogenesis" in our recent work on these models).The annihilation produces a distinctive peaked spectrum, at a frequency corresponding to the inverse of the cosmic horizon at annihilation.Thus, this peak frequency is related to the annihilation temperature.
We require ALPs to decay into Standard Model (SM) products which thermalize much before BBN.In particular, we have shown that ALPs decaying into two photons in the region of masses and couplings necessary to explain the signal can evade existing observational limits, the most relevant of which are derived from supernova data (see Fig. 1), for ALP masses from about 300 MeV to 1.8 GeV.The model closest to ours that NANOGrav fitted to their signal is that of domain walls decaying into SM products (DW-SM).Our model is very similar to this one if the ALP decay happens very shortly after string-wall system annihilation, which we showed is possible.Thus we use the NANOGrav fits to this model to select a range of annihilation temperatures and thus peak frequencies.
We have found a range of c f values (as defined in Eq. ( 1)) which corresponds to the range of peak frequencies from f peak = 2.9 × 10 −8 Hz to f peak = 1.5 × 10 −7 Hz.This corresponds to annihilation temperatures in the range 200 MeV to 1 GeV.This temperature range overlaps with the upper portion of the 68% credible interval (which goes to 275 MeV) and the 95% credible interval (which goes to 505 MeV) quoted by NANOGrav [2] if their DW-SM model is the sole origin of the signal.Considering their fit done with the addition of a SMBHB contribution, our temperature range overlaps with a larger portion of both the 68% credible interval (which goes to 309 MeV) and the 95% credible interval (which goes to 843 MeV) quoted in the text, and is included within the red region in the lower left corner of Fig. 12 of [2] (for its DW-SM + SMBHB fit).
The upper portion of the region of ALP-photon coupling and mass necessary to explain the NANOGrav signal (shown in Fig. 1), for couplings above 10 −7 GeV −1 , could be tested in the future by DarkQuest, HIKE-dump and SHiP, as shown e.g. in Fig. 133 of [66] (see references therein).The lower portion would be tested if a new supernova is observed by the Supernova Early Warning System (SNEWS) [86], which would allow to extend considerably all the limits derived from SN1987A.

Appendix A
Before annihilation starts, the energy density of the walls in the scaling regime is ρ walls ≃ σ/t ≫ V bias .The annihilation of the string-wall system starts when the bias volume energy density, or magnitude of volume pressure, V bias becomes of the same order as the energy density, or surface tension, of the walls σ/t (t ann ≃ σ/V bias ), after which V bias dominates and accelerates walls towards the higher-energy vacuum adjacent to each wall.If PBHs do not form at annihilation, i.e. p(t ann ) < 1, the energy contained in a closed wall will need to increase with time for PBHs to form later, at a time t ⋆ such that p(t ⋆ ) = 1.Since the energy density V bias is constant, this requires that the dimensions of the closed walls keep growing for t > t ann .In fact, if the characteristic linear dimension L of walls continues being close to t, L ≃ t, then p(t > t ann ) ∼ V bias L 3 /L grows with time as t 2 and eventually reaches 1.However, if L decreases with time and never becomes larger than t ann , the figure of merit p(t) decreases after annihilation starts and never reaches 1.
Based on the simple power-law parameterization we used in our previous recent work [4,5] for the evolution of the energy density after annihilation starts, for T < T ann , ρ walls (T ) ρ walls (T ann ) with a real positive power α that needs to be extracted from simulations of the annihilation process, we can make a naive estimate of how the characteristic linear wall size L within a Hubble volume t 3 evolves with time.It is easy to do it either assuming that walls dominate the energy density, or that volume density dominates.In both cases we find the same condition on α for L to continue growing with time, i.e.L > t ann for t > t ann .Therefore it is reasonable to assume that the same condition holds in the transition period, when both volume and walls contribute significantly to the energy density of the annihilating string-wall system.If the energy in walls still dominates and requiring L/t ann > 1 for t > t ann , means that (t/t ann ) (6−α)/4 > 1, i.e. (6 − α)/4 > 0, thus α < 6.
In both cases we find the condition α < 6 for L to become larger than t ann after annihilation stars.However the only simulations available to estimate values of α [85] lead to α ≥ 7 (see [4,5] for details).In this case, with our naive estimates the linear size of walls would decrease with time after annihilation starts and PBHs would not form if p(t ann ) < 1.

Figure 2 .
Figure2.Approximate spectra which could account for the NANOGrav signal in our catastrogenesis model, with peak amplitude Ω GW h 2 peak = 10 −8 and peak frequencies f peak = 2.9 × 10 −8 Hz and 1.5 × 10 −7 Hz, which are the minimum and maximum frequencies we found (see text).Allowed spectra would have peak frequencies in between these two.Also shown are the approximate NANOGrav 15 yr signal[2] (in purple) and limits (solid line boundaries) or reach (dashed line boundaries) of other GW detectors: the European Pulsar Timing Array (EPTA)[72] and the Square Kilometre Array (SKA)[73] in purple; the space-based experiments TianQin[74], Taiji[75], and the Laser Interferometer Space Antenna (LISA)[76] in green; the Atom Interferometer Observatory and Network (AION)[77], the Atomic Experiment for Dark Matter and Gravity Exploration in Space (AEDGE)[78], the Deci-hertz Interferometer Gravitational wave Observatory (DECIGO)[79], and the Big Bang Observer (BBO)[80] in blue; the ground-based Einstein Telescope (ET) in red[81]; and the Laser Interferometer Gravitational-wave Observatory (LIGO) in gray[82].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[83], which imposes Ω GW h 2 < 10 −6 .