Gravitational waves from cosmic strings associated with pseudo-Nambu-Goldstone dark matter

We study stochastic gravitational waves from cosmic strings generated in an ultraviolet-complete model for pseudo-Nambu-Goldstone dark matter with a hidden $\mathrm{U(1)}$ gauge symmetry. The dark matter candidate in this model can naturally evade direct detection bounds and easily satisfy other phenomenological constraints. The bound on the dark matter lifetime implies an ultraviolet scale higher than $10^9~\mathrm{GeV}$. The spontaneous $\mathrm{U(1)}$ symmetry breaking at such a high scale would induce cosmic strings with high tension, resulting in a stochastic gravitational wave background with a high energy density. We investigate the constraints from current gravitational wave experiments as well as the future sensitivity. We find that most viable parameter points can be well studied in future gravitational wave experiments.


I. INTRODUCTION
Null experimental results from the direct detection of dark matter (DM) have put strong constraints on the DM-nucleon scattering cross section [1][2][3], implying that DM particles might not have a weak interaction strength, which is required by the conventional freeze-out mechanism of DM production [4][5][6]. Nevertheless, it is possible to suppress DM-nucleon scattering but leave the DM annihilation cross section at the freeze-out epoch unaffected. This can be gracefully realized by assuming that the DM particle is a pseudo-Nambu-Goldstone boson (pNGB) whose scattering off nucleons is extremely suppressed by low momenta , evading the constraints from direct detection experiments.
Such a pNGB DM framework typically requires a proper ultraviolet (UV) completion with an extra gauged U(1) symmetry [7,18,19,34]. The experimental bound on the DM lifetime demands that the UV scale where this gauge symmetry breaks down could be higher than O(10 10 ) GeV [18,19,34]. The spontaneous breaking of the U(1) gauge symmetry in the early universe may results in cosmic strings [39,40], whose tension, i.e., the energy per unit length, would be rather large because of the high UV scale. The network of cosmic strings are expected to radiate gravitational waves (GWs) [41,42], which would form a stochastic background remaining as a relic in the present universe.
The stochastic GW background (SGWB) from cosmic strings covers an extremely broad range of GW frequencies [43]. Therefore, such a background is an interesting target for various types of GW experiments in different frequency bands, including pulsar timing arrays (PTAs) in 10 −9 -10 −7 Hz, ground interferometers in 10-10 3 Hz, and future space interferometers such as LISA [44], TianQin [45,46], and Taiji [47] in 10 −4 -10 −1 Hz. As a result, searching for an SGWB provides a unique way to test many new physics theories beyond the standard model (SM) that would induce cosmic strings in the early universe .
In this work, we will study the stochastic GWs radiated by cosmic strings originating from the UV-complete theory of pNGB DM. As mentioned above, the required high UV scale leads to cosmic strings with a rather high tension, which results in a high energy density of stochastic GWs. Thus, GW experiments could be more sensitive to this theory than other types of experiments, which would typically lose sensitivity for high UV scales. More specifically, we will discuss a simple model with a hidden U(1) gauge symmetry [34] as an illuminating example and explore the constraints from previous GW experiments and the sensitivities of future GW experiments.
The remainder of this paper is organized as follows. In Section II, we briefly introduce the UVcomplete model for pNGB DM with a hidden U(1) gauge symmetry. In Section III, we study the phenomenological constraints and perform a random scan in the parameter space. In Section IV, we discuss the SGWB arising from cosmic strings, which are induced by the spontaneous breaking of the U(1) gauge symmetry. In Section V, we investigate the constraints from current GW experiments and estimate the sensitivity in future GW experiments. In Section VI, we summarize the paper.

II. PNGB DM MODEL
In this section, we briefly discuss the UV-complete model of pNGB DM that extends the SM with a hidden U(1) X gauge symmetry. Details can be found in Ref. [34]. "Hidden" means that all the SM fields do not have U(1) X charges. Two complex scalar fields S and Φ carrying U(1) X charges 1 and 2 are introduced, and they are singlets under the SM gauge group SU(3) C × SU(2) L × U(1) Y . We denote the SM Higgs field as H. The related Lagrangian is where are the field strengths of the SU(2) L , U(1) Y , and U(1) X gauge fields W a,µ , B µ , and X µ . The parameter s ε ≡ sin ε induces the kinetic mixing between B µ and X µ . The covariant derivatives are given by For realizing the pNGB DM framework, S and Φ gains nonzero vacuum expectation values (VEVs) v S and v Φ , which satisfy v S v Φ . Thus, v Φ is a UV scale below which the U(1) X gauge symmetry spontaneously broken to a U(1) X global symmetry, which is approximate, because the µ SΦ term softly breaks it. v S is a lower scale where the U(1) X global symmetry is spontaneously broken, leading to a pNGB with a mass arising from µ SΦ . Such a pNGB is a DM candidate whose scattering off nucleons is extremely suppressed [7].
We decompose the scalar fields as where the SM Higgs VEV is v = 246.22 GeV. There are mass mixing terms among the CP -even scalars (h, s, φ) and between the CP -odd scalars (η S , η Φ ). After diagonalizing the mass-squared matrices, we obtain the physical scalar mass terms The Higgs bosons h 1 , h 2 , and h 3 are linear combinations of h, s, and φ. We further require that h 1 is the SM-like Higgs boson, whose mass is measured to be m h 1 = 125.25 ± 0.17 GeV [79], and that h 2 and h 3 are the s-like and φ-like exotic Higgs bosons, respectively. χ is a linear combination of η S and η Φ with a mass given by This is the pNGB DM candidate we desire. Another linear combination of η S and η Φ is a massless Nambu-Goldstone boson eaten by a gauge boson.
After the spontaneous breaking of the gauge symmetry, the SU(2) L × U(1) Y × U(1) X gauge fields acquire the mass terms Taking into account the B µ -X µ kinetic mixing and the B µ -W 3,µ mass mixing, the physical neutral gauge fields (A µ , Z µ , Z µ ) can be derived through linear combinations of (B µ , W 3,µ , X µ ) [80]. The electromagnetic field A µ is massless, while the masses for the Z and Z bosons are given by [81] where t ε ≡ tan ε,ŝ W ≡ sinθ W ,ĉ W ≡ cosθ W , andθ W ≡ tan −1 (g /g). We define r ≡ m 2 Z /m 2 Z , and t ξ ≡ tan ξ can be expressed as [82] Note that Z is an exotic neutral vector boson.
The interactions in this model have been explicitly discussed in Ref. [34]. Here, we briefly summarize the phenomenologically important couplings: • The h i χχ and h i h j h k couplings come from the scalar interactions.
• The h i f f Yukawa couplings for any SM fermion f arise from the SM Yukawa couplings and the mixings among the Higgs bosons.
• The h i W W , h i ZZ, h i Z Z , and h i ZZ couplings come from the covariant kinetic terms of the scalar fields.
• The Zf f , Z f f , Zχχ, and Z χχ neutral current couplings are induced by the gauge interactions and the kinetic and mass mixings of the gauge fields.
• The Zχh i and Z χh i couplings originating from the neutral current gauge interactions could lead to χ decays. These couplings vanish in the v Φ → ∞ limit, where χ is stable. Thus, they must be greatly suppressed by a high UV scale v Φ to give a sufficiently long lifetime for the DM candidate χ.

III. PHENOMENOLOGICAL CONSTRAINTS
Ten free parameters in this pNGB DM model can be chosen as s ε , λ HS , λ HΦ , λ SΦ , v S , v Φ , m χ , m h 2 , m h 3 , and m Z . As a UV completion of pNGB DM, we are interested in the parameter According to effective field theory [83], the tree-level of the spin-independent (SI) χ-nucleon scattering cross section is obtained as [34] σ SI where m N represents the nucleon mass and f N u,d,s are nucleon form factors [84]. Note that σ SI χN is highly suppressed by v −4 Φ . For v Φ 10 5 GeV, σ SI χN is typically below 10 −50 cm 2 [34]. Therefore, direct detection experiments would hardly probe the pNGB DM candidate χ.
As mentioned above, the Zχh i and Z χh i couplings result in χ decays. For m χ m h 3 ∼ m Z , the pNGB DM candidate χ could decay via χ → h where h 3 and Z are off shell, but Z, h 1 , and h 2 can be either on or off shell. These processes are all suppressed by the UV scale v Φ . In order to explain dark matter in the present universe, χ must have a very long lifetime. A conservative bound on the DM lifetime, i.e., τ χ 10 27 s, has been given by Fermi-LAT γ-ray observations of nearby dwarf spheroidal galaxies [85]. This puts a strong constraint on v Φ .
In this model, χχ annihilation channels involve h i h j (i, j = 1, 2), W + W − , ZZ, and ff . Contrary to χ-nucleon scattering, χχ annihilation has no particular suppression; thus, the observed DM relic abundance can be easily obtained via the freeze-out mechanism. The χχ annihilation processes at the present would induce γ rays, which can be probed by indirect detection experiments.
Moreover, the couplings of the SM-like Higgs boson h 1 to W , Z, and the SM fermions deviate from the SM. Exotic decays h 1 → χχ, h 1 → χZ, and h 1 → h 2 h 2 may occur for m h 1 > 2m χ , m h 1 > m χ + m Z , and m h 1 > 2m h 2 , respectively. Therefore, the LHC Higgs measurements have put some constraints on the parameter space of this model. In addition, the exotic Higgs boson h 2 could be directly produced at high energy hadron colliders. The h 2 production processes include gluon-gluon fusion gg → h 2 , vector boson fusion W W/ZZ → h 2 , and h 2 production associated with W , Z, tt, or bb, while the h 2 decay channels involve tt, bb, W + W − , ZZ, γγ, Zγ, etc. These processes lead to detectable signals of h 2 at the LHC and the Tevatron. Nonetheless, the constraints from LHC and Tevatron direct searches can be evaded if the h component in h 2 is small and/or h 2 is heavy.
In order to study the phenomenological constraints on the pNGB DM model, we use FeynRules 2 [86] to implement the model and generate model files for the numerical package micrOMEGAs 5.2 [87]. Utilizing micrOMEGAs, we calculate the χ lifetime τ χ , the DM relic abundance Ω χ h 2 , and the χχ annihilation cross section σ ann v . HiggsSignals 2 [88] is applied to test whether the properties of the h 1 boson are consistent with the LHC measurements, while HiggsBounds 5 [89] is employed to constrain the h 2 boson by LHC and Tevatron direct searches.
A random parameter scan is carried out within the following ranges of the free parameters: The parameter points are randomly generated by uniform distributions on the logarithmic scale. We further demand the induced couplings g X , λ H , λ S , and λ Φ lying from 10 −2 to 1. Because we have required v Φ > 10 8 GeV, the direct detection constraints [1][2][3] are totally irrelevant.
The following criteria are used to select the parameter points.
1. In order to guarantee the vacuum stability, the following conditions for the scalar potential from copositivity criteria [90] should be satisfied: 2. The lifetime of the pNGB DM particle χ should satisfy the Fermi-LAT bound τ χ 10 27 s [85].
4. The total χχ annihilation cross section σ ann v should not be excluded by the upper limits at the 95% confidence level (C.L.) given by the combined Fermi-LAT and MAGIC γ-ray observations of dwarf spheroidal galaxies in the bb channel [92].
5. The signal strengths of the SM-like Higgs boson h 1 should be consistent with the LHC Higgs measurements at 95% C.L. based on the HiggsSignals calculation.
6. The exotic Higgs boson h 2 should not be excluded at 95% C.L. by the direct searches at the LHC and the Tevatron according to HiggsBounds.
In Fig. 1, the selected parameter points are projected onto the m χ -σ ann v and v Φ -λ Φ planes. In order to achieve the correct DM relic abundance, most of the parameter points have a canonical annihilation cross section σ ann v 2 × 10 −26 cm 3 /s, as shown in Fig. 1(a). Deviations from this value are mainly due to the h 1 and h 2 resonance effects. Such resonances could greatly increase the χχ annihilation cross section at the freeze-out epoch, but the present σ ann v would be rather small, because the center-of-mass energy of χχ at low velocities today would not be close to the resonance centers. Obviously, the deep valley around m h 1 /2 63 GeV in Fig. 1(a) corresponds to the h 1 resonance. The h 2 mass is not fixed, and the h 2 resonance effect is probably responsible for other parameter points with σ ann v that deviate from the canonical value. From Fig. 1(a), we find that in order to cover all the parameter points, the sensitivity of indirect detection experiments needs to be improved by three orders of magnitude. The color axis in Fig. 1(a) indicates the tree-level SI DM-nucleon scattering cross section σ SI χN , whose values are far too low to be reached by direct detection experiments. Fig. 1(b) shows that the UV scale v Φ should be higher than 10 9 GeV 1 , because of the bound on the DM lifetime. Since a smaller v Φ tends to give a shorter lifetime of χ, the parameter points with smaller v Φ have higher probabilities to be rejected by the DM lifetime bound than those with larger v Φ . Considering the DM detection experiments and the measurements and searches at colliders, GW experiments provide a complementary way to test this model. In the following sections, we will discuss the potential GW signals associated with the selected parameter points.

IV. GRAVITATIONAL WAVES GENERATED BY COSMIC STRINGS
The spontaneous breaking of the U(1) X gauge symmetry at the high UV scale v Φ is expected to induce cosmic strings, which are one-dimensional topological defects concentrated with energies of the scalar and gauge fields [93]. The key quantity of a cosmic string is its tension µ, which is the energy per unit length. According to the analysis based on the Abelian Higgs model [94], the tension in the pNGB DM model can be estimated as follows: , a high UV scale v Φ suggested by the Fermi-LAT bound on the χ lifetime would lead to cosmic strings with high tension. As we will see below, this would result in an SGWB with a high energy density ∝ Gµ 2 .
The dimensionless quantity Gµ, where G is the Newtonian constant of gravitation, is commonly used to describe the tension of cosmic strings. According to Eq. (13), we calculate Gµ for the selected parameter points, and the result is shown on the v Φ -λ φ plane in Fig. 1(b) with the color axis ranging from 10 −19 to 10 −7 . The positive correlation between Gµ and v Φ is clearly demonstrated.
A network of cosmic strings would be formed in the early universe after the spontaneous breaking of the U(1) X gauge symmetry. According to the analysis of string dynamics [43], the intersections of long strings could produce closed loops, whose size is smaller than the Hubble radius. Cosmic string loops could further fragment into smaller loops or reconnect to long strings. The relativistic oscillations of the loops due to their tension emit GWs, and the loops would shrink because of energy loss. Moreover, the loops typically have special features called cusps and kinks, which could produce GW bursts [95]. Consequently, the energy of the cosmic string network is converted into the energy of GWs, and an SGWB is formed due to the incoherent superposition of GWs.
At the emission time t e , a cosmic string loop of length L emits GWs with frequencies where the integer n = 1, 2, 3, · · · indicates the harmonic modes of the loop oscillation. Denoting P n as the power of gravitational radiation for a mode n in units of Gµ 2 , the total power is given by [96] P = Gµ 2 n P n .
According to the numerical simulation of smoothed cosmic string loops, the averaged power spectra P n for loops in the radiation-and matter-dominated eras are obtained in Ref. [97], including the contributions from cusps. The total power can be characterized by the dimensionless quantity which is estimated to be ∼ 50 in the two eras [97]. We will adopt this result for P n and Γ = 50 in the calculation of the GW spectrum below.
Defining n(L, t) dL as the number density per physical volume of cosmic string loops with length L at cosmic time t in length interval dL, the GW energy density ρ GW induced by the cosmic string network per unit time at the emission time t e can be expressed as Utilizing Eq. (14), we derive Because of the cosmological redshift effect, the GW frequency f observed at the present time t 0 differs from the emission frequency f e . They are related by where a(t) is the cosmological scale factor and z represents the redshift. Setting a(t 0 ) = 1, the GW energy at the present is given by ρ GW (t 0 ) = ρ GW (t e )a 4 (t e ), and we obtain where t * represents the time when the GW emissions start. According to dt e = −[H(z)(1 + z)] −1 dz, where H(z) represents the Hubble expansion rate, we can rephrase Eq. (20) as The dimensionless quantity commonly used to characterize the SGWB is the spectrum of the GW energy density per logarithmic frequency interval divided by the critical density ρ c = 3H 2 0 /(8πG), i.e. , The Hubble constant is usually expressed as H 0 = 100h km s −1 Mpc −1 with h = 0.674 ± 0.005 [91]. In order to avoid the uncertainty on the Hubble constant, one prefers to use Ω GW (f )h 2 . For calculating Eq. (21), we need to know H(z). In a flat ΛCDM universe, the Hubble rate is given by [98] H where Ω m = 0.315 ± 0.007 [91], Ω r = 1.68 × (5.38 ± 0.15) × 10 −5 [79], and Ω Λ = 1 − Ω r − Ω m represent the energy densities of matter, radiation, and dark energy relative to the critical density at the present, respectively. The function takes account of the change in the number of radiation degrees of freedom between redshift z and the present, where g * and g * s represent the effective numbers of relativistic degrees of freedom for the energy and entropy densities. Considering electron-positron annihilation at z 10 9 and the QCD phase transition at z 2 × 10 12 , we have [98] 0.83, 10 9 < z < 2 × 10 12 , 0.39, z > 2 × 10 12 .
The last gradient that we need to compute the energy density spectrum of the SGWB is the loop number density distribution n(L, t). There are various approaches for modeling n(L, t), which can lead to significant differences in the GW spectrum. Here, we discuss two typical models for loop production.
The first model introduced by Blanco-Pillado, Olum, and Shlaer (BOS) [99] makes use of the scaling nature of the cosmic string network and takes the horizon distance to be the only kinematic scale. By extrapolating the loop production function found in numerical simulations of Nambu-Goto strings [100], n(L, t) is derived for any given cosmic time. Then, the loop number density produced in the radiation era can be approximated as n r (L, t) 0.18 t 3/2 (L + ΓGµt) 5/2 θ(0.1t − L), where θ(x) is the Heaviside step function. The cosmic string loops produced in the matter era give a contribution of Moreover, there would be cosmic string loops produced in the radiation era and still surviving in the matter era. Their number density is given by where t eq = 51.1 ± 0.8 kyr [79] represents the cosmic time at the matter-radiation equality. Therefore, we have n(L, t) = n r (L, t) for t < t eq and n(L, t) = n m (L, t) + n r→m (L, t) for t > t eq .
The second model given by Lorenz, Ringeval, and Sakellariadou (LRS) [101,102] follows the Polchinski-Rocha approach [103] and assumes that the number density distribution of produced loops per unit time P(L, t) is given by a power law in the form of where γ ≡ L/t is a dimensionless variable and χ is a parameter characterizing the power law. Furthermore, they took into account the gravitational backreaction effect, which prevents loop production below a certain scale, by considering a different power law for γ < γ c with χ changed to another parameter χ c . Here, γ c 20(Gµ) 1+2χ [104] characterizes the length scale of gravitational backreaction. The resulting loop number density distribution in scaling is found to be insensitive to the value of χ c and can be expressed as n(L, t) = t −4 N (γ), with the function N (γ) approximately given by [101] N (γ) Here, γ d = −dL/dt ΓGµ represents the shrinking rate of cosmic string loops. ν is the exponent in the relation between the scale factor and the cosmic time, i.e., a(t) ∝ t ν , with ν = 1/2 in the radiation era and ν = 2/3 in the matter era. The parameters C and χ can be expressed as where C • and p are extracted by the distribution of scaling loops from Nambu-Goto simulations [105]: Now, we can evaluate the energy density spectrum of the SGWB induced by the cosmic string network arising from the U(1) X gauge symmetry breaking in the pNGB DM model. From the selected parameter points, we choose four benchmark points (BPs), whose Gµ varies from O(10 −17 ) to O(10 −11 ), to demonstrate the results. Table I presents the detailed information for the four BPs. For each BP, we calculate the corresponding Ω GW h 2 assuming the BOS and LRS loop production models, as plotted in Fig. 2(a) with solid and dashed lines, respectively. In the BOS model, the high frequency behavior is controlled by the GW emissions during the radiation era [99]. As the integral contributed by n r (L, t) tends to be a constant at high frequencies, Ω GW h 2 becomes quite flat in the high frequency regime, as indicated by the solid lines in Fig. 2(a). For smaller Gµ, the total GW emission power is smaller, and cosmic string loops could survive longer, leading to more smaller loops radiating at higher frequencies. Consequently, the GW spectrum moves downward and rightward as Gµ decreases [97]. Compared with the BOS model, the LRS model leads to GW spectra with far higher amplitudes, especially at high frequencies. This is because the LRS model gives a very high number density of small loops in the γ < γ c regime, which significantly contribute to high frequency GWs [106,107].

V. CONSTRAINTS AND SENSITIVITY OF GW EXPERIMENTS
As shown in Fig. 2(a), the SGWB arising from cosmic strings extends over a very broad frequency range. Thus, various GW experiments can probe it, and their sensitive bands are f (Hz) denoted by shaded regions in Fig. 2(a). In this section, we investigate the constraints from current GW experiments on the pNGB DM model and the sensitivity in future experiments.
In Fig. 2(b), we demonstrate the sensitivity curves of several GW experiments and compare them with the GW spectra of BP2 estimated by the BOS and LRS models. For PTA experiments aiming at 10 −9 -10 −7 Hz, we show the 95% C.L. upper limits from the European Pulsar Timing Array (EPTA) [108], the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) [109], and the Parkes Pulsar Timing Array (PPTA) [110], as well as the projected strain noise spectra expressed in terms of the GW energy density [111] (denoted as Ω n h 2 below) for the International Pulsar Timing Array (IPTA) [112] and the Square Kilometer Array (SKA) [113]. For ground-based laser interferometers, whose sensitive band is 10 1 -10 3 Hz, we plot Ω n h 2 for the O5 runs of the Laser Interferometer Gravitational Wave Observatory (LIGO) [114] and for the Cosmic Explorer (CE) [115]. For future space-borne laser interferometers searching for GWs at 10 −4 -10 −1 Hz, we display Ω n h 2 for LISA [44], TianQin [116], and Taiji [117]. For TianQin, Ω n h 2 is calculated from the detector noise power spectrum density and response function of the A/E channel given in Ref. [116].
There are constraints on Gµ from previous GW experiments, depending on the assumptions of the loop production models. According to the LIGO-Virgo O3 data, assuming that the average numbers of cusps and kinks per loop oscillation are both 1, the constraint for the BOS model at 95% C.L. is Gµ < 1.96×10 −8 , while that for the LRS model is Gµ < 4.83×10 −15 [118]. Compared with the BOS model, the LRS model predicts a higher amplitude in the frequency band of ground interferometers, leading to a far stronger constraint on Gµ. The 11-year NANOGrav data yield a 95% C.L. upper limit of Gµ < 5.3 × 10 −11 for the BOS model [109]. Utilizing the PPTA data v Φ (GeV) over 15 years assuming a non-auto Hellings-Downs correlation [119], the 95% C.L. constraints are estimated to be Gµ < 2.88×10 −11 for the BOS model 2 and Gµ < 9.12×10 −12 for the LRS model. These constraints and the selected parameter points are displayed on the v Φ -Gµ plane for the BOS and LRS models in Figs. 3(a) and 3(b), respectively. The color axes indicate the parameter b = 8g 2 X /λ Φ , which relates v Φ to Gµ through Eq. (13). Because of the strong correlation between Gµ and v Φ , the parameter points align around a straight line, and the dependence on b scatters the points. We find that the parameter points with v Φ 5 × 10 13 (7 × 10 11 ) GeV in the pNGB DM model have been excluded assuming the BOS (LRS) model for loop production.
Below, we evaluate the sensitivity of future GW experiments according to their expected Ω n h 2 shown in Fig. 2(b). For LISA, Taiji, CE, or SKA with an idealized auto-correlation measurement in an accessible frequency range f min ≤ f ≤ f max , the signal-to-noise ratio (SNR) ρ can be estimated as [111,121] where T represents the practical observation time. If the SNR ρ turns out to be larger than an appropriate threshold, which will be set as ρ thr = 10 in the following, the predicted GW signal is likely to be detected.
For estimating the SNR of TianQin, we follow the strategy in Ref. [116] with the null channel method, where the T channel is constructed to highly suppress the SGWB signal and the A and E channels are sensitive to the signal. Assuming an ideal symmetric scenario where the A and Because two channels are used, this expression has an additional factor √ 2 compared with Eq. (34).
We assume that the observation time for LISA, Taiji, TianQin, or CE is T = 1 yr, and that for SKA is T = 10 yr [122]. Then, the SNR ρ in these experiments is evaluated as a function of Gµ for the BOS and LRS models, as shown in Figs. 4(a) and 4(b), respectively. The estimated SNRs for the four BPs of the pNGB DM model are presented in Table I.
Taking the SNR threshold to be ρ thr = 10, the expected upper limits on Gµ for LISA, Taiji, TianQin, CE, and SKA are presented in Table II. These upper limits are plotted in Fig. 3 for comparison with previous experimental results. For the BOS model, the most sensitive experiments are LISA, Taiji, and CE, which have comparable sensitivities and could probe v Φ down to ∼ 2 × 10 10 GeV. For the LRS model, CE has the highest sensitivity, and the parameter points with v Φ down to ∼ 5 × 10 9 GeV could be detected. In order to clearly demonstrate the sensitivity of future GW experiments to the pNGB DM model, we also project the selected parameter points onto the v φ -ρ LISA plane for the BOS model and onto the v φ -ρ CE plane for the v Φ (GeV)   Finally, we discuss a related experimental anomaly. According to the 12.5-yr data set, the NANOGrav collaboration reported strong evidence of a stochastic common-spectrum process, but they found no significant evidence of the quadrupolar spatial correlations, which are necessary to claim the detection of an SGWB [123]. If this signal is genuine, it can be explained by the SGWB induced by cosmic strings, and the 68% (95%) confidence interval of the string tension implied by the data is 4 × 10 −11 < Gµ < 10 −10 (2 × 10 −11 < Gµ < 3 × 10 −10 ) [124]. In Fig. 6, we compare the NANOGrav confidence intervals with the selected parameter points on the v Φ -Gµ plane and find that it is possible to explain such a suspicious signal at 95% C.L. by 2 × 10 13 GeV v Φ 2 × 10 14 GeV in the pNGB DM model. Of course, more data from future GW experiments are needed to reveal the nature of this anomaly.

VI. SUMMARY AND OUTLOOK
We studied the stochastic GW signals from cosmic strings induced by the spontaneous breaking of the hidden U(1) X gauge symmetry from the UV-complete model for pNGB DM. Because of the pNGB nature of the DM candidate χ, the tree-level χ-nucleon scattering cross section is highly suppressed by the UV scale v Φ , which characterizes the breaking scale of the U(1) X gauge symmetry. Thus, DM direct detection experiments would not be able to probe χ for v Φ 10 5 GeV. Meanwhile, χχ annihilation processes are not suppressed, and the observed DM relic abundance can be easily achieved by the conventional freeze-out mechanism. Therefore, the pNGB DM model can naturally explain dark matter in the universe and satisfy current experimental constraints.
In order to investigate the phenomenological constraints on the pNGB DM model, we carried out a random scan in the 10-dimensional parameter space. The parameter points simultaneously satisfying the constraints from the DM lifetime, the DM relic abundance, indirect detection experiments, measurements of the SM-like Higgs boson h 1 , and collider searches for the exotic Higgs boson h 2 have been selected.
The spontaneous U(1) X symmetry breaking at the UV scale v Φ is expected to generate cosmic strings in the early universe. The intersections among the cosmic strings would lead to closed loops, which could emit GWs via relativistic oscillations. The incoherent superposition of GWs leads to an SGWB, which is an important target for GW experiments. We evaluated the corresponding GW spectra for the viable parameter points according to the BOS and LRS models for loop production.
Moreover, we studied the constraints from current GW experiments including NANOGrav, PPTA, and LIGO-Virgo. These constraints excluded the parameter points with v Φ 5 × 10 13 (7 × 10 11 ) GeV for the BOS (LRS) model. Furthermore, we estimated the sensitivity of future GW experiments LISA, Taiji, TianQin, CE, and SKA. For the BOS (LRS) model, these future experiments could explore the parameter points with v Φ down to ∼ 2 × 10 10 (5 × 10 9 ) GeV. Note that the bound on the DM lifetime has basically excluded a UV scale v Φ lower than 10 9 GeV. Therefore, almost all the viable parameter points of the pNGB DM model can be well studied in the future.
Typical experiments in particle physics would lose their sensitivity if the related energy scale is too high. In this study, however, a higher UV scale v Φ would lead to a higher tension of cosmic strings, resulting in an SGWB with higher energy density, which would be more easily discovered in GW experiments. Remarkably, we have demonstrated that GW experiments could be complementary to other types of experiments in exploring new physics beyond the SM.
In the calculation above, the value of v Φ is assumed to be a constant. Nonetheless, if the temperature corrections to the effective potential of the scalar fields are considered, v Φ would depend on the temperature at the epoch when the spontaneous breaking of the U(1) X gauge symmetry occurs. This could affect the tension of the cosmic strings formed and hence the SGWB spectrum. The treatment in this paper should be regarded as an approximation without the temperature effect. Such a effect is worthy to be studied in the future.