^†^†institutetext: Institute for Fundamental Theory, Department of Physics, University of Florida,
Gainesville, FL 32611, USA

Gauged Global Strings

Xuce Niu Wei Xue and Fengwei Yang xuce.niu@ufl.edu weixue@ufl.edu fengwei.yang@ufl.edu

Abstract

We investigate the string solutions and cosmological implications of the gauge ${\rm U(1)_{Z}}\,\times$ global ${\rm U(1)_{PQ}}$ model. With two hierarchical symmetry-breaking scales, the model exhibits three distinct string solutions: a conventional global string, a global string with a heavy core, and a gauge string as a bound state of the two global strings. This model reveals rich phenomenological implications in cosmology. During the evolution of the universe, these three types of strings can form a Y-junction configuration. Intriguingly, when incorporating this model with the QCD axion framework, the heavy-core global strings emit more axion particles compared to conventional axion cosmic strings due to their higher tension. This radiation significantly enhances the QCD axion dark matter abundance, thereby opening up the QCD axion mass window. Consequently, axions with masses exceeding $\sim 10^{-5}\,{\rm eV}$ have the potential to constitute the whole dark matter abundance. Furthermore, in contrast to conventional gauge strings, the gauge strings in this model exhibit a distinctive behavior by radiating axions.

1 Introduction

Cosmic strings, soliton solutions in field theory, arise when loops in the vacuum manifold cannot be contracted to a point [1, 2]. In the context of symmetry-breaking, where a symmetry group $G$ is spontaneously broken to a subgroup $H$ , $G\to H$ , the vacuum manifold is a quotient space, $G/H$ . Cosmic string solutions are connected to the non-trivial first homotopy group $\pi_{1}(G/H)$ . The simplest cosmic string originates from a complex scalar field, denoted as $\Phi$ , with a global $\rm U(1)$ symmetry. Its vacuum expectation value (VEV), $\langle\Phi\rangle=\frac{1}{\sqrt{2}}f_{a}$ spontaneously breaks its global $\rm U(1)$ symmetry. A cosmic string along the $z$ -direction exhibits a scalar field configuration in polar coordinates $(r,\theta)$ ,

\Phi(r,\theta)=\frac{1}{\sqrt{2}}f_{a}{\rm e}^{i\theta}\,,\quad r\to\infty

(1)

with a winding number $1$ . The energy per unit length of the string is estimated by integrating radially from the inverse of the scalar mass $m^{-1}$ to a long distance cutoff $L$ ,

\mu\simeq 2\pi\int_{m^{-1}}^{L}{\rm d}r\frac{1}{r}|{\partial_{\theta}}\Phi(r,\theta)|^{2}=\pi f_{a}^{2}\ln(mL)

(2)

The tension of a global string is predominantly contributed by the gradient term outside of the string core. Alternatively, considering an Abelian Higgs model where the $\rm U(1)$ symmetry is a gauge symmetry, the cosmic strings, known as gauge strings, exhibit finite tension concentrated inside the string core. While the gauge strings can have similar scalar configurations as global strings, the energy from the gradient term outside the core regime is eliminated by a gauge configuration.

Cosmic strings generically form from a phase transition in the early universe through the Kibble mechanism [1]. The broken symmetry is restored at the high temperature of the universe, $T\gtrsim f_{a}$ . A phase transition occurs when the temperature falls around $f_{a}$ . During this transition, the vacuum expectation value $\langle\Phi\rangle$ turns on, and the symmetry is spontaneously broken. The phase of $\langle\Phi\rangle$ in the vacuum manifold is chosen at random beyond the correlation length of the phase transition, resulting in the formation of cosmic strings. Subsequently, the interactions between cosmic strings lead to a few long strings per Hubble volume, entering a string scaling regime.

Cosmic string networks in the universe provide intriguing signatures, and their detection is an exciting direction to probe UV physics. The existence of cosmic strings influences the large-scale structure of the universe [3, 4, 5, 6, 7]. The current constraint on the cosmic string tension arises from analyzing the angular power spectrum of cosmic microwave background (CMB) [8, 9, 10, 11]. Furthermore, a cosmic string loop cannot survive in the universe forever. Gauge strings emit gravitational radiations from loop oscillations, detectable by current and future gravitational wave detectors [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22]. Additionally, considering interactions between cosmic string and particles, such as photons, opens new possibilities for detection in the CMB or astrophysical observations[23, 24, 25, 26, 27, 28, 29, 30].

One well-motivated cosmic string is the QCD axion string, as a global string. The QCD axions [31, 32, 33, 34, 35, 36, 37, 38, 39] provide a physically intriguing solution, solving the strong CP problem [40, 41], and serving as a dark matter candidate [37, 38, 39]. Axion string emission to axions can be a dominant contributor to dark matter abundance, though the emitted axion energy spectrum is still an unsettled question [42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59]. Recently, numerical simulations [60, 61, 62, 63, 64] make efforts to address axion dark matter abundance from the QCD axion string radiation, predicting an axion mass of $\mathcal{O}(10-100)\,\mu{\rm eV}$ to explain the full dark matter abundance. The QCD axions motivate world-wide efforts for their search, not only focusing on a specific mass range but also planning to cover a broad mass range in the future experiments [65, 66, 67, 68, 69, 70, 71, 72, 73]. Also, various mechanisms, including parametric resonance [74, 75, 76], anharmonicity effect [77, 78, 79], domain walls decays [80, 81, 82, 58, 83, 84, 85, 86], axion production during a kination era [87], and kinetic misalignment [88], allow for different axion masses to explain the dark matter abundance. Here, we provide another mechanism from cosmic string decays.

In this paper, we introduce a simple model having two complex fields, $\Phi_{1}$ and $\Phi_{2}$ , with a symmetry of gauge $\rm U(1)_{Z}\,\times$ global $\rm U(1)_{PQ}$ . From this model, we obtain two kinds of global string solutions and one gauge string solution from the spontaneous symmetry breaking, as the three most energetically favorable string configurations. The two global strings arise from the winding around $\Phi_{1}$ and $\Phi_{2}$ , respectively. Assuming that the vacuum expectation values of the two fields are hierarchical, we observe that one global string is heavier than another. The lighter global string tension is close to the conventional QCD axion string. The third string is a bound gauge string formed by combining the two global strings together. With the hybrid string solutions, we investigate their cosmological implications. The system displays unique cosmological dynamics and signatures. This is evident in the formation and evolution of string networks, as well as in the radiation emitted by string loops. In the context of QCD axion physics, heavier global strings emit more axions due to their larger tension, contributing more to the axion dark matter abundance. The large tension has been used to emulate the behavior of axion strings in a cosmological simulation [89]. Additionally, for a gauge string from the Abelian Higgs model, gravitational radiation is the dominant channel for the comic strings to lose energy. However, the gauge string in this model couples to massless Goldstone modes, raising the intriguing question of which particles the gauge strings radiate dominantly.

The structure of this paper is organized as follows. Section˜2 introduces the ${\rm U(1)_{Z}}$ $\times\,{\rm U(1)_{PQ}}$ model, and then we embed the model into the QCD axion framework. Following this, the section provides the string solutions for this model, which are further validated through numerical analysis in section˜3. The cosmological implications of the model are explored in section˜4, and our findings are summarized in section˜5.

2 ${\rm U}(1)_{\rm Z}\times{\rm U}(1)_{\rm PQ}$

In this section, we introduce a model having a gauge symmetry ${\rm U}(1)_{\rm Z}$ and a global symmetry ${\rm U}(1)_{\rm PQ}$ . The dynamics of this model are driven by two scalar fields, $\Phi_{1}$ and $\Phi_{2}$ , which will break the symmetry sequentially at distinct energy scales in the cosmic evolution. Furthermore, we consider the intriguing possibility of integrating this model with the QCD axion framework, so that the full model merges the Abelian gauge symmetry with the original QCD axion models [34, 33, 36, 35, 90].

2.1 The model

We consider the gauge symmetry and the global symmetry, ${\rm U}(1)_{\rm Z}\times{\rm U}(1)_{\rm PQ}$ , both of which are broken by two complex scalar $\Phi_{1}$ and $\Phi_{2}$ . The dynamics are described by their gauge-invariant and renormalizable Lagrangian density, which takes the form

{\cal L}=-\frac{1}{4}Z_{\mu\nu}Z^{\mu\nu}+D_{\mu}\Phi_{1}^{\dagger}D^{\mu}\Phi_{1}-\frac{\lambda_{1}}{4}\left(|\Phi_{1}|^{2}-\frac{v_{1}^{2}}{2}\right)^{2}+D_{\mu}\Phi_{2}^{\dagger}D^{\mu}\Phi_{2}-\frac{\lambda_{2}}{4}\left(|\Phi_{2}|^{2}-\frac{v_{2}^{2}}{2}\right)^{2}\,.

(3)

$Z^{\mu\nu}$ represents the field strength of the $\rm U(1)_{Z}$ gauge boson $Z_{\mu}$ , defined as $Z^{\mu\nu}=\partial^{\mu}Z^{\nu}-\partial^{\nu}Z^{\mu}$ . The scalar fields interact with the gauge boson $Z_{\mu}$ with the gauge coupling $e$ through the covariant derivative terms, as expressed by

D_{\mu}\Phi_{i}=(\partial_{\mu}-ieZ_{\mu})\Phi_{i},~~i=1,2\,.

(4)

We assign that $\Phi_{1}$ carries charges $(+1,+1)$ under ${\rm U}(1)_{\rm Z}\times{\rm U}(1)_{\rm PQ}$ , while $\Phi_{2}$ has charges $(+1,-1)$ . It is important to note that the charges cannot be entirely determined by the above Lagrangian density. They should be determined by the scalars’ coupling to other particles or a UV theory. The charges of $\Phi_{1}$ and $\Phi_{2}$ can be extended to other values, and the implication for cosmic strings will be discussed in sections˜2.3 and 4.1. For simplicity, we consider two independent Mexican-hat potentials for the scalar fields with self-couplings $\lambda_{1}$ and $\lambda_{2}$ while assuming that interactions between the two scalars, namely $|\Phi_{1}|^{2}|\Phi_{2}|^{2}$ , are either negligible or zero. In the vacuum, $\Phi_{1}$ and $\Phi_{2}$ acquire non-zero expectation values, denoted as $v_{1}$ and $v_{2}$ ,

\langle\Phi_{1}\rangle=\frac{v_{1}}{\sqrt{2}},~\langle\Phi_{2}\rangle=\frac{v_{2}}{\sqrt{2}}\,.

(5)

This spontaneous symmetry breaking leads to the sequential breaking of the $\rm U(1)_{Z}\times U(1)_{PQ}$ symmetry in the cosmic evolution, particularly when $v_{1}\gg v_{2}$ . The first phase transition results in a non-zero VEV of $\langle\Phi_{1}\rangle$ , leading to the breaking of $\rm U(1)_{Z}$ and giving the gauge boson $Z^{\mu}$ mass, $m_{Z}=ev_{1}$ . After the second phase transition, the non-zero VEV of $\langle\Phi_{1}\rangle$ and $\langle\Phi_{2}\rangle$ further breaks the global symmetry $\rm U(1)_{PQ}$ , with the gauge boson mass increasing, $m_{Z}=e\sqrt{v_{1}^{2}+v_{2}^{2}}$ .

We conveniently parametrize the perturbations of $\Phi_{1}$ and $\Phi_{2}$ using real scalar fields $\tilde{\phi}_{1}(x)$ , $\tilde{\phi}_{2}(x)$ , $\pi_{1}(x)$ and $\pi_{2}(x)$ ,

\begin{split}\Phi_{1}(x)&=\frac{1}{\sqrt{2}}\phi_{1}(x)\,e^{i\alpha_{1}}=\frac{1}{\sqrt{2}}(v_{1}+\tilde{\phi}_{1}(x))\,e^{i\pi_{1}(x)/v_{1}}\,,\\ \Phi_{2}(x)&=\frac{1}{\sqrt{2}}\phi_{2}(x)\,e^{i\alpha_{2}}=\frac{1}{\sqrt{2}}(v_{2}+\tilde{\phi}_{2}(x))\,e^{i\pi_{2}(x)/v_{2}}\,,\end{split}

(6)

where $\alpha_{1}$ and $\alpha_{2}$ represent the rotation angle of $\Phi_{1}$ and $\Phi_{2}$ , respectively. We identify the axion, denoted as $a(x)$ , by ensuring that it is orthogonal to the would-be Goldstone boson $\pi_{z}$ or the longitudinal mode of the gauge boson $Z^{\mu}$ . Using the expression of the $\rm U(1)_{Z}$ current $J_{z}^{\mu}$ and $\pi_{z}$ -to-vacuum matrix element, $\langle{\rm vac}|J_{z}^{\mu}(0)|\pi_{z}(p)\rangle=ip^{\mu}v$ , we identify the would-be Goldstone boson of $\rm U(1)_{Z}$ , which takes the form

\pi_{z}(x)=\frac{1}{v}\left(v_{1}\,\pi_{1}+v_{2}\,\pi_{2}\right)\,,

(7)

where $v=\sqrt{v_{1}^{2}+v_{2}^{2}}$ . Requiring orthogonality between the axion field and $\pi_{z}$ yields the expression for $a(x)$

a(x)=\frac{1}{v}\left(v_{2}\,\pi_{1}-v_{1}\,\pi_{2}\right)\,.

(8)

According to the symmetry of $\Phi_{1}$ and $\Phi_{2}$ , we can express $\pi_{1}$ and $\pi_{2}$ in terms of the rotation angles of $\rm U(1)_{Z}$ and $\rm U(1)_{PQ}$ in the vacuum manifold, specifically

\pi_{1}/v_{1}=\alpha_{Z}+\alpha_{\rm PQ}\,,\quad\pi_{2}/v_{2}=\alpha_{Z}-\alpha_{\rm PQ}\,.

(9)

Here, $\alpha_{Z}$ and $\alpha_{\rm PQ}$ represent the rotation angles of $\rm U(1)_{Z}$ and $\rm U(1)_{PQ}$ , respectively. Consequently, we can rewrite the axion fields in terms of the $\rm U(1)_{PQ}$ rotation angle $\alpha_{\rm PQ}$ ,

a(x)=v_{a}\alpha_{\rm PQ}\,,\quad v_{a}=\frac{2v_{1}v_{2}}{\sqrt{v_{1}^{2}+v_{2}^{2}}}\,.

(10)

The parameter $v_{a}$ tells us the magnitude of the vacuum expectation value that spontaneously breaks $\rm U(1)_{PQ}$ .

2.2 QCD axion

The ${\rm U}(1)_{\rm Z}\times{\rm U}(1)_{\rm PQ}$ model provides an elegant framework for realization, and it becomes particularly intriguing when integrated into QCD axion models.

One possible approach, based on the KSVZ model, involves introducing a vector-like heavy fermion that can be decomposed into its left-handed and right-handed components, ${\cal Q}={\cal Q}_{L}+{\cal Q}_{R}$ . This fermion resides in the fundamental representation of the Standard Model $\rm SU(3)_{c}$ color symmetry and is a singlet under the $\rm U(1)_{Z}$ symmetry and other Standard Model symmetries. Moreover, the fermion carries a chiral charge under the PQ symmetry, with ${\cal Q}_{L}$ having a $+1$ charge and ${\cal Q}_{R}$ a $-1$ charge. Since it is a singlet in the gauge $\rm U(1)_{Z}$ symmetry, the fermion does not introduce any additional $\rm U(1)_{Z}$ anomaly, preserving the gauge symmetry anomaly-free.

Taking into account the symmetry of the fermion $\cal Q$ , we construct the Lagrangian

{\cal L}=i\bar{\cal Q}\not{D}{\cal Q}-\frac{y}{\Lambda}\left(\Phi_{1}\Phi_{2}^{*}\bar{\cal Q}_{L}{\cal Q}_{R}+h.c.\right)\,,

(11)

where $\Lambda$ represents the UV cutoff of the theory. The expectation values of $\Phi_{1}$ and $\Phi_{2}$ yield the fermion mass, $m_{\cal Q}=yv_{1}v_{2}/(2\Lambda)$ . By integrating out the heavy scalars and employing eqs.˜6 and 8, we derive the effective Lagrangian for axion couplings to fermions

{\cal L}\supset-m_{\cal Q}e^{i2a/v_{a}}\bar{\cal Q}_{L}{\cal Q}_{R}+h.c.\,.

(12)

Subsequently, we can deduce axion interactions with the Standard Model particles, including axion couplings to gluons. The derivation and results parallel those of the KSVZ model. Through a field-dependent axial transformation,

\mathcal{Q}_{L}\rightarrow e^{ia/v_{a}}\mathcal{Q}_{L}\,,\quad\mathcal{Q}_{R}\rightarrow e^{-ia/v_{a}}\mathcal{Q}_{R},

(13)

the heavy fermion becomes disentangled from the axion, introducing an axion-gluon coupling term

{\cal L}\supset\frac{g_{s}^{2}}{16\pi^{2}}\frac{a}{v_{a}}G^{a}_{\mu\nu}\tilde{G}^{a}_{\mu\nu}\,=\frac{g_{s}^{2}}{32\pi^{2}}\frac{a}{f_{a}}G^{a}_{\mu\nu}\tilde{G}^{a}_{\mu\nu}\,,

(14)

where $g_{s}$ is the coupling constant of $\rm SU(3)_{c}$ . The axion-gluon coupling defines the axion decay constant $f_{a}$ . Notably, in this model, $v_{a}=2f_{a}$ . The ratio of $v_{a}$ and $f_{a}$ , often referred to as the domain wall number, is represented by $\frac{v_{a}}{f_{a}}$ . However, in this case, the domain wall number for cosmic string solutions is $N=1$ . This arises because, in the vacuum manifold of $\rm U(1)_{PQ}$ , the minimal winding is achieved by choosing the angle $\alpha_{\rm PQ}$ from $0$ to $\pi$ . $\alpha_{\rm PQ}=0$ and $\alpha_{\rm PQ}=\pi$ are gauge-equivalent, and $0$ can return to $\pi$ through the gauge group $\rm U(1)_{Z}$ (see fig.˜1). A similar method of counting domain wall number is observed in the PQWW model [40, 31, 32], which has three domain walls in an axion string, but $\frac{v_{a}}{f_{a}}=6$ .

Refer to caption — Figure 1: The cross-section of the vacuum manifold of ${\rm U(1)_{Z}\times U(1)_{PQ}}$ , where the two black dots represent the angle $\alpha_{\rm PQ}=0$ and $\alpha_{\rm PQ}=\pi$ .

Alternatively, we can construct a model by introducing two sets of quarks, denoted as ${\cal Q}_{1}$ and ${\cal Q}_{2}$ , as discussed by Barr and Seckel [91]. These quarks are color-triplets under $\rm SU(3)_{c}$ and are assigned the following chiral charges under $\rm U(1)_{Z}\times U(1)_{PQ}$ : ${\cal Q}_{1L}$ has charges $(1/2,1/2)$ , ${\cal Q}_{1R}$ has charges $(-1/2,-1/2)$ , ${\cal Q}_{2L}$ has charges $(-1/2,-1/2)$ , and ${\cal Q}_{2R}$ has charges $(1/2,1/2)$ . This charge assignment ensures that $\rm U(1)_{Z}$ remains anomaly-free. Furthermore, it leads to ${\cal Q}_{1}$ interacting with $\Phi_{1}$ and ${\cal Q}_{2}$ interacting with $\Phi_{2}$ though the Lagrangian density,

{\cal L}=\Phi_{1}\bar{\cal Q}_{1L}{\cal Q}_{1R}+\Phi_{2}\bar{\cal Q}_{2L}{\cal Q}_{2R}+h.c.\,.

(15)

Following a similar procedure, we arrive at the same axion gluon coupling shown in eq.˜14.

Let us comment on other extensions of these models. In both realizations, we can introduce $N_{f}$ flavors of the heavy quarks, thereby suppressing the axion decay constant $f_{a}$ as a factor of $1/N_{f}$ . This increases domain wall numbers in the cosmic string solutions. Additionally, in the Standard Model, $\rm U(1)_{B-L}$ is anomaly-free and could serve as a gauge symmetry. Another intriguing possibility is that $\rm U(1)_{Z}$ symmetry can be identified with the $\rm U(1)_{B-L}$ [92].

2.3 String solutions

In this section, we study the string solutions of ${\rm U(1)_{Z}}\times{\rm U(1)_{PQ}}$ analytically. Initially, we analyze the string solutions beyond their core regions to pinpoint the three most energetically favorable string configurations. Then, we extend our analysis to include string solutions with arbitrary winding numbers $(j,k)$ . Finally, we extrapolate our findings to encompass generic gauge charges.

The gradient energy of $(1,0)$ , $(0,1)$ and $(1,1)$ strings

Before delving into the total energy per unit length of strings, we focus on the gradient term for $(1,0)$ , $(0,1)$ and $(1,1)$ strings. Examining the gradient energy allows us to distinguish the $(1,0)$ , $(0,1)$ strings are global strings, while $(1,1)$ strings are gauge strings. Also, these three string configurations are the lightest ones, playing a major role in cosmology.

For the winding number $(1,0)$ , the classical configurations of the scalar and gauge fields outside the string cores take the form

\Phi_{1}=\frac{1}{\sqrt{2}}v_{1}\,e^{i\theta}\,,\quad\Phi_{2}=\frac{1}{\sqrt{2}}v_{2}\,,\quad Z_{\mu}=c\,\partial_{\mu}\theta\,,\quad r\to\infty

(16)

Here, $c$ is a normalization of the gauge field, determined by minimizing the gradient energy of strings. The gradient energy per unit length is then calculated by integrating the 2D cross-section of the string, yielding

\begin{split}\mu_{k,(1,0)}&=\int_{0}^{2\pi}{\rm d}\theta\int_{\delta}^{L}{\rm d}r\,r\left(|(\frac{1}{r}\partial_{\theta}-ieZ_{\theta})\Phi_{1}|^{2}+|(-ieZ_{\theta})\Phi_{2}|^{2}\right)\\ &=\pi\ln(\frac{L}{\delta})\left[v_{1}^{2}(1-ec)^{2}+v_{2}^{2}(ec)^{2}\right]\end{split}

(17)

where we take the core size as $\delta$ . Since $m_{Z}$ is typically lighter than the scalar field mass $m_{1}$ and $m_{2}$ , we can choose the core size $\delta\simeq m_{Z}^{-1}$ .¹¹1 We consider the gradient term outside the core region. There are also gradient corrections from the radius of $m_{1}^{-1}$ to $m_{Z}^{-1}$ , being neglected here but included in the full tension calculation. The contribution is about $\pi v_{1}^{2}\ln(m_{1}/m_{Z})$ . It can be dominant when $v_{1}\gg v_{2}$ . The parameter $c$ is determined by minimizing the tension,

c=\frac{1}{e}\frac{v_{1}^{2}}{v_{1}^{2}+v_{2}^{2}}\,.

(18)

Therefore, the $(1,0)$ tension outside the core is

\mu_{k,(1,0)}=\pi\frac{v_{1}^{2}v_{2}^{2}}{v_{1}^{2}+v_{2}^{2}}\ln(\frac{L}{\delta})=\pi f_{a}^{2}\ln(\frac{L}{\delta})

(19)

Next, for the $(0,1)$ string tension outside the string core, the configurations take the form $\Phi_{1}=\frac{1}{\sqrt{2}}v_{1}$ , $\Phi_{2}=\frac{1}{\sqrt{2}}v_{2}e^{i\theta}$ and $Z_{\mu}=c\,\partial_{\mu}\theta$ . $c=\frac{1}{e}\frac{v_{2}^{2}}{v_{1}^{2}+v_{2}^{2}}$ minimizes the energy of the $(0,1)$ string. The result should be the same as $(1,0)$ string by switching $v_{1}$ and $v_{2}$ . The tension of $(0,1)$ string is the same as the $(1,0)$ string,

\mu_{k,(0,1)}=\pi f_{a}^{2}\ln(\frac{L}{\delta})

(20)

There is an intuitive way to understand that the $(1,0)$ and $(0,1)$ strings share the same tension outside the core regimes. The $(1,0)$ string configuration outside core is equivalent to a $(0,-1)$ string through a gauge transformation, $\alpha_{Z}\to\alpha_{Z}+\theta$ . Consequently, $(1,0)$ and $(0,1)$ are string and anti-string.²²2We thank Pierre Sikivie for the discussion on this point.

The $(1,1)$ string results from combing a $(1,0)$ string with its anti-string, $(0,1)$ string, implying zero tension outside the core. This result aligns with our study of field configurations. Using the same energy minimization procedure, we find that the $Z^{\mu}$ configuration with $c=1/e$ cancels the $\Phi_{1}$ and $\Phi_{2}$ gradient terms simultaneously. Hence, we conclude that the tension outside the core of $(1,1)$ string is zero, $\mu_{k,(1,1)}=0$ . Consequently, the $(1,1)$ string is identified as a gauge string.

The tension of $(j,k)$ strings

By considering the time-independent stable field configuration and choosing the gauge $Z_{0}=0$ , the full string tension is obtained by the two-dimensional spatial integral of the Hamiltonian density,

\mu=\int d^{2}x\left[\frac{B_{z}^{2}}{2}+|D_{i}\Phi_{1}|^{2}+|D_{i}\Phi_{2}|^{2}+\frac{\lambda_{1}}{4}\left(|\Phi_{1}|^{2}-\frac{v_{1}^{2}}{2}\right)^{2}+\frac{\lambda_{2}}{4}\left(|\Phi_{2}|^{2}-\frac{v_{2}^{2}}{2}\right)^{2}\right],

(21)

where $B_{z}=\frac{1}{r}\frac{d(rZ_{\theta}(r))}{dr}$ is the $z$ -component of the “magnetic” field ${\bf B}=\nabla\times{\bf Z}$ .

For a $(j,k)$ string, we choose the Ansatz for the fields:

$\displaystyle\Phi_{1}(r,\theta)$	$\displaystyle=$	$\displaystyle\frac{v_{1}}{\sqrt{2}}f_{1}(r)e^{ij\theta},$
$\displaystyle\Phi_{2}(r,\theta)$	$\displaystyle=$	$\displaystyle\frac{v_{2}}{\sqrt{2}}f_{2}(r)e^{ik\theta},$
$\displaystyle Z_{\theta}(r)$	$\displaystyle=$	$\displaystyle c\frac{g(r)}{r},$	(22)

where $f_{1}(r),~f_{2}(r)$ , and $g(r)$ are the profile functions and where we have $c\equiv\frac{1}{e}\frac{jv_{1}^{2}+kv_{2}^{2}}{v_{1}^{2}+v_{2}^{2}}$ to minimize the gradient energy (deviation in appendix A). They need to satisfy the boundary conditions when $j\neq 0$ and $k\neq 0$ :

			$\displaystyle f_{1,2}(r\rightarrow 0)\rightarrow 0,~f_{1,2}(r\rightarrow\infty)\rightarrow 1,$		(23)
			$\displaystyle g(r\rightarrow 0)\rightarrow 0,~g(r\rightarrow\infty)\rightarrow 1.$		(23)

Note that the boundary conditions of the profile function $f_{1,2}(r)$ for $r\rightarrow 0$ are not necessarily equal to 0 when the corresponding winding number is 0. Substitute the above Ansatz into eq.˜21, we find the string tension in terms of $f_{1},~f_{2},~g$ :

\begin{split}\mu=\int d^{2}x\bigg{[}&\frac{c^{2}g^{\prime 2}}{2r^{2}}+\frac{1}{2}f_{1}^{2}v_{1}^{2}\left(\frac{(j-ceg)^{2}}{r^{2}}+\frac{f_{1}^{\prime 2}}{f_{1}^{2}}\right)+\frac{1}{2}f_{2}^{2}v_{2}^{2}\left(\frac{(k-ceg)^{2}}{r^{2}}+\frac{f_{2}^{\prime 2}}{f_{2}^{2}}\right)\\ &+\frac{\lambda_{1}}{16}(f_{1}^{2}-1)^{2}+\frac{\lambda_{2}}{16}(f_{2}^{2}-1)^{2}\bigg{]},\end{split}

(24)

where ^′ represents the derivative with respect to the radius $r$ .

Here, we introduce an assumption about the profile functions, allowing us to derive analytical estimates for string tensions. More precise string profile functions are determined through numerical solutions to the equations of motion, as presented in section˜3. We define critical radii, denoted as $r_{1,c}$ , $r_{2,c}$ , where the string profiles $f_{1}(r)$ and $f_{2}(r)$ respectively reach asymptotic values at large radii. Also, we identify $r_{c}$ as the radius of the magnetic flux. A Heaviside function $\Theta(r)$ is introduced to simplify the profile functions, taking the form as

\displaystyle f_{1}(r)

\displaystyle=

\displaystyle\Theta(r-r_{1,c}),\quad f_{2}(r)=\Theta(r-r_{2,c}),\quad g(r)=\Theta(r-r_{c}).

(25)

Considering naturalness of the scalar mass, we take $r_{c}>r_{1,c},\,r_{2,c}$ . We further assume a uniform magnetic field $B_{z}(r)=B_{0}$ and evaluate the magnetic flux at the radius of $r_{c}$ ,

B_{0}\times\pi r_{c}^{2}=\int_{0}^{2\pi}rd\theta Z_{\theta}=2\pi c\,,

(26)

yielding

B_{0}=\frac{2c}{r_{c}^{2}}\,.

(27)

Substituting the above expressions into eq.˜24, we get the string tension with the winding number $(j,k)$ as

	$\displaystyle\mu$	$\displaystyle=$	$\displaystyle\frac{2\pi c^{2}}{r_{c}^{2}}+\pi v_{1}^{2}j^{2}\ln\left(\frac{r_{c}}{r_{1,c}}\right)+\pi v_{2}^{2}k^{2}\ln\left(\frac{r_{c}}{r_{2,c}}\right)+\pi\left(j-k\right)^{2}\frac{v_{1}^{2}v_{2}^{2}}{v_{1}^{2}+v_{2}^{2}}\ln\left(\frac{L}{r_{c}}\right)$		(28)
		$\displaystyle+$	$\displaystyle\frac{\pi}{16}\left(\lambda_{1}r_{1,c}^{2}v_{1}^{4}(1-\delta_{j0})+\lambda_{2}r_{2,c}^{2}v_{2}^{4}(1-\delta_{k0})\right)\,.$		(28)

We use $\int_{0}^{L}=\int_{0}^{r_{1(2),c}}+\int_{r_{1(2),c}}^{r_{c}}+\int_{r_{c}}^{L}$ when performing the radial integral. We introduce the Kronecker delta functions, $\delta_{j0}$ , $\delta_{k0}$ , so that our result is valid for generic $j$ and $k$ , including either $j,k=0$ . We can find the relation between the core sizes and the mass of the fields by minimizing the string tension in eq.˜28 with respect to $r_{c}$ , $r_{1,c}$ , and $r_{2,c}$ individually. With the two Higgs masses $m_{i}=\sqrt{\lambda_{i}/2}v_{i}~(i=1,2)$ and the gauge boson mass $m_{Z}=e\sqrt{v_{1}^{2}+v_{2}^{2}}$ , the radii are given as

\displaystyle r_{1,c}=\frac{2j}{m_{1}},\quad r_{2,c}=\frac{2k}{m_{2}},\quad r_{c}=\frac{2}{m_{Z}}.

(29)

Hence, the string tension can be written in terms of the mass of the massive fields

	$\displaystyle\mu$	$\displaystyle=$	$\displaystyle\frac{\pi e^{2}(jv_{1}^{2}+kv_{2}^{2})^{2}}{2m_{Z}^{2}}+\pi v_{1}^{2}j^{2}\ln\left(\frac{m_{1}}{jm_{Z}}\right)+\pi v_{2}^{2}k^{2}\ln\left(\frac{m_{2}}{km_{Z}}\right)$		(30)
		$\displaystyle+$	$\displaystyle\pi\left(j-k\right)^{2}\frac{e^{2}v_{1}^{2}v_{2}^{2}}{m_{Z}^{2}}\ln\left(\frac{m_{Z}L}{2}\right)+\frac{\pi}{2}\left(j^{2}v_{1}^{2}+k^{2}v_{2}^{2}\right)\,.$		(30)

We define a string global charge $q_{\rm global}\equiv j-k$ . When $q_{\rm global}=0$ , the IR logarithmic divergence of the string tension $\ln(m_{Z}L/2)$ vanishes, implying gauge string solutions.

The full tensions of the $(1,0)$ , $(0,1)$ , and $(1,1)$ strings take the form

\mu_{(1,0)}=\frac{\pi v_{1}^{4}}{2(v_{1}^{2}+v_{2}^{2})}+\pi v_{1}^{2}\ln\left(\frac{m_{1}}{m_{Z}}\right)+\pi\frac{v_{1}^{2}v_{2}^{2}}{v_{1}^{2}+v_{2}^{2}}\ln\left(\frac{m_{Z}L}{2}\right)+\frac{\pi}{2}v_{1}^{2},

(31)

\mu_{(0,1)}=\frac{\pi v_{2}^{4}}{2(v_{1}^{2}+v_{2}^{2})}+\pi v_{2}^{2}\ln\left(\frac{m_{2}}{m_{Z}}\right)+\pi\frac{v_{1}^{2}v_{2}^{2}}{v_{1}^{2}+v_{2}^{2}}\ln\left(\frac{m_{Z}L}{2}\right)+\frac{\pi}{2}v_{2}^{2},

(32)

\mu_{(1,1)}=\frac{\pi(v_{1}^{2}+v_{2}^{2})}{2}+\pi v_{1}^{2}\ln\left(\frac{m_{1}}{m_{Z}}\right)+\pi v_{2}^{2}\ln\left(\frac{m_{2}}{m_{Z}}\right)+\frac{\pi}{2}(v_{1}^{2}+v_{2}^{2}).

(33)

Through the above calculation, we confirm the disappearance of the IR divergent gradient energy term for a $(1,1)$ string, indicating that it is a gauge string. In contrast, $(1,0)$ and $(0,1)$ strings exhibit global-like characteristics, since they have non-zero $|q_{\rm global}|=1$ . The tension difference between gauge strings and the two global strings gives rise to the binding energy, represented by the expression

\mu_{(1,0)}+\mu_{(0,1)}-\mu_{(1,1)}=\frac{\pi v_{1}^{2}v_{2}^{2}}{v_{1}^{2}+v_{2}^{2}}\left[2\ln\left(\frac{m_{Z}L}{2}\right)-1\right].

(34)

For $m_{Z}L\gtrsim 3.3$ , when the $(1,1)$ string string tension is less than the combined tension of $(1,0)$ and $(0,1)$ strings, it indicates an attractive force between a $(1,0)$ string and a $(0,1)$ string. Hence, a $(1,1)$ string is more stable. Moreover, the hierarchical symmetry breaking scale $v_{1}>v_{2}$ implies $\mu_{(1,0)}>\mu_{(0,1)}$ , due to a larger energy stored inside the core of a $(1,0)$ string compared to the one of a $(0,1)$ string.

Generic gauge charges $q_{1}$ and $q_{2}$

$\Phi_{1}$ and $\Phi_{2}$ can carry generic gauge charges $q_{1}$ and $q_{2}$ , different from $+1$ that we assign before. The covariant derivative of $\Phi_{1}$ and $\Phi_{2}$ is, therefore, $D_{\mu}\Phi_{1}=(\partial_{\mu}-iq_{1}eZ_{\mu})\Phi_{1}$ and $D_{\mu}\Phi_{2}=(\partial_{\mu}-iq_{2}eZ_{\mu})\Phi_{2}$ . Similarly, we deduce the string tension, taking the form

	$\displaystyle\mu$	$\displaystyle=$	$\displaystyle\frac{\pi e^{2}(jq_{1}v_{1}^{2}+kq_{2}v_{2}^{2})^{2}}{2m_{Z}^{2}}+\pi v_{1}^{2}j^{2}\ln\left(\frac{m_{1}}{jm_{Z}}\right)+\pi v_{2}^{2}k^{2}\ln\left(\frac{m_{2}}{km_{Z}}\right)$		(35)
		$\displaystyle+$	$\displaystyle\pi\left(jq_{2}-kq_{1}\right)^{2}\frac{e^{2}v_{1}^{2}v_{2}^{2}}{m_{Z}^{2}}\ln\left(\frac{m_{Z}L}{2}\right)+\frac{\pi}{2}\left(j^{2}v_{1}^{2}+k^{2}v_{2}^{2}\right),$		(35)

where the string global charge $q_{\rm global}\equiv jq_{2}-kq_{1}$ .

3 Numerical Study of String Solutions

In this section, we present the results of our numerical investigation into the string profiles and tensions with three string solutions: $(1,0)$ and $(0,1)$ global strings, and $(1,1)$ gauge string. We confirm that inside the core, the string tension of $(1,0)$ , $(0,1)$ global strings scale with the square of the symmetry-breaking scale. Outside of the core, it exhibits logarithmic divergence. In contrast, the tension of the $(1,1)$ gauge string predominantly originates from its core region.

String profile

We employ the multiparameter shooting method to solve the profile functions for different string configurations with winding numbers $(j,k)$ . Starting with the Ansatz, as shown in eq. (22), the equations of motion of $\Phi_{1}$ , $\Phi_{2}$ and $Z^{\mu}$ leads to the following differential equations:

	$\displaystyle f_{1}^{\prime\prime}+\frac{f_{1}^{\prime}}{r}-\frac{f_{1}}{r^{2}}(j-ceq_{1}g)^{2}+\frac{\lambda_{1}}{4}v_{1}^{2}(f_{1}^{2}-1)f_{1}=0,$
	$\displaystyle f_{2}^{\prime\prime}+\frac{f_{2}^{\prime}}{r}-\frac{f_{2}}{r^{2}}(k-ceq_{2}g)^{2}+\frac{\lambda_{2}}{4}v_{2}^{2}(f_{2}^{2}-1)f_{2}=0,$
	$\displaystyle cg^{\prime\prime}-\frac{cg^{\prime}}{r}-e^{2}v_{1}^{2}f_{1}^{2}(cq_{1}g-j)-e^{2}v_{2}^{2}f_{2}^{2}(cq_{2}g-k)=0,$		(36)

where ^′ denotes the derivative with respect to $r$ . These non-linear differential equations can be solved numerically by imposing boundary conditions at the origin, while the shooting method gives $f_{1}\rightarrow 1$ , $f_{2}\rightarrow 1$ , $g\rightarrow 1$ at large radius. The profile functions would approach zero at the origin for non-zero winding numbers due to the absence of singularity. However, for zero winding numbers, the corresponding profile functions can be non-zero values at the origin. A more detailed discussion of the profile functions for the three configurations is provided in appendix˜B. We consider a hierarchy in scale with $v_{1}=3v_{2}$ . The remaining free parameters are chosen as $e=0.1,\lambda_{1}=10e^{2},\lambda_{2}=40e^{2},q_{1}=q_{2}=1$ . To simplify the analysis, we introduce the dimensionless parameter $\rho=ev_{1}r$ . The profile functions for $(1,0)$ , $(0,1)$ , and $(1,1)$ strings are shown in fig.˜2. More details of the numerical solutions are given in table˜1 of appendix˜B. For zero winding numbers, we find that the profile functions of the corresponding scalar fields are non-zero at the origin.³³3 The study in [93] exhibits a distinct density profile when the winding number is zero, where the scalar field configuration approaches zero at the origin. Discrepancies between these two results may arise from variations in the boundary conditions near $r=0$ or model parameters.

String tension

We numerically compute the string tension, confirming that the $(1,0)$ and $(0,1)$ strings are global strings, while the $(1,1)$ string is a local (gauge) string with no tension increase outside the core regime. Furthermore, we compare the string tension among the three configurations. At small scales (near the string cores), the $(1,0)$ and $(1,1)$ strings exhibit larger tension than the $(0,1)$ string with $v_{1}=3v_{2}$ . At large scales ( $r>\frac{2}{m_{Z}}$ ), the tensions of $(1,0)$ and $(0,1)$ strings show logarithmic divergence due to the contribution from the kinetic energy.

We plot the string tension and its components explicitly for the three configurations in fig.˜3. The components are kinetic energy per unit length $\mu_{k}=\int d^{2}x(|D_{i}\Phi_{1}|^{2}+|D_{i}\Phi_{2}|^{2})$ , energy from gauge field $\mu_{Z}=\int d^{2}x\frac{B_{z}^{2}}{2}$ , and potential energy $\mu_{U}=\int d^{2}xV(\Phi_{1},\Phi_{2})$ . Notably, the $(1,0)$ and $(0,1)$ string tension exhibit infrared logarithm divergence $\pi f_{a}^{2}\ln(m_{Z}L/2)$ , while the $(1,1)$ string tension converges. The string separation $L$ becomes crucial in determining which string is heavier between the $(1,0)$ and the $(1,1)$ . For small separations, the $(1,0)$ string has a lower tension compared to the $(1,1)$ string, while for larger separations, the $(1,0)$ string becomes heavier, as shown in the last plot in fig.˜3.

The numerical results of string tension align with the analytic results presented in section˜2.3. The tension contribution inside the string core is of order $\pi v_{1}^{2}\ln(\frac{m_{1}}{m_{Z}})$ for the $(1,0)$ string and $\pi v_{2}^{2}\ln(\frac{m_{2}}{m_{Z}})$ for the $(0,1)$ string. We parameterize the $(1,0)$ and $(0,1)$ string tension as a combination of string core tension plus a string tension outside the core, introducing free parameters $x_{10}$ and $x_{01}$ ,

	$\displaystyle\mu_{(1,0)}(r)$	$\displaystyle=$	$\displaystyle\mu_{(1,0)}(r<\delta)+\mu_{(1,0)}(r>\delta)$		(37)
		$\displaystyle\simeq$	$\displaystyle x_{10}\,\pi v_{1}^{2}\ln(\frac{m_{1}}{m_{Z}})+\pi f_{a}^{2}\ln(\frac{rm_{Z}}{2}),$		(37)

	$\displaystyle\mu_{(0,1)}(r)$	$\displaystyle=$	$\displaystyle\mu_{(0,1)}(r<\delta)+\mu_{(0,1)}(r>\delta)$		(38)
		$\displaystyle\simeq$	$\displaystyle x_{01}\,\pi v_{2}^{2}\ln(\frac{m_{2}}{m_{Z}})+\pi f_{a}^{2}\ln(\frac{rm_{Z}}{2}),$		(38)

where we determine $x_{10}=1.85$ and $x_{01}=3.19$ by fitting the string tension outside of the strings, as shown in fig.˜4. Since $x_{10}$ and $x_{01}$ are of order $\mathcal{O}(1)$ , we conclude that eq.˜31 and eq.˜32 provide the correct order of string tensions.

4 Cosmological Implication

The $\rm U(1)_{Z}\times U(1)_{PQ}$ model allows for the existence of two hierarchical symmetry-breaking scales, denoted as $v_{1}$ and $v_{2}$ ( $v_{1}>v_{2}$ ). The occurrence of the two $\rm U(1)$ symmetry-breaking in the thermal history of the universe has profound implications, giving rise to the formation of various comic strings. These cosmic strings consist of $(1,0)$ , $(0,1)$ global strings, as well as $(1,1)$ gauge strings. In this section, we delve into the complexities of their formation, evolution, and radiation, which are generic and applicable to the $\rm U(1)_{Z}\times U(1)_{PQ}$ model and its extensions. Moreover, we scrutinize their significance in the context of the QCD axion. The presence of this new string network can substantially impact axion cosmology and the calculation of axion relic abundance when the $\rm U(1)_{Z}\times U(1)_{PQ}$ symmetry is incorporated into QCD axion framework. Finally, we examine the radiation emitted by $(1,1)$ gauge strings, addressing the interesting question of whether they predominantly emit axions or gravitational waves.

4.1 Formation of string network

Two distinct phase transitions occurred during the evolution of the universe, sequentially breaking the $\rm U(1)_{Z}\times U(1)_{PQ}$ . We assume that the phase transitions are second-order, eliminating the possibility of bubble nucleation during the transition. The first phase transition leads to the formation of gauge string, denoted as the $\pi_{1}$ string, in contrast to those generated after the second phase transition when the symmetry is completely broken. In the second phase transition, the $\pi_{1}$ strings undergo modifications due to the scalar field $\Phi_{2}(x)$ configuration, and, in addition, $(0,1)$ strings form. The interaction between these stings will establish a network of strings as the universe evolves.

During the first phase transition, the complex scalar field $\Phi_{1}(x)$ acquires a non-zero VEV, $\langle\Phi_{1}(x)\rangle=v_{1}/\sqrt{2}$ , while the second scalar field $\Phi_{2}(x)$ remains at its potential minimum, $\Phi_{2}(x)=0$ . As a second-order phase transition, this occurs around the temperature $T_{c,1}\sim\sqrt{\frac{6\lambda_{1}}{\lambda_{1}+3e^{2}}}v_{1}$ . The $\pi_{1}$ string with a winding number $1$ is therefore formed by the Kibble mechanism [1]. The $\pi_{1}$ strings with higher winding numbers, $|j|>1$ , are rarely formed. Since $\Phi_{2}(x)=0$ remains at its origin, $\pi_{1}$ strings can be regarded as $(1,n)$ string, with the integer $n$ determined by the second phase transition when $\Phi_{2}$ acquires a VEV. To minimize the energy, the gauge field $Z^{\mu}$ compensates for the gradient of the phase of $\Phi_{1}(x)$ , and thus the $\pi_{1}$ string is a $\rm U(1)$ gauge string. Following their production, these strings inevitably collide, either passing through each other or breaking and reconnecting with other strings. Due to these interactions, after a few Hubble times, the $\pi_{1}$ string network enters a scaling regime where a few long $\pi_{1}$ strings persist, and the correlation length of the long strings becomes comparable to the horizon size $H^{-1}$ .

When the temperature of the universe continuously drops below $T_{c,2}\sim\sqrt{\frac{6\lambda_{2}}{\lambda_{2}+3e^{2}}}v_{2}$ , another second-order phase transition occurs. During this transition, the complex scalar field $\Phi_{2}(x)$ acquires a VEV $\langle\Phi_{2}(x)\rangle=v_{2}/\sqrt{2}$ . The correlation length of the $\Phi_{2}$ at the beginning of the phase transition can be estimated by the Ginzburg length [94],

\xi_{2}\sim(\lambda_{2}T_{c,2})^{-1}.

(39)

Since this correlation length is much shorter than the separation of the $\pi_{1}$ long strings, $H^{-1}$ , within the correlation length, $\Phi_{1}(x)$ can be considered homogeneous. Consequently, we expect the formation of $(0,1)$ string independently of the pre-existence of the $\pi_{1}$ strings. In the case of $\pi_{1}$ string, within a distance $\sim\xi_{2}$ , the presence of non-trivial winding number ( $j=1$ ) in the $\Phi_{1}$ fields influences the $\Phi_{2}$ and gauge field configurations. We assume that the field configurations adjust themselves with a slowly changing vacuum during the second-order the phase transition to energetically favorable solutions. In this scenario, the $\Phi_{2}$ and gauge field reach new configurations, and the integer $n$ in the previous $(1,n)$ string is determined by minimizing the energy. Therefore, we compare the tension $(1,0)$ to $(1,1)$ strings. The difference is estimated using eq.˜31 and eq.˜33,

\mu_{(1,1)}-\mu_{(1,0)}\simeq\pi v_{2}^{2}\ln\left(\frac{m_{2}}{m_{Z}}\right)-\pi v_{2}^{2}\ln\left(\frac{m_{Z}L}{2}\right)=\pi v_{2}^{2}\ln\left(\frac{2m_{2}}{m_{Z}^{2}L}\right)\,.

(40)

The $(1,0)$ strings often is lighter than $(1,1)$ strings, particularly considering $v_{2}L\sim{\cal O}(1)$ during the second phase transition (see fig.˜5). Therefore, we expect the predominant formation of $(1,0)$ strings. However, it is worth noting that the phase transition may exhibit more complex dynamics, or the adiabatic argument may fail, potentially leading to the simultaneous production of $(1,0)$ , $(0,1)$ and gauge string.

The evolution of $(1,0)$ and $(0,1)$ strings is more complicated than the evolution of conventional global string networks. The complexity arises due to the formation of $(1,1)$ gauge string segments, known as Y-junctions, when the two types of strings encounter and intercommute (see fig.˜6). The formation of Y-junctions can be understood as the result of attractive forces between the two strings or $\mu_{(1,0)}+\mu_{(0,1)}-\mu_{(1,1)}>0$ , as shown in eq.˜34. Although the analysis of Y-junctions requires string simulations, these Y-junctions cannot transform the whole $(1,0)$ and $(0,1)$ strings into a $(1,1)$ gauge string. First, the intersecting probability of $(1,0)$ and $(0,1)$ strings is not frequent in an expanding universe. Even when intersecting, they are easily unzipped by the high velocities of the strings. Furthermore, due to the long-range interaction between a global string and its anti-string, an attractive force from the other side of the string loop balances the Y-junction.

These Y-junctions also occur in cosmic super-strings, with lattice simulations to analyze their evolution. These simulations have explored two local $\rm U(1)$ strings [95, 96, 97, 98], global $\rm SU(2)/Z_{3}$ strings [99], and others [100, 101, 102]. They provide insights into how Y-junctions influence the string network. They consistently observe that the evolution of the string network tends towards a scaling regime. We anticipate similar behavior in the $\rm U(1)_{Z}\times U(1)_{PQ}$ model, but cosmological simulations for this model is encouraged to validate this conclusion.

In summary, we reveal that $(1,0)$ and $(0,1)$ are produced during the second phase transition by considering the adiabatic condition. The $(1,0)$ and $(0,1)$ strings form Y-junctions during the evolution of the string network, and is expected to enter a scaling regime according to [95, 96, 97, 98, 99]. They remain in the network with some fraction of the Y-junctions.

As a caveat, more complicated dynamics may happen, requiring dedicated cosmological simulations. In the case that $(1,1)$ strings are generated in the second phase transition, or afterward, the evolution of $(1,1)$ strings is independent of $(1,0)$ and $(0,1)$ strings. In the string network, the $(1,1)$ string independently enters a scaling regime. Since they have the smallest tension compared to other strings with higher winding numbers, the $(1,1)$ strings do not bind with the $(1,0)$ or $(0,1)$ strings. Moreover, having a generic gauge charge $q_{1}$ and $q_{2}$ may lead to gauge string formation during the evolution and the string network dynamics are more complicated. The cosmological simulation conducted in [103] revealed that the fraction of gauge strings could be on the same order of magnitude as that of global strings with $q_{1}=1,q_{2}=4$ and $v_{1}=4v_{2}$ . The dynamics of the string network, such as the reconnection of string bundles, are highly intricate. The long-range force between the $(1,0)$ and $(0,1)$ strings plays an important role in forming the bound state of gauge strings when the charge ratio $q_{2}/q_{1}$ are large.

4.2 QCD axion abundance

In this section, we explore the interesting possibility of the Goldstone modes in the $\rm U(1)_{Z}\times U(1)_{PQ}$ being the well-motivated QCD axion. We consider that the symmetry breaking of $\rm U(1)_{Z}\times\rm U(1)_{PQ}$ occurs after inflation. In this new scenario, we calculate the axion production. The contribution from axion vacuum realignment is the same as the predictions from other post-inflationary scenarios of the QCD axions. However, the axion production resulting from the decay of the defect network, including both string and domain wall decay, can considerably alter the axion abundance projection.

Let us present the axion production from cosmic strings and domain walls sequentially:

Cosmic strings

Following [104], we analyze the emission of axions from both $(1,0)$ and $(0,1)$ strings. We exclude the contribution from $(1,1)$ strings for two primary reasons. First, according to section˜4.1, we do not anticipate that the $(1,1)$ string number in the network is dominant. Second, even considering the presence of $(1,1)$ strings, the axion radiation from $(1,1)$ string loops is less than the radiation from $(1,0)$ . This is because, during the QCD phase transition, the string tension of $(1,1)$ strings is smaller than that of $(1,0)$ global strings. This can be verified using eq.˜40 by setting $L\sim H^{-1}(T=\Lambda_{\rm QCD})$ . The $(1,1)$ string radiation is further discussed in section˜4.3. While a more complicated string network may arise, introducing uncertainties in the number density of $(1,1)$ strings, the findings in this section remain valid under the assumption that the scaling regime is attained.

When the string loops collapse, they convert their entire energy into axion particles. Consequently, the tension of the strings and the energy spectrum of the emitted axions play crucial roles in determining the axion number density. The tension of $(1,0)$ strings surpasses that of $(0,1)$ strings due to their heavier cores, resulting in a greater abundance of axion dark matter. The final axion density resulting from string decay also relies on the characteristics of the energy spectrum. Some attempts have been made to address this question through numerical simulations [47, 52, 61, 105, 63, 64], but as far as we know, this issue has not been definitively resolved. Furthermore, the presence of a heavier core in the $(1,0)$ string could potentially impact the energy spectrum. Therefore, we take an agnostic approach regarding the energy spectrum and consider two scenarios: one where axion emission is equally important across all energy scale (Scenario A), and another where emitted axions are dominated in the infrared modes (Scenario B).

As $(1,0)$ and $(0,1)$ strings enter the scaling regime, the correlation length becomes causality-limited $L\sim t$ , leading to the string tension of a long string taking the form

\mu_{(1,0)}(t)\simeq\pi v_{1}^{2}\ln\left(\frac{m_{1}}{m_{Z}}\right)+\pi f_{a}^{2}\ln\left(\frac{m_{Z}t}{2}\right),

(41)

\mu_{(0,1)}(t)\simeq\pi v_{2}^{2}\ln\left(\frac{m_{2}}{m_{Z}}\right)+\pi f_{a}^{2}\ln\left(\frac{m_{Z}t}{2}\right)\,.

(42)

Here, we neglect the contributions of magnetic energy and potential energy since they are subdominant and confirmed by our numerical study as shown in eqs. (37), (38). Notably, the first term in eq.˜41 introduces a substantial correction to the string tension, owing to $v_{1}^{2}>f_{a}^{2}=\frac{v_{1}^{2}v_{2}^{2}}{v_{1}^{2}+v_{2}^{2}}$ . In the case of $(0,1)$ strings, their tension closely resembles that of the standard QCD axion strings, as $v_{2}\simeq f_{a}$ when $v_{1}>v_{2}$ .

Considering that the string network enters a scaling regime, the number of long strings within one Hubble patch is of the order of ${\cal O}(1)$ . The energy density of the long string within one Hubble volume at time $t$ can thus be expressed as

\rho_{\rm str,i}(t)=N_{i}\frac{\mu_{i}(t)}{t^{2}},

(43)

where the subscript $i$ represents the two types of strings, $(1,0)$ or $(0,1)$ , and $N_{i}$ is the number of long strings in a Hubble patch, which is roughly $\mathcal{O}(1)$ .

To consider the decay of strings into axions, we track the evolution of the density of long strings using the following equation

\frac{d\rho_{\rm str,i}}{dt}+3H\rho_{\rm str,i}-H\rho_{\rm str,i}=-\frac{d\rho_{{\rm str,i}\rightarrow a}}{dt}\,,

(44)

where $\frac{d\rho_{{\rm str,i}\rightarrow a}}{dt}$ is the rate at which the energy density of strings gets converted into axions at time $t$ . The term $3H$ accounts for the dilution due to the expansion of the universe, while $-H$ arises from the stretching of the strings. The string decay enhances the number density of axions, resulting in the following form for the number density of axions from string decays

\frac{dn_{a}^{\rm str}}{dt}+3Hn_{a}^{\rm str}=\sum_{i}\frac{1}{\bar{\omega}_{i}(t)}\frac{d\rho_{{\rm str,i}\rightarrow a}}{dt}\simeq\sum_{i}\frac{1}{\bar{\omega}_{i}(t)}\frac{N_{i}\,\mu_{i}(t)}{t^{3}}\,.

(45)

To achieve the last approximate form, we use the string density evolution function, eq.˜44, and neglect the time dependence in $N_{i}\mu_{i}(t)$ . The average energy of axions, $\frac{1}{\bar{\omega}_{i}(t)}$ , is calculated through the energy spectrum of the emitted axions, $\frac{dE_{i}}{d\omega}$ ,

\frac{1}{\bar{\omega}_{i}(t)}=\frac{\int\frac{1}{\omega}\frac{dE_{i}}{d\omega}{d\omega}}{\int\frac{dE_{i}}{d\omega}d\omega}\,.

(46)

The spectrum of emitted axions has some theoretical uncertainties. To address these uncertainties, we introduce two distinct scenarios as mentioned. Further, we assume that the $(1,0)$ and $(0,1)$ share the same spectrum shape to encompass these uncertainties.

In Scenario A, the spectrum peaks in the infrared region at $\omega=2\pi/t$ , i.e., $\frac{dE}{d\omega}\propto\delta(\omega-2\pi t^{-1})$ . In Scenario B, the spectrum takes the form $\frac{dE}{d\omega}\propto\frac{1}{\omega}$ , with the frequency range $\omega\in[2\pi t^{-1},2\pi\delta^{-1}]$ , where $\delta$ represents the UV cutoff. For a $(1,0)$ string, $\delta\simeq 1/m_{1}$ ; for a $(0,1)$ string, $\delta\simeq 1/m_{2}$ . Although we choose a delta function in the spectrum in Scenario A, the analysis should encompass the case of $\frac{dE}{d\omega}\propto\frac{1}{\omega^{q}}$ with $q>1$ , corresponding to an IR-dominant spectrum. Substituting $H=\frac{1}{2t}$ , we solve eq.˜45 and find that the number density of axions from strings can be expressed as

n_{a}^{\rm str}(t)\simeq\frac{1}{t^{3/2}}\sum_{i}\int_{t_{*}}^{t}dt^{\prime}\frac{1}{\bar{\omega}_{i}(t^{\prime})}\frac{N_{i}\,\mu_{i}(t^{\prime})}{t^{\prime 3/2}}\,.

(47)

The lower limit of the time integration, $t_{*}$ , is taken as the second phase transition time, $t(T_{c,2})$ , but the final result is not sensitive to the initial time.

In Scenario A, taking $\frac{1}{\bar{\omega}}=\frac{t}{2\pi}$ and neglecting the ${\cal O}(1)$ coefficients, we estimate the axion number density as follows

n_{a}^{\rm str}(t)\sim\frac{t^{1/2}}{t^{3/2}}\sum_{i}\frac{N_{i}\,\mu_{i}(t)}{2\pi}\,.

(48)

In Scenario B, where $\frac{1}{\bar{\omega}}=\frac{t}{2\pi}\left[\ln\left(\frac{t}{\delta}\right)\right]^{-1}$ , the axion number density is estimated as

n_{a}^{\rm str}(t)\sim\frac{t^{1/2}}{t^{3/2}}\sum_{i}\frac{N_{i}\,\mu_{i}(t)}{2\pi\ln\left(\frac{t}{\delta_{i}}\right)}\,.

(49)

In the two equations above, the factor of $\frac{1}{t^{3/2}}$ accounts for the dilution due to the expansion of the universe. Following axion production, their number density scales as $\frac{1}{t^{3/2}}$ in a radiation-dominant universe. The factor of $t^{1/2}$ indicates that axions produced from strings at later times dominate the axion population. Note that these formulas are valid only before the time $t_{1}$ , representing the horizon crossing time for the axion field, $3H(t_{1})=m_{a}(t_{1})$ . Beyond this crossing time, the IR cutoff needs to be replaced by the axion mass, resulting in the dominant axion contribution around $t_{1}$ . Therefore, to compute the axion density originating from the string decays, we need to evaluate the number density in eq.˜48 or eq.˜49 at the time $t_{1}$ . As expected, the number density depends on the tension of strings and the energy spectrum.

Domain walls

The QCD axion potential exhibits degenerate minima as the QCD axion becomes massive, leading to kinks between neighboring potential minima and the presence of domain walls. Here, we examine the QCD axion model discussed in section˜2.2, which has strings connected by a single domain wall with $N=1$ . When the axion mass turns on, each string is connected to a nearby anti-string via a domain wall. The domain wall solution is present in $(1,0)$ and $(0,1)$ strings. However, for $(1,1)$ strings, the PQ symmetry rotation angle $\alpha_{PQ}=\alpha_{1}-\alpha_{2}=0$ , resulting in no domain wall associated with the $(1,1)$ string.

The tension of domain walls when connected to $(1,0)$ and $(0,1)$ strings is given by

\sigma(t)=8m_{a}(t)f_{a}^{2},

(50)

The axion mass $m_{a}(t)$ is time- or temperature-dependent. At high temperature ( $T\gtrsim 1\,{\rm GeV}$ ), the mass can be approximated as

m_{a}(T)\simeq 4\times 10^{-9}\,{\rm eV}\left(\frac{10^{12}\,{\rm GeV}}{f_{a}}\right)\left(\frac{{\rm GeV}}{T}\right)^{4},

(51)

with the dilute instanton gas approximation [106], and at temperatures below 100 MeV

m_{a}\simeq 6\,\mu{\rm eV}\left(\frac{10^{12}\,\rm GeV}{f_{a}}\right)\,.

(52)

Once the axion mass turns on at time $t_{1}$ , the strings bound to the wall have an IR cutoff $L\sim m_{a}^{-1}$ , comparable to the wall thickness. Hence the energy stored in a $(1,0)$ or $(0,1)$ string per unit length can be expressed as

\begin{split}\mu_{(1,0)}(t)&\simeq\pi v_{1}^{2}\ln\left(\frac{m_{1}}{m_{Z}}\right)+\pi f_{a}^{2}\ln\left(\frac{m_{Z}}{2m_{a}(t)}\right),\\ \mu_{(0,1)}(t)&\simeq\pi v_{2}^{2}\ln\left(\frac{m_{2}}{m_{Z}}\right)+\pi f_{a}^{2}\ln\left(\frac{m_{Z}}{2m_{a}(t)}\right).\end{split}

(53)

After the axion mass turns on at $t_{1}$ , the domain walls become bound to the strings. The evolution of the domain walls goes through several stages. Initially, the dynamics of the wall-string system are governed by the strings until time $t_{2}$ , and the scaling solution of long strings remains; subsequently, starting from time $t_{2}$ , the domain walls begin pulling the strings, and the dynamics of the system become dominated by the walls. Finally, at time $t_{3}$ , the domain walls and strings collapse, emitting axion particles.

The transition time $t_{2}$ is determined when the energy stored in the strings is comparable to the energy in the walls,

\frac{E_{\sigma}(t_{2})}{E_{\mu}(t_{2})}\simeq\frac{8m_{a}(t_{2})f_{a}^{2}t_{2}^{2}}{\mu(t_{2})t_{2}}=1.

(54)

The $(1,0)$ and $(0,1)$ strings have different transition times due to their distinct string tension $\mu(t)$ , as shown in eq.˜53. For walls connected to the light $(0,1)$ strings, the situation aligns with the standard scenario [47], which approximates $t_{2}\sim t_{1}$ . However, when the walls are connected to the heavier $(1,0)$ strings, where $\mu_{(1,0)}>\mu_{(0,1)}$ , it takes longer to achieve the energy balance between the wall and the string, giving us $t_{2}>t_{1}$ .

Beyond $t_{2}$ , the energy stored in the wall surpasses that stored in the strings bound to the wall due to their different scaling with time. The walls pull the strings and accelerate them. The strings eventually unzip the wall into several smaller walls until the wall’s size is comparable to $1/m_{a}$ . Finally, the walls collapse, emitting axions at $t_{3}$ .

The energy density of the walls at time $t_{1}$ was estimated to be approximately $0.7$ per horizon volume [47]. This energy density scales as $\sigma(t)/t$ when the dynamics are governed by the strings. After $t_{2}$ , the domain wall area does not change significantly within one Hubble volume, so the average energy density is scaled by the volume or the inverse cubic power of the cosmological scale factor $t^{-3/2}$ ,

\rho_{\text{DW}}(t)\sim 0.7\frac{\sigma(t)}{t_{2}}\left(\frac{t_{2}}{t}\right)^{3/2},\quad t_{2}<t<t_{3}\,.

(55)

After time $t_{3}$ , the axions produced during the collapse of the wall-string system are boosted. We define an average Lorentz factor $\gamma$ as the ratio of the energy density to the axion mass at $t_{3}$ , denoted as $\gamma\equiv\frac{\langle\omega\rangle}{m_{a}(t_{3})}$ . The number density of the produced axions is scaled by the volume

n_{a}^{\rm DW}(t)\sim\frac{\rho_{\rm DW}(t_{3})}{\langle\omega\rangle}\left(\frac{t_{3}}{t}\right)^{3/2}\sim\frac{6}{\gamma}\frac{f^{2}_{a}}{t_{2}}\left(\frac{t_{2}}{t}\right)^{3/2}\,.

(56)

In the second equality, substituting eq.˜55, we find that the dependence on time $t_{3}$ drops out of the estimate of $n_{a}^{\rm DW}(t)$ . According to simulation results in Ref. [47], $\gamma\sim 60$ , though there is significant uncertainty on this value.

We consider the $(1,0)$ string and its wall-string system here since it has large tension and contributes significantly to the axion density, but in the numerical analysis, we include contributions from both global strings. We compare the axion number density produced by the $(1,0)$ string bound to the wall using eq.˜56 with the one produced by the $(1,0)$ string decay in Scenario A and B. The domain wall contribution to axion number density does not depend on $t_{3}$ , as shown in eq.˜56, allowing us to compare the string decay contribution with the string-wall collapse contribution at $t_{2}$ directly. In Scenario A, we compare $n_{a}^{\rm str}$ with $n_{a}^{\rm DW}$ at $t_{2}$ , and set the number of $(1,0)$ long strings is $N_{(1,0)}\sim 1$ ,

\begin{split}n_{a}^{\rm str}(t_{2})&=n_{a}^{\rm str}(t_{1})\left(\frac{t_{1}}{t_{2}}\right)^{3/2}\sim\frac{\mu}{m_{a}(t_{1})t_{2}^{2}}\left(\frac{t_{2}}{t_{1}}\right)^{1/2}\\ &>\frac{\mu}{m_{a}(t_{1})t_{2}^{2}}\sim\gamma n_{a}^{\rm DW}(t_{2}).\end{split}

(57)

In the last line, we use eq.˜54 to replace $\mu(t_{2})$ . Considering $t2>t_{1}$ and $\gamma\gtrsim 1$ , the axion production from the wall decay is sub-dominant in Scenario A. In Scenario B, the axion spectrum is harder, leading to a smaller number density of axion from strings. Numerically, we find that the domain wall contribution remains sub-dominant in most regions of the model’s parameter space even when $\gamma=1$ .

There is one more complication stemming from Y-junctions. The bound $(1,1)$ strings do not connect with the domain walls. The collapse of global strings with the Y-junctions is expected to close the $(1,1)$ strings and leave $(1,1)$ string loops in the universe. However, the $(1,1)$ strings, even when they eventually emit axions, represent a sub-dominant contribution to the axion dark matter abundance, and thus, they are neglected in our analysis.

Results

The total axion energy density at present, $t=t_{0}$ , is given by the sum of the three contributions: misalignment, string decays, and domain wall decays,

\rho_{a,0}=\rho^{\rm vac}_{a}(t_{0})+m_{a}n_{a}^{\rm str}(t_{0})+m_{a}n_{a}^{\rm DW}(t_{0})\,.

(58)

$n_{a}^{\rm str}$ is the number density of axions from both $(1,0)$ and $(0,1)$ string decay. There is an uncertainty of the number of long strings $N_{i}$ in a Hubble patch. Here, we just set $N_{i}=1$ as an estimation. $n_{a}^{\rm DW}$ is the number density of axions from domain wall decay bounded by $(1,0)$ and $(0,1)$ strings, and $\rho^{\rm vac}_{a}=\Omega_{a}^{\rm vac}\rho_{c,0}$ is the energy density of axions produced by the misalignment mechanism, with $\rho_{c,0}$ the current critical density of the universe, for $f_{a}<2\times 10^{15}~$ GeV,

\Omega_{a}^{\rm vac}h^{2}\sim 2\times 10^{4}\left(\frac{f_{a}}{10^{16}\text{ GeV}}\right)^{7/6}\langle\theta_{a}^{2}\rangle,

(59)

where $\langle\theta_{a}^{2}\rangle=\pi^{2}/3$ is the average value of the square of the misalignment angle.

If one considers axion to explain the 100% of the relic abundance of dark matter observed today, this will require

\Omega_{a}h^{2}\equiv\frac{\rho_{a,0}}{\rho_{c,0}}h^{2}=\Omega_{\rm DM}h^{2}\sim 0.12.

(60)

Figure˜7 summarizes the three contributions to axion energy density as a function of $v_{2}$ for both Scenario A and B. It confirms that the decay from the $(1,0)$ string with the heavy core is the dominant contribution to dark matter relic abundance for most regions of $v_{2}$ .

In the KSVZ model, the axion-photon coupling is linked to the axion mass through the relation $g_{a\gamma\gamma}=2.0\times 10^{-16}C_{a\gamma}(m_{a}/\mu{\rm eV})\,{\rm GeV}^{-1}$ , with $C_{a\gamma}=-1.92$ . Given the results obtained above for making QCD axions to explain the 100% of the dark matter relic abundance, we show the axion-photon coupling versus the axion mass predicted by our model in fig.˜8 (Top), see the magenta (Scenario A) and cyan (Scenario B) bands. The current ADMX experiments [107, 108, 109, 110] exclude the QCD axion mass range $2-4\,\mu{\rm eV}$ , while the globular cluster observations [111] put an upper bound on the axion mass around 0.1 eV.

We do not analyze mass below $\mu{\rm eV}$ , as they are associated with the pre-inflation scenario, where strings do not form. According to a recent axion string simulation [64], the axion mass is found to be within the range of $m_{a}\in(40,180)\,\mu{\rm eV}$ . In contrast, our proposed mechanism opens a large mass window for QCD axions as 100% cold dark matter. This offers motivation for the upcoming haloscope experiments [65, 66, 67, 68, 69, 70, 71, 72, 73] to probe axions across a broad mass spectrum.

4.3 $(1,1)$ gauge string radiation: gravitons or axions?

For macroscopic gauge strings, the production of massive particles is significantly suppressed, while the emission of massless gravitons becomes the dominant channel for these strings to lose energy. However, the $(1,1)$ gauge string not only couple to massless gravitons but also to axions. This section aims to address whether the $(1,1)$ gauge strings predominantly decay into axions or gravitons.

To comprehend the interaction between strings and axions, we start with the Lagrangian density of the two scalar fields, as shown in eq.˜3, and investigate the axion interactions with classical configurations of $Z^{\mu}(x)$ , $\Phi_{1}(x)$ and $\Phi_{2}(x)$ . Subsequently, we perform a gauge transformation, $Z^{\mu}\to Z^{\mu}+\frac{1}{e}\partial_{\mu}\alpha_{Z}$ . This transformation eliminates the quadratic term of $Z^{\mu}\partial_{\mu}a$ and $Z^{\mu}\partial_{\mu}\pi_{z}$ . The Lagrangian density thus takes the form

{\cal L}=-\frac{1}{4}Z_{\mu\nu}Z^{\mu\nu}+\frac{1}{2}e^{2}(\phi_{1}^{2}+\phi_{2}^{2})Z_{\mu}^{2}-g(\phi_{1},\phi_{2})\,e\,Z^{\mu}\partial_{\mu}a+\frac{1}{2}\,f(\phi_{1},\phi_{2})\,(\partial_{\mu}a)^{2}\,.

(61)

In this equation, we introduce the functions $g(\phi_{1},\phi_{2})$ and $f(\phi_{1},\phi_{2})$ to simplify the notation,

g(\phi_{1},\phi_{2})=f_{a}\frac{\phi_{1}^{2}}{v_{1}^{2}}-f_{a}\frac{\phi_{2}^{2}}{v_{2}^{2}}\,,

(62)

f(\phi_{1},\phi_{2})=\frac{g^{2}(\phi_{1},\phi_{2})+\frac{1}{f_{a}^{2}}\phi_{1}^{2}\phi_{2}^{2}}{\phi_{1}^{2}+\phi_{2}^{2}}\,.

(63)

When considering the vacuum expectation values of $\phi_{1}$ and $\phi_{2}$ , $g(v_{1},v_{2})=0$ and $f(v_{1},v_{2})=1$ . This implies that there is no interaction between axion and classical fields outside string cores. Equation˜61 already reveals that $(1,1)$ strings interact with axions through the term $g(\phi_{1},\phi_{2})\,e\,Z^{\mu}\partial_{\mu}a$ , owing to the non-trivial field configurations of $Z^{\mu}(x)$ , $\phi_{1}(x)$ , $\phi_{2}(x)$ . One example of these field configurations of $(1,1)$ string is shown in fig.˜3, where the classical fields gradually change in the core size of $\sim 1/m_{Z}$ .

To calculate the radiation power of $(1,1)$ strings to axions, we can employ the Kalb-Ramond field $B^{\mu\nu}$ [112, 113, 114]. A duality relationship exists between real massless field, axion, and the two-form antisymmetry gauge field, $B^{\mu\nu}$ , given by

\partial_{\mu}a=\frac{1}{2}\epsilon_{\mu\nu\alpha\beta}\partial^{\nu}B^{\alpha\beta}\,

(64)

This relation is satisfied beyond the string core. However, in the presence of $\phi_{1}(x)$ , $\phi_{2}(x)$ and $Z^{\mu}(x)$ , this duality relationship is modified. Based on the field equation of $a(x)$ derived from eq.˜61,

\partial^{\mu}(f(\phi_{1},\phi_{2})\,\partial_{\mu}a-g(\phi_{1},\phi_{2})\,eZ_{\mu})=0\,,

(65)

the relation between $a$ and $B^{\mu\nu}$ should be replaced with

f(\phi_{1},\phi_{2})\,\partial_{\mu}a-g(\phi_{1},\phi_{2})\,eZ_{\mu}=\frac{1}{2}\epsilon_{\mu\nu\alpha\beta}\partial^{\nu}B^{\alpha\beta}\,

(66)

Substituting axion with $B^{\alpha\beta}$ or its gauge-invariant field strength $H_{\mu\nu\lambda}$ in the Lagrangian density, we obtain the following expression,

\begin{split}{\cal L}&=-\frac{1}{4}Z_{\mu\nu}Z^{\mu\nu}+\frac{1}{2}e^{2}(\phi_{1}^{2}+\phi_{2}^{2})Z_{\mu}^{2}-\frac{1}{2}e^{2}\frac{g(\phi_{1},\phi_{2})^{2}}{f(\phi_{1},\phi_{2})}Z_{\mu}^{2}\\ &+\frac{1}{12}\frac{1}{f(\phi_{1},\phi_{2})}H_{\mu\alpha\beta}H^{\mu\alpha\beta}+\frac{1}{2}\epsilon_{\mu\nu\alpha\beta}B^{\alpha\beta}\partial^{\mu}\left(\frac{g(\phi_{1},\phi_{2})}{f(\phi_{1},\phi_{2})}eZ^{\nu}\right)\,.\end{split}

(67)

where the field strength of $B^{\mu\nu}$ is defined as,

H_{\mu\alpha\beta}=\partial_{\mu}B_{\alpha\beta}+\partial_{\alpha}B_{\beta\mu}+\partial_{\beta}B_{\mu\alpha}\,.

(68)

The last term in each line of eq.˜67 is included to ensure that the field equations for $Z^{\mu}$ and $a$ remain consistent with the original Lagrangian. The field equation for $H^{\nu\alpha\beta}$ is given by

\partial^{\nu}\frac{H_{\nu\alpha\beta}}{f(\phi_{1},\phi_{2})}+\epsilon_{\mu\nu\alpha\beta}\,\partial^{\mu}\left(\frac{g(\phi_{1},\phi_{2})}{f(\phi_{1},\phi_{2})}eZ^{\nu}\right)=\epsilon_{\mu\nu\alpha\beta}\partial^{\mu}\partial^{\nu}a

(69)

The right-hand side of this equation is zero in the vacuum or the background of $(1,1)$ string but is connected to the winding number in the global strings, leading to nontrivial coupling between the global strings and axions. Although the gauge $(1,1)$ string lacks coupling between the winding number flux and $B^{\mu\nu}$ , the coupling to $Z^{\mu}$ term as shown in the last term in eq.˜67 leads to the radiation of axion from the string. The action for the interaction takes the form

S_{j}=\int d^{4}x\,\frac{1}{2}B^{\mu\nu}j_{\mu\nu}=\int d^{4}x\,\frac{1}{2}\epsilon_{\mu\nu\alpha\beta}B^{\alpha\beta}\partial^{\mu}\left(\frac{g(\phi_{1},\phi_{2})}{f(\phi_{1},\phi_{2})}eZ^{\nu}\right)\,.

(70)

Analogous to the derivation of the gravitational wave radiation power [115], the radiation power per solid angle at a frequency $\omega$ in a direction ${\bf k}$ is proportional to the amplitude square of the Fourier transformation of $j^{\mu\nu}$ ,

\frac{{\rm d}P_{a}}{{\rm d}\Omega}=\frac{\omega^{2}}{32\pi^{2}}j^{\mu\nu*}(\omega,{\bf k})j_{\mu\nu}(\omega,{\bf k})

(71)

where $|{\bf k}|=\omega$ for massless axions, and

j_{\mu\nu}(\omega,{\bf k})=\int d^{4}x\,e^{-ik\cdot x}j_{\mu\nu}(x)\,.

(72)

We then estimate the radiation power as

\frac{{\rm d}P_{a}}{{\rm d}\Omega}\sim e^{2}f_{a}^{2}\,.

(73)

This estimation considers that the classical field $Z^{\mu}$ varies on the order of $m_{Z}$ within a core of size on the order of $1/m_{Z}$ . The radiation power of global strings is approximately $\sim f_{a}^{2}$ , and the radiation power of gauge strings to graviton is $\sim Gf_{a}^{4}$ , with $G$ the Newtonian gravitational constant. In contrast, the radiation to axions from $(1,1)$ strings is smaller than that from global strings due to the gauge coupling suppression, $e^{2}$ . Nevertheless, the radiation of axions is dominant over that of gravitons, which is suppressed by the Planck mass unless an exceedingly smaller gauge coupling is considered. For the case when the gauge coupling is highly suppressed, we estimate that the gravitational wave spectrum plateau produced by $(1,1)$ string loop oscillations is $\Omega_{\rm GW}(f)\sim O(10^{-4})\sqrt{{G\mu}/{\Gamma}}$ , with the constant efficiency of gravitational wave emission $\Gamma\sim 50$ .

5 Conclusion

In conclusion, we investigate the string solutions and cosmological implications within the ${\rm U}(1)_{\rm Z}\times{\rm U(1)_{PQ}}$ model. This model has two hierarchical symmetry-breaking scales given by two complex scalar fields. By assuming that both $\rm U(1)$ symmetries spontaneously broke after the inflation, we find that the ${\rm U}(1)_{\rm Z}\times{\rm U(1)_{PQ}}$ model can give a rich feature of string networks in cosmology, produce more axions as cold dark matter when embedding the model into the QCD axion framework, and, further, have a new decay product from the gauge string.

We identify three distinct types of string solutions: $(1,0)$ global string with a heavy core, $(0,1)$ global string with the tension close to QCD axion strings, and $(1,1)$ gauge string as a bound state of the former two global strings. Numerical studies confirm the existence of these string solutions, providing a solid foundation for our cosmological considerations.

In the early universe, the formation of string networks during the second phase transition predominantly yield $(1,0)$ strings instead of $(1,1)$ gauge strings, assuming an adiabatic process to generate these strings. Due to the attractive interaction between $(1,0)$ and $(0,1)$ strings, it allows the global strings to form a bound state of $(1,1)$ gauge strings during the evolution of the string network, such that Y-junctions of three types of strings are expected to be present in this string network.

Furthermore, we introduce a KSVZ-like model and make the Nambu-Goldstone mode become the QCD axion. We find that the decay of $(1,0)$ strings with heavy cores dominates axion radiation, providing a potential explanation for the observed dark matter relic abundance for the masses exceeding $10^{5}\,{\rm eV}$ . Additionally, the $(1,1)$ gauge string in this model has a coupling with axions due to the spatial dependent profile of the gauge field $Z^{\mu}$ within the string core. This presents a novel decay channel of gauge strings into axions.

While our work initiates the exploration of this hybrid cosmic strings network and its QCD axion implication, it raises several interesting questions requiring string simulations. It includes the evolution of Y-junctions, the existence of scaling solutions, and the energy spectrum of QCD axion radiated by $(1,0)$ strings. Also, the preference of $(1,1)$ gauge strings to radiate axions or gravitons awaits confirmation through numerical studies. Simulations of the $\rm U(1)_{Z}\times U(1)_{PQ}$ model with a large gauge charge ratio [103], as well as studies that replicate the axion cosmic strings with heavy core global strings [89] have been performed, but many of the questions raised above remain unanswered.

The rich phenomenology of this model offers testable predictions. Upon improving the detection sensitivity of future haloscope experiments, the parameter space in $g_{a\gamma\gamma}-m_{a}$ predicted by this model in the QCD axion framework can be probed. Additionally, if the gauge coupling is small enough, the gravitational radiation is predominantly produced by the $(1,1)$ gauge string. The existence of Y-junctions in the string network can modify the loop distribution function compared to conventional cosmic string scenarios [116]. Therefore, the gauge strings could yield a distinctive gravitational wave signal from this model.

Acknowledgments

We thank Pierre Sikivie, Robert Brandenberger, Jeff Dror, Keisuke Harigaya, Sungwoo Hong, Rachel Houtz, Junwu Huang, Subir Sarkar, Sergey Sibiryakov, and Tanmay Vachaspati for valuable discussions. This work was supported in part by the U.S. Department of Energy under grant DE-SC0022148 at the University of Florida.

Appendix A Asymptotic value of the gauge field configuration

The asymptotic value of $Z_{\theta}(r)=\frac{c}{r}$ at large $r$ can be obtained by requiring $\frac{\partial\mu}{\partial c}=0$ . At large $r$ , we write down the string tension with only $c$ -dependent terms, i.e., the kinetic terms,

$\displaystyle\mu$	$\displaystyle=$	$\displaystyle\int d^{2}x\frac{1}{2}\left[f_{1}^{2}v_{1}^{2}\left(\frac{(j-q_{1}ec)^{2}}{r^{2}}+\frac{f_{1}^{\prime 2}}{f_{1}^{2}}\right)+f_{2}^{2}v_{2}^{2}\left(\frac{(k-q_{2}ec)^{2}}{r^{2}}+\frac{f_{2}^{\prime 2}}{f_{2}^{2}}\right)\right]$	(74)
	$\displaystyle=$	$\displaystyle\pi\int drr\left[v_{1}^{2}\left(\frac{(j-q_{1}ec)^{2}}{r^{2}}\right)+v_{2}^{2}\left(\frac{(k-q_{2}ec)^{2}}{r^{2}}\right)\right]$
	$\displaystyle=$	$\displaystyle\pi\left[v_{1}^{2}(j-q_{1}ec)^{2}+v_{2}^{2}(k-q_{2}ec)^{2}\right]\int\frac{dr}{r}\,.$

In the second equality, we use $f_{1,2}\to 1$ and $f_{1,2}^{\prime}\to 0$ . By minimizing the string tension with respect to $c$ , we can find the value of $c$ :

	$\displaystyle\frac{\partial\mu}{\partial c}$	$\displaystyle=$	$\displaystyle 2\pi\left[(v_{1}^{2}q_{1}^{2}+v_{2}^{2}q_{2}^{2})ec-(jv_{1}^{2}q_{1}+kv_{2}^{2}q_{2})\right]\int\frac{dr}{r}=0$		(75)
	$\displaystyle\Rightarrow$	$\displaystyle c$	$\displaystyle=\frac{1}{e}\frac{jv_{1}^{2}q_{1}+kv_{2}^{2}q_{2}}{v_{1}^{2}q_{1}^{2}+v_{2}^{2}q_{2}^{2}}.$		(76)

Appendix B Numerical solutions to the string profile functions

We take large- $\rho$ behaviors $f_{1}\rightarrow 1$ , $f_{2}\rightarrow 1$ , $g\rightarrow 1$ are satisfied at $\rho=8$ in our numerical study, and fig.˜2) confirms that the profile functions capture the asymptotic values in all configurations. The small $\rho$ behaviors can be found by taking $\rho\rightarrow 0$ in the equation of motions eq.˜36. For $(1,0)$ string, $f_{1}\rightarrow 0$ and $g\rightarrow 0$ must be satisfied at the origin because of no singularity. However, $f_{2}$ doesn’t need to be so. Explicitly, the boundary conditions for a $(1,0)$ string takes the form

\begin{split}\lim_{\rho\to\epsilon}f_{1}(\rho)&=c_{1}+\beta_{1}\rho+O(\rho^{2})\\ \lim_{\rho\to\epsilon}f_{2}(\rho)&=c_{2}+\beta_{2}\rho^{2}+O(\rho^{3})\\ \lim_{\rho\to\epsilon}g(\rho)&=m\rho^{2}+O(\rho^{3})\,,\end{split}

(77)

so that we take six boundary conditions from the values at $\rho=\epsilon$ and their derivatives. Here $c_{1}=0$ , $\beta_{2}=\frac{\lambda_{2}}{16}(c_{2}^{2}-1)c_{2}$ , and we take $\epsilon=10^{-3}$ . The left three parameters $\beta_{1}$ , $c_{2}$ , and $m$ are determined by the shooting method at $\rho=8$ . Similarly, the $(0,1)$ string shares the same large $\rho$ asymptotic behaviors as the $(1,0)$ string. The boundary conditions for $(0,1)$ string at $\rho\rightarrow 0$ are

\begin{split}\lim_{\rho\to\epsilon}f_{1}(\rho)&=c_{1}+\beta_{1}\rho^{2}+O(\rho^{3})\\ \lim_{\rho\to\epsilon}f_{2}(\rho)&=c_{2}+\beta_{2}\rho+O(\rho^{2})\\ \lim_{\rho\to\epsilon}g(\rho)&=m\rho^{2}+O(\rho^{3})\end{split}

(78)

where $c_{2}=0$ and $\beta_{1}=\frac{\lambda_{1}}{16}(c_{1}^{2}-1)c_{1}$ . The profile functions of the $(1,1)$ string must satisfy the smoothness at the origin:

\begin{split}\lim_{\rho\to\epsilon}f_{1}(\rho)&=c_{1}+\beta_{1}\rho+O(\rho^{2})\\ \lim_{\rho\to\epsilon}f_{2}(\rho)&=c_{2}+\beta_{2}\rho+O(\rho^{2})\\ \lim_{\rho\to\epsilon}g(\rho)&=m\rho^{2}+O(\rho^{3})\,,\end{split}

(79)

where $c_{1}=c_{2}=0$ due to no singularity at the origin. Under the parameter space mentioned in section˜3, we find the numerical result in table˜1.

Table 1: Numerical result of string profile functions.

configurations	$\beta_{1}$	$\beta_{2}$	$c_{1}$	$c_{2}$	$m$
$(1,0)$ string	1.1025	-0.0012	0	0.9741	0.4248
$(0,1)$ string	$-1.24\times 10^{-6}$	0.6302	0.9999	0	0.1992
$(1,1)$ string	$1.12143$	0.847589	0	0	0.428086

References

[1] T.W.B. Kibble, Topology of Cosmic Domains and Strings, J. Phys. A 9 (1976) 1387.
[2] A. Vilenkin, Cosmic Strings and Domain Walls, Phys. Rept. 121 (1985) 263.
[3] Y.B. Zeldovich, Cosmological fluctuations produced near a singularity, Mon. Not. Roy. Astron. Soc. 192 (1980) 663.
[4] A. Vilenkin, Cosmological Density Fluctuations Produced by Vacuum Strings, Phys. Rev. Lett. 46 (1981) 1169.
[5] T.W.B. Kibble, Some Implications of a Cosmological Phase Transition, Phys. Rept. 67 (1980) 183.
[6] U.-L. Pen, U. Seljak and N. Turok, Power spectra in global defect theories of cosmic structure formation, Phys. Rev. Lett. 79 (1997) 1611 [astro-ph/9704165].
[7] J. Magueijo, A. Albrecht, D. Coulson and P. Ferreira, Doppler peaks from active perturbations, Phys. Rev. Lett. 76 (1996) 2617 [astro-ph/9511042].
[8] Planck collaboration, Planck 2013 results. XXV. Searches for cosmic strings and other topological defects, Astron. Astrophys. 571 (2014) A25 [1303.5085].
[9] J. Urrestilla, N. Bevis, M. Hindmarsh and M. Kunz, Cosmic string parameter constraints and model analysis using small scale Cosmic Microwave Background data, JCAP 12 (2011) 021 [1108.2730].
[10] J. Lizarraga, J. Urrestilla, D. Daverio, M. Hindmarsh, M. Kunz and A.R. Liddle, Constraining topological defects with temperature and polarization anisotropies, Phys. Rev. D 90 (2014) 103504 [1408.4126].
[11] A. Lazanu, E.P.S. Shellard and M. Landriau, CMB power spectrum of Nambu-Goto cosmic strings, Phys. Rev. D 91 (2015) 083519 [1410.4860].
[12] LIGO Scientific, Virgo, KAGRA collaboration, Constraints on Cosmic Strings Using Data from the Third Advanced LIGO–Virgo Observing Run, Phys. Rev. Lett. 126 (2021) 241102 [2101.12248].
[13] J. Ellis and M. Lewicki, Cosmic String Interpretation of NANOGrav Pulsar Timing Data, Phys. Rev. Lett. 126 (2021) 041304 [2009.06555].
[14] J.J. Blanco-Pillado, K.D. Olum and J.M. Wachter, Comparison of cosmic string and superstring models to NANOGrav 12.5-year results, Phys. Rev. D 103 (2021) 103512 [2102.08194].
[15] P. Auclair et al., Probing the gravitational wave background from cosmic strings with LISA, JCAP 04 (2020) 034 [1909.00819].
[16] G. Boileau, A.C. Jenkins, M. Sakellariadou, R. Meyer and N. Christensen, Ability of LISA to detect a gravitational-wave background of cosmological origin: The cosmic string case, Phys. Rev. D 105 (2022) 023510 [2109.06552].
[17] M. Punturo et al., The Einstein Telescope: A third-generation gravitational wave observatory, Class. Quant. Grav. 27 (2010) 194002.
[18] K. Yagi and N. Seto, Detector configuration of DECIGO/BBO and identification of cosmological neutron-star binaries, Phys. Rev. D 83 (2011) 044011 [1101.3940].
[19] AEDGE collaboration, AEDGE: Atomic Experiment for Dark Matter and Gravity Exploration in Space, EPJ Quant. Technol. 7 (2020) 6 [1908.00802].
[20] S. Hild et al., Sensitivity Studies for Third-Generation Gravitational Wave Observatories, Class. Quant. Grav. 28 (2011) 094013 [1012.0908].
[21] A. Sesana et al., Unveiling the gravitational universe at $\mu$ -Hz frequencies, Exper. Astron. 51 (2021) 1333 [1908.11391].
[22] The Theia Collaboration, C. Boehm, A. Krone-Martins, A. Amorim, G. Anglada-Escude, A. Brandeker et al., Theia: Faint objects in motion or the new astrometry frontier, arXiv e-prints (2017) arXiv:1707.01348 [1707.01348].
[23] P. Agrawal, A. Hook and J. Huang, A CMB Millikan experiment with cosmic axiverse strings, JHEP 07 (2020) 138 [1912.02823].
[24] P. Agrawal, A. Hook, J. Huang and G. Marques-Tavares, Axion string signatures: a cosmological plasma collider, JHEP 01 (2022) 103 [2010.15848].
[25] M. Jain, A.J. Long and M.A. Amin, CMB birefringence from ultralight-axion string networks, JCAP 05 (2021) 055 [2103.10962].
[26] W.W. Yin, L. Dai and S. Ferraro, Probing cosmic strings by reconstructing polarization rotation of the cosmic microwave background, JCAP 06 (2022) 033 [2111.12741].
[27] M. Jain, R. Hagimoto, A.J. Long and M.A. Amin, Searching for axion-like particles through CMB birefringence from string-wall networks, JCAP 10 (2022) 090 [2208.08391].
[28] W.W. Yin, L. Dai and S. Ferraro, Testing charge quantization with axion string-induced cosmic birefringence, JCAP 07 (2023) 052 [2305.02318].
[29] A. Hook, G. Marques-Tavares and C. Ristow, CMB Spectral Distortions from an Axion-Dark Photon-Photon Interaction, 2306.13135.
[30] R. Hagimoto and A.J. Long, Measures of non-Gaussianity in axion-string-induced CMB birefringence, JCAP 09 (2023) 024 [2306.07351].
[31] S. Weinberg, A New Light Boson?, Phys. Rev. Lett. 40 (1978) 223.
[32] F. Wilczek, Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons, Phys. Rev. Lett. 40 (1978) 279.
[33] M.A. Shifman, A.I. Vainshtein and V.I. Zakharov, Can Confinement Ensure Natural CP Invariance of Strong Interactions?, Nucl. Phys. B166 (1980) 493.
[34] J.E. Kim, Weak Interaction Singlet and Strong CP Invariance, Phys. Rev. Lett. 43 (1979) 103.
[35] A.R. Zhitnitsky, On Possible Suppression of the Axion Hadron Interactions. (In Russian), Sov. J. Nucl. Phys. 31 (1980) 260.
[36] M. Dine, W. Fischler and M. Srednicki, A Simple Solution to the Strong CP Problem with a Harmless Axion, Phys. Lett. 104B (1981) 199.
[37] J. Preskill, M.B. Wise and F. Wilczek, Cosmology of the Invisible Axion, Phys. Lett. 120B (1983) 127.
[38] L.F. Abbott and P. Sikivie, A Cosmological Bound on the Invisible Axion, Phys. Lett. 120B (1983) 133.
[39] M. Dine and W. Fischler, The Not So Harmless Axion, Phys. Lett. 120B (1983) 137.
[40] R.D. Peccei and H.R. Quinn, CP Conservation in the Presence of Instantons, Phys. Rev. Lett. 38 (1977) 1440.
[41] R.D. Peccei and H.R. Quinn, Constraints Imposed by CP Conservation in the Presence of Instantons, Phys. Rev. D16 (1977) 1791.
[42] R.L. Davis and E.P.S. Shellard, DO AXIONS NEED INFLATION?, Nucl. Phys. B 324 (1989) 167.
[43] A. Dabholkar and J.M. Quashnock, Pinning Down the Axion, Nucl. Phys. B 333 (1990) 815.
[44] C. Hagmann and P. Sikivie, Computer simulations of the motion and decay of global strings, Nucl. Phys. B 363 (1991) 247.
[45] R.A. Battye and E.P.S. Shellard, Global string radiation, Nucl. Phys. B 423 (1994) 260 [astro-ph/9311017].
[46] R.A. Battye and E.P.S. Shellard, Radiative back reaction on global strings, Phys. Rev. D 53 (1996) 1811 [hep-ph/9508301].
[47] S. Chang, C. Hagmann and P. Sikivie, Studies of the motion and decay of axion walls bounded by strings, Phys. Rev. D 59 (1999) 023505 [hep-ph/9807374].
[48] M. Yamaguchi, J. Yokoyama and M. Kawasaki, Numerical analysis of formation and evolution of global strings in ( 2+1)-dimensions, Prog. Theor. Phys. 100 (1998) 535 [hep-ph/9808326].
[49] M. Yamaguchi, M. Kawasaki and J. Yokoyama, Evolution of axionic strings and spectrum of axions radiated from them, Phys. Rev. Lett. 82 (1999) 4578 [hep-ph/9811311].
[50] M. Yamaguchi, Scaling property of the global string in the radiation dominated universe, Phys. Rev. D 60 (1999) 103511 [hep-ph/9907506].
[51] M. Yamaguchi, J. Yokoyama and M. Kawasaki, Evolution of a global string network in a matter dominated universe, Phys. Rev. D 61 (2000) 061301 [hep-ph/9910352].
[52] C. Hagmann, S. Chang and P. Sikivie, Axion radiation from strings, Phys. Rev. D 63 (2001) 125018 [hep-ph/0012361].
[53] C.J.A.P. Martins, J.N. Moore and E.P.S. Shellard, A Unified model for vortex string network evolution, Phys. Rev. Lett. 92 (2004) 251601 [hep-ph/0310255].
[54] O. Wantz and E.P.S. Shellard, Axion Cosmology Revisited, Phys. Rev. D 82 (2010) 123508 [0910.1066].
[55] T. Hiramatsu, M. Kawasaki and K. Saikawa, Evolution of String-Wall Networks and Axionic Domain Wall Problem, JCAP 08 (2011) 030 [1012.4558].
[56] T. Hiramatsu, M. Kawasaki, T. Sekiguchi, M. Yamaguchi and J. Yokoyama, Improved estimation of radiated axions from cosmological axionic strings, Phys. Rev. D 83 (2011) 123531 [1012.5502].
[57] T. Hiramatsu, M. Kawasaki, K. Saikawa and T. Sekiguchi, Production of dark matter axions from collapse of string-wall systems, Phys. Rev. D 85 (2012) 105020 [1202.5851].
[58] T. Hiramatsu, M. Kawasaki, K. Saikawa and T. Sekiguchi, Axion cosmology with long-lived domain walls, JCAP 01 (2013) 001 [1207.3166].
[59] M. Kawasaki, K. Saikawa and T. Sekiguchi, Axion dark matter from topological defects, Phys. Rev. D 91 (2015) 065014 [1412.0789].
[60] V.B.. Klaer and G.D. Moore, The dark-matter axion mass, JCAP 11 (2017) 049 [1708.07521].
[61] M. Gorghetto, E. Hardy and G. Villadoro, Axions from Strings: the Attractive Solution, JHEP 07 (2018) 151 [1806.04677].
[62] M. Buschmann, J.W. Foster and B.R. Safdi, Early-Universe Simulations of the Cosmological Axion, Phys. Rev. Lett. 124 (2020) 161103 [1906.00967].
[63] M. Gorghetto, E. Hardy and G. Villadoro, More axions from strings, SciPost Phys. 10 (2021) 050 [2007.04990].
[64] M. Buschmann, J.W. Foster, A. Hook, A. Peterson, D.E. Willcox, W. Zhang et al., Dark matter from axion strings with adaptive mesh refinement, Nature Commun. 13 (2022) 1049 [2108.05368].
[65] I. Stern, ADMX Status, PoS ICHEP2016 (2016) 198 [1612.08296].
[66] D. Alesini et al., Search for invisible axion dark matter of mass m ${}_{a}=43~\mu$ eV with the QUAX– $a\gamma$ experiment, Phys. Rev. D 103 (2021) 102004 [2012.09498].
[67] J. De Miguel and J.F. Hernández-Cabrera, Discovery prospects with the Dark-photons & Axion-Like particles Interferometer part I, 2303.03997.
[68] M. Lawson, A.J. Millar, M. Pancaldi, E. Vitagliano and F. Wilczek, Tunable axion plasma haloscopes, Phys. Rev. Lett. 123 (2019) 141802 [1904.11872].
[69] S. Beurthey et al., MADMAX Status Report, 2003.10894.
[70] B.T. McAllister, G. Flower, J. Kruger, E.N. Ivanov, M. Goryachev, J. Bourhill et al., The ORGAN Experiment: An axion haloscope above 15 GHz, Phys. Dark Univ. 18 (2017) 67 [1706.00209].
[71] B. Aja et al., The Canfranc Axion Detection Experiment (CADEx): search for axions at 90 GHz with Kinetic Inductance Detectors, JCAP 11 (2022) 044 [2206.02980].
[72] BRASS: Broadband Radiometric Axion SearcheS collaboration. https://www.physik.uni-hamburg.de/iexp/gruppe-horns/forschung/brass.html.
[73] BREAD collaboration, Broadband Solenoidal Haloscope for Terahertz Axion Detection, Phys. Rev. Lett. 128 (2022) 131801 [2111.12103].
[74] P. Agrawal, G. Marques-Tavares and W. Xue, Opening up the QCD axion window, JHEP 03 (2018) 049 [1708.05008].
[75] R.T. Co, L.J. Hall and K. Harigaya, Qcd axion dark matter with a small decay constant, Phys. Rev. Lett. 120 (2018) 211602.
[76] K. Harigaya and J.M. Leedom, QCD Axion Dark Matter from a Late Time Phase Transition, JHEP 06 (2020) 034 [1910.04163].
[77] M.S. Turner, Cosmic and local mass density of “invisible” axions, Phys. Rev. D 33 (1986) 889.
[78] D.H. Lyth, Axions and inflation: Vacuum fluctuations, Phys. Rev. D 45 (1992) 3394.
[79] L. Visinelli and P. Gondolo, Dark matter axions revisited, Phys. Rev. D 80 (2009) 035024.
[80] P. Sikivie, Axions, domain walls, and the early universe, Phys. Rev. Lett. 48 (1982) 1156.
[81] S. Chang, C. Hagmann and P. Sikivie, Studies of the motion and decay of axion walls bounded by strings, Phys. Rev. D 59 (1998) 023505.
[82] T. Hiramatsu, M. Kawasaki and K. Saikawa, Evolution of string-wall networks and axionic domain wall problem, Journal of Cosmology and Astroparticle Physics 2011 (2011) 030.
[83] M. Kawasaki, K. Saikawa and T. Sekiguchi, Axion dark matter from topological defects, Phys. Rev. D 91 (2015) 065014.
[84] A. Ringwald and K. Saikawa, Publisher’s note: Axion dark matter in the post-inflationary peccei-quinn symmetry breaking scenario [phys. rev. d 93, 085031 (2016)], Phys. Rev. D 94 (2016) 049908.
[85] K. Harigaya and M. Kawasaki, Qcd axion dark matter from long-lived domain walls during matter domination, Physics Letters B 782 (2018) 1.
[86] A. Caputo and M. Reig, Cosmic implications of a low-scale solution to the axion domain wall problem, Phys. Rev. D 100 (2019) 063530.
[87] L. Visinelli and P. Gondolo, Axion cold dark matter in nonstandard cosmologies, Phys. Rev. D 81 (2010) 063508.
[88] R.T. Co, L.J. Hall and K. Harigaya, Axion Kinetic Misalignment Mechanism, Phys. Rev. Lett. 124 (2020) 251802 [1910.14152].
[89] V.B. Klaer and G.D. Moore, How to simulate global cosmic strings with large string tension, JCAP 10 (2017) 043 [1707.05566].
[90] L. Di Luzio, M. Giannotti, E. Nardi and L. Visinelli, The landscape of QCD axion models, Phys. Rept. 870 (2020) 1 [2003.01100].
[91] S.M. Barr and D. Seckel, Planck scale corrections to axion models, Phys. Rev. D 46 (1992) 539.
[92] M. Ibe, M. Suzuki and T.T. Yanagida, $B-L$ as a Gauged Peccei-Quinn Symmetry, JHEP 08 (2018) 049 [1805.10029].
[93] T. Hiramatsu, M. Ibe and M. Suzuki, New Type of String Solutions with Long Range Forces, JHEP 02 (2020) 058 [1910.14321].
[94] T.W.B. Kibble and A. Vilenkin, Density of strings formed at a second order cosmological phase transition, hep-ph/9501207.
[95] J. Urrestilla and A. Vilenkin, Evolution of cosmic superstring networks: A Numerical simulation, JHEP 02 (2008) 037 [0712.1146].
[96] N. Bevis and P.M. Saffin, Cosmic string Y-junctions: A Comparison between field theoretic and Nambu-Goto dynamics, Phys. Rev. D 78 (2008) 023503 [0804.0200].
[97] J. Lizarraga and J. Urrestilla, Survival of pq-superstrings in field theory simulations, JCAP 04 (2016) 053 [1602.08014].
[98] J.R.C.C.C. Correia and C.J.A.P. Martins, Multitension strings in high-resolution U(1)×U(1) simulations, Phys. Rev. D 106 (2022) 043521 [2208.01525].
[99] A. Rajantie, M. Sakellariadou and H. Stoica, Numerical experiments with p F- and q D-strings: The Formation of (p,q) bound states, JCAP 11 (2007) 021 [0706.3662].
[100] E.J. Copeland and P.M. Saffin, On the evolution of cosmic-superstring networks, JHEP 11 (2005) 023 [hep-th/0505110].
[101] M. Sakellariadou and H. Stoica, Dynamics of F/D networks: The Role of bound states, JCAP 08 (2008) 038 [0806.3219].
[102] A. Avgoustidis and E.J. Copeland, The effect of kinematic constraints on multi-tension string network evolution, Phys. Rev. D 81 (2010) 063517 [0912.4004].
[103] T. Hiramatsu, M. Ibe and M. Suzuki, Cosmic string in Abelian-Higgs model with enhanced symmetry — Implication to the axion domain-wall problem, JHEP 09 (2020) 054 [2005.10421].
[104] P. Sikivie, Axion Cosmology, Lect. Notes Phys. 741 (2008) 19 [astro-ph/0610440].
[105] A. Saurabh, T. Vachaspati and L. Pogosian, Decay of Cosmic Global String Loops, Phys. Rev. D 101 (2020) 083522 [2001.01030].
[106] D.J. Gross, R.D. Pisarski and L.G. Yaffe, QCD and Instantons at Finite Temperature, Rev. Mod. Phys. 53 (1981) 43.
[107] ADMX collaboration, A SQUID-based microwave cavity search for dark-matter axions, Phys. Rev. Lett. 104 (2010) 041301 [0910.5914].
[108] ADMX collaboration, A Search for Invisible Axion Dark Matter with the Axion Dark Matter Experiment, Phys. Rev. Lett. 120 (2018) 151301 [1804.05750].
[109] ADMX collaboration, Extended Search for the Invisible Axion with the Axion Dark Matter Experiment, Phys. Rev. Lett. 124 (2020) 101303 [1910.08638].
[110] ADMX collaboration, Search for Invisible Axion Dark Matter in the 3.3–4.2 $\mu$ eV Mass Range, Phys. Rev. Lett. 127 (2021) 261803 [2110.06096].
[111] M.J. Dolan, F.J. Hiskens and R.R. Volkas, Advancing globular cluster constraints on the axion-photon coupling, JCAP 10 (2022) 096 [2207.03102].
[112] M. Kalb and P. Ramond, Classical direct interstring action, Phys. Rev. D 9 (1974) 2273.
[113] A. Vilenkin and T. Vachaspati, Radiation of Goldstone Bosons From Cosmic Strings, Phys. Rev. D 35 (1987) 1138.
[114] R.L. Davis and E.P.S. Shellard, Antisymmetric Tensors and Spontaneous Symmetry Breaking, Phys. Lett. B 214 (1988) 219.
[115] S. Weinberg, Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity, John Wiley and Sons, New York (1972).
[116] R. Brandenberger, H. Firouzjahi, J. Karouby and S. Khosravi, Gravitational Radiation by Cosmic Strings in a Junction, JCAP 01 (2009) 008 [0810.4521].