††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<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### 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$.

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}$.111 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.222We 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.$

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)\,.$

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)\,.$

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),$

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.333 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.

Figure 2: Profile functions of three string configurations as a dimensionless
parameter $\rho=ev_{1}r$.

Figure 3: String tension and its components vs $\rho$ for winding numbers
$(1,0)$, $(0,1)$, and $(1,1)$. The long-dashed lines represent the kinetic
energy of complex scalar fields, $\mu_{k}$, while The dashed lines represent
the summation of gauge field and potential energy, $\mu_{Z}+\mu_{U}$. These
three components sum up to be the total tension $\mu$, as shown by solid
lines.

#### 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}),$
$\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}),$

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.

Figure 4: String tension and fitting curves in eq. (37) and eq. (38) as a
dimensionless parameter $\rho$.

##  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.

Figure 5: The difference between $(1,0)$ and $(1,1)$ sting tension in unit of
$\pi v_{2}^{2}$ as a function of $v_{2}L$. We choose $e=4\times 10^{-5}$, and
$\lambda_{1}=\lambda_{2}=1$. $\lambda_{2}=1$ implies $v_{2}L\sim L/\xi_{2}$.

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.

Figure 6: The sketch of a Y-junction.

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.

Figure 7: The three contributions to the dark matter relic abundance as a
function of $v_{2}$. The red solid lines show the string decay contribution
from $(1,0)$ and $(0,1)$ strings. The blue solid line shows the misalignment
mechanism contribution. The orange dot-dashed lines show the domain wall decay
contribution from the walls bound by $(1,0)$ and $(0,1)$ strings. The total
axion abundance is shown by a horizontal black dashed line. We set the gauge
coupling $e=4\times 10^{-5}$, $\lambda_{1}=\lambda_{2}=1$, the Lorentz factor
$\gamma=60$ for domain wall decay.

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}$.

Figure 8: The bound on the QCD axion-photon coupling as a function of the
axion mass (Top). The existing constraints are shown by gray-shaded regions
from ADMX experiments [107, 108, 109, 110] and globular clusters [111]. The
combined projection using future haloscope experiments is shown by a blue
dashed line [65, 66, 67, 68, 69, 70, 71, 72, 73]. The gray solid line denotes
the power-law dependence of the axion-photon coupling on QCD axion mass in the
KSVZ model. Considering the QCD axion to explain 100% dark matter relic
abundance, the magenta and cyan bands show the mass region opened up by the
gauge global string model in Scenarios A and B. The bottom panel shows the
value of $v_{1}$ and $v_{2}$ to reproduce the dark matter relic abundance. The
lines are truncated at $v_{1}=v_{2}$ since we assume $v_{1}\gtrsim 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].