Stochastic gravitational-wave background from metastable cosmic strings

A metastable cosmic-string network is a generic consequence of many grand unified theories (GUTs) when combined with cosmic inflation. Metastable cosmic strings are not topologically stable, but decay on cosmic time scales due to pair production of GUT monopoles. This leads to a network consisting of metastable long strings on superhorizon scales as well as of string loops and segments on subhorizon scales. We compute for the first time the complete stochastic gravitational-wave background (SGWB) arising from all these network constituents, including several technical improvements to both the derivation of the loop and segment contributions. We find that the gravitational waves emitted by string loops provide the main contribution to the gravitational-wave spectrum in the relevant parameter space. The resulting spectrum is consistent with the tentative signal observed by the NANOGrav and Parkes pulsar timing collaborations for a string tension of G μ ∼ 10-11…-7 and has ample discovery space for ground- and space-based detectors. For GUT-scale string tensions, G μ ∼ 10-8…-7, metastable strings predict a SGWB in the LIGO-Virgo-KAGRA band that could be discovered in the near future.


Introduction
The formation of cosmic defects is a generic feature of cosmological phase transitions [1].Defects such as monopoles and domain walls can easily overclose the universe and must therefore be avoided.Cosmic strings, on the other hand, evolve towards a scaling regime in which their relative contribution to the total energy density of the Universe remains constant.Cosmic strings have characteristic signatures in gravitational lensing, the cosmic microwave background (CMB), and the stochastic gravitational-wave background (SGWB) and are therefore a potentially very interesting messenger from the early universe (for reviews and references, see, e.g., Refs.[2,3]).
During the past two decades, much progress has been made to describe the time evolution of a cosmic-string network (for a recent review and references, see [4]).Since after an initial transient period, the characteristic width of cosmic strings is much smaller than the horizon, cosmic strings are often described by the Nambu-Goto (NG) action.The cosmic-string network consists of "long" superhorizon strings and "short" subhorizon loops, which are formed in intercommutation events of long strings and which decay slowly by emitting gravitational radiation.The approach to the scaling regime can be understood analytically in the velocity-dependent one-scale (VOS) model [5,6], and it has also been established by large simulations of NG string networks [7,8].
Despite two decades of research, predictions of the SGWB signal from cosmic strings still have large uncertainties.In a cosmic-string network, gravitational waves (GWs) are primarily emitted by oscillating string loops as well as in the form of GW bursts emitted by sharp features propagating on string loops, so-called cusps and kinks [9,10].To compute the SGWB signal, one has to know the number density of non-self-interacting loops per unit string length, which can be determined by NG string simulations [7,8], as well as the average power radiated off in GWs by each loop.Following essentially the same strategy, different groups have nevertheless obtained significantly different results [4,[11][12][13].The main differences concern the number density of small loops and the treatment of gravitational backreaction, which can smooth out string singularities.Moreover, the entire picture of NG string loops decaying by gravitational radiation has been challenged by fieldtheoretic simulations suggesting a much faster decay of the network whose origin, however, remains mysterious [4,14].In this paper, we follow the approach in Ref. [13], relying on the evidence for long-string dominance in recent large simulations [12] and assuming suppression of GW radiation from kinks after gravitational backreaction.
We consider cosmic strings associated with the spontaneous breaking of a local U(1) symmetry embedded in a grand unified theory (GUT) [2,15,16], a prominent example being the breaking of B −L, the difference of baryon and lepton number [17].GUT-scale strings have a string tension in the range Gµ 10 −8 . . . 10 −6 , which seems excluded by the SGWB bound set by pulsar timing array (PTA) experiments [18][19][20], which constrain topologically stable cosmic strings to Gµ < 1.5× 10 −11 [21] for a standard loop size parameter α ∼ 0.1 (see below).However, it was recently pointed out that this bound can be avoided for metastable cosmic strings, which opens a new window for a SGWB signal close to the current upper limit in the LIGO-Virgo-KAGRA (LVK) frequency band that is consistent with the PTA bounds [22].Metastable cosmic strings decay by quantum tunneling into string segments connecting monopole-antimonopole pairs.In the semiclassical approximation, the decay rate per string unit length is given by [23][24][25][26] where m is the monopole mass and µ is the string tension.This suppresses the GW spectrum at low frequencies, rendering large string tensions compatible with PTA bounds.This opens up a new window to explore GUT-scale physics with gravitational waves [22,[27][28][29][30][31][32], which has received considerable attention since the recent report by the NANOGrav collaboration of evidence for a stochastic common-spectrum process at nanohertz frequencies [33], which has been interpreted as a SGWB in a large number of recent papers.Beyond the astrophysical interpretation in terms of supermassive black-hole binaries [34], possible cosmological interpretations include stable [31,35,36] as well as metastable strings [28].In a first calculation, string tensions in the range 10 −10 Gµ 10 −6 were shown to be consistent with the NANOGrav data, with a monopole mass to string tension ratio of around √ κ 8 [28].Such values can indeed be obtained in typical GUT models [37].Alternatively, quasi-stable strings with √ κ 1 associated with intermediate scales, which occur in symmetry breaking chains of GUT models [38], may provide a fit to the NANOGrav data while predicting a suppressed SGWB contribution at LIGO scales.Note that, also independent of grand unification, a SGWB signal from a cosmic-string network is a wellmotivated signature for physics beyond the Standard Model (SM) [39][40][41].
In this paper, we will analyze the GW spectrum produced by a metastable string network in detail.For metastable strings, both long superhorizon strings and short subhorizon loops decay into string segments.There are then two qualitatively different possibilities.In the first case, all the monopole magnetic flux is confined to the string, whereas in the second case, only some of the flux is confined, while the remaining flux is unconfined (for a review and references, see, for example, [42]).The pattern of GW radiation is very different in the two cases.In the first one, both loops and segments radiate GWs, with loops typically yielding the dominant contribution.In the second one, only loops radiate GWs, whereas segments loose energy much more rapidly by radiating gauge quanta corresponding to the unconfined flux [43,44].Symmetry breaking in GUTs generically leads to monopoles with partially unconfined flux (for recent examples, see [27,37,45]).
The GW spectrum of oscillating string segments connecting a monopole-antimonopole pair was first calculated by Martin and Vilenkin in a straight-string approximation [46].Subsequently, Leblond, Shlaer, and Siemens computed the GW spectrum from bursts and the SGWB of metastable strings based on string segments [25].A crucial point in this analysis is the matching of an early scaling regime, where metastable strings behave like stable strings, to a decay regime, where new loops are no longer produced, at a time t s = 1/Γ 1/2 d , where the scaling regime ends.In our work, we follow this approach, including also segments from decaying loops in the analysis.Our analysis of the GW contribution emitted by loops largely follows the analysis in Ref. [22], adding a refined treatment of the time scale of the loop decay.
The paper is organized as follows.In Sec. 2, we outline the derivation of the GW spectrum emitted by cosmic-string loops and segments, providing simple analytical expressions for all relevant quantities for GW emission in the radiation era.The more technical components of this analysis are deferred to the appendix.In Sec. 3, we perform a numerical evaluation of the GW spectrum, demonstrating the detection prospects for PTAs and ground-based GW detectors.Sec. 4 contains our conclusions.

Gravitational waves from string loops and segments
The time evolution of a network of stable cosmic strings emitting GWs has been extensively studied.After an initial transient period, the network reaches a "scaling regime", where the number densities of long superhorizon strings and subhorizon loops per Hubble volume are preserved.The long strings loose their energy mostly by loop formation, whereas the oscillating loops loose their energy by gravitational radiation.This picture is strongly supported by analytical models as well as large numerical simulations.
Much less is known for metastable cosmic-string networks characterized by a decay time It is expected that they behave similarly to stable networks at early times t < t s .However, they start decaying immediately after the phase transition in the course of which they are initially formed; and at t s , typical string segments enter the horizon.These segments then start oscillating under their own tension, radiating off GWs, and loop production stops.At present, numerical simulations for metastable string networks do not exist.We therefore follow the approach of Leblond, Shlaer, and Siemens [25] and match at the time t s , which marks the end of the scaling regime, two different periods to each other: Loops, short subhorizon segments from long strings and loops Loops and segments The second period ends at a time t e , which marks the end of GW emission, that is, the time when segments and loops have emitted all their energy in GWs and the string network disappears.In the following, we shall briefly recall the main ingredients in the computation of the GW signal from loops and segments.We shall also give explicit formulas for all relevant quantities in the radiation era, which is where most of the GW emission occurs for a large part of the parameter space.For a detailed discussion of the matter era, we refer the reader to the appendix.
Let us first consider stable cosmic strings, where the approach to scaling can be understood analytically in the velocity-dependent one-scale (VOS) model [5,6].The name of the VOS model derives from the fact that it assumes several length scales in the network to coincide up to constant factors: the inter-string separation L = (µ/ρ ∞ ) 1/2 , the string correlation length, and the string curvature radius.Here, µ is the energy density per unit of string length and ρ ∞ is the energy density of the long strings.In addition, the VOS model assumes the initial loop size at the time of formation to be controlled by the universal scale L ∝ H −1 ∝ t times a constant factor.
The energy density of long strings, ρ ∞ , is diluted in consequence of the Hubble expansion and depleted by loop production [2], where H = ȧ/a is the Hubble rate and v = v 2 ∞ 1/2 denotes the root-mean-square velocity of long strings, which tends to a constant value in the scaling regime and which affects the redshift behavior of the network.f ( , t) is the loop production function, which gives the number density of non-self-interacting loops produced per unit time and unit string length.The model accounts for scaling as a fixed point of differential equations for L (t) and v (t); and in the radiation era, one finds L (t) → ξt, ξ = 0.271 and v (t) → 0.662.Using Eq. ( 2), the loop production function then takes the form [6], where the constant c parametrizes the efficiency of loop-chopping.As discussed in the appendix, this loop production function yields the loop number density during radiation domination.Here, ΓGµ 2 is the total power radiated by a loop, G is Newton's constant, and A = 0.54 is obtained from a fit to numerical simulations.An important quantity is the ratio of the energy densities in loops and long strings.Using ρ ∞ = µ/L 2 → µ/ (tξ) 2 , Eq. (4) yields for α ΓGµ, which increases with decreasing radiated power Γ as well as with decreasing string tension µ.
Large numerical simulations of cosmic-string networks have led to the number density of the Blanco-Pillado-Olum-Shlaer (BOS) model [12], which is similar to Eq. ( 4).Compared to the VOS model, also the parameter α is determined as α 0.1.The BOS number density describing the loop population in the radiation era is given by [4] • n ( , t) = B t 3/2 ( + ΓGµt) 5/2 Θ (αt − ) Θ (t eq − t) with B 0.18 and t eq denoting the time of matter-radiation equality.A more general discussion of loop number densities, accounting also for the loop population after matter-radiation equality, as well as corresponding references can be found in Ref. [4].Number densities for stable as well as decaying loops and string segments satisfy kinetic equations.Their general form and solutions are described in the appendix.The number densities for decaying loops and segments can be obtained by matching the early scaling regime to a decay regime at t s = 1/Γ 1/2 d .This procedure leads to the result for decaying loops in Eq. (A.13), where denotes the length of a loop at time t s that evolves to the length at time t due to the emission of GWs.The number density differs from the one for stable loops by two damping terms, which become effective for t > 1/Γ d and t > √ 2 (ΓGµ) −1/2 t s ≡ t e , respectively.The first Heaviside theta function reflects the fact that only loops produced before t s contribute to the number density; the second theta function indicates that this expression is only valid during the radiation era.
The number density of segments, ñ, can be obtained in a similar way.It receives contributions from long strings decaying into segments as well as from loops decaying into segments, leading to the kinetic equation given in Eq. (A.29).The former has the analytic solution [25], where Γ parametrizes the GW emission of the segments and correspondingly te = √ 2 ΓGµ −1/2 t s is the time after which segments have disappeared.Since in the scaling regime, t < t s , superhorizon segments behave like stable long strings, the normalization factor C can be determined by matching the energy density of segments to the energy density ρ ∞ at time which yields C = 4/ξ 2 .For the contribution to segments from decaying loops, ñ(l) , the kinetic equation obtains the form of a partial integro-differential equation.We provide an exact analytical solution to this equation in terms of an infinite series in the appendix.We, moreover, find that this exact result can be reproduced to good approximation, as far as the consequences for the GW spectrum are concerned, by multiplying the first term in the series by an overall numerical ("fudge") factor σ, where σ 5.The total segment number density in the radiation era for t > t s thus reads Here and below, we follow Ref. [25] and set Γ Γ 50, for simplicity.The corresponding number densities in the matter era, which enter our numerical results in Sec. 3, are derived in the appendix.
In view of their contributions to GWs, it is interesting to compare the energy densities at t e .From Eqs. ( 7) and ( 12), one obtains For typical numbers, B = 0.18, σ = 5, ξ = 0.271, Γ = 50, α = 0.1, and large tensions, Gµ ∼ 10 −7 , all contributions are roughly of similar size.For smaller string tensions, the energy density stored in segments sourced by cosmic-string loops dominates over the other two contributions at t e .

Stable loops
We first consider stable loops and the corresponding GW spectrum from a network in the scaling regime following the discussion in Ref. [13].The GW energy density relative to the critical density per logarithmic frequency unit at cosmic time t is given by Here, H(t) is the Hubble rate and ρ gw (t, f ) is the energy density in GWs per frequency unit, which is obtained from the redshifted power density in GWs integrated from some initial time t i to the time of observation t, with f = (1 + z (t )) f ≡ (1 + z ) f .Loops of length oscillating in their kth harmonic emit GWs with frequency f = 2k/ .Hence, the power density per unit frequency is related to a loop number density • n ( , t ) and the power Gµ 2 P k per unit length as4 In the following, we focus on the contribution from cusps, which corresponds to Integrating from t i to t and changing variables, dt = −dz / (H (z ) (1 + z )), Eq. ( 15) yields with Note that the z integral depends only very weakly on the upper limit z i .
The loop density is diluted by the expansion of the universe and sourced by the interactions of long strings, which is encoded in the loop production function f ( , t), where a (t) is the scale factor and where the function ¯ has been introduced in Eq. ( 8), ¯ (t ) = + ΓGµ (t − t ).In the VOS model, and approximately also in the BOS model, strings are formed with a fixed fraction α of the horizon, ¯ (t ) ≈ αt .For a particular choice of and t, the formation time t is fixed by Eq. ( 8), which then determines the loop number density (20) via the value of the scale factor at the formation time, yielding Eqs. ( 4) and ( 6), respectively.
In the following, we discuss the resulting GW spectrum observed today, at t = t 0 , resulting from loops emitting GWs during the radiation era.The computation of the GW contribution emitted during the matter era is fully analogous and can be performed numerically by inserting the corresponding loop and segment number densities derived in the appendix into Eq.(19).For pedagogical reasons, we focus on the radiation epoch in this section; however, the numerical results shown in Sec. 3 will contain also the contributions from the matter era.
During radiation domination, the Hubble rate and cosmic time are given by with h 2 Ω r = 4.15 × 10 −5 , h = 0.68 [47], where the increase of the effective number of degrees of freedom with z has been neglected.Using the BOS loop number density in Eq. ( 6), this yields for the coefficient functions defined in Eq. ( 19), where z eq occurs as lower integration limit in Eq. ( 19) for GWs produced in the radiation era.
From this expression, one reads off the main qualitative features of the GW spectrum today.For 5 only the first of the two terms in the square bracket contributes, and the behavior of this term as a function of frequency depends on whether f is larger or smaller than f eq ≡ f p (z eq ).For small frequencies, i.e., f < f eq , the GW spectrum increases as Ω gw ∝ f 3/2 , as one reads off from Eqs. ( 14) and (22).Starting at around f eq , the GW spectrum then begins to approach a flat plateau, with the turnover frequency f eq being inversely proportional to Gµ. 6 In between f eq and k max f eq , the behavior of the spectrum is controlled by the sum over the harmonic oscillation modes, where k eq is the largest integer that is smaller than f /f eq and H is the kth harmonic number of order 4/3.For f eq < f < k max f eq , the deviation of the spectrum from a flat plateau therefore decays like f −1/3 .Once the top of the plateau is reached around f k max f eq , we have which, together with Eqs. ( 14), (18), and ( 22), yields the familiar result [4], The summation over all modes also plays an important role at high frequencies.For f high < f < k max f high , all modes between k high and k max contribute with a flat plateau to the GW spectrum, where k high is the smallest integer that is larger than f /f high , while all modes n < f /f high decay like f −1 (see below).The sum over k can then be expressed in terms of the Hurwitz zeta function [48], indicating that the spectrum decreases like f −1/3 in the interval f high < f < k max f high .At even higher frequencies, the two terms in the square bracket of Eq. ( 22) are always of similar size, for all values of k.Expanding C k at large f then shows that the spectrum only receives contributions falling off like f −1 , resulting in a total spectrum summed over all modes also falling off like f −1 .

Decaying loops
Loops are produced in the early scaling regime, t < t s .At later times, the produced loops shrink by emitting GWs.Their number density • n > satisfies a kinetic equation discussed in the appendix.Matching , t s , one obtains the number density in Eq. (A.13), which corresponds to the number density in Eq. (7).
Knowing the loop number densities, we can evaluate the coefficient functions C k (t, f ), see Eq. (19), which determine the GW spectrum.For GWs generated in the radiation era, one obtains The exponential factor yields frequency-independent and frequency-dependent cutoffs z e and z f , so that the range z < z m = max {z e , z f } does not contribute7 to the integral for C k (t 0 , f ).In this section, we will for simplicity focus on the regime z e > z eq , ensuring that the GW production is limited to the radiation-dominated regime.This is the case for [see Eq. ( 1)], κ 81 + 0.32 ln (Gµ) .
Cutting off the first integral in Eq. ( 29) at z m , the coefficients C k are approximately given by Compared to Eq. ( 22) for stable loops, the redshift z eq has been replaced by the cutoff z m .Above the frequency 2 f low , with one has z m = z e , and for f > f low , the first term in the square brackets in Eq. ( 32) approaches a constant.Hence, the GW spectrum today approximately features a plateau for f low < f < f high .On the other hand, at smaller frequencies, one has z m = z f > z e , and consequently the GW spectrum falls off as f 2 for f < f low .For frequencies above k max f high , the GW spectrum falls off again like 1/f , where for reference Eq. ( 23) can be expressed as with a redshift z i ∼ 10 23 corresponding to a reheating temperature of T rh ∼ 10 10 GeV.[20] and the LIGO/VIRGO collaboration [49], the lighter shaded regions show the prospective reach of SKA [50], LISA [51], LIGO and the Einstein Telescope (ET) [52].The orange shaded region indicates the region preferred by the NANOGrav hint [33].For simplicity, we fix the number of SM degrees of freedom to its high-temperature value in this figure, g * = 106.75.
Note that the decaying loops have to be created before t s , i.e., the argument of the first theta function in the number density in Eq. ( 7) has to be positive.This implies For the parameters given in Eq. ( 33), one finds k Θ (f ) ∼ 100 f /f low .To obtain the GW spectrum, one has to sum over all modes.We focus on the contribution from cusps given in Eq. ( 17).From Eqs. ( 14), ( 18), (32), and (35), one then obtains For large enough frequencies, i.e., k Θ (f ) 1, and f > f low , the GW spectrum reaches a plateau as for stable loops, where now k Θ (f ) states contribute to the total power Γ.
The final GW spectrum sourced by loops decaying during radiation domination is shown in the left panel of Fig. 1 for Gµ = 10 −11 . . . 10 −7 with √ κ = 8 (solid) and √ κ = 7 (dashed).These results are obtained by inserting the loop number density in Eq. ( 7) into Eq.( 19), leading to the GW energy density in Eq. (18).It suffices to take into account the loop number density at t > t s , since the GW spectrum is sourced largely at t ∼ t e t s .For easier comparison with the analytical results, we have fixed the number of SM degrees of freedom to its high-temperature value in this figure, g * = 106.75.The resulting agreement with the analytical expressions is very good.

Decaying segments
At early times, t < t s , the superhorizon strings decay and loose energy by chopping off loops.As discussed in the appendix, the segment density ñ(s) < sourced by long strings satisfies the kinetic equation [25] where u ( , t) = 3H (t) − 2 /t, and where the decay of segments acts as a source term for smaller segments.The rate for producing a segment with length between and + d in the decay of a segment of length is Γ d d .The segment with length can be chopped off at either side, hence the factor of 2. In the case = 2 , one breaking produces two segments of length .The solution exhibits the expected scaling behaviour ρ cs (t) ∼ µ d ñ(s) < ( , t) ∼ µ/t 2 .The normalization constant C = 4/ξ 2 is determined by the scaling solution.
At t = t s , typical segments enter the horizon, loop production terminates and GW radiation begins.Now the relevant kinetic equation reads and the solution of this integro-differential equation satisfying the initial condition ñ(s) < ¯ (t s ) , t s is given by Eq. ( 9).In addition, the segment number density obtains a second contribution, ñ(l) , from loop decays.The exact solution of the full kinetic equation (A.29) is given in the appendix.To good approximation their contribution is described by Eq. ( 11), so that the full segment number density is then given by the sum of both contributions, see Eq. (12).Given the number density of decaying segments, we can evaluate the GW spectrum as in the previous section.Using Eq. ( 19) but replacing the loop number density with the segment number densities given in Eqs. ( 9) and (11), one obtains for the coefficient functions Now the exponential yields the frequency-independent cutoff ze , and the range z < ze does not contribute to the integral for Ck (t, f ).
The power in mode k of the oscillating segment is quasi-constant [46], up to a very large maximum value determined by the Lorentz factor of the oscillating monopole, k max ∼ γ 2 0 , beyond which P k decreases like 1/k 2 .The total power is given by and Γ ∼ 50 would imply k max ∼ 10 5 .To obtain the GW spectrum, one has to sum over all modes and integrate over z .For ñ(s) > , this sum extends to k max , whereas for ñ(l) > introduced above, it extends to k Θ (f ) given in Eq. (35).Because of the quasi-constant behaviour of k P k , it is convenient to perform the summation over k first, approximated as an integral.From Eqs. ( 14) and ( 18) one finds for the GW spectrum, From Eq. ( 45), the qualitative features of the GW spectrum from segments are easily understood.For small frequencies, f < f low , the first term in the bracket behaves exactly as the contribution from decaying loops in Eq. ( 36), i.e., the spectrum increases as Ω gw ∝ f 2 , and above f low it approaches a plateau.However, contrary to decaying loops, the spectrum from loop segments is cut off as k Θ (f ) approaches k max around fhigh = k max f low .Above fhigh , the constant parts of the first and second terms cancel, and the spectrum falls off as Ω gw ∝ 1/f .A similar cancellation takes place around fhigh between the third and fourth terms arising from long-string segments, where we have again assumed Γ = Γ, i.e. ze = z e .For f < f low , one has z m = z f .The z integral in Eq. ( 45) yields a factor (1 + z s ) / (1 + z f ) 5 , which implies Ω gw ∝ f 5/3 .Finally, one finds for the height of the plateau, Comparing the result with Eq. ( 37), one observes that, for Γ 50 and σ 5, the GW spectrum from loop segments is suppressed by about one order of magnitude compared to the GW spectrum from decaying loops.The contributions from loop segments and long-string segments are about equal at Gµ = 10 −8 .Hence, the total contribution from segments to the GW spectrum is always subdominant with respect to the one from decaying loops.A detailed comparison of the GW spectrum from decaying loops, segments from loop decays, and segments from long-string segments is given in Fig. 5 in the appendix.The final GW spectrum sourced by segments decaying during radiation domination is shown in the right panel of Fig. 1 for Gµ = 10 −11 . . . 10 −7 with √ κ = 8 (solid) and √ κ = 7 (dashed).Again we find very good agreement with our analytical expression for the GW spectrum.
The extension of the plateau from the segment contribution depends on the number of contributing modes k max , which is determined by the Lorentz factor γ 2 0 ∼ µ 2 2 /m 2 = M 2 Pl Gµ/κ.The calculation in Ref. [46] obtained γ 2 0 ∼ 300, corresponding to Γ ∼ 25.In Ref. [25], Γ ∼ 50 has been used, corresponding to k max ∼ γ 2 0 ∼ 10 5 .For the parameters given above, extending the plateau up to the LVK band around 100 Hz would require k max ∼ 10 10 with Γ ∼ 100.The calculation of the radiated power has thus far been restricted to a straight string.The extension to realistic configurations, where the two monopoles can pass each other without forming a black hole, remains a problem for future research. 8It is expected that, above some critical mode number k c , the radiated power P k falls off exponentially [46] .

Detection prospects
Our final results for the GW spectrum produced by metastable cosmic strings are shown in Figs. 2  to 4. Starting from the analytical expressions for the loop number density in Eqs. ( 7) and (A.16) as well as (for Fig. 4) the segment number density in Eqs. ( 12), (A.28), (A.39), and (A.40), we numerically compute the coefficient functions C k (see Eq. ( 19)) and the resulting GW spectrum (see Eqs. ( 14) and ( 18)).Here, we include the changes in the number of degrees of freedom in the SM thermal plasma, leading to deviations from the analytical prediction of a perfectly flat spectrum, particularly visible in the left panel of Fig. 2. The summation over higher harmonics k is performed up to the maximum relevant mode, as discussed in the previous section.We vary the string tension from Gµ = 10 −11 to 10 −7 , covering the entire range of interest for the existing pulsar timing and ground-based interferometer experiments.We consider values of √ κ all the way from quasi-stable cosmic strings, which have a life time comparable to the age of the Universe, √ κ 9, to metastable cosmic strings with a strongly suppressed spectrum in the pulsar timing array band, √ κ 7.5.For this entire range of κ, the cosmic strings are stable enough to give a large signal in the LVK band.We contrast these predictions with existing bounds, the expected sensitivity 8 The GW spectrum in Eq. ( 45) can be expressed as an integral over the redshift z and the segment length instead of the mode number k [25].The integral over has the lower bound min ∼ 1/((1 + z )f ) and the upper bound max ∼ ts, i.e., the horizon at the beginning of the "short-string period".Moreover, the frequency satisfies the upper bound Pl Gµ/κ.The integral over z is dominated by the contribution close to the lower limit z ∼ ze.From these inequalities, one obtains lower and upper bounds on f , which for large κ read: log 10 (f low ) ∼ βκ, β = − log 10 (e) π/4 −0.34 and log 10 (f high ) ∼ γκ, γ = 3 log 10 (e) π/4 1.02.These straight lines represent the boundaries of the plateau in figure 7 of Ref. [25].At κ ∼ 60, the trans-Planckian upper bound is f high ∼ 10 70 Hz, corresponding to kmax ∼ 10 78 and Γ ∼ 800. Figure 2: GW spectrum from metastable cosmic strings for monopoles with unconfined fluxes.The experimental constraints depicted in the left panel are as in Fig. 1.The black dotted curves indicated the spectra obtained for topologically stable cosmic strings for the corresponding value of Gµ.In the right panel, we show predictions for the frequency range of pulsar timing arrays, together with the bound from the Parkes Pulsar Timing Array (PPTA) [20] published in 2015 (gray) and the more recently reported preferred regions of NANOGrav [33] (orange) and PPTA [53] (black).These include contributions from the loops decaying during the matter era.
of upcoming experiments, as well as with the tentative GW signal reported by the NANOGrav collaboration [33].

Monopoles with unconfined fluxes
If the monopoles feature unconfined fluxes, any cosmic string segments formed from long strings or loops will rapidly decay radiating massless gauge bosons as the monopoles and antimonopoles oscillate and finally annihilate.In this case, the resulting GW spectrum is dominated by the GW emission from cosmic string loops.The left panel of Fig. 2 shows the resulting GW spectrum for different values of the model parameters Gµ and κ.This extends the result shown in the left panel of Fig. 1 by including the change of degrees of freedom in the SM thermal bath as well as the GW emission during the matter domination era. 9 The latter is relevant only for t e > t eq and leads to an enhancement at low frequencies for quasi-stable strings (indicated by the dotted black curves).We recall that the origin of this enhancement can be traced back to the scaling behaviour.During radiation domination, the loop production and subsequent GW emission have to be efficient enough to compensate the T 4 decrease of the energy in the SM thermal bath.During matter domination, scaling dictates a reduced GW emission.In this sense, loops generated during the radiation era but surviving until the matter era radiate a disproportionate amount of energy, leading to an enhancement of the GW spectrum.
In the right panel of Fig. 2, we perform a more detailed comparison with the existing pulsar timing results for the case of monopoles with unconfined fluxes.Parameterizing the GW power spectrum as Ω gw = Ω gw (f PTA )•(f /f PTA ) nt , we determine the amplitude and tilt of this power law in 9 Both of these affect the cosmological expansion history, which is encoded in the Hubble parameter H(z) = s (z)), Ωm = 0.308, ΩΛ = 0.702, and g * (z) (gs(z)) denoting the effective number of degrees of freedom relevant for the energy (entropy) density of the SM thermal bath.
the range of [2 . . .4] nHz around the peak sensitivity (10 yr) −1 3 nHz of current PTA experiments.Note that, once depicted at this frequency instead of the more conventional reference frequency of 32 nHz = 1/year, it becomes more transparent that the NANOGrav measurement is essentially a measurement of the amplitude with the tilt still subject to a large uncertainty and largely uncorrelated with the amplitude.We show the predictions for metastable cosmic strings for the entire range of model parameters considered, with the solid lines denoting contours of constant Gµ and the dotted lines indicating contours of constant κ (with the maximal value, √ κ = 9.2, indicated by a white dot).The orange region indicates the region suggested by interpreting the recent NANOGrav result as a GW signal [33], with the solid (dashed) contours showing the 68 % and 95 % regions reported by NANOGrav when performing a fit to the first five frequency bins (when performing a fit with a broken power law).The black solid lines show the preferred region reported by PPTA when performing a similar analysis based on the first five frequency bins [53].This region is in tension with previous results from the PPTA [20] and NANOGrav [18] collaborations, see Refs.[33,53,54] for a discussion.
The analysis presented here updates and largely justifies the simpler analyses performed in Refs.[22,28].The main difference is the inclusion of the decay of the cosmic string network, now encoded in the exponential function in Eq. (7).Contrary to the previous analysis, this allows for loops of different length to decay at different times, as expected from the different probability of forming a monopole pair along the loop.The main consequence of this is that the fall-off of the GW spectrum at small frequencies is described by an f 2 power law instead of f 3/2 as found in Ref. [22].The resulting overall shift to larger values of n t for small κ mildly reduces the overlap with the preferred NANOGrav region, implying in particular that there is now barely any overlap with the 1σ NANOGrav region, while within 2σ there is good agreement.We note, however, that the very recent results reported by the PPTA collaboration prefer a larger spectral tilt, yielding significantly better agreement with our predictions.In Ref. [28], the range of string tensions Gµ = 10 −10 . . . 10 −6 was considered.Given our calculation of the GW spectrum, the LIGO O3 upper bound on Ω gw [49] restricts Gµ to values below 2 × 10 −7 (see left panel of Fig. 3).In this paper, we therefore focus on the range of string tensions Gµ = 10 −11 . . . 10 −7 .
We stress that the current significant uncertainties both in the interpretation of the PTA data as a GW signal and in the modelling of the cosmic-string network force us to take any model-todata comparison with some grain of salt.It is nevertheless instructive to contrast the tentative NANOGrav signal with other existing and upcoming GW observations, in particular by LIGO, see Fig. 3.The left panel shows the NANOGrav signal (95 % C. L. region of the broken power law fit) in the model parameter plane, together with the PPTA exclusion limit [20], the LIGO O3 bound on stochastic backgrounds [49], 10 and the design sensitivities of LISA and LIGO.Note that 10 In Ref. [55], upper bounds on the string tension were derived for different NG models.In a model (A) for the loop number density [12], the obtained GW spectrum essentially shows a plateau between the nHz and the LVK band, similar to Ref. [13].This is not the case for the loop number density of model (B) [7], which leads to a GW spectrum [11] that differs from model (A) by up to four orders of magnitude.Correspondingly, the derived bounds on the string tension are very different.For model (A), the upper bound varies in the range Gµ 10 −8 . . . 10 −6 , whereas for model (B), the upper bound is Gµ (4.0 . . .6.3) × 10 −15 [55].A range of two orders of magnitude for model (A) appears reasonable in view of current theoretical uncertainties.The difference by seven orders of Figure 3: Parameter space and GW detection prospects for monopoles with unconfined fluxes.The orange region and the contour lines in dark gray show the tentative signal reported by NANOGrav [33] and PPTA [53], which lies within the exclusion region of the previous PTA bounds (grey region in left panel).The blue shaded regions indicate bounds and prospects reported by the LIGO/VIRGO collaboration [49].The entire parameter space shown in the left panel can be probed by LISA.The right panel focuses on two likely future GW observables: the tilt nt of the SGWB at PTA frequencies and the amplitude at LIGO frequencies.Metastable cosmic strings can explain a signal anywhere in this plain outside the grey shaded region.
Ref. [49] quotes both a more conservative bound, Ω GW < 1.7 × 10 −8 at 95 % C. L. (labeled "O3" in Fig. 3), and a more aggressive bound, Ω GW < 5.8 × 10 −9 (labeled "O3 log prior"), depending on the choice of priors.We conclude that, within the framework of metastable cosmic strings decaying through the production of monopoles with unconfined fluxes, the current NANOGrav data are compatible with 2 × 10 −11 Gµ 2 × 10 −7 , with current (and possibly upcoming) LIGO data pushing this to lower values, towards the regime of quasi-stable cosmic strings.We note, however, that large values of Gµ are more sensitive to the formation time of the cosmic-string network and / or the reheating temperature of the Universe, which in our analysis we have taken to be at very high redshift.Lowering this can suppress the GW spectrum at LIGO scales while leaving the predictions in the pulsar timing array band untouched.The entire parameter space compatible with the NANOGrav signal will be finally probed by LISA.
In the right panel of Fig. 3, we focus on the most likely observables of the near future: the tilt measured in pulsar timing arrays (horizontal axis) and the amplitude measured in groundbased interferometers (vertical axis).The entire white region can be reached by varying the model parameters Gµ and κ, with the orange region indicating the preferred NANOGrav region.On the contrary, a GW signal in the gray region could not be explained within this setup.

Monopoles with no unconfined fluxes
If on the contrary the monopoles do not feature any unconfined fluxes, the channel of energy loss for the cosmic-string segments is gravitational radiation.In this case, the total GW spectrum receives an additional contribution from decaying segments, which can originate both from long strings or from loops.The corresponding number densities are derived in the appendix and given magnitude between models (A) and (B) is largely due to the fact that the evidence for large-loop dominance from large numerical simulations [8] is not taken into account in model (B).by Eqs. ( 12), (A.28), (A.39), and (A.40).Fig. 4 shows the resulting GW spectrum as well as the predictions for pulsar timing arrays.Comparing with Fig. 2, we note that the GW contribution from cosmic-string segments gives at most a minor correction to the GW contribution from cosmic-string loops only.The two contributions become comparable only for large values of Gµ and f < f low .This in particular entails a flatter slope at frequencies below the onset of the plateau in radiation domination, noticeable in the right panel of Fig. 4 by the marginally more limited range of the tilt n t in the pulsar timing regime.Once κ becomes sufficiently large for the string network to survive into the matter era, the additional boost in the number density of the cosmic-string loops makes them completely dominate over the segment contribution.For a more detailed comparison of the different contributions, see Fig. 5 in the appendix.In summary, we conclude that the inclusion of the GW emission from cosmic-string segments is at most a minor correction, in particular given the overall uncertainties in modelling GW emission by a cosmic-string network.We note, however, that this conclusion is based on some model assumptions, in particular on the GW emission rate of loops and segments ( Γ Γ = 50) and on the loop size (α = 0.1), which may need to be revisited as our understanding of the dynamics of cosmic-string networks improves.

Other observables
The GW spectrum has to satisfy constraints imposed by big-bang nucleosynthesis (BBN) and the cosmic microwave background.During BBN, the expansion rate of the universe is tightly constrained, which limits the contribution of GWs to the energy density to the contribution of about one relativistic neutrino [56,57].With T BBN ∼ 0.05 MeV, one has z BBN ∼ 5 × 10 8 and H(z BBN ) ∼ 10 −3 Hz.The contribution of one relativistic neutrino to Ω today is ∆Ω 1ν 7/43 Ω r /(1 + z eq ) 5 × 10 −5 .The GW spectrum produced until t BBN has to be integrated over all subhorizon GWs present at the time of BBN.From Eq. ( 23), one obtains for the frequency f BBN = f p (z BBN ) 20 Hz (10 −7 /(Gµ)) H(z BBN ).Hence, all frequencies of GWs produced before t BBN fit into the BBN Hubble horizon.This conclusion also immediately follows from the fact that string loops and segments represent causal GW sources on subhorizon scales in a decel-erating expansion background.If the network decays after BBN, i.e. t e > t BBN , one obtains the contribution of GWs to Ω during BBN by integrating the plateau in Eq. ( 26) from f BBN to f high .If the network decays before t BBN , this integral yields an upper bound on the contribution of GWs.In this way, one obtains ∼ 10 −8 Gµ 10 −7 1/2 32 + ln T rh 10 10 GeV .
For the considered parameter values, ∆Ω BBN gw is smaller than ∆Ω 1ν by at least three orders of magnitude.
Precision measurements of the CMB constrain cosmic-string networks in several ways.Temperature anisotropies yield an upper bound on the tension of quasi-stable strings, Gµ < 10 −7 [58], which is the largest string tension that we consider.Other interesting observables are spectral distortions.Current bounds are not yet very stringent, but future experiments may indeed be able to probe metastable strings [59].In principle, also monopole annihilation from string segments could lead to interesting signatures [60,61], which requires further investigations.

Conclusions
The formation of cosmic strings that are not topologically stable is a rather common feature in GUT models [23,27,42].If the symmetry breaking step responsible for monopole production is separated from the symmetry breaking step generating cosmic strings by a phase of cosmic inflation, we generically obtain a network of metastable cosmic strings.Their decay is triggered by pair production of monopoles along the cores of the cosmic strings.This process is exponentially suppressed by the ratio of the monopole mass to the cosmic string tension, κ = m 2 /µ, leading to a cosmological life time.For mass ratios in the range of √ κ ∼ 7 . . .8, this leads to a strong suppression of the GW at low frequencies compared to the signal expected from topologically stable cosmic strings, while allowing for a large signal in the Hz regime.Consequently, a large scale-invariant SGWB at LIGO frequencies would be perfectly compatible with a null detection at pulsar timing arrays.This in particular demonstrates the significant discovery space for GWs from cosmic strings that ground-based interferometers are currently probing and which will be further significantly enlarged by the space-based interferometer LISA.For GUT-scale string tensions, Gµ ∼ 10 −8 . . . 10 −7 , metastable strings predict a SGWB in the LVK band that could be discovered in the very near future.
A significant theoretical distinction in the computation of the GW spectrum is the existence of unconfined fluxes in pair-produced GUT monopoles.Monopoles featuring unconfined fluxes are produced as monopole-antimonopole pairs, and the resulting cosmic-string segments decay rapidly under the emission of massless gauge bosons.On the other hand, if the monopoles do not feature any unconfined fluxes, the cosmic-string segments decay only due to GW emission, leading to an additional contribution to the GW spectrum.In the present paper, we computed for the first time all contributions in a systematic way, allowing us to compare the different contributions in both scenarios.In conclusion, we find that the GW contribution from cosmic-string loops is the dominant contribution in essentially the entire parameter space of interest, though if present, the contribution from cosmic-string segments has the potential to mildly influence the slope of the spectrum at PTA frequencies.
At the technical level, we improve the estimation of the GW spectrum from metastable cosmicstring loops first given in Ref. [22] to allow for cosmic-string loops of different size decaying at different times, which changes the estimate of the low-frequency slope from 3/2 to 2. For the contribution from metastable cosmic-string segments, our main finding with respect to Ref. [25] is that the frequency range of the plateau in the GW signal is limited in the ultraviolet by the number of modes contributing.We, moreover, provide analytical formulas for the loop and segment number densities in all relevant epochs of cosmic history that take into account recent progress in the modelling topologically stable cosmic-string networks [12].
Metastable cosmic strings provide a possible explanation for the tentative SGWB signal reported by the NANOGrav [33] and PPTA [53] collaborations.For 2 × 10 −11 < Gµ < 2 × 10 −7 and √ κ > 8, the metastable cosmic-string signal is compatible with the NANOGrav 2σ region and the PPTA 1σ region, respectively.Upcoming, more precise determinations of the spectral tilt of this signal will be decisive in distinguishing between this explanation and other astrophysical or cosmological sources.This, moreover, demonstrates the great potential of future multi-band GW observations involving PTAs, space-based, and ground-based interferometers.

Note added
Shortly after submitting this work to the arXiv, the PPTA collaboration [53] presented an analysis of their latest data set [19], reporting results in agreement with NANOGrav [33] when performing a similar analysis.For further tests and a discussion of possible interpretations of these results, see Ref. [53].We included these new results in Figs. 2 to 4 in the current version, demonstrating that they fit our predictions very well.
where S is a source term, u describes the change of the string length due to Hubble stretching and energy loss of the network, H is the Hubble parameter and Γ d is the decay rate.Eq. (A.1) can be derived by considering the changes in n ( , t) ∆ , i.e., the number density of loops whose lengths lie in the interval [ , + ∆ ], after some infinitesimally small time step from t to t + ∆t, n ( + u ∆t, t + ∆t) ∆ = S ( , t) ∆t ∆ + a (t) a (t + ∆t)

A.1 Loops
This subsection is dedicated to the derivation of the loop number density of metastable cosmic strings throughout different cosmological epochs.We begin with a detailed derivation of the loop number density during radiation domination, which underlies the analytical estimates presented in the main body of the text.We then proceed to include the matter era, which we include in our numerical computations of the GW spectra.We follow a similar procedure for the segment number density in Secs.A.2 and A.3.For a comparison of the various contributions both in radiation and matter domination, see Figs. 5 and 6.The left panel of Fig. 6, moreover, illustrates the relevant time scales.

Figure 1 :
Figure 1: GW spectrum from metastable cosmic-string loops (left) and segments (right) for different values of Gµ and κ.The dark gray-shaded regions indicate existing bounds from pulsar timing arrays[20] and the LIGO/VIRGO collaboration[49], the lighter shaded regions show the prospective reach of SKA[50], LISA[51], LIGO and the Einstein Telescope (ET)[52].The orange shaded region indicates the region preferred by the NANOGrav hint[33].For simplicity, we fix the number of SM degrees of freedom to its high-temperature value in this figure, g * = 106.75.

− 3 ,
where ζ m = 1 for z m = z e and ζ m = 4/3 for z m = z f , and zm = max (z e , z f ).

Figure 4 :
Figure 4: GW spectrum from metastable cosmic strings for monopoles with no unconfined fluxes, i.e., including the additional contribution from cosmic-string segments.Color coding as in Fig. 2.

3 n
( , t) ∆ − Γ d n ( , t) ∆t ∆ , (A.2) where ∆ = ∆ +∂ u ∆t ∆ accounts for the change in the interval length ∆ during ∆t.Expanding all terms in Eq. (A.2) up to first order in ∆t and collecting all terms of order ∆t ∆ on both sides reproduces Eq. (A.1).In standard form, the partial differential equation (A.1) reads[u ( , t) ∂ + w ( , t) ∂ t ] n ( , t) = F ( , t, n ( , t)) , w = 1 , F = S − (3H + Γ d + ∂ u) n , (A.3)which can be solved by integrating the ordinary differential equations for the three characteristic curves l( t ), t(t ) and n(t ) as functions of an auxiliary parameter t ∈ [t i , t], Imposing the boundary conditions ¯ (t) = , t (t) = t, and n ¯ (t i ) , t i = n i ¯ (t i ) , one obtainsn ( , t) = n ¯ (t) , t (t) = exp − t t i dt 3H ( t) + Γ d ¯ + ∂¯ u ¯ , t × n i ¯ (t i ) + t t i dt S ¯ , t exp t t i dt 3H ( t) + Γ d ¯ + ∂¯ u ¯ , t = exp − t t i dt Γ d ¯ t + ∂¯ u ¯ t , t t × a (t i ) a (t) 3 n i ¯ (t i ) + 3 S ¯ t ,t exp t t i dt Γ d ¯ t + ∂¯ u ¯ t , t t .(A.5)