Gravitational Waves from Global Cosmic Strings and Cosmic Archaeology

Global cosmic strings are predicted in many motivated extensions to the Standard Model of particle physics, with close connections to axion dark matter physics. Recent studies suggest that, although subdominant relative to Goldstone emission, gravitational wave (GW) signals from global strings can be detectable with current and planned GW detectors such as LIGO, LISA, DECIGO/BBO, ET/CE and AEDGE/AION, as well as pulsar timing arrays such as PPTA, NANOGrav and SKA. This work is an extensive, updated study on GWs from a global cosmic string network, taking into account of the most recent developments related to the subject. The main analysis is based on the analytical Velocity-dependent One-Scale (VOS) model calibrated with recent simulation results, which provides a generic protocol for such calculations with details given. We also demonstrate how the GW signal can be influenced with variations to the baseline model: this includes considering the uncertainties of model parameters and the potential deviation from the conventional VOS model prediction (i.e. the scaling behavior) as suggested by some of the recent simulation results. Furthermore, we investigated in detail the effect of a non-standard cosmology (e.g. early matter domination or kination) or new particle species on the GW signals from global strings. We demonstrate that the frequency spectrum of GW background from global cosmic strings can be used to probe the cosmic history prior to the Big Bang nucleosynthesis (BBN) (i.e. the primordial dark age) up to a temperature of $T\sim 10^8$ GeV.

Among the known cosmological sources of SGWB (see review [32]), cosmic strings stand out as one that can yield strong signals over a wide frequency range due to continuous emission throughout a long period of time. Cosmic strings are one-dimensional, topologically stable objects that are generically predicted by many theoretical extensions of the Standard Model of particle physics, e.g., field theories with a spontaneously broken U (1) symmetry (gauge or global) [33][34][35][36][37][38][39], and the fundamental and/or composite strings in superstring theory [40][41][42][43][44]. After formation, the strings quickly evolve towards a scaling regime where the string network consists of a few Hubble-length long strings per horizon volume, along with more copious loops formed by long string intersections. The loops then oscillate and radiate energy in the form of GWs and/or other particles until they decay away. Most literature on GW signatures from cosmic strings have been focused on those sourced by local strings or superstrings which typically can be described by Nambu-Goto (NG) action. In contrast, a global string network as a potential source of GWs has been largely ignored since by naive estimate GW radiation would be overwhelmed by Goldstone emission which occurs with a much larger rate. Very recently, inspired by its intimate connections to axion dark matter physics, significant progress has been made in simulating global topological defects and on the GW signals originated from it [45][46][47][48][49]. With a semianalytical approach based on the Velocity-dependent One-Scale (VOS) model, our earlier work [24] demonstrated that the GW signal from global strings, albeit notably smaller than that from its NG string counterpart, can be within reach of future GW experiments such as LISA [6,7], AEDGE [50], DECIGO and BBO [51]. Such a positive prospect of detection has been confirmed by simulation-based work [48,49], although details differ which will be addressed in this work.
The frequency spectrum of the SGWB from a cosmic string network can also serve as a powerful tool to probe the very early cosmic history that is not accessible by existing means. The ΛCDM cosmology was established based on precise measurements of electromagnetic radiation over different frequency ranges with a variety of experiments. A simple extrapolation of ΛCDM cosmology back in time suggests that the Universe is radiation dominated from the recent matter-radiation equality all the way back to the end of inflation. This paradigm is supported by observing cosmic microwave background (CMB), the relic photons that started free traveling when the radiation temperature was about 0.3 eV. The success of BBN theory in predicting primordial abundances of light elements also provides evidence for a radiation dominated era up to T ∼ 5 MeV. However, the hypothesis of radiation domination (RD) for epochs prior to BBN or at radiation temperature higher than ∼ 5 MeV is yet to be experimentally tested. On the other hand, possibilities of non-standard pre-BBN cosmologies are well motivated by many grounds, such as dark matter [52,53], axion physics [54,55], baryogenesis [56,57], non-minimal inflation/reheating [58,59], and string compactification [60,61]. In particular, recently there has been an increased interest in the impact of non-standard cosmology on dark matter physics [62][63][64]. The discovery of GWs leads to unprecedented opportunities to shed light on this mysterious pre-BBN primordial dark age [65][66][67]. GWs are the only cosmic messengers that can travel freely throughout space-time since the Big Bang. They carry unique information about the earliest phases of the Universe's evolution, beyond what can be assessed by observing EM radiations. Due to the continuous, potentially strong GW emissions from a string network throughout a long era of cosmic history, the SGWB frequency spectrum from cosmic strings is particularly appealing as a tool for looking back in time or cosmic archaeology [21][22][23]. The application of this idea in the context of NG strings was recently proposed and studied in [21,22], based on a frequency-time (temperature) correspondence. Cosmic archaeology with global string induced GWs was only briefly discussed in [24], which we will explore in great detail in this update.
In this work, we aim at an extensive study of SGWB signals originated from a global string network, and a comprehensive investigation into the potential new physics imprints in the pre-BBN Universe that can be detected with such a GW spectrum. Greater technical details are given, which may serve as a handy reference for future studies. Our primary approach is to use the analytic VOS model calibrated with simulation results (directly obtained for early times). Due to technical difficulties of simultaneously capturing physics at hierarchical scales, current simulations can only cover the evolution history of a global string network up to a few e-folds of Hubble expansion after the formation time. Thus, whether it is reliable to make a direct extrapolation of simulation results to late times (most relevant for observations today) requires further investigation. On the other hand, while VOS model for global strings are still being tested and needs to be calibrated with simulation data, the prediction for late times by the VOS model is obtained by solving the evolution equation incorporating the known physics effects instead of simple extrapolation. Therefore, such a semi-analytical approach is highly complementary to the simulation efforts and the two approaches can lead to insights to help improve each other. We significantly updated and expanded the related studies initiated in our earlier paper [24], taking into consideration the very recent developments since then. For instance, [26,30] show that the inclusion of the very high oscillation modes can drastically change the shape of the GW spectrum from NG strings in (early-)matter dominated era, which was neglected in earlier literature. We included the contribution from these high modes in this updated study, which leads to substantial modifications to the GW spectrum at low f for standard thermal history as well as at high f with the presence of an early matter domination epoch. We also discuss the consequence for the prediction of SGWB if the non-scaling behavior found in some simulation results for early evolution sustains in the late-time evolution of a global string network, compare with the results found in [49,68], and suggest potential modifications to the VOS model to accommodate such a feature. We will dive into the time-frequency correspondence for global strings, which is the guiding principle for testing standard cosmology. We conduct an extensive study on probing a potentially existing non-standard equation of state of the pre-BBN Universe such as early matter domination (EMD) or kination, where we also include a concrete example for a finite duration of a kination epoch. In addition, we study the effects on the GW spectrum with the presence of new massive degrees of freedom. Furthermore, a detailed discussion is given to address uncertainties such as loop size distribution, radiation parameters, and distinguishing from other SGWB sources. Our results directly apply to pure global strings associated with massless Goldstone. The application to the axion case where the Goldstone acquires a mass at a QCD(-like) phase transition is more complex and requires treatments of the axion domain walls in addition to the strings, see for example [69]. We reserve a dedicated study on the axion case for future work. We also comment on the prospect of addressing the recent NANOGrav result with global strings.
The rest of this article is organized as follows. In Section 2.1 we will present our methodology based on the analytical Velocity-dependent One-Scale (VOS) model for global strings calibrated with recent simulation results. In Section 2.2 we derive the GW frequency spectrum from a global string network in the context of standard thermal history. In Section 3 we illustrate the relation between the frequency of a GW signal observed today and its emission time in the early Universe. With several benchmark examples, we show how this relation can be used to test standard cosmology and detect potential new physics. Related experimental constraints and sensitivities are also demonstrated in Section 3 and 4. In Section 5 we will address various uncertainty factors that may affect the results, as well as how to distinguish global string induced SGWB from other potential SGWB sources. We make our conclusions in Section 6.
2 Evolution of a Global Cosmic String Network 2.1 Velocity-dependent One-Scale (VOS) model for global strings grow with the scale factor a due to the energy loss from the decay of the loops. While GWs constitute the leading radiation by the NG strings, they are irreducible but subdominant mode for global strings for which the emission of Goldstone particles is more important 1 . The energy density of the global string network (mainly stored in long strings) is where the dimensionless parameter ξ(t) is defined as the number of long strings per horizon volume. µ(t) is the time-dependent tension (i.e. energy per unit length) of the global strings (µ is a constant for NG or local strings), where δ is the width of the string core, λ is the coupling in φ 4 theory and m φ sets the mass of the Higgs-like complex φ whose VEV breaks the global U (1), and we have defined the time-dependent parameter N which will be used in later discussions. The temperature T dependent thermal mass contribution is negligible well after the symmetry breaking phase transition (T η), and thus we ignore it in our analysis. We consider λ ∼ 1 such that m φ and η are comparable. The evolution equation for the correlation length L is [36,76,77,81] which couples to the evolution equation for the average long string velocityv ∞ : where k v is the momentum parameter. While we will investigate the detectability of GW signal, we left out the GW radiation term in these evolution equations because its contribution here is sufficiently suppressed [36]. The terms on the RHS of Eq. 2.4 represent, in order, the dilution effect from the expansion of the Universe, thermal friction effect with characteristic scale f ∝ µT −3 , loop chopping rate parameterc, and the back-reaction due to Goldstone boson emission [77,81]. The thermal friction is negligible as the Universe cools down such that T η.
In the following analysis we consider various possibilities of background cosmology parametrized by n, defined as ρ ∝ a −n , a(t) ∝ t 2/n (2.6) 1 We neglect the emission of radial mode which is shown to decouple soon after the network formation [47,49] and may be generally suppressed when the loop size is larger than ∼ 1/η [71].
where a(t) is the expansion parameter as a function of time t, ρ is the background cosmic energy density. n = 3, 4 correspond to the cases of matter (MD) and radiation (RD) domination, respectively. We focus on the range of 2 < n ≤ 6 (n = 6 corresponds to kination epoch which we will discuss more in Sec. 4).
In the scaling regime, the parameters ξ andv ∞ are approximately time-independent. For a specific n, the solution to the evolution equations in the VOS model can be expressed as [81] (2.8) The Goldstone particle radiation term sv 6 ∞ /N is treated as a perturbation (valid when ∆ 1), and v 0 is the solution tov ∞ in the limit where the Goldstone emission term is set to 0. The model parameters {c, k v , σ} can be extracted by calibrating with current simulation results, as we will discuss.
Although the presence of a scaling regime with a constant ξ as in the VOS model has been confirmed by simulations for NG or gauge strings [85][86][87][88][89][90], the situation is not yet clear for global strings. Some of the recent simulation studies such as [47] suggest a time-dependent ξ(t) that grows linearly with N in small 4 N 7 region, where β is a constant bearing large uncertainty related to initial condition. The linear increase of ξ in N is also found in some other simulation studies [46,68,70,75,91,92], but is in conflict with other groups' simulation results, which predict a nearly constant ξ [48,72,73,78,89,[93][94][95][96][97]. This discrepancy is an intriguing puzzle, and requires further investigation with higher resolution simulations. Given the uncertainty, while we mainly focus on the application of the VOS model, here and in Sec. 5.3 we also carefully considered the effect of potential deviation from scaling and suggest modification/extension to the current VOS model. In order to calibrate the parameters {c, k v , σ} for the VOS model, we fit data extracted from simulation results in [45,73,91], as summarized in Table. 1. The error bars are visually estimated from the plots in [45,73,91], as we are doing a simplified statistical analysis as in [81]. Our best fitting result for the VOS model parameters are as follows: {c, k v , σ} {0.497, 0.284, 5.827}, (2.10) and the fitting quality is about 3.3-σ significance (p-value < 0.001). Such a fitting quality reflects moderate tensions among simulation data listed in Table.  different simulation methods as well as the different ways of counting the number of strings that are employed in the literature. We will assume that these same parameters apply for different scenarios of cosmological background, e.g. radiation domination (RD) or matter domination (MD) 2 . As an example, for N = 70, we obtain the number of strings per Hubble volume ξ ∼ 4.0 andv ∞ ∼ 0.57 in RD, and ξ ∼ 3.55 withv ∞ ∼ 0.40 in MD.
In Fig. 1 we show the evolution of ξ andv ∞ as functions of N using the VOS model with the fitting model parameters listed above. Given the recent findings suggesting deviation from scaling (Eq.(2.9)), in the sub-figure of Fig. 1 where the small N region is zoomed in, we also show the 1-σ area of Eq.(2.9) as the yellow band (the error bar is given in [47]), in comparison with the VOS model prediction (sold curve). We found that in the region of small N 7 the VOS model prediction is consistent with a linear growth of ξ in N , provided that β is not too small, e.g. β ∼ 0.20 is taken as an example in our analysis. The late-time evolution in the scaling regime is insensitive to the exact value of β which depends on initial conditions. Nevertheless, as can be seen in Fig. 1, at large N VOS model prediction approaches scaling, i.e. a nearly constant ξ. For our later analysis of the GW signals, the late-time evolution in the region of N 50 region is most relevant, yet is beyond the reach of most of current simulations. In contrast, VOS model provides a reasonable prediction for the entire time range of interest, after calibrating with low N simulation data.
The data points we used and listed in Table.1 were also applied in [81]. We initially considered using a larger data set for the fitting by including more simulation results, but they are in some way in conflict with the data in Table.1. In order to have meaningful results, we decided to leave out the data sets that fit VOS model poorly. The discrepancies among different simulation results could be in part because these simulations are done with very different methods, covering different ranges of N , and the number of strings and the velocities are counted by different numerical algorithms [81]. Further investigations and developments are certainly required to reach a convergence among different simulation results. To fairly consider these other data sets, in the following, we further discuss their implications and why we left them out of our analysis.
First, note that the VOS model is only valid once the string network enters the scaling region N ∼ 6 (e.g. see Fig. 3 of [45]). The evolution in the very early stage of N 5 is sensitive to initial condition. Therefore, for our fitting, we exclude simulation data points with very low N 's such as in [92] (N = 2−4). The result from [72] is not included because we found that its large velocityv = 0.609 ± 0.014 leads to a poor χ 2 fit with other simulation data 3 . Consequently, with these parameters the GW production would be reduced to about 44% compare to the GWs prediction with parameter set Eq.(2.10) as we will see in later section. Ref. [68] simulated the global string network with cosmological background parameter n ≤ 3, without a data point simulated with a radiation dominated background, thus cannot be analyzed with the results included in Table. 1. Another reason we did not include data from [68] is that their results suggest a time-dependent loop chopping rate, which does not match the VOS model. [46] suggests another pattern of deviation from scaling that is inconsistent with Eq.(2.9). While these suggested non-scaling behaviors only directly apply to low N range and do not converge among literature, it is intriguing to consider their potential effects on GW signals (if the non-scaling persist till large N ) and how VOS model would need to be revised accordingly. We leave more discussion on this topic in Sec. 5.3.
In this study, we simply keep the velocity parameter k v as a constant as in the conventional VOS model. Nevertheless, some studies suggested the possibility of velocitydependent momentum parameter k v = k v (v) [78,79,98,99] k where q 2.3, β 1.5, and k 0 1.37 [79,99]. We found that in RD background Eq.(2.11) gives a numerical value of velocity parameter k v (v =v ∞ ) ∼ 0.3 at high N 10 which is consistent with our fitting result Eq.(2.10). In addition, there is a debate about whether the chopping parameterc is time-independent: e.g. [68] suggests thatc decreases with N , while [72] fits a constant valuec = 0.843 ± 0.039 in radiation background. We will not elaborate on these two particular types of uncertainty.  Table 1. Results from recent global string network simulations (in a radiation dominated background) for the number of strings per Hubble volume ξ and the average velocity of long stringsv in radiation dominated background. These data points were also applied in [81]. In the main text, we explain why some other recent simulation results were left out of this table (thus our analysis) and their implications.

Dynamics of global string loops: formation and radiation into GWs and Goldstones
A global cosmic string network forms during the phase transition around T ∼ η. The dynamics of the very early stage of evolution is sensitive to initial conditions. However, the string network would soon evolve towards an initial condition independent scaling regime [45,47,68,75,91], namely, ξ ∼ constant (or with potential deviation from scaling suggested by some recent work, see earlier discussion and later in Sec. 5.3 ). The horizonsized long strings randomly intersect each other and lose energy via forming sub-horizon sized loops, which subsequently oscillate and radiate energy until they decay away. The loop size distribution at formation time can be parameterized by: a distribution function F α and the fraction of energy stored in loops that can be released as radiation (GWs or Goldstones), F α . In this work, we consider two representative scenarios in detail, both inspired by simulation results: (1) a nearly monochromatic loop size at formation i ∼ αt i with α ∼ 0.1, such that F α∼0.1 ∼ 1, while ∼ 90% fraction of loop's energy is in the form of kinetic energy which would eventually redshift away without contributing to GWs, thus F α ∼ 0.1, as inspired by [82,83]; and (2) a flatter, log-uniform distribution of loop size as suggested in [45]. In this section we focus on the simpler first case, and the second scenario will be discussed in Sec. 5.1. In Sec. 5.1 we also comment on other loop distribution possibilities to account for the related uncertainties [100].
By energy conservation, for a specific α the formation rate of string loops in a scaling string network is given by where the chopping rate parameterc is given in Eq.(2.10), and ρ o denotes the energy density of string loops. The number density of loops with length = αt is then The right panel of Fig. 2 illustrates this point with numerical results. With our calibrated parameters, t η as defined corresponds to N ∼ 6 − 7, which implies that the perturbative VOS model [81] has large uncertainties in such a low N range. After formation, a global string loop would rapidly oscillate and emit energy in the form of GWs and Goldstones by the following energy loss rates until the loop disappears completely [101]: Note that the parameter Γ (a) only depends on the loop trajectory [36,102], thus we expect that the Goldstone radiation constant Γ a should be close to the value of the GW radiation constant Γ ∼ Γ a [102] which is also determined by the loop shape. In the following, we assume benchmark values Γ ∼ 50 [82,83,103,104], and Γ a ∼ 65 [36,102]. We will discuss the effect of varying Γ a , Γ in Sec. 5.2 to account for the potential uncertainty on the radiation parameters.
The size of a loop with initial length i = αt i therefore decreases as where t i is the loop formation time. The radiation of GW and Goldstone from a loop can be decomposed into a set of normal-mode oscillations with frequenciesf k = 2k/˜ , where mode numbers k = 1, 2, 3 · · · , and˜ ≡ (t) is the instantaneous size of the loop when it radiates att. We can rewrite the radiation parameters in a decomposed form We have assumed that the cusps are the dominating source of GW and Goldstone emissions as found in NG string simulations [105][106][107]. The contributions from kinks and kink-kink collisions follow different power laws: Γ (k) ∝ k −5/3 and k −2 for kinks and kink-kink collisions, respectively [108][109][110]. As shown in Eq.(2.16), relative to Goldstone emission, GW radiation is suppressed by a factor of ∼ η 2 /m 2 p , where m p is Planck scale. Nevertheless, the suppression factor becomes less severe as the symmetry breaking scale η gets closer to m p .
Our main analysis results shown in Sec. 3 are obtained by focusing on the simple, motivated assumptions made in this section. Nevertheless, we acknowledge other possibilities of C eff and Γ (a) that were suggested in literature. We further discuss the effects on phenomenology in light of possible deviations from our assumptions on these factors in Sec. 5.2 and Sec. 5.3.

SGWB Spectrum from Global Cosmic Strings (Standard Cosmology)
In this section we will first show the derivation and numerical results of SGWB frequency spectrum from a global cosmic string network assuming a standard cosmic history (Sec. 3.1, 3.2). Then in order to give more physics explanation and insights, in Sec. 3.3 we provide parametric estimates for the relic densities of Goldstones and GWs emitted from global strings, and compare with GW signals from NG strings. Comparison with related results in other literature is also given.

Derivation of GW spectrum from global strings
The generic form of the relic energy density of a SGWB is given by where ρ GW is the energy density of GWs, and ρ c = 3H 2 0 /8πG is the critical density. String loops emit GWs from normal mode oscillations with frequenciesf k = 2k , where k ∈ Z + , is the loop size at emission timet [47,111]. Taking into account of redshift effects, the observed frequencies today are then The relic GW background is obtained by summing over all harmonic modes Using Eq.(2.13) and Eq.(2.17) that we derived earlier and integrating over emission timet, we can derive the contribution Ω where t F is the formation time of the global string network, t 0 is the current time, and the causality and energy conversation conditions are imposed by With Eq.(2.17) we can derive that a loop that emits GW at timet leading to an observed frequency f was formed at the time Note that to consider the radiation of Goldstones, we may define Ω Gold (f ) in analogy to Eq.(3.4) with the simple replacements: Γ → Γ a , and ΓGµ 2 → Γ a η 2 . We will apply this prescription in later discussions involving Goldstone radiation (e.g. Sec. 3.3).
Earlier studies based on radiation dominated background [102,112,113] found that the first few k-modes dominate the GW radiation from loops. However, recent work (in the context of NG strings) showed that a large value of k 10 5 may be needed to converge, depending on the background cosmology [26,30,114]. For instance, including higher k modes changes the power-law index of Ω GW (f ) from −1 to −1/3 in a MD epoch. In this work we investigated the importance of high k modes in the context of global strings and found similar results. We found that to reach a converging result for Ω GW (f ) up to ∼ 100 Hz, i.e. within the frequency range relevant for current and near-future GW detections, up to k ∼ 10 8 modes need to be included. Higher f range requires more k modes to converge. To draw the full spectrum shown in Fig. 3 we took into account of up to k ∼ 10 15 modes in our analysis. Fig. 3 demonstrates the SGWB spectrum calculated numerically with the method we outlined. The corresponding results for NG strings are also shown in comparison, with more explanation given in Sec. 3.3. While the very high f range f 100 Hz is well beyond the reach of any foreseeable experiment, we keep it in Fig. 3 for theoretical completeness by capturing physics at times as early as the formation time of the string network. As can be seen, towards high f the spectrum falls most significantly starting around the frequency the falling spectrum at frequencies f > f η has uncertainties and is sensitive to the initial condition and the very early stage of string network evolution, which is not captured by the VOS model. Then over a wide range of f the spectrum gradually declines towards higher f (∼ ln 3 (1/f ), see Eq.(3.7) below), corresponding to the emissions during the RD era. Note that this feature of the SGWB spectrum from global strings is in contrast to a nearly flat plateau as in its NG string counterpart. Starting at f 0 2 αt 0 ∼ 3.6 × 10 −16 Hz, the spectrum behaves as f −1/3 until f eq ∼ 1.8 × 10 −7 Hz, which is due to the transition to the late MD era. f eq indicates the frequency corresponding to the emission around the matter-radiation equality time. We will elaborate the f -T or f -t correspondence later in Sec. 4.1. Note that the f −1/3 behavior was obtained by summing up to high oscillation modes k 10 5 which was shown to be important for a MD background [26,30]. By only summing up to low k modes (k 10 5 ) it would be f −1 . The low end of the frequency spectrum has a cutoff corresponding to emission at the present time t 0 , with a maximum point shortly before the ending of the spectrum at f 0 .
By combining Eq.(3.13) and Eq.(4.11) we derive the following analytical approximation for global string SGWB spectrum in different f regions, which shows the parametric dependence: where t eq and z eq denote the time and redshift at the matter-radiation equality, respectively. ∆ R (f ) accounts for the effect of varying the number of relativistic degrees of freedom, g * and g * S , over time: where g * (f ) and g * S (f ), are obtained by applying the f -T relation which will be introduced later in Eq.(4.1), and the superscript 0 indicates the values today. Note that here we focus on global strings associated with massless Goldstones, thus they are stable until the current time. In the case of axion strings with massive Goldstones, the string network would turn to domain walls and finally disintegrate around the transition time when Goldstones acquire masses. In that case, the GW spectrum would beget a cut with potentially distinct structure around a characteristic f cut that is larger than f 0 .

GW frequency spectrum and experimental sensitivities
Fig . 4 illustrates the SGWB signal originated from global strings based on our numerical results. We also include related experimental sensitivities: current constraints (solid lines) from LIGO [3,4,116,117] and European Pulsar Timing Array (EPTA) [118], Parkes Pulsar Timing Array (PPTA) [84,119,120]; the projected future sensitivities (dashed lines) with LIGO A+ [121], LISA [7], DECIGO/BBO [51], AEDGE/AION [50,122], Einstein Telescope (ET) [123,124], Cosmic Explorer (CE) [1], and Square Kilometer Array (SKA) [125]; as well as the region corresponding to the recent NANOGrav excess [9,10]. We can see that as expected from Eq.(3.13) the global string GW spectrum is sensitive to the symmetry breaking scale η. Experiments such as LISA, BBO and SKA can probe η 10 14 GeV. Among the existing searches, PPTA gives the strongest constraint of η 2 × 10 15 GeV. These sensitivities/constraints on η may be improved/relieved with non-standard cosmology and alternative modelings, see discussions in Sec. 4 and Sec. 5. Various intriguing interpretations of the recent NANOGrav excess as a SGWB signal have been considered [11][12][13][14][15][16][17][18]126]. In particular, [49] and [17] investigated the possibility of fitting the NANOGrav signal with GWs from QCD axion strings or general ALP strings. The former [49] found that the GW amplitude hinted by the NANOGrav data requires f a 10 15 GeV which is in conflict with bound on ∆N eff from BBN and CMB data, given that the axions are emitted as radiation from the strings. Nevertheless, the latter suggests that a non-standard cosmological history may improve the fit [17]. In our independent check by including high k modes in the summation, we find that the GW frequency spectrum follows a power-law f −1/3 in the range of f ≤ f eq (defined before Eq.(3.7)), and with 4.3 × 10 15 GeV ≤ η ≤ 6.1 × 10 15 GeV, global strings can lead to a good 1-σ fit to the NANOGrav 12.5-year data [9]. However, as also discussed in [9,11], such a spectrum with a gentle slope is in tension with previous bounds from PPTA [119,120], EPTA [118], and NANOGrav 11-year data [8]. Such a tension may be eased by re-analyzing the data sets using different choices of the red noise model [127] which is being investigated. Variations to the standard theoretical assumptions may allow a viable interpretation of the NANOGrav signal as originated from a global/axion string network, consistent with PPTA data and ∆N eff constraints, which we will explore in future study. In Sec. 3.3, we will discuss the ∆N eff bound on Goldstone and GW emissions in detail. Other relevant constraints on the global U (1) breaking scale η include inflation scale and CMB anisotropy bound, which were discussed in [24,93], also pointing to η O(10 15 ) GeV. CMB polarization data potentially yields stronger bound on GW in the frequency range of 10 −17 − 10 −14 Hz [119,128,129]. Nevertheless, in Sec. 4.1 we will demonstrate that this latter constraint does not apply to our case following the introduction of the f -T relation.
Comparison with literature: GW signals from a global string network have also been recently investigated by simulation approaches based on a Nambu-Goto effective theory [49] or field theory for global defects [48]. Our results agree with others' on some general features such as Ω GW ∝ η 4 , but differ in details. [49] simulated the global string network in a radiation background during a very early stage of evolution, i.e. N 7-8, and extrapolated the linear growth of ξ ∝ N to high N when computing the GW spectrum. They agree with our finding that the global strings can lead to detectable GW signals, but found that the GW spectrum scales as Ω GW ∝ η 4 N 4 , instead of η 4 N 3 as found in our analysis (see Eq.(3.13) in Sec. 3.3). The N 3 dependence we found results from the prediction of the conventional scaling VOS model. The difference may be resolved if the loop emission factor C eff in the VOS model is not (nearly) a constant but evolves as C eff ∝ N (see Eq.(2.14). We further discuss the effect of such a non-scaling behavior or deviation from the conventional VOS in Sec. 5.3. On the other hand, [48] found that the GW spectrum asymptotes to an exact scale invariant form, and the amplitude of the signal is below the prediction by both our method and by [49]. The possible explanations for this discrepancy was suggested in [24,49], while further investigation is certainly needed to fully resolve this issue.  [130] or d(lnf )∆Ω rad h 2 8.1 × 10 −7 [130,131], which requires η 3.5 × 10 15 GeV.

Comparison with GWs from NG strings, relic densities of GWs and (massless) Goldstones
In this subsection, we give a simple estimate for the relic density of GWs from global strings which captures key parametric dependence, and compare it with that for NG strings. This can help us gain insights about the detectability of the GW signal from global strings. In addition, while in this work we focus on GW radiation from global strings, it is important to better understand Goldstone emission which is the dominant radiation mode in this case. We thus also present a parametric estimate for the relic density of the emitted Goldstones and compare it with GW emission. As shown in Fig. 5, with these analyses we can find the constraint on η considering the upper limit on extra radiation energy density ∆N eff from BBN/CMB data. The emitted GWs can potentially affect CMB observables in other ways (e.g. CMB polarization) and lead to other constraints, which we will discuss later in Sec. 4.1. As mentioned earlier, in this work we focus on the simple case with massless Goldstones and our discussion about Goldstone emission is illustrative and concise. Nevertheless, some key insights can be applied to axion strings where the Goldstones are massive as potential dark matter candidates. We leave a detailed discussion on Goldstone radiation and its impact on axion DM physics for future work.
A key difference between the dynamics of a global and a NG string network is that the global string loops are rather short-lived due to the strong Goldstone emission rate. We consider a loop formed at time t i which decays away at time t r ≡ γ r (t i )t i , where we have adopted a unified notation for the cases of NG and global string loops for an easy comparison: r = {NG, global}. Using Eq.(2.17) we find the following expression for estimating the lifetime parameter γ(t i ) for the two cases: (3.10) Our ansatz of α ∼ 0.1 ΓGµ is applied to derive the final results. The lifetime of a loop formed at time t i with an initial length of αt i can then be estimated as τ r = (γ r (t i ) − 1)t i . Recent simulations support our estimates of global string loop's lifetime [47,71]. Due to the time dependence of global string tension (i.e. the N -dependence), κ varies in the range of 0.6 1/κ 10 throughout the expansion history of universe. Therefore, the global strings are short-lived and are expected to decay in about one Hubble time after its formation (but the lifetime is still sufficient to yield detectable GWs with large η). In contrast, as can be seen from Eq.(3.9), the NG string loops generally survive a much longer time after formation. Due to this drastic difference in loop lifetime, with the same parameters such as η and loop distribution function, GWs from a global string network on average experience a larger redshift effect after emission, which contributes to a suppressed GW amplitude (along with the suppression effect due to the Goldstone dominance) and shifts the spectrum towards lower frequencies. These considerations help us understand the numerical results as shown in Fig. 3.
We now estimate the relic densities of GW and Goldstone emitted from global strings. As mentioned in Sec. 3.1 the formulation for GW calculation given in Eq.(3.4) and Eq. (3.5) can be applied to the Goldstone case with the replacements of Γ → Γ a , and ΓGµ 2 → Γ a η 2 . (based on Eq.(2.16)). We can then express the total relic densities (integrated over f ) of GWs and Goldstones from global string radiation in the following unified form: Our numerical results of the relic energy densities are illustrated in Fig. 5 as functions of symmetry breaking scale η, along with Ω GW from NG strings for comparison. The upper limit on the total relic radiation energy density from the CMB data is also shown [130][131][132]. One can see that the constraint is dominantly driven by the emission of radiation-like Goldstones, which requires η 3.5×10 15 GeV, while for GWs alone the constraint is relaxed to η 9 × 10 15 GeV. In comparison, with a non-scaling solution as suggested in Eq.(2.9) this CMB constraint on η would be tighter: η 9 × 10 14 GeV [49], as the total energy of the string network would increase relative to the scaling scenario (see Sec. 5.3 for more related discussion).
Next we further discuss the parametric dependence of Ω GW and Ω Gold for global strings and compare with Ω GW for NG string. For a fair comparison, we assume that the symmetry breaking scale η and the string network evolution parameters such as the long string number density ξ and loop size α are the same for the NG and global string network under consideration. Then we consider Ω Global GW , Ω Global Gold and Ω NG GW as observed at a time parametrized by N ≡ ln(Lη). Based on simple analytic estimates checked with numerical fitting, we find the following relations: where α ΓGµ is the lifetime parameter for NG string, γ NG r , as defined in Eq.(3.9), which accounts for the aforementioned difference in redshift effects between global and NG case, and the square-root of ΓGµ α is due to the redshift of the GW energy ∝ a(t) ∝ t 1/2 ; the N factors account for the log enhanced string tension for global strings; ΓGµ Γa/(2πN ) represents the different energy loss rates to GWs vs. to Goldstones. We also find the following key parametric dependencies (focusing on η and N ) for each of these Ω's: The η dependence of GWs from NG strings that we found agrees with earlier literature [21,22,36,83,84,103,108,[133][134][135], and Ω Global GW ∝ η 4 agrees with two most recent independent simulations [49] and [48]. Nevertheless, the N -dependence of the scaling solution of long string number density ξ in the VOS model (see Eq.(2.7)) disagrees with some of the simulation results which suggest a logarithmic increase in ξ based on low N data [49]. The effect of a non-scaling ξ persisting till late times, e.g. N > 70, will be discussed in Sec. 5.3, including a comparison with the result in [49].

Probing the Early Universe with SGWB Spectrum from Global Strings
In this section we investigate how the SGWB spectrum from a global string network would alter if the cosmic history and particle content of the early Universe differ from the standard scenario which we assumed in Sec. 3. This in turn allows us to use such GW signals to test the standard paradigms and probe the dynamics of the early Universe well before BBN. Such an idea of using GWs for cosmic archaeology was proposed and developed in the context of NG strings [21,22]. The situation with global strings bear similarities with that of NG strings, yet with significant differences. In the following, we will demonstrate our findings and make comparison with NG string results.

The connection between the observed GW frequencies and emission times
In the context of NG strings, the frequency-temperature (f -T ) correspondence during a RD era was derived in [22], and serves as the foundation of cosmic archaeology with the f spectrum of GWs from strings. The analogous relation for global strings can be derived following the same method. Nevertheless, the derivation can be greatly simplified in this case. As explained in Sec. 3.3 (Eq.(3.13)), a key difference between NG and global string loop dynamics is that, global string loops decay away within ∼ 1 Hubble time after formation due to the strong Goldstone emission rate. Therefore, the timescale when the GW emission from a loop occurs is approximately the same as the loop's formation time, i.e.t ∼ t i (Eq. (2.17)). For an estimate, it suffices to focus on the k = 1 mode which we find to be the dominant one in the cases of interest. With this understanding and following the calculation in Sec. 3.1, we find that a specific f ∆ band observed today relates to a particular emission temperature T ∆ in the following way: where the loop size at the emission time (t) αt i ≡ αt ∆ (see Eq.(2.17)), z eq 3387 is the redshift at the matter-radiation equality, and t eq , T eq are the corresponding time and temperature, respectively. Note that f ∆ linearly depends on T ∆ , but is insensitive to the symmetry breaking scale η, unlike in the case of local strings. Eq. (4.1) applies to RD era, while f -T relation varies with background cosmology, which we will discuss in Sec. 4.2. A departure from the standard cosmology at T ∆ would thus imprint itself in the GW spectrum around the corresponding f ∆ .
In Fig. 6 we illustrate the f -T relation derived for SGWB spectrum from global strings, in comparison with the recent results for NG strings [21,22]. There are two major differences between the two cases: f -T correspondence for NG strings has η-dependence while for global strings it is almost independent of η which makes it more robust in a way; for the same f the corresponding emission T is much earlier for global strings than for NG strings. Both these differences originate from the aforementioned fact that global string loops decay shortly after formation and their resultant GW signal observed in a certain f band has undergone longer period of redshift after emission (relative to its NG counterpart). Note that due to the current bounds from LIGO and PPTA, η for NG strings is constrained as η 1.89 × 10 13 GeV [22,84]. If the recent NANOGrav excess indeed originates from NG cosmic strings, it favors η 3 − 5 × 10 13 GeV [11,12]. According to Fig. 6 these constraints/potential signal implies that GW spectrum from NG strings can reach up to T ∼ 10 4 GeV (with ET and CE). In contrast, as shown in Fig. 6 global strings can probe much earlier cosmic history, up to T ∼ 10 8 GeV. As discussed in Sec. 3.2, η 10 14 GeV is needed to be within experimental sensitivity reach in terms of Ω GW , while other constraints require η O(10 15 ) GeV. Fig. 7 illustrates the f -T relation for global strings in a different manner where the sensitivity to η is explicitly shown. As demonstrated in Fig. 7, global string GWs can trace the cosmic history over a rather wide range in time: up to T ∼ 10 8 GeV (with ET and CE) and down to T ∼ 10 −4 GeV (with PPTA and SKA) which intriguingly corresponds to the beginning of the BBN era.
The f -T relation as we have elaborated can also help us understand why the global string scenario safely evades the potentially strong bound on Ω GW in the range of f ∼ 10 −17 − 10 −14 Hz by the CMB polarization data [119,128,129,132]. The f -T relation in Eq.(4.1), together with the observation that global string loops decay in one Hubble  time, indicate that the SGWB signal below a certain f range could not be generated until after a certain time or below a certain T . In Fig. 8 we illustrate the constraints from CMB polarization data, and the decomposed contributions to a SGWB induced by global strings: the signal in the low f range of f ∼ 10 −17 −10 −14 Hz in fact is not populated until after the photon decoupling, thus is not present at the CMB epoch to be subject to the constraint. One can also simply estimate f corresponding to the photon decoupling T γ ∼ 0.3 eV using Eq.

Probing new phases of cosmological evolution
According to the standard thermal history, the Universe is radiation dominated starting from the end of inflation all the way down to the matter-radiation equality at z eq ∼ 3000. Nevertheless, so far there is no data evidence to support this assumption for the epoch prior to the BBN time, i.e. the primordial dark age. On the other hand, recently there has been substantial interest to consider well-motivated non-standard cosmology scenarios, where the standard RD era transits to a different phase at some point in the early Universe, such as EMD or kination. An EMD era can be due to the temporary domination of a long-lived massive particle or oscillations of a scalar moduli field [136]. More generic possibilities arise in models where a scalar field φ oscillates in a polynomial potential V (φ) ∝ φ N , characterized by an averaged equation of state w = (N − 2)/(N + 2). In the limit N → ∞, we have n = 6 in Eq.(2.6) which is called kination phase, as the kinetic energy of the scalar dominates. Kination can generally arise in inflation [137], quintessence, dark energy [138], and axion-like particle (ALP) models with varying power of sin-Gordon potential [139] or with a non-zero initial field velocity [140,141]. In order to retain the successful predictions of BBN theory, for all these scenarios the Universe needs to settle to RD before the BBN time T ∆ ∼ 5 MeV.
It is thus intriguing to see how the SGWB from global strings would alter in a nonstandard cosmology and the related implication for detections. From another perspective, similar to the finding in the context of NG strings, SGWB from global strings thus opens up the possibility of probing the early Universe during the primordial dark age that may not be directly accessible otherwise. This allows us to test the standard assumption about cosmology while uncovering potential deviations. The base of this method lies in the f -T relation during RD (Eq.(4.1)) which allows us to relate a deviation from the standard prediction for the SGWB frequency spectrum to a time point in history where RD transits to a new (earlier) phase. To calculate Ω GW (f ) with a non-standard cosmology background, we follow the method given in [22,24]: we assume that the Universe transits from RD to a new equation state parametrized by n (Eq.(2.6): ρ ∝ a −n ) at T ∆ , and match the energy density at T ∆ for a smooth transition. Ω GW (f ) can then be calculated using Eq.(3.4) with the input of a non-standard evolution of a(t).
We therefore expect that a non-standard cosmology leads to a modified GW spectrum in the high frequency region starting round f ∆ corresponding to the transition in cosmic history occurring around T ∆ (see Eq.(4.1)). Numerically, we found that Ω GW (f ) can be parametrized in the following way in large f region for general cosmologies (parametrized by n) : where t ∆ is the time corresponding to the temperature T ∆ , and we have assumed α = 0.1. Eq. (4.2) shows that Ω GW (f f ∆ ) ∝ f +1 for kination (n = 6) and ∝ f −1/3 for MD (n = 3). Note that the validity of the VOS model approach requires n > 2, otherwise both ξ andv ∞ would me imaginary valued according to Eq. 2.7. Another caveat is that, at sufficiently large f f η such that log(...) ∼ 1 or N ∼ O(1) (corresponding to the very early stage after the string network formation), Ω GW would universally fall as Ω GW ∝ f −1/3 , for different background cosmologies.
In Fig. 9 we show our numerical results for benchmark examples of GW spectrum from a scaling global cosmic string network with a non-standard cosmology background such as kination or EMD, contrasted by the standard prediction shown in solid black line. We can see that compared to standard cosmology, with the presence of an EMD phase Ω GW (f ) falls faster towards higher f , making it harder to observe in that f range. On the other hand, the spectrum rises above the standard prediction at high f in case of an early kination phase, leading to a stronger signal. The kination case thus is more subject to existing constraint from LIGO. In general, the LIGO O3 constraint on T ∆ can be expressed as where we have dropped the logarithmic dependence from Eq. concrete example with two stages of transitions where kination is preceded by an earlier RD era and investigate the constraints on the model due to the CMB/BBN bound on ∆N eff (both Goldstone and GW).

A two-stage transition scenario with kination:
We assume that at T = T ∆2 > T ∆1 kination transits to an early RD era (note: not the later standard RD era). Such a scenario can be realized if, for instance, a dominating radiation species decays to kination particles around T ∆2 . Other possibilities of exiting kination at high T exist, e.g. by a vacuum energy dominated phase such as inflation. However, a long period of vacuum energy domination would dilute the overall GW signal significantly [26,142]. An alternative is to have a short duration of vacuum energy domination (mini-inflation) which then transits to an early RD. Kination can also be preceded by a dominating matter-like species that decay to kination particles. Here we choose to consider the simple scenario of RD-kination-RD for illustration.
Examples of GW spectrum of such a two-stage transition scenario are shown in the leftpanel of Fig. 10: Ω GW linearly rises with f in a finite frequency range of f ∆2 > f > f ∆1 due to kination, then restores the logarithmically decreasing behavior in the range of f > f ∆2 (Eq.(4.2)) during the early RD. The three benchmark cases shown satisfy both LIGO O3 bound and the CMB ∆N eff bound that we will discuss next. The characteristic frequency f ∆1 corresponding to the later stage of transition at T ∆1 can be estimated by Eq. (4.1). Similarly, based on Eq.(3.2), the frequency corresponding to the earlier transition at T ∆2 can be estimated as (applying ρ ∝ a −6 for kination) Now we consider the implication of the CMB ∆N eff bound on additional relic radiation for this kination example. We assume the equation of state of the Goldstones emitted from global strings is radiation-like. As suggested by the sharp rising of GW spectrum shown in Fig. 10 in the presence of a kination phase, the relic radiation energy densities of Goldstone and GW from the string network are dominated by the emission during the kination epoch, T ∆1 < T < T ∆2 , which can be roughly estimated as (4.5) where Ω {GW,Gold} are defined in Eq. (3.11), and numerically computed/fitted based on Eq. (3.4). We also defined the reference values with η = 10 15 GeV assuming the standard cosmology: As discussed in Sec. 3.3, the emission of radiation-like Goldstone dominates the bound. The right panel of Fig. 10 illustrates three viable benchmark scenarios assuming η = 10 15 GeV, parametrized by T ∆1 , T ∆2 : with T ∆2 = 25 GeV, 180 GeV, and 3150 GeV, the CMB ∆N eff bound requires T ∆1 1 GeV, 10 GeV, and 100 GeV, respectively.

Probing new degrees of freedom
Many BSM theories involve new particles that are relativistic and in thermal equilibrium in the early Universe, e.g. in many potential solutions to the electroweak hierarchy problem [143][144][145][146], and theories of dark sectors [147][148][149][150][151][152][153][154][155][156][157]. These particles would contribute to the effective number of relativistic degrees of freedom (DOFs) in energy, g * , and in entropy, g * S , in the high T Universe, but can generally be out of reach of available probes such as by the LHC or CMB experiments due to heavy masses or feeble interactions with the SM. The methodology for calculating the effect of new DOFs on the SGWB spectrum of NG strings was introduced in [22]. In this work, we briefly review the method and apply it to obtain results in the case of global strings.
We illustrate the effect of new massive DOFs on the string GW spectrum without referring to the details of the underlying theory. We model the change in the number of DOF with the following assumption where g * rapidly decreases as T drops below a mass threshold T ∆g [22]: To numerically demonstrate the effect, we choose the well-motivated scenario, where T ∆g is of weak scale, which may be motivated from solutions to the Hierarchy Problem. In particular, in Fig. 11 we choose the benchmark values of T ∆g = 200 GeV, ∆g * = 0, 10 2 , 10 3 , and assume g * g * S . As can be seen, relative to the prediction with SM DOFs only, the spectrum falls towards higher f starting from a frequency f ∆g that agrees with the prediction by the f -T relation in the RD era (see Eq.(4.1)).
Such an effect can be understood by analytical estimates following [22]. Deep in the RD regime, the Hubble rate and the corresponding time depend on g * in the following way: where H 0 is the current Hubble constant, and Ω R is the radiation energy relic density observed today. Note that ∆ R (a) is simply a variational form of ∆ R (f ) as defined in Eq. (3.8). Applying this simplification in Eq.(3.4), we have where Ω SM GW (f ) indicates the amplitude with SM DOFs only. Eq.(4.11) clearly shows that the overall amplitude of the high f tail (f > f ∆g ) of Ω GW decreases with the presence of additional DOFs, agreeing with numerical findings.  Figure 11. Modification to the GW spectrum from a global string network due to an increase in the number of relativistic degrees of freedom above T ∆g = 200 GeV. In the example shown, η = 10 15 GeV, α = 0.1, and ∆g * = 0, 10 2 , 10 3 (shown in black, red, and blue, respectively). The relevant experimental sensitivities are also shown.

Sensitivity to the loop size parameter α and its distribution
As introduced in Sec. 2.2, throughout our work we have used α 0.1 as the peak value of loop sizes at their formation time, which is inspired by results from NG string simulations [82,83]. However, there are still uncertainties about loop distribution for global strings.
To investigate how such uncertainties may impact the predicted GW spectrum, in this subsection we consider two alternative scenarios of loop distribution: 1. varying α for the peak value, and 2. a log uniform distribution of loops as suggested in [45]. Alternative-1: varying α for the peak value.
The analysis with different α values is straightforward with our formulations in Sec. 3. In the left panel of Fig. 12 we show the α dependence of Ω GW (f ) for specific f 's normalized by the prediction with α = 0.1 (the benchmark choice used in earlier sections) with different background cosmologies, assuming η = 10 15 GeV. As shown, RD, MD and kination dominated eras have different dependencies on α, which are insensitive to f for the cases of MD and kination. The α dependence can be discussed in two distinct regions. Firstly, in the range of α < κ ∼ 0.11, the loop lifetime is shorter than a Hubble time, and thus the analysis and discussion in the previous sections can apply: the 2nd line in Eq.(3.7) explains the result for RD. As the loops decay quickly after formation, for a given GW frequency observed today, on average scenarios with smaller α would be associated with larger N at the emission time in order to allow for a longer period of redshifting. Larger N corresponds to higher string tension and thus larger string energy density available for GW production. In addition, with smaller α the loops emit GW with higher frequency and thus the spectrum blue shifts as shown in the right panel of Fig. 12. Given these two effects, the GW spectrum appears to increase in amplitude as α decreases in this regime. We find that in a kination epoch the spectrum linearly increases with α, while in matter  domination Ω GW (f ) ∝ α −1/3 . In the other region of α > κ, the loops are long-lived, and thus the spectrum in RD agrees with the NG string case, which gives Ω GW (f ) ∝ α 1/2 [22]. In this large α region, Ω GW (f ) still linearly increases with α in kination, while becomes approximately α independent in MD. Alternative-2: A log uniform distribution. Fig. 5 of the recent global string simulation [45] suggests a logarithmic uniform distribution of the size of string loops at formation time, which is very different from the nearly monotonous α that we have assumed inspired by NG string simulation. While this hint of log uniform distribution is yet to be further tested, we consider how this variation can impact the prediction for GW signals. A log uniform distribution indicates that at formation time the loop number density dn( )/d at size follows dn( )/d(log ) ∼ const, or dn( )/d ∝ 1/ . Our method of calculating SGWB signal with a monotonous loop formation size α as shown in Eq. (3.4) can be adapted to this alternative distribution by replacing F α /α... in Eq.(3.4) (the "..." part represent other parts in the formula for computing GWs) with a sum over thinly sliced loop sizes in the range of π/η < < π/H: For the cases with kination or EMD, the departure from standard cosmology is assumed to occur at T ∆ = 1 GeV. and taken the continuous limit to get the second equality. The parameters y, δ, α 0 , α 1 are introduced to rewrite the integration limits in more convenient forms: max ∼ π/H ≡ α 0 t ≡ e y t, and min ∼ π/η ≡ α 1 t ≡ e −δ+y t. However, in our numerical calculation we found that including small loops down to the scale of π/η leads to very large Ω GW (f ) 1 in certain f range as a consequence of energy conservation. Therefore, we assume a lower cutoff of α at α ∼ 10 −4 . The exact value of small scale cutoff on α is not essential for our study here, as our purpose is to simply show an example of how a log uniform distribution can alter the GW spectrum.
In Fig. 13 we show the GW spectrum predicted with the assumed log uniform distribution for different cosmology scenarios. We find that by summing over the loop sizes in the range of 10 −4 ≤ α ≤ 2π, the GW amplitude is generally increased over many decades in the frequency range except around the cutoff around f 0 ∼ 10 −16 Hz. Due to the inclusion of larger loops up to α = 2π in the distribution, the low frequency cutoff extends to ∼ 2/(2πt 0 ).

Sensitivity to the loop radiation parameter Γ and Γ a
While we chose motivated benchmark values of loop radiation parameters Γ and Γ a in our main studies, we acknowledge that there are still uncertainties around these values. Here we investigate how the GW signal would change by varying Γ and Γ a . Considering energy conservation law and energy loss rates in Eq.(2.16), naively, we expect the GW density to depend on Γ, Γ a simply as Ω GW ∝ ΓGµ 2 Γaη 2 . However, such dependencies can be more complex as the redshift-related factors a(t) and t merical results (fixing f = 1 Hz and η = 10 15 GeV and Γ = 50 for example). Γ dependence is simpler, linear as naively expected, unless GW becomes the dominant radiation mode (Γ a Γ). We will show the Γ dependence explicitly in the following formulae/discussion. As can be seen in Fig. 14 there are three distinct regions in the Ω GW (f ) − Γ a relation, which we can understand analytically as follows: • Large Γ a , such that loops decay within a Hubble time after formation, driven by strong Goldstone emission. In this region α < κ, where κ ≡ Γ a /(2πN ) (Eq.(3.6)), and we can estimate with N ∼ 70 for relevant observations. The Γ a term thus dominates both numerator and denominator of Eq.(3.6), which implies that both a(t) and t • Medium size Γ a , such that loops survive beyond a Hubble time after formation, while Goldstone radiation still dominates over GWs. In this region, α > κ > ΓGµ, and thus redshift factors a(t) and t (k) i in Eq.(3.4) depend on Γ a . By fitting numerical results we find the following relations which depend on background cosmologies: for RD and EMD, Γ Γ a , for Kination.
As discussed in Sec. 4.2, the GW frequency spectrum with an EMD is dominated by loop radiation during the later radiation domination era. Consequently, the Ω GW (f )-Γ a relation is approximately the same as RD for the benchmark frequency f = 1 Hz.
• Small Γ a , such that Γ a ΓGµ 2 /η 2 (i.e. κ < ΓGµ, and the Goldstone emission term in Eq.(2.16) becomes negligible relative to the GW radiation). Given the hierarchy between the Planck mass and the viable η value considering the relevant constraints, this scenario is only possible for very small Γ a Γ. In this case, GW radiation would become the dominant energy loss mechanism and Ω GW (f ) would increase as, Γ −1/2 which agrees with the related result for NG strings [22].

Non-scaling solution
In this subsection, we consider the impact of possible non-scaling solutions on the GW signals. The violation of the scaling properties in the case of global strings were found in some recent simulation studies [45-47, 68, 70, 75, 91, 92]. This suggests that the attractor solution of the average number of strings per Hubble patch, ξ, logarithmic growing with N . Note that in most of these studies the non-scaling behavior is found in the low N regime which is within direct reach of current simulations, and whether such a behavior can apply to large N still needs to be investigated. For example, Ref. [158] pointed out that the discrepancy in literature could be due to the different interpretations of the initial growth of ξ by different groups: the data set in the earlier study Ref. [45] may be too close to the string formation, and thus is sensitive to initial condition. They further showed that the solution should converge to the constant estimation of ξ 1.19 as given in [72,73]. We found that by assuming the same VOS parameters, such a constant ξ would decrease GW amplitude to about 15% due to C eff ∝ ξ 3/2 dependence (see Eq.(2.14)).
As earlier shown in Fig. 1, the VOS model can be consistent with the non-scaling solution Eq.(2.9) (or Eq.(5.4) below) within the range of low N , 3 N 7, then predicts ξ ∼ const. for larger N . Nevertheless, it is intriguing to see how the GW spectrum would change if such a behavior does sustain throughout the evolution history of the string network, and whether/how a variation to the original analytical VOS model may match this behavior. We will focus on the following two examples and then comment on other possibilities, and in both cases we adopt the non-scaling solution as suggested in simulations [47,49,68] ξ = 0.24(2)N + 0.2, (5.4) where N ≡ ln(η/H(t)). We consider ξ taking the above non-scaling form in both examples that we will discuss next, and adopt the relevant parameters from the VOS model for the GW calculations (Eqs.(2.12,2.13,2.14)).
In the first possibility we consider, in addition to Eq. 5.4, we apply the following benchmark parameters: a constant average velocity of long stringsv ∞ 0.50 ± 0.04, and a loop chopping parameterc = 0.497, which we obtained in Sec. 2.1 based on fitting simulation results (Table. 1). With Eq.(2.13), primarily derived based on energy conservation, we find the prediction for effective loop formation parameter C eff ∝ N 3/2 . With this C eff as an input for Eq.(3.4) we computed the GW spectrum, and found that the amplitude is larger than the prediction in [49] by a factor of O(10 − 100), depending on frequencies. This discrepancy motivated us to introduce the second scenario which is found to lead to a good agreement with [47]: while still assuming Eq. (5.4), this example involves a time-dependent cv ∞ , and consequently a different form of C eff : Based on the analysis method given in Sec. 3.3, the GW spectrum with the non-scaling solution Eq.(5.5) in a RD background (the spectrum would be cut off at lower frequencies by a QCD-like phase transition as shown in [49]) can be estimated as The notable difference between this result and that based on the scaling VOS model solution (Eq.(3.7)) is the power law index of the log term (i.e. log 4 vs. log 3 ), which enhances the GW amplitude by O(10) for this non-scaling example. The enhancement is due to the increase in loop number density (Eq. 2.13). Fig. 15 illustrates the GW spectrum predicted with a non-scaling solution where ξ ∝ N , including a comparison between our results based on a variation to the VOS model (Eq. (5.5)) and the result in [49] based on extrapolating simulation results to large N . A good agreement between our second scenario (Eq.(5.5)) and that in [49] can be seen in Fig. 15. In particular ,our analytic fit for the GW spectrum (Eq.(5.6)) captures the key log 4 dependence that agrees with [49]. This agreement suggests that the extrapolation of the non-scaling solution to large N may be reproduced in a variation to the original VOS model, where the relationscv ∞ ∝ N −1/2 and ξ ∝ N are realized. This hint may be helpful for future investigations. As a supplemental discussion, in Fig. 16 we illustrate and compare the different predictions of GW spectrum based on the two aforementioned non-scaling scenarios, with various background cosmologies. As shown, in the second scenario (Eq.(5.5)) the GW spectrum amplitude is lowered by O(10) relative to the first scenario (Eq.(5.4)), which is due to the different predictions for loop number density (C eff ∝ N v.s. C eff ∝ N 3/2 ). In the frequency range of our interest, the result is insensitive to the initial condition dependent parameter β: as discussed in [45,47], the linearly growing term in Eq.(5.4) would quickly dominate the string network evolution.
A very different form of non-scaling solution was suggested in another simulation work [46]: where T PQ ∼ η is the temperature when the PQ symmetry breaking occurs. The prediction for ξ as in Eq.(5.7) is significantly larger than non-scaling results from other groups' simulations [45,47,68,91,92]. As suggested by the authors of [46], the prediction of Eq.(5.7) only provides a rough counting for cosmic strings, which may address the discrepancy, while further investigations are needed. We attempted to fit Eq.  Figure 15. GW spectrum in the radiation dominated epoch assuming a non-scaling solution (ξ = 0.24(2)N + 0.2): a comparison between the result with our assumption/method and that obtained in the recent simulation work [49]. The data points with error bars are taken from [49]. The blue dashed curve is based on our analytical estimate Eq.(5.6). The red curves with shadowed uncertainty band is based on our numerical calculation of Eq.(3.4) with linear growth of C eff ∝ N . Further details are given in the main text. VOS model, but found a rather poor VOS model fit for the 13 data points provided in [46] due to the large value ξ given in Eq.(5.7), which is inconsistent with other simulation results. Assuming the non-scaling behavior as in form of Eq.(5.7) sustains till late times, we expect the GW amplitude to be amplified by a factor of O(10 − 100) relative to the scaling case due to the larger loop density implied (similar to the case inspired by [45]).

Distinguish from other SGWB sources
In this subsection, we discuss potential challenges for detecting a SGWB signal from global strings in practice, including astrophysical background and a comparison with other cosmological sources of SGWB.
Upon detection of a cosmogenic SGWB signal, it is important to analyze and identify the nature of the underlying physics. Global cosmic string is among many motivated new physics sources that can give rise to a SGWB [32,178,179], for example, primordial inflation [180,181] and black hole [182], preheating [183][184][185][186], first-order phase transitions [187][188][189][190][191], and other types of topological defects [134,192] including local/NG strings [35,103,[193][194][195][196]. A key to distinguishing the various cosmological sources lies in the GW spectral information. For instance, SGWB from a first-order phase transition features a peaky spectrum in frequency associated with specific split power laws, which results from the fact that the GWs were emitted during a specific epoch in the early Universe. In contrast, SGWBs from cosmic strings (both global and NG) feature a rather long (nearly) flat plateau towards high frequencies, due to the continuous emission throughout the cosmic history. We refer to [22] for more detail regarding the general comparison of SGWB originated from cosmic strings with other cosmological sources. Here we highlight the prospect of distinguishing SGWB from global strings vs. that from NG strings. As seen in Sec. 3.1 and Fig. 3 the GW spectrum from global strings has a long tail which logarithmically declines towards high frequencies, whereas the spectrum from NG strings is very close to simple flatness (except for the mild steps due to the change in g * ). A main cause of such a difference is the logarithmic time-dependence of the global string tension, Eq.(5.4). The difference would be further amplified if the non-scaling behavior as discussed in Sec. 5.3 is confirmed to last till late times. In practice, we therefore expect that for global strings, GW searches at lower frequencies such as SKA in general have a better prospect of detection than those at higher frequencies such as LIGO (the prospect also depends on the experimental sensitivities). In summary, while challenges for experimentally detecting a global string sourced SGWB are present, potentially promising solutions exist and will be further developed in coming years. Using frequency band information is a common potential solution for disentangling a global string signal from both astrophysical background and other cosmological sources, which will be strengthened with a multi-band GW experimental program [119,164,173,176].

Conclusion
Global or axion cosmic strings are well motivated sources of SGWB, and have attracted growing interest in the past few years. In this work, we applied the analytical VOS model in solving the evolution of a global string network over the course of the cosmic history, and illustrated the procedure of calculating the resultant GW signals with great detail. We demonstrated how our VOS model parameters were calibrated by simulation data which are most reliable for the early time of evolution N 7, and commented on the compatibility between VOS model prediction and simulation results found by various groups. We found that the deviation from the scaling property as found by some simulation studies can be consistent with conventional VOS model prediction in the early regime of 3 N 7, but the simple extrapolation of such a non-scaling behavior to large N or late times contradicts the conventional VOS model. Nevertheless, we also investigated how the SGWB signal would alter if the non-scaling does persist to late times, and suggested a possible revision to the VOS model that addresses this difference which can lead to a GW signal prediction consistent with that given in [49] based on simulation. While it will take time to resolve the discrepancies among different simulation data sets as well between VOS model prediction and some simulation results, our methodology of analysis and related discussions are timely complements to the literature, and can be further improved/updated in light of future developments.
Our main results are presented following the standard VOS model and analytical calculation of SGWB by summing over harmonic modes and taking into account the significant effect of Goldstone emissions. In light of the recent findings on the important effect of high k modes for NG strings, we summed over 10 5 modes in all our analysis, leading to updated spectra relative to our earlier results in [24]. We first demonstrated the results assuming standard cosmology, and then considered the possible presence of a non-standard equation of state before BBN, e.g. an EMD or kination, which would lead to a drastic departure from the standard prediction at high f ranges. Since an indefinitely long kination period is subject to strong constraints on additional relic radiation energy density from CMB/BBN data due to a blue-tilted spectrum, we also considered an example where the kination epoch has a finite window of span and is preceded by an early stage of RD. We further demonstrated how the presence of new relativistic degrees of freedom in the early Universe can alter the GW spectrum. We showed the current and projected future sensitivities of GW detectors in detecting global string signals, and found that a detectable signal requires the corresponding spontaneous symmetry breaking scale η > 2 × 10 14 GeV. Different from NG strings, GW amplitude from global strings is very sensitive to η (Ω GW ∝ η 4 ). The frequency-time (temperature) relation, which is the foundation for the method of GW cosmic archaeology, takes a very different form for global strings relative to its NG string counterpart. In particular, there is no Gµ dependence in the f -T relation and the same f band corresponds to a much earlier emission time for global strings, which enables us to test the standard radiation era up to T ∼ 10 8 GeV. We explained the physics behind the notable differences between SGWB from global strings and from NG strings, where the strong rate of Goldstone emission and the consequent short lifetime of global string loops play an important role.
We further considered how the GW signal based on our baseline assumptions and model choices could vary with alternative possibilities. We studied the effects of different loop distribution patterns, the uncertainty in the radiation parameters (Γ, Γ a ), as well as a persisting non-scaling regime during the string network evolution. For example, by adopting the suggested non-scaling solution with ξ ∝ N while assumingcv ∞ ∝ N −1/2 , we found that the predicted GW frequency spectrum (including the log 4 relation) can be consistent with the simulation-based finding in [49]. We also briefly discussed the prospect of distinguishing a SGWB sourced by global strings from other cosmogenic sources or astrophysical background.
It is worth noting the importance of studying GWs from global/axion strings in light of its connection to axion physics (for QCD axion or general axion-like particles (ALPs)). Axion strings are indispensable companions of axion particles when the U (1) PQ symmetry breaking occurs after inflation. The detection of axion particles is being actively pursued, but the prospect is model-dependent due to the uncertainty of the interaction between the SM and the axions. The prospect would be particularly dim in the case of the hidden axion scenario (e.g. motivated by the string axiverse [197]) where non-gravitational coupling is absent. Thus, the universal GW signal from axion strings could be the smoking-gun for the underlying axion physics. While our current work focused on the simpler case of pure global strings with massless Goldstones, the methodology and results are relevant for the massive axion case. The GW spectrum from axion topological defects at high frequencies is expected to be dominated by cosmic strings, while at a low frequency corresponding to a QCD(-like) phase transition when the domain wall forms, the signal is expected to change form and die down. The recent work based on extrapolating simulation results to late times shows such a pattern [49]. It is intriguing to apply our analytical approaches to the more complex axion scenario including domain wall contributions, which will be pursued in future work.
completed. The work is supported in part by the US Department of Energy under award number DE-SC0008541.