High Quality QCD Axion at Gravitational Wave Observatories

The axion solution to the strong CP problem is delicately sensitive to Peccei-Quinn breaking contributions that are misaligned with respect to QCD instantons. Heavy QCD axion models are appealing because they avoid this so-called"quality problem". We show that generic realizations of this framework can be probed by the LIGO-Virgo-KAGRA interferometers, through the stochastic gravitational wave (GW) signal sourced by the long-lived axionic string-domain wall network, and by upcoming measurements of the neutron and proton Electric Dipole Moments. Additionally, we provide predictions for searches at future GW observatories, which will further explore the parameter space of heavy QCD axion models.

Introduction-A great amount of experimental effort has been aiming at discovering the QCD axion [1,2], the pseudo-Goldstone boson of a spontaneously broken axial U (1) Peccei-Quinn (PQ) symmetry [3,4] that explains the smallness of CP violation in strong interactions.
While attractive, the PQ mechanism is vulnerable to possible additional sources of symmetry breaking, generically misaligned with respect to the axion potential from QCD instantons. This quality problem (originally formulated with various perspectives in [5][6][7][8][9][10][11]) is alleviated in heavy axion models (see [6] and [12][13][14] for earlier related work), where a "heavy QCD" sector provides a larger contribution to the axion potential, aligned with that from QCD instantons.
Existing realizations of this idea rely on: the QCD coupling becoming strong at high energies [6,[15][16][17][18][19], see also [20] for a 5D model; a separate confining gauge group, whose alignment is ensured by unification at high scales [14,[21][22][23] or by a softly-broken Z 2 symmetry [24][25][26]. When the strong coupling scale Λ H of the heavy sector is above the QCD scale Λ QCD , the axion mass is larger than in the standard window, and the cosmological evolution of the axion field in the early Universe is shifted to higher energy scales. Despite its appeal, it is not immediately clear what the signatures of such a scenario are since generically the axion can be very heavy, e.g. above the electroweak scale, while its interactions remain very weak. 1 Furthermore, in contrast to the standard case, a heavy QCD axion can easily decay in the early Universe, and thus leaves no detectable relic dark matter today.
GWs are indeed radiated [34] by the network of axionic topological defects (domain walls attached to strings) [35] (see also [36]), which are abundant in the early Universe if the PQ symmetry is broken after inflation.
In standard QCD axion models the network necessarily annihilates while making up only a very tiny fraction of the energy density of the Universe, and therefore the GW signal is too weak to be detectable [37]. In contrast, the heavy QCD axion network can carry much more energy because of its larger domain wall tension. Furthermore, in generic realizations (e.g. DFSZ [38,39] and simple generalizations of KSVZ models [40,41]), the network can be long-lived while still avoiding the overproduction of relics, since radiated axion quanta are unstable. Annihilation of the network can be triggered by the misaligned PQ breaking effects that motivate the scenario in the first place (see [42]). These also induce a small but potentially observable shift of the QCD vacuum angle.
Our work points out a new source of observable gravitational waves from the dynamics of the QCD axion (see e.g. [43][44][45][46] for previous work, unrelated to the axion quality problem and [47][48][49][50] for related scenarios with ALPs).
The Heavy QCD Axion-Heavy QCD axion models are characterized by an extra contribution to the axion potential that is larger than and aligned with the contribution from QCD. The zero temperature potential is For our discussion, it is important to recall that gauge instantons generically break the original U (1) PQ symmetry to a discrete Z NDW subgroup, where N DW is a model-dependent integer number related to the axion coupling to gluons. Therefore the periodicity 2πf induced by the potential (1) can be smaller than the fundamental axion field range 2πf N DW and the potential V a can feature N DW degenerate minima. In writing (1), we assumed that the periodicity induced by the heavy sector coincides with that of QCD instantons. This appears to be linked to the requirement of alignment between the two sectors, as is evidently the case in constructions with a Z 2 symmetry [26] and in simple unification frameworks where SM and heavy sector fermions descend from the same fundamental representation of a higher-rank gauge group [18,22]. This feature implies that the low-energy QCD-induced potential does not lift the degeneracy of the N DW minima.
Generically, however, we may expect further contributions to the axion potential, misaligned with V a . Independently of its specific origin, such a contribution can be written as where N b defines the subgroup Z N b of the PQ symmetry which is preserved by (3) and δ is a CP violating phase.
In the absence of tuning, this offset is naturally O(1) and µ b Λ H is required to solve the strong CP problem. The low temperature potential is V = V a + V b and when N b = 1 or is co-prime with N DW , the degeneracy of the N DW minima is lifted. In particular, the vacuum energy difference between the global CP preserving minimum and its nearest neighbor is of the order ∆V µ 4 b [1 − cos(2πN b /N DW )] (provided that δ is not too close to π/N DW ). Broadly speaking, (3) can originate at a scale Λ b , such that f Λ H can arise from UV physics via: non-perturbative effects, κ b ∼ e −S/2 (see e.g. [51][52][53]); higher-dimensional operators when the axion is the phase of a complex scalar field (see e.g. [9]); another gauge sector with confinement scale Λ b and a light fermion of mass m q , κ b ∼ (m q /Λ b ) 1/2 . Λ b Λ H can also arise from a confining gauge sector. Further details are provided in Appendix A.
Despite its smallness, a contribution from (3) can lead to potentially observable CP violation. In particular, at low temperatures one finds: where r ≡ µ b /Λ H . Current bounds from nEDM measurements [54] require ∆θ 10 −10 . Clearly, (4) shows that Λ H Λ QCD makes the PQ mechanism more robust against misaligned contributions.
In the early Universe, the mass m a and the scale µ b are generally temperature-dependent, for instance in the standard QCD axion case m a (T ) m a (T 0 /T ) 4 for T ≥ T 0 134 MeV and m a (T ) = m a otherwise [55]. Nonetheless, our results are mostly independent of the detailed temperature dependence.
Axionic Defects-Let us now move to the cosmological evolution of topological defects, whose history begins at the PQ symmetry breaking scale ∼ N DW f . Our investigation concerns scenarios where this occurs during radiation domination after inflation, which will be generic for the values of f considered in this work.
Axionic strings form at T N DW f and continuously radiate axion quanta and gravitational waves. In the absence of significant friction due to the plasma, they quickly achieve a scaling regime [56,57] (see also [58] and [59][60][61][62] for recent updates), with energy density scaling as ρ s = λµH 2 , with λ a O(1) parameter and µ ∼ N 2 DW f 2 the string tension. This behavior is altered once 3H m a (T ). This occurs at a temperature T osc Λ H (see Appendix B) when the axion field, with average initial value a i /(N DW f ) ∼ O(1), starts oscillating in its potential V a and domain walls form, attached to the strings, with a tension σ 8m a (T )f 2 . At this epoch, two possibilities arise: i) when N DW = 1, the network of topological defects is rapidly annihilated by string-wall interactions (see e.g. [58]); ii) when N DW > 1, the network persists because multiple domain walls pull each string in different directions.
In both cases, the tension of the walls is larger than in the standard QCD axion case by a factor Λ H /Λ QCD 1. For N DW > 1, in the absence of significant friction from the plasma (we show in Appendix D that this has a minor impact on our conclusions), the network rapidly achieves a scaling regime, with its energy density dominated by domain walls, ρ dw cσH, where c is a O(1) numerical prefactor (in this regime V a is normally already temperature-independent). This scales slower than matter and radiation and thus the network is potentially dangerous for cosmology [63]. However, domain wall domination can be generically avoided in the heavy axion scenario, thanks to the misaligned potential contribution V b . The resulting vacuum pressure causes the contraction of the false vacuum regions and the collapse of the network [36] at a temperature T ann , which can be estimated by imposing ρ dw ∆V and more precisely determined via numerical simulations [64]. Here we focus on the case where V b is temperature-independent below Λ H , as occurs generically when PQ breaking is due to physics above Λ H , see Appendix B for the temperature-dependent case. To set ideas and simplify expressions, in the following we set N b = 1, N DW = 6 as example values, fix numerical prefactors according to the simulations of [64] and also fix the number of (entropy) relativistic degrees of freedom at T ann to the SM value at high temperatures (g * s,ann ) g * ,ann = 106.75 (see also Appendix B), although our results are only mildly affected by these precise choices. We then find showing that for r 1 network annihilation is significant delayed.
At (5) the network collapses and its energy density is transferred mostly to mildly relativistic axion quanta (see e.g. [58] and [64]). In contrast to the standard QCD axion case, in the heavy axion scenario these relics can efficiently decay to SM gluons, above the QCD Phase Transition (PT) and to photons and/or fermions, above and below the QCD PT, depending on the specific axion model. Focusing on the decay to gluons, since m a GeV in most of the parameter space of interest [65], we find that decay is efficient below the temperature obtained by setting Γ a→gg H (see Appendix B). This temperature can be larger than T ann for r 0.001 and/or f 10 12 GeV. Therefore, axion relics from the network will in general decay immediately. Crucially, however, the string-wall network can source a significant relic abundance of gravitational waves [70][71][72][73][74][75][76]. The simple quadrupole estimate for their energy density ρ gw (T ann ) ∼ c 2 σ 2 /(32πM 2 p ) has been confirmed by numerical simulations [76] (see also [37]). Assuming a standard radiation-dominated cosmological history after domain wall annihilation, one finds that the relic abundance of gravitational waves today is where˜ 0.1 − 1 is a numerical efficiency factor [76] and ρ rad and Ω rad h 2 4 · 10 −5 are the energy density and relic abundance of radiation today respectively. The formula above shows that when the network makes up O(5%) fraction of the energy density of the Universe at annihilation, its gravitational wave signal is detectable by present interferometers, i.e. Ω gw h 2 ∼ 10 −9 . This fraction at the annihilation temperature reads Constraints are also shown, as dark-shaded regions, from: domain wall domination (lower right corner), nEDM [54] (upper part), LIGO-Virgo O3 run [66] (dark-blue shaded). Dashed contours bound regions probed by LVK at design sensitivity and ET (sensitivity curves taken from [67]). The gray shaded region will be also probed by neutron [68] (dot-dashed line) and proton EDM [69] measurements.
The GW signal is peaked at a frequency corresponding to H at annihilation (see e.g. [37]). Redshifted to today: According to (7), (8) and (9), the signal from a heavy axion with f 10 11 GeV, Λ H 10 10 GeV and r 10 −3 sits right in the reach of the LVK interferometers [77,78]. The GW spectrum away from the peak frequency [37] decreases as ω 3 for ω < ω peak , whereas for ω > ω peak it behaves as ∼ ω −1 , until a cutoff frequency corresponding to the domain wall width. However, further numerical simulations are required to understand the precise behavior of the spectrum around the peak frequency.
Predictions-Although the N DW = 1 case does not leave observable GW signals (see Appendix E) due to the quick decay of the network, the situation is radically different for N DW > 1, where network annihilation is delayed. To simplify the presentation, we fix N b = 1, N DW = 6, κ H = 1 and g * ,ann = g * s,ann = 106.75 (see Appendix B for the case κ H 1) and δ = 0.3 and present results varying Λ H , r ≡ µ b /Λ H and f . According to (4), r can then be traded for ∆θ.
We first present results for large values of Λ H , which maximally reduce the sensitivity to misaligned contributions. We consider V b to be temperature-independent at T Λ H and Λ H to have a QCD-like temperature behavior (see also Appendix C).
Fixing Λ H = 10 10 GeV as a representative example, we show values of r and f that can be probed by gravitational wave observatories (dashed contours), together with constraints (solid contours), in Fig. 1. In the lower right half, the string-wall network dominates before annihilation. While this region might not be completely ruled out, annihilation of the network in this case would require a dedicated study. In the upper part of the parameter space the PQ solution is spoiled, i.e. ∆θ 10 −10 [54]. In the dark blue shaded region the GW signal is incompatible with the latest 2σ upper bound from LIGO-Virgo (LV) [66], this corresponds to the region close to DW domination. The dashed blue contours bound regions where the GW signal is detectable at the design sensitivity of LVK and Einstein Telescope (ET) respectively. The change of slope in the GW regions arises because of an intermediate phase of matter domination driven by the axions produced by the string-wall annihilation. This occurs at small decay temperatures (6), corresponding to the right half of the figure.
Very interestingly, we find that a significant fraction of these GW-observable regions also predicts a detectable nEDM (and/or pEDM) in the near future [32,68,69], i.e. ∆θ 10 −12 (∆θ 10 −14 ) (above the dot-dashed and dotted gray contours respectively). Motivated by the exciting possibility of a combined heavy axion discovery via nEDM experiments and GW observatories, we fix ∆θ = 8 · 10 −13 and broaden our analysis to different values of Λ H in Fig. 2. We find that any Λ H 10 6 GeV leads to a GW signal in the foreseen reach of future ground (design LVK, ET, CE) and space based interferometers (LISA [79,80], DECIGO [81,82]). The lower right corner in the figure is strongly constrained by the slow decay of axion quanta and a phase of matter domination which spoils BBN. As in Fig. 1, the change of slope in the GW contours is due to an intermediate phase of matter domination, in the lower right part of the figure.
Values of Λ H smaller than 10 6 GeV can also lead to viable cosmologies, observable GWs and detectable nEDM and/or pEDM, if the potential V b is temperaturedependent below Λ H , see Appendix C.
Finally, let us mention that for Λ H 10 10 GeV, LVK are expected to probe the high frequency tail of the GW signal (∼ ω −1 ), ET can investigate the peak and LISA can probe the low frequency tail (∼ ω 3 ). Full GW spectra for some representative choices of parameters are shown in Fig. 3.
A caveat is in order before our conclusions: the parameter space shown in Fig. 1 and Fig. 2 is further constrained if V b arises from dimension-five (and to a less relevant extent from dimension-six) operators with large coefficients. While we discuss this quantitatively in Appendix A, we note here that a large region of parameter space remains unaffected if such operators originate from non-perturbative effects (as can be expected if they are due to gravity, see e.g. [7,52,83]).
Conclusions-Heavy QCD axion models can feature a long-lived network of topological defects. The main finding of this Letter is that these models predict: i) a stochastic gravitational wave signal measurable by the design LVK interferometers in a large region of parameter space (further broadened by ET and CE, with the possibility of correlated signals also at LISA and DECIGO); ii) a nEDM and/or pEDM measurable in the near future; when: a) the new "heavy QCD" scale is large, i.e. Λ H 10 10 GeV, thus making the PQ mechanism more robust; and b) misaligned PQ breaking terms that motivate these models in the first place are not strongly suppressed.
Furthermore, we showed that combined GW (at LISA and DECIGO) and nEDM/pEDM signatures also arise for 10 6 GeV Λ H 10 10 GeV.
Our results do not strongly depend on the specific heavy QCD axion model, as long as its potential has approximately degenerate minima.
We necessarily left several interesting points for future work. First, in order to precisely characterize the GW signal, numerical simulations of axionic stringwall networks beyond the current literature, possibly including friction and plasma effects, are crucial [64].
Second, we left unspecified the particle content and properties of the "heavy QCD" and of the misaligned PQ breaking sectors. However, these sectors may contain dark matter candidates (see e.g. [84]) or light states that can contribute to the number of extra relativistic degrees of freedom, ∆N eff [85,86].
Furthermore, it is interesting to understand whether the collapse of the network of topological defects may also lead to a significant fraction of Primordial Black Holes, in a scaled-up version of the mechanism proposed in [87].
Finally, a more complete exploration of the parameter space of heavy axion models may lead to further interesting signatures. For instance, GWs of very low frequency may arise for misaligned sectors lighter than QCD, and may provide a new explanation for recent NANOGrav observations [88]. An investigation of all these aspects is ongoing [89].
Our work is relevant for the ongoing effort to probe well-motivated regions in the parameter space of the PQ mechanism. Guided by the theoretical pursuit of "higher quality" models, we suggest that gravitational wave interferometers and nEDM/pEDM experiments may be the right laboratories to discover the heavy QCD axion.
Acknowledgments-We thank Jaume Garriga, Alex Pomarol and Alex Vilenkin for useful discussions. The work of FR is supported in part by National Science In this section, we critically review well-motivated origins for misaligned contributions of the form Such contributions can be generated at a scale Λ b such that Let us first consider misaligned terms generated at Λ b f Λ H . The following possibilities can then be envisioned: 1 -Non-perturbative corrections to the axion potential, among them: stringy and/or gravitational instantons, in which case κ b ∼ e −S/2 , where S is the instanton action of interest (see e.g. the discussions in [52] for stringy instantons and [53] for gravitational instantons); instantons of another gauge sector with confinement scale Λ b and a light fermion of mass m q , with ∆ > 4 in scenarios where the axion is the phase of a complex scalar field Φ, whose renormalizable Lagrangian respects the PQ symmetry, giving The size of the coefficients c ∆ depends on the specific origin of the higher dimensional operators.
The contributions above are often discussed as a source of the quality problem, since they induce the shift of the QCD vacuum angle given in eq. (4) of the main text and reported here for the reader's convenience: ∆θ r 4 (N b /N DW )(sin δ/κ H ), with r ≡ µ b /Λ H . Imposing that the QCD axion solution to the strong CP problem is kept, the following constraints on the two contributions above arise: 1 -For contributions with κ b ∼ e −S/2 , one finds: whereas for contributions from confining gauge sectors with κ b ∼ (m q /Λ b ) 1/2 , one finds the following constraint on the light fermion mass: (A4) In both cases, it is evident that taking Λ H as large as possible provides the most robust implementation of the PQ mechanism. In particular, the lower bound on the instanton action (A3) is reduced by approximately a factor of two for Λ H 10 10 GeV with respect to the standard QCD axion. This can importantly relax the quality problem in certain setups [7]. , for ∆ = 6. (A6) The lower bounds above should be complemented with the upper bound Λ H ≤ f , which arises from consistency of the axion EFT. We have normalized Λ b to the Planck scale, to reflect the common assumption in the literature that the operators O ∆ may be induced by gravitational interactions, see e.g. [9]. The relevance of these constraints depends crucially on the size of the coefficients c ∆ . It is sometimes assumed that c ∆ ∼ O(1), which may be justified if UV physics, in particular gravity, breaks the PQ symmetry perturbatively. Under this assumption, (A5) imposes a very stringent constraint on Λ H , which can be satisfied only for 5 TeV Λ H f ∼ 10 TeV. This would severely limit the interesting region of the parameter space shown in Figs. 1 and 2 of our main text.
On the other hand, arguments against global symmetries in gravitational theories suggest that the latter may break the PQ symmetry only via nonperturbative effects (see e.g. [7,52,83,90]). In this case, higher dimensional operators may still be generated, albeit the coefficients c ∆ would then be expected to be exponentially suppressed by the corresponding instanton action, as in case 1 above. Typical values S inst ∼ O(100) (see e.g. [52,91]) render Planck-suppressed operators irrelevant for the heavy QCD axion scenario considered in this work. Nonetheless, we may consider much smaller values of the action (as suggested by ) to understand what suppression would be required to affect our parameter space, as we do in Fig. 4 (left) for dimension ∆ = 5 operators and κ H = 1. We note that even for c 5 10 −7 , corresponding to a very small instanton action S inst 15 (this roughly coincides with S inst ∼ ln M p /f for f 10 11 GeV, which is the smallest estimate of the action from non-perturbative gravitational contributions [83]), there exists a region of available parameter space where the GW signal is detectable at LIGO-Virgo-KAGRA (LVK), and that for c 5 10 −15 , corresponding to S inst 34 the LVK region is entirely available.
It is also conceivable that dimension-5 operators may be forbidden in a concrete model. In this case, the most dangerous operators would be those of dimension 6. Constraints from (A6) are shown in Fig. 4 (right). We note that in this case even the choice c 6 ∼ 1 leaves an available region of parameter space which overlaps with the LVK reach, and that for c 6 ∼ 10 −7 the entire region is available. This is dramatically different from the situation in the standard (light) QCD axion case, where all operators up to dimension 9 (or 12 depending on the value of f ) have to be forbidden in order for the axion to solve the strong CP problem.

Scenarios with Λ b
Λ H can also arise from a confining gauge sector. Then κ b ∼ (m q /Λ b ) 1/2 as in case 1 above. The quality problem in this case can then be evaded by taking Λ b sufficiently below Λ H . Nonetheless, in contrast to the heavy QCD axion case, this is difficult in the standard QCD axion case where Λ b is bounded from below by the overproduction of axion dark matter from the decay of the axion string-wall network (see e.g. [87]).

Appendix B: Network Evolution
Here we report full formulae for quantities defined in the main text, including the dependence on all parameters. We start with the temperature-dependent axion mass m a, H . We used the parametrization where T 0,H is a critical temperature, analogous to the critical temperature in QCD. Usually, T 0,H Λ H , for instance in QCD T 0 134 MeV 0.4 Λ QCD [55]. The exponent γ can be estimated analytically at T Λ H in the Dilute Instanton Gas Approximation (DIGA). For G = SU (N c ) with N f light flavors, one finds [92]: γ 11N c /6 + N f /6 − 2. Lattice calculations are required to compute γ at T Λ H . However, in the case of QCD, available lattice results agree well with the DIGA estimate [55]. In order to produce the blue shaded region in our figures, we have fixed γ = 4, T 0,H = Λ H .
The parameter κ H suppresses the heavy axion mass in models where G contains at least one light fermion ψ, in which case κ H = m ψ /Λ H . Alternatively, κ H can also be much smaller than unity in models where m a arises from small instantons of QCD itself, when the latter becomes strongly coupled in the UV, above the EW phase transition [6,[15][16][17]20]. In this case SM quarks are massless and κ H is proportional to powers of their Yukawa couplings, see e.g. [20].
The temperature at which the heavy axion starts its oscillations is determined by imposing the condition 3H = m a (T ): T osc c γ · 10 11 GeV κ .
(B2) where c γ varies by a O(1) factor with γ. For γ = 4 and f = 10 12 GeV, one finds T osc 10 11 GeV. For certain values of f, r ≡ µ b /Λ H and Λ H , the heavy axion starts oscillating in the potential V b at temperatures above (B2). Obviously, this can only happen if V b is larger than V a at high temperatures, for instance when V b is temperature-independent at temperatures slightly above Λ H (the temperature-dependent case is discussed in Appendix C). In this case, the temperature at which the axion would start oscillating around one of the minima of V b is given by In this case the scenario may still be cosmologically viable (if it satisfies other constraints), but the first domain walls to form are those of V b , rather than of V a . Their tension is significantly smaller than that of the V a walls, since r 1. Furthermore, if N b = 1, then the V b network annihilates at high temperatures and there are no walls formed at T osc . If N b is co-prime with N DW , annihilation occurs close to Λ H . Thus in this case the gravitational wave signal is suppressed compared to the case analyzed in the main text. However, we find that this occurs only in a tiny region of the parameter space, located in the upper left corner of Fig. 1 of the main text.
When T osc, b T osc , the network forms as described in the text. We assume that it quickly reaches the scaling regime, where ρ dw = cσH(T ). The annihilation temperature is then found by imposing the condition dσH = ∆V , where d is a numerical parameter. One finds: where we kept the explicit dependence on all parameters, while in the text we fixed N b = 1, N DW = 6, κ H = 1 as representative example. The parameters c and d are determined by numerical simulations [64] and depend on N DW (and presumably N b ). The relations to the parameters extracted from numerical simulations in [64] are c = 2A, d = 2C d . For instance, for our example choice N DW = 6 the numerical simulations of [64] find central values A = 2.24, C d = 3.14 (we follow the 1% criterion in [64]), so that c = 4.48, d = 6.28.
The string-wall network may dominate the energy density of the Universe, if annihilation is too slow. The temperature corresponding to network domination is T dw-dom 6 · 10 6 GeV √ cκ H 106.75 g * ,dw-dom In the shaded regions denoted with "DW domination" in Figures 1 and 2 in the main text, T dw-dom > T ann . In this case it is in general problematic to gracefully exit the domain wall dominated phase, at least in the axion scenario considered in this work. The heavy axion decay rate into SM gluons (photons) is given by: where c g = 1 and c γ = E/N DW − 1.92, with E being the model-dependent electromagnetic anomaly (see e.g. [93] for a review in standard DFSZ and KSVZ constructions) and the −1.92 arising only below the QCD PT. The decay is efficient once Γ H, which gives: T a→gg(γγ) 9 · 10 6 GeV α s(EM) κ In models where E/N = 0, the decay channel to photons is in principle available above the EWPT, but it is less efficient than the channel to gluons (we use α s = 0.1 in our Figures). Below the QCD PT, the latter remains available when m a 1 GeV [65], which happens in most of the parameter space of interest, otherwise the decay temperature should be computed according to the a → γγ rate. In the shaded blue region denoted by "Life-Time" in Fig. 2 of the main text the decay temperature is below 10 MeV. In this case products of axion decay can be dangerous for BBN. When the decay temperatures above are larger than the annihilation temperature of the network, the latter rapidly transfers its energy to the SM at T ann . The annihilation of the string-wall network leads to a population of mildly relativistic axions, which soon contributes to the dark matter abundance. When the decay rates (B6) are not fast enough, these axions can dominate the energy density of the Universe and lead to a temporary phase of matter domination. The temperature at which the matter dominated phase starts is found by imposing The temperature above should be compared with the decay temperatures (B7): whenever the former is 1 3 , which shows that the frequency can be much smaller than in the RD case, when T MD T a→gg(γγ) . The amplitude of the signal is also suppressed by the duration of the MD epoch, i.e. .
When the suppression of the amplitude due to MD is not too strong, the gravitational signal can be probed by interferometers which are sensitive at low frequencies, such as DECIGO, LISA and possibly PTAs. The spectrum of GW waves is characterized by a peak at a frequency ω peak roughly corresponding to the horizon at the time of the annihilation of the network, as explained in the main text. For smaller frequencies the spectrum decay as ω 3 , as dictated by causality, whereas for larger frequencies numerical simulations have shown a 1/ω dependence which is cut off at a frequency corresponding to the axion mass [37]. The details of the spectrum near the peak are subject to more uncertainties and motivate further numerical studies. In this work we consider the peak to have a O(1)ω peak width, as suggested by numerical simulations, and use the following simple formula for the spectrum to produce Fig. 3 of the main text: Smaller widths would not significantly change the predictions.
If the peak width is much larger than O(1)ω peak the detectability of the signal would significantly improve.
Additionally, we present results for a representative example with κ H 1 in Figures 5 and 6. In this case, we take √ κ H f ≤ f as cutoff of the axion effective theory.
Correspondingly, LIGO-Virgo-KAGRA can now probe much larger scales Λ H 10 15 GeV than in the case κ H = 1.
In the main text we focused on the case where the misaligned contribution to the axion potential is temperature independent. This is typically the case when the source of PQ breaking comes from energy scales above f . Instead, if there are infrared contributions, e.g. coming from sectors which confine at scales µ b f , the axion potential will be temperature dependent and saturate at temperatures T b µ b . In this Section we present the results for the latter case and consider a bias potential given by T n , 1 . When n = 0 we recover the case described in the main text. On the other hand, for large n the potential is strongly temperature dependent, the potential approaches a step function and the predictions become almost independent of n. The main consequence of a temperature dependent bias is that the annihilation of the network will happen at a later temperature T T b thus causing an enhancement of the relative energy in the network, ρ dw /ρ tot , and so on the GW signals. This allows to probe a larger region of parameter space, in particular smaller values of f , simultaneously with nEDM and GW experiments.

Appendix D: Friction
Particles in the primordial plasma that interact with the domain wall can exert thermal friction [56,63,94]. This is usually negligible for the QCD axion, since the axion mass is very small compared to the temperature of the plasma, but it can be relevant in the heavy axion scenario. Quantifying these effects is a modeldependent task, since axion interactions are not fixed by the requirement that the axion solves the strong CP problem. Here, we consider a conservative case, where we assume one particle in the plasma to have a probability of order one to be reflected by the axion wall, and show that even in this case friction is only important in a region of parameter space relevant for future interferometers. Even so, we find only a minor impact on the GW signal in this case.
The thermal pressure exerted on the wall by relativistic particle in thermal equilibrium is given by [70] F = P r n ∆p (D1) where n = g fr ζ(3)T 3 /π 2 is the density of relativistic particles that are reflected by the wall with a probability P r , g fr is the number of degrees of freedom and ∆p = vT the momentum transfer (for non-relativistic walls and after averaging over the two sides of the wall [70]). If the friction caused by the thermal pressure, F/σ, overcomes the Hubble friction, H, the wall will reach an attractor regime where the self-acceleration due to its own tension balances the friction force [70] σ/R F .
where R is the typical curvature of the network. In this attractor, the velocity of the wall increases with time as v(t) 8.7 1 g fr P r where we used eq. D1 and R = vt. The friction domination regime ends when F/σ H or, equivalently, when the wall reaches the relativistic limit. Interestingly, when the background is radiation dominated, this always happen at a temperature T fr close to the time of domain wall domination T fr 0.43 √ g fr P r T dw-dom . Another consequence of thermal friction is the delay of the network annihilation. In the absence of friction, the wall annihilates when the bias pressure overcomes the domain wall tension, V b ρ dw . However, because R is now smaller due to friction, ρ dw is enhanced and so it will take more time for the bias pressure to dominate and cause the collapse of the network. For a temperature independent bias in a radiation dominated era, that happens when H fr ann = 1.2 g * ,ann g fr P r in gravitational waves at the time of annihilation as [95] P fr GW = G where we used I ∼ σR 4 for the quadrupole moment of the wall. There are three main changes compared to the frictionless case: i) the amount of gravitational waves is suppressed by v 3 RH 2 , ii) the network annihilates later thus increasing the relative energy in the GWs (when compared to the background), iii) the peak frequency of the spectrum at the time of annihilation is expected to move from H fr ann to 1/R fr ann . What remains to assess is the value of the reflection probability P r .
To do that we need to specify the interactions of the axion with the plasma. At temperatures above 1 GeV, the focus of this work, there is one unavoidable interaction due to the axion coupling to gluons. Other interactions are also possible, for example with SM fermions, which can be sizeable in DFSZ-like constructions, or with particles in the dark sector (both in the bias or in Λ H ), but those are model dependent.
Here we will take a conservative approach and assume the reflection probability to be maximal, P r = 1. We will present a more detailed analysis of friction due to gluons in [89]. In Figs. 9 and 10 we show the region of parameter space where friction is important, assuming g fr = 1, and the respective prediction for the GW spectrum. In order to perform a fair comparison with the frictionless case, we also included the numerical factors c, d, β in the expressions for friction. The results show that friction is only relevant in a small region of parameter space which is bounded by the condition that the network annihilates during the friction domination regime T fr < T fr ann and by the requirement that T fr ann < m a , otherwise the reflection probability will be suppressed. Even in the corner where friction is relevant, the predictions are only mildly different from those without friction (c.f. Figs.1 and 2 of the main text).