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

Null experimental results from the direct detection of dark matter (DM) have put strong constraints on the DM-nucleon scattering cross section [1−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−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 [7−38], 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 $ \mathcal{O}(10^{10})\; \text{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}\text{–}10^{-7}\; \text{Hz} $, ground interferometers in $ 10\text{–}10^{3}\; \text{Hz} $, and future space interferometers such as LISA [44], TianQin [45, 46], and Taiji [47] in $ 10^{-4}\text{–}10^{-1}\; \text{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 [48−78].

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 Sec. II, we briefly introduce the UV-complete model for pNGB DM with a hidden ${U}(1)$ gauge symmetry. In Sec. III, we study the phenomenological constraints and perform a random scan in the parameter space. In Sec. IV, we discuss the SGWB arising from cosmic strings, which are induced by the spontaneous breaking of the ${U}(1)$ gauge symmetry. In Sec. V, we investigate the constraints from current GW experiments and estimate the sensitivity in future GW experiments. In Sec. VI, we summarize the paper.

In this section, we briefly discuss the UV-complete model of pNGB DM that extends the SM with a hidden ${U}(1)_\mathrm{X}$ gauge symmetry. Details can be found in Ref. [34]. "Hidden" means that all the SM fields do not have ${U}(1)_\mathrm{X}$ charges. Two complex scalar fields S and Φ carrying ${U}(1)_\mathrm{X}$ charges $ 1 $ and $ 2 $ are introduced, and they are singlets under the SM gauge group ${S U}(3)_\mathrm{C} \times {S U}(2)_\mathrm{L} \times $ $ {U}(1)_\mathrm{Y}$. We denote the SM Higgs field as H. The related Lagrangian is

where $ W^a_{\mu\nu} = \partial_{\mu}W^a_\nu - \partial_{\nu}W^a_\mu + g\varepsilon^{abc}W^b_\mu W^c_\nu $, $ B_{\mu\nu} \equiv \partial_\mu B_\nu - \partial_\nu B_\mu $, and $ X_{\mu\nu} \equiv \partial_\mu X_\nu - \partial_\nu X_\mu $ are the field strengths of the ${S U}(2)_\mathrm{L}$, ${U}(1)_\mathrm{Y}$, and ${U}(1)_\mathrm{X}$ gauge fields $ W^{a,\mu} $, $ B^\mu $, and $ X^\mu $. The parameter $ s_\varepsilon \equiv \sin\varepsilon $ induces the kinetic mixing between $ B^\mu $ and $ X^\mu $. The covariant derivatives are given by $D_{\mu}H = (\partial_{\mu}-{\rm i}g' B_{\mu}/2 - \mathrm{i}g W_\mu^a \sigma^a/2)H$, $ D_{\mu}S = (\partial_{\mu}- \mathrm{i} g_{X}X_{\mu})S $, and $ D_{\mu}\Phi = (\partial_{\mu}- 2\mathrm{i} g_{X}X_{\mu})\Phi $, with $ g_X $ denoting the ${U}(1)_\mathrm{X}$ gauge coupling.

For realizing the pNGB DM framework, S and Φ gains nonzero vacuum expectation values (VEVs) $ v_S $ and $ v_\Phi $, which satisfy $ v_S\ll v_{\Phi} $. Thus, $ v_\Phi $ is a UV scale below which the ${U}(1)_\mathrm{X}$ gauge symmetry spontaneously broken to a ${U}(1)_\mathrm{X}$ global symmetry, which is approximate, because the $ \mu_{S\Phi} $ term softly breaks it. $ v_S $ is a lower scale where the ${U}(1)_\mathrm{X}$ global symmetry is spontaneously broken, leading to a pNGB with a mass arising from $ \mu_{S\Phi} $. 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\; \text{GeV} $. There are mass mixing terms among the $ CP $-even scalars $ (h, s, \phi) $ and between the $ CP $-odd scalars $ (\eta_S, \eta_\Phi) $. 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 \pm 0.17\; \text{GeV} $ [79], and that $ h_2 $ and $ h_3 $ are the s-like and ϕ-like exotic Higgs bosons, respectively. χ is a linear combination of $ \eta_S $ and $ \eta_\Phi $ with a mass given by

This is the pNGB DM candidate we desire. Another linear combination of $ \eta_S $ and $ \eta_\Phi $ is a massless Nambu-Goldstone boson eaten by a gauge boson.

After the spontaneous breaking of the gauge symmetry, the ${S U}(2)_\mathrm{L} \times {U}(1)_\mathrm{Y} \times {U}(1)_\mathrm{X}$ gauge fields acquire the mass terms

Taking into account the $ B^\mu $-$ X^\mu $ kinetic mixing and the $ B^\mu $-$ W^{3,\mu} $ mass mixing, the physical neutral gauge fields $ (A^\mu, Z^\mu, Z'^\mu) $ can be derived through linear combinations of $ (B^\mu, W^{3,\mu}, X^\mu) $ [80]. The electromagnetic field $ A^\mu $ is massless, while the masses for the Z and $ Z' $ bosons are given by [81]

where $ t_\varepsilon \equiv \tan \varepsilon $, $ \hat{s}_\mathrm{W} \equiv \sin \hat{\theta}_\mathrm{W} $, $ \hat{c}_\mathrm{W} \equiv \cos \hat{\theta}_\mathrm{W} $, and $ \hat{\theta}_\mathrm{W} \equiv \tan^{-1} (g'/g) $. We define $ r \equiv m_{Z'}^2/m_Z^2 $, and $ t_\xi \equiv \tan \xi $ 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 \chi\chi $ 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 $ Z f f $, $ Z' f f $, $ Z \chi\chi $, and $ Z' \chi\chi $ neutral current couplings are induced by the gauge interactions and the kinetic and mass mixings of the gauge fields.

● The $ Z\chi h_i $ and $ Z'\chi h_i $ couplings originating from the neutral current gauge interactions could lead to χ decays. These couplings vanish in the $ v_\Phi \to \infty $ limit, where χ is stable. Thus, they must be greatly suppressed by a high UV scale $ v_\Phi $ to give a sufficiently long lifetime for the DM candidate χ.

Ten free parameters in this pNGB DM model can be chosen as $ s_{\varepsilon} $, $ \lambda_{HS} $, $ \lambda_{H\Phi} $, $ \lambda_{S\Phi} $, $ v_{S} $, $ v_{\Phi} $, $ m_{\chi} $, $ m_{h_{2}} $, $ m_{h_{3}} $, and $ m_{Z'} $. As a UV completion of pNGB DM, we are interested in the parameter regions with $ v \sim v_S \ll v_\Phi $, implying a mass hierarchy $ m_{h_1} \sim m_{h_2} \ll m_{h_3} \sim m_{Z'} $. According to effective field theory [83], the tree-level of the spin-independent (SI) χ-nucleon scattering cross section is obtained as [34]

where $ m_N $ represents the nucleon mass and $ f_{u,d,s}^{N} $ are nucleon form factors [84]. Note that $ \sigma_{\chi N}^\mathrm{SI} $ is highly suppressed by $ v_\Phi^{-4} $. For $ v_\Phi \lesssim 10^5\; {\text{GeV}} $, $ \sigma_{\chi N}^\mathrm{SI} $ is typically below $ 10^{-50}\; {{\rm{cm}}^2} $ [34]. Therefore, direct detection experiments would hardly probe the pNGB DM candidate χ.

As mentioned above, the $ Z\chi h_i $ and $ Z'\chi h_i $ couplings result in χ decays. For $ m_\chi \ll m_{h_3} \sim m_{Z'} $, the pNGB DM candidate χ could decay via $ \chi \to h_i^{(*)} Z^{(*)} $ and $ \chi \to h_i^{(*)} Z'^{*} $, 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_\Phi $. 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., $ \tau_{\chi} \gtrsim 10^{27}\; \text{s} $, has been given by Fermi-LAT γ-ray observations of nearby dwarf spheroidal galaxies [85]. This puts a strong constraint on $ v_\Phi $.

In this model, $ \chi \chi $ annihilation channels involve $ h_i h_j $ ($ i,j = 1,2 $), $ W^+W^- $, $ ZZ $, and $ f\bar{f} $. Contrary to χ-nucleon scattering, $ \chi \chi $ annihilation has no particular suppression; thus, the observed DM relic abundance can be easily obtained via the freeze-out mechanism. The $ \chi \chi $ 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 \to \chi\chi $, $ h_1 \to \chi Z $, and $ h_1 \to h_2 h_2 $ may occur for $ m_{h_1} > 2m_\chi $, $ m_{h_1} > m_\chi + 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\to h_2 $, vector boson fusion $ WW/ZZ\to h_2 $, and $ h_2 $ production associated with W, Z, $ t\bar{t} $, or $ b\bar{b} $, while the $ h_2 $ decay channels involve $ t\bar{t} $, $ b\bar{b} $, $ W^+W^- $, $ ZZ $, $ \gamma\gamma $, $ Z\gamma $, 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 $\mathrm{FeynRules 2}$ [86] to implement the model and generate model files for the numerical package $\mathrm{micrOMEGAs 5.2}$ [87]. Utilizing $\mathrm{micrOMEGAs}$, we calculate the χ lifetime $ \tau_\chi $, the DM relic abundance $ \Omega_\chi h^2 $, and the $ \chi\chi $ annihilation cross section $ \langle {\sigma}_{\mathrm{ann}}v\rangle $. $\mathrm{HiggsSignals 2}$ [88] is applied to test whether the properties of the $ h_1 $ boson are consistent with the LHC measurements, while $\mathrm{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 $, $ \lambda_H $, $ \lambda_S $, and $ \lambda_\Phi $ lying from $ 10^{-2} $ to $ 1 $. Because we have required $ v_\Phi > 10^8\; \text{GeV} $, the direct detection constraints [1−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 $ \tau_{\chi} \gtrsim 10^{27}\; \text{s} $ [85].

3. The DM relic abundance $ \Omega_\chi h^2 $ should be in the $ 3\sigma $ range of the Planck value $ \Omega_\mathrm{DM} h^2 = 0.1200\pm 0.0012 $ [91].

4. The total $ \chi\chi $ annihilation cross section $ \langle {\sigma}_{\mathrm{ann}}v\rangle $ 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 $ b\bar{b} $ 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 $\mathrm{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 $\mathrm{HiggsBounds}$.

In Fig. 1, the selected parameter points are projected onto the $ m_{\chi} $-$ \langle {\sigma}_{\mathrm{ann}}v\rangle $ and $ v_\Phi $-$ \lambda_{\Phi} $ planes. In order to achieve the correct DM relic abundance, most of the parameter points have a canonical annihilation cross section $ \langle {\sigma}_{\mathrm{ann}}v\rangle \simeq 2\times 10^{-26}\; {{\rm{cm}}^3/{\rm 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 $ \chi\chi $ annihilation cross section at the freeze-out epoch, but the present $ \langle {\sigma}_{\mathrm{ann}}v\rangle $ would be rather small, because the center-of-mass energy of $ \chi\chi $ at low velocities today would not be close to the resonance centers. Obviously, the deep valley around $ m_{h_1}/2 \simeq 63\; \text{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 $ \langle {\sigma}_{\mathrm{ann}}v\rangle $ 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 $ \sigma_{\chi N}^\mathrm{SI} $, whose values are far too low to be reached by direct detection experiments.

Figure 1.  (color online) Selected parameter points projected onto the $ m_{\chi} $-$ \langle {\sigma}_{\mathrm{ann}}v\rangle $ (a) and $ v_{\Phi} $-$ \lambda_{\phi} $ (b) planes, with colors denoting $ \sigma_{\chi N}^\mathrm{SI} $ and $ G\mu $, respectively.

Figure 1 (b) shows that the UV scale $ v_\Phi $ should be higher than $ 10^{9}\; \text{GeV} $ ¹, because of the bound on the DM lifetime. Since a smaller $ v_{\Phi} $ tends to give a shorter lifetime of χ, the parameter points with smaller $ v_{\Phi} $ have higher probabilities to be rejected by the DM lifetime bound than those with larger $ v_{\Phi} $. 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.

The spontaneous breaking of the $ {U}(1)_\mathrm{X} $ gauge symmetry at the high UV scale $ v_\Phi $ 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:

with $ b \equiv 8 g_X^2 / \lambda_{\Phi} $. Because $ \mu \propto v_\Phi^2 $, a high UV scale $ v_\Phi $ 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 $ \propto G \mu^2 $.

The dimensionless quantity $ G\mu $, 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\mu $ for the selected parameter points, and the result is shown on the $ v_{\Phi} $-$ \lambda_{\phi} $ plane in Fig. 1 (b) with the color axis ranging from $ 10^{-19} $ to $ 10^{-7} $. The positive correlation between $ G\mu $ and $ v_\Phi $ is clearly demonstrated.

A network of cosmic strings would be formed in the early universe after the spontaneous breaking of the ${U}(1)_\mathrm{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_\mathrm{e} $, a cosmic string loop of length L emits GWs with frequencies

where the integer $ n=1,2,3,\cdots $ 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\mu^2 $, the total power is given by [96]

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 $ \sim 50 $ in the two eras [97]. We will adopt this result for $ P_{n} $ and $ \Gamma = 50 $ in the calculation of the GW spectrum below.

Defining $ \mathsf{n}(L,t) \,\mathrm{d}L $ as the number density per physical volume of cosmic string loops with length L at cosmic time t in length interval $ \mathrm{d} L $, the GW energy density $ \rho_\mathrm{GW} $ induced by the cosmic string network per unit time at the emission time $ t_\mathrm{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_\mathrm{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 $ \rho_{\mathrm{GW}}(t_0) = \rho_{\mathrm{GW}}(t_\mathrm{e}) a^4(t_\mathrm{e}) $, and we obtain

where $ t_* $ represents the time when the GW emissions start. According to $ \mathrm{d}t_\mathrm{e} = - [H(z)(1+z)]^{-1} \mathrm{d}z $, 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 $ \rho_\mathrm{c} = 3H_0^2/(8\pi G) $, i.e.,

The Hubble constant is usually expressed as $ H_0 = 100h \mathrm{km}\; \mathrm{s}^{-1}\; \mathrm{Mpc}^{-1} $ with $ h=0.674 \pm 0.005 $ [91]. In order to avoid the uncertainty on the Hubble constant, one prefers to use $ \Omega_{\mathrm{GW}}(f)h^2 $. For solving Eq. (21), we need to know $ H(z) $. In a flat $ \Lambda\mathrm{CDM} $ universe, the Hubble rate is given by [98]

where $ \Omega_\mathrm{m} = 0.315 \pm 0.007 $ [91], $ \Omega_\mathrm{r} = 1.68 \times (5.38\pm 0.15) \times 10^{-5} $ [79], and $ \Omega_{\Lambda} = 1 - \Omega_\mathrm{r} - \Omega_\mathrm{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 \simeq 10^9 $ and the QCD phase transition at $ z \simeq 2\times 10^{12} $, we have [98]

The last gradient that we need to compute the energy density spectrum of the SGWB is the loop number density distribution $ \mathsf{n}(L,t) $. There are various approaches for modeling $ \mathsf{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], $ \mathsf{n}(L,t) $ is derived for any given cosmic time. Then, the loop number density produced in the radiation era can be approximated as

where $ \theta(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_{\mathrm{eq}}=51.1 \pm 0.8\; \mathrm{kyr} $ [79] represents the cosmic time at the matter-radiation equality. Therefore, we have $\mathsf{n}(L,t) = \mathsf{n}_\mathrm{r}(L, t)$ for $ t < t_\mathrm{eq} $ and $\mathsf{n}(L,t) = \mathsf{n}_\mathrm{m}(L,t) + \mathsf{n}_{\mathrm{r}\to\mathrm{m}}(L,t)$ for $ t > t_\mathrm{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 $ \mathcal{P}(L,t) $ is given by a power law in the form of

where $ \gamma \equiv 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 $ \gamma < \gamma_\mathrm{c} $ with χ changed to another parameter $ \chi_\mathrm{c} $. Here, $ \gamma_\mathrm{c} \simeq 20 (G\mu)^{1+2\chi} $ [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 $ \chi_\mathrm{c} $ and can be expressed as $ \mathsf{n}(L,t)=t^{-4}N(\gamma) $, with the function $ N(\gamma) $ approximately given by [101]

Here, $ \gamma_\mathrm{d} = -\mathrm{d}L/\mathrm{d}t \simeq \Gamma G\mu $ 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)\propto t^{\nu} $, with $ \nu=1/2 $ in the radiation era and $ \nu=2/3 $ in the matter era. The parameters C and χ can be expressed as

where $ C_\circ $ 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)_\mathrm{X}$ gauge symmetry breaking in the pNGB DM model. From the selected parameter points, we choose four benchmark points (BPs), whose $ G\mu $ varies from $ \mathcal{O}(10^{-17}) $ to $ \mathcal{O}(10^{-11}) $, to demonstrate the results. Table 1 presents the detailed information for the four BPs. For each BP, we calculate the corresponding $ \Omega_{\mathrm{GW}} h^2 $ assuming the BOS and LRS loop production models, as plotted in Fig. 2 (a) with solid and dashed lines, respectively.

	BP1	BP2	BP3	BP4
$v_S /\mathrm{GeV}$	1953	2101	548.5	1388
$v_{\Phi} /\mathrm{GeV}$	$ 1.335\times10^{13} $	$ 1.939\times10^{12} $	$ 1.969\times10^{11} $	$ 3.179\times10^{10} $
$m_{\chi} /\mathrm{GeV}$	199.8	56.26	98.16	123.1
$m_{h_2} /\mathrm{GeV}$	986.7	627.7	484.3	362.6
$m_{h_3} /\mathrm{GeV}$	$ 8.403\times10^{12} $	$ 1.469\times10^{12} $	$ 1.893\times10^{11} $	$ 8.312\times10^{9} $
$m_{Z'} / \mathrm{GeV}$	$ 7.255\times10^{11} $	$ 5.929\times10^{11} $	$ 9.661\times10^{10} $	$ 4.979\times10^{10} $
$ \lambda_{H \Phi} $	$ -6.330\times10^{-2} $	$ -3.786\times10^{-1} $	$ -1.278\times10^{-2} $	$ -6.114\times10^{-2} $
$ \lambda_{S \Phi} $	$ -2.870\times10^{-1} $	$ -5.416\times10^{-2} $	$ 2.813\times10^{-1} $	$ 3.188\times10^{-2} $
$ \lambda_{HS} $	$ 3.259\times10^{-1} $	$ 1.189\times10^{-1} $	$ -1.750\times10^{-1} $	$ 1.819\times10^{-2} $
$ s_{\varepsilon} $	$ 4.840\times10^{-3} $	$ 3.222\times10^{-1} $	$ 7.161\times10^{-2} $	$ 1.929\times10^{-3} $
$ \lambda_{\Phi} $	$ 3.97 \times10^{-1} $	$ 5.74\times10^{-1} $	$ 9.24\times10^{-1} $	$ 6.84\times10^{-2} $
$ g_X $	$ 2.72 \times10^{-2} $	$ 1.45\times10^{-1} $	$ 2.45\times10^{-1} $	$ 7.83\times10^{-1} $
$ G\mu $	$ 1.01 \times10^{-11} $	$ 1.20\times10^{-13} $	$ 1.11\times10^{-15} $	$ 1.10\times10^{-17} $
$ {\Omega}_{\chi}h^2 $	$ 0.118 $	$ 0.121 $	$ 0.120 $	$ 0.119 $
${\sigma}_{\chi N}^{\mathrm{SI} }\ /{{\rm{cm}}}^2$	$ 1.38\times10^{-86} $	$ 1.62\times10^{-86} $	$ 1.59\times10^{-82} $	$ 8.45\times10^{-77} $
$\langle {\sigma}_{\mathrm{ann} }v\rangle\ ({{\rm{cm}}^3/{\rm s}})$	$ 2.00\times10^{-26} $	$ 2.87\times10^{-29} $	$ 2.01\times10^{-26} $	$ 1.71\times10^{-26} $
$ \rho_{\mathrm{LISA}} $ (BOS)	$ 1.15\times10^{4} $	$ 1.48\times10^{3} $	$ 2.00\times10^{2} $	$ 3.97 $
$ \rho_{\mathrm{Taiji}} $ (BOS)	$ 7.26\times10^{3} $	$ 9.37\times10^{2} $	$ 1.26\times10^{2} $	$ 2.45 $
$ \rho_{\mathrm{TianQin}} $ (BOS)	$ 9.25\times10^{2} $	$ 1.15\times10^{2} $	$ 1.59\times10^{1} $	$ 5.28\times10^{-1} $
$ \rho_{\mathrm{CE}} $ (BOS)	$ 3.49\times10^{3} $	$ 4.33\times10^{2} $	$ 4.42\times10^{1} $	$ 5.48 $
$ \rho_{\mathrm{SKA}} $ (BOS)	$ 2.45\times10^{4} $	$ 4.61 $	$ 3.93\times10^{-4} $	$ 3.90\times10^{-8} $
$ \rho_{\mathrm{LISA}} $ (LRS)	$ 1.15\times10^{7} $	$ 1.38\times10^{5} $	$ 1.28\times10^{3} $	$ 4.93 $
$ \rho_{\mathrm{Taiji}} $ (LRS)	$ 7.19\times10^{6} $	$ 8.57\times10^{4} $	$ 7.95\times10^{2} $	$ 3.05 $
$ \rho_{\mathrm{TianQin}} $ (LRS)	$ 1.20\times10^{6} $	$ 1.42\times10^{4} $	$ 1.36\times10^{2} $	$ 6.48\times10^{-1} $
$ \rho_{\mathrm{CE}} $ (LRS)	$ 4.36\times10^{6} $	$ 2.18\times10^{6} $	$ 2.02\times10^{4} $	$ 2.11\times10^{2} $
$ \rho_{\mathrm{SKA}} $ (LRS)	$ 1.13\times10^{5} $	$ 2.19\times10^{1} $	$ 1.87\times10^{-3} $	$ 1.86\times10^{-7} $

Table 1.  Information for four benchmark points.

Figure 2.  (color online) Estimated GW spectra of four BPs in the pNGB DM model and sensitivity curves of various experiments. In the left panel, the spectra of the SGWB estimated for the BOS (LRS) loop production model are indicated by solid (dashed) lines. The sensitive frequency bands of PTAs, space interferometers, and ground interferometers are denoted by light blue, light yellow, and light green shaded bands, respectively. In the right panel, the upper limits from current GW experiments and the sensitivity curves of future GW experiments are indicated by solid and dashed lines, respectively, and the GW spectra of BP2 are shown for comparison.

In the BOS model, the high frequency behavior is controlled by the GW emissions during the radiation era [99]. As the integral contributed by $ \mathsf{n}_\mathrm{r}(L,t) $ tends to be a constant at high frequencies, $ \Omega_{\mathrm{GW}}h^2 $ becomes quite flat in the high frequency regime, as indicated by the solid lines in Fig. 2 (a). For smaller $ G\mu $, 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\mu $ 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 $ \gamma < \gamma_c $ regime, which significantly contribute to high frequency GWs [106, 107].

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 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}\text{–} 10^{-7}\; \mathrm{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 $ \Omega_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\text{–} 10^{3}\; \mathrm{Hz} $, we plot $ \Omega_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}\text{–} 10^{-1}\; \mathrm{Hz} $, we display $ \Omega_n h^2 $ for LISA [44], TianQin [116], and Taiji [117]. For TianQin, $ \Omega_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\mu $ 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\mu < 1.96 \times 10^{-8} $, while that for the LRS model is $ G\mu < 4.83\times 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\mu $. The 11-year NANOGrav data yield a 95% C.L. upper limit of $ G\mu < 5.3 \times 10^{-11} $ for the BOS model [109]. Utilizing the PPTA data over 15 years assuming a non-auto Hellings-Downs correlation [119], the 95% C.L. constraints are estimated to be $ G\mu < 2.88\times 10^{-11} $ for the BOS model ² and $ G\mu < 9.12\times 10^{-12} $ for the LRS model. These constraints and the selected parameter points are displayed on the $ v_{\Phi} $-$ G\mu $ plane for the BOS and LRS models in Figs. 3 (a) and 3 (b), respectively. The color axes indicate the parameter $ b = 8 g_X^2 / \lambda_{\Phi} $, which relates $ v_\Phi $ to $ G\mu $ through Eq. (13). Because of the strong correlation between $ G\mu $ and $ v_\Phi $, the parameter points align around a straight line, and the dependence on b scatters the points. We find that the parameter points with $ v_\Phi \gtrsim 5 \times 10^{13}\left(7 \times 10^{11}\right)\; \text{GeV} $ in the pNGB DM model have been excluded assuming the BOS (LRS) model for loop production.

Figure 3.  (color online) Selected parameter points projected onto the $ v_{\Phi} $-$ G\mu $ plane for the BOS (a) and LRS (b) models, with color axes corresponding to $ b = 8 g_X^2 / \lambda_{\Phi} $. The constraints on $ G\mu $ from current GW experiments are indicated by solid lines, and the expected upper limits in future experiments are denoted by dashed lines.

Below, we evaluate the sensitivity of future GW experiments according to their expected $ \Omega_{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} \leq f \leq f_{\max} $, the signal-to-noise ratio (SNR) ρ can be estimated as [111, 121]

where $ \mathcal{T} $ represents the practical observation time. If the SNR ρ turns out to be larger than an appropriate threshold, which will be set as $ \rho_{\mathrm{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 E channels have the same $ \Omega_{n}(f) $, the SNR can be evaluated as follows:

Because two channels are used, this expression has an additional factor $ \sqrt{2} $ compared with Eq. (34).

We assume that the observation time for LISA, Taiji, TianQin, or CE is $ \mathcal{T} = 1\; \text{yr} $, and that for SKA is $ \mathcal{T} = 10\; \text{yr} $ [122]. Then, the SNR ρ in these experiments is evaluated as a function of $ G\mu $ 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 1.

Figure 4.  (color online) Estimated SNRs (ρ) with respect to $ G\mu $ in future GW experiments LISA, Taiji, TianQin, CE, and SKA for the BOS (a) and LRS (b) models. The horizontal dashed lines correspond to the SNR threshold $ \rho_\mathrm{thr} = 10 $. The vertical dashed lines indicate the values of $ G\mu $ for the four BPs.

Taking the SNR threshold to be $ \rho_\mathrm{thr} = 10 $, the expected upper limits on $ G\mu $ for LISA, Taiji, TianQin, CE, and SKA are presented in Table 2. 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_\Phi $ down to $\sim 2 \times 10^{10} \text{GeV}$. For the LRS model, CE has the highest sensitivity, and the parameter points with $ v_\Phi $ down to $ \sim 5 \times 10^{9}\; \text{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_\phi $-$ \rho_\mathrm{LISA} $ plane for the BOS model and onto the $ v_\phi $-$ \rho_\mathrm{CE} $ plane for the LRS model in Fig. 5, where $ \rho_\mathrm{LISA} $ and $ \rho_\mathrm{CE} $ represent the SNRs of LISA and CE, respectively.

	LISA	Taiji	TianQin	CE	SKA
BOS model	$ 2.21\times10^{-17} $	$ 3.34\times10^{-17} $	$ 4.28\times10^{-16} $	$ 4.54\times10^{-17} $	$ 1.77\times 10^{-13} $
LRS model	$ 1.79\times10^{-17} $	$ 2.51\times10^{-17} $	$ 9.67\times10^{-17} $	$ 4.66\times10^{-19} $	$ 8.09\times 10^{-14} $

Table 2.  Expected upper limits on $ G\mu $ corresponding to $ \rho_\mathrm{thr} = 10 $.

Figure 5.  (color online) Selected parameter points projected onto the $ v_{\Phi} $-$ \rho_{\mathrm{LISA}} $ plane for the BOS model (a) and onto the $ v_{\Phi} $-$ \rho_{\mathrm{CE}} $ plane for the LRS model (b), with color axes corresponding to $ G\mu $. The horizontal dashed lines indicate the SNR threshold $ \rho_\mathrm{thr} = 10 $.

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 \times 10^{-11} < G\mu < 10^{-10} $ ($ 2 \times 10^{-11} < G\mu < 3 \times 10^{-10} $) [124]. In Fig. 6, we compare the NANOGrav confidence intervals with the selected parameter points on the $ v_\Phi $-$ G\mu $ plane and find that it is possible to explain such a suspicious signal at 95% C.L. by $2 \times 10^{13}\; \text{GeV} \lesssim v_\Phi \lesssim 2 \times 10^{14}\; \text{GeV}$ in the pNGB DM model. Of course, more data from future GW experiments are needed to reveal the nature of this anomaly.

Figure 6.  (color online) 68% (gray) and 95% (light gray) confidence intervals of $ G\mu $ for the NANOGrav anomaly compared with the selected parameter points in the pNGB DM model.

We studied the stochastic GW signals from cosmic strings induced by the spontaneous breaking of the hidden ${U}(1)_\mathrm{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_\Phi $, which characterizes the breaking scale of the $ \mathrm{U}(1)_\mathrm{X} $ gauge symmetry. Thus, DM direct detection experiments would not be able to probe χ for $ v_\Phi \gtrsim 10^5\; \text{GeV} $. Meanwhile, $ \chi\chi $ 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)_\mathrm{X} $ symmetry breaking at the UV scale $ v_\Phi $ 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_\Phi \gtrsim 5 \times 10^{13}\; ( 7 \times 10^{11})\; \text{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_\Phi $ down to $ \sim 2 \times 10^{10}\; ( 5 \times 10^{9})\; \text{GeV} $. Note that the bound on the DM lifetime has basically excluded a UV scale $ v_\Phi $ lower than $ 10^9\; \text{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_\Phi $ 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_{\Phi} $ is assumed to be a constant. Nonetheless, if the temperature corrections to the effective potential of the scalar fields are considered, $ v_{\Phi} $ would depend on the temperature at the epoch when the spontaneous breaking of the $ {U}(1)_\mathrm{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.

We thank Shi-Qi Ling for helpful discussions.