Gauged Global Strings

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.


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 → H, the vacuum manifold is a quotient space, G/H.Cosmic string solutions are connected to the non-trivial first homotopy group π 1 (G/H).The simplest cosmic string originates from a complex scalar field, denoted as Φ, with a global U(1) symmetry.Its vacuum expectation value (VEV), ⟨Φ⟩ = 1 √ 2 f a spontaneously breaks its global U(1) symmetry.A cosmic string along the z-direction exhibits a scalar field configuration in polar coordinates (r, θ), 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, µ ≃ 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 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 ≳ f a .A phase transition occurs when the temperature falls around f a .During this transition, the vacuum expectation value ⟨Φ⟩ turns on, and the symmetry is spontaneously broken.The phase of ⟨Φ⟩ 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.
In this paper, we introduce a simple model having two complex fields, Φ 1 and Φ 2 , with a symmetry of gauge U(1) Z × global 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 Φ 1 and Φ 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 U(1) Z × 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.

U(1) Z × U(1) PQ
In this section, we introduce a model having a gauge symmetry U(1) Z and a global symmetry U(1) PQ .The dynamics of this model are driven by two scalar fields, Φ 1 and Φ 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 [33][34][35][36]90].

The model
We consider the gauge symmetry and the global symmetry, U(1) Z × U(1) PQ , both of which are broken by two complex scalar Φ 1 and Φ 2 .The dynamics are described by their gaugeinvariant and renormalizable Lagrangian density, which takes the form (2.1) Z µν represents the field strength of the U(1) Z gauge boson Z µ , defined as The scalar fields interact with the gauge boson Z µ with the gauge coupling e through the covariant derivative terms, as expressed by (2.2) We assign that Φ 1 carries charges (+1, +1) under U(1) Z × U(1) PQ , while Φ 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 Φ 1 and Φ 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 λ 1 and λ 2 while assuming that interactions between the two scalars, namely , are either negligible or zero.In the vacuum, Φ 1 and Φ 2 acquire non-zero expectation values, denoted as v 1 and v 2 , This spontaneous symmetry breaking leads to the sequential breaking of the U(1) Z ×U(1) PQ symmetry in the cosmic evolution, particularly when v 1 ≫ v 2 .The first phase transition results in a non-zero VEV of ⟨Φ 1 ⟩, leading to the breaking of U(1) Z and giving the gauge boson Z µ mass, m Z = ev 1 .After the second phase transition, the non-zero VEV of ⟨Φ 1 ⟩ and ⟨Φ 2 ⟩ further breaks the global symmetry U(1) PQ , with the gauge boson mass increasing, We conveniently parametrize the perturbations of Φ 1 and Φ 2 using real scalar fields φ1 (x), φ2 (x), π 1 (x) and π 2 (x), where α 1 and α 2 represent the rotation angle of Φ 1 and Φ 2 , respectively.We identify the axion, denoted as a(x), by ensuring that it is orthogonal to the would-be Goldstone boson π z or the longitudinal mode of the gauge boson Z µ .Using the expression of the U(1) Z current J µ z and π z -to-vacuum matrix element, ⟨vac|J µ z (0)|π z (p)⟩ = ip µ v, we identify the would-be Goldstone boson of U(1) Z , which takes the form ) Requiring orthogonality between the axion field and π z yields the expression for a(x) According to the symmetry of Φ 1 and Φ 2 , we can express π 1 and π 2 in terms of the rotation angles of U(1) Z and U(1) PQ in the vacuum manifold, specifically Here, α Z and α PQ represent the rotation angles of U(1) Z and U(1) PQ , respectively.Consequently, we can rewrite the axion fields in terms of the U(1) PQ rotation angle α PQ , (2.8) The parameter v a tells us the magnitude of the vacuum expectation value that spontaneously breaks U(1) PQ .

QCD axion
The U(1) Z × U(1) 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, This fermion resides in the fundamental representation of the Standard Model SU(3) c color symmetry and is a singlet under the U(1) Z symmetry and other Standard Model symmetries.Moreover, the fermion carries a chiral charge under the PQ symmetry, with Q L having a +1 charge and Q R a −1 charge.Since it is a singlet in the gauge U(1) Z symmetry, the fermion does not introduce any additional U(1) Z anomaly, preserving the gauge symmetry anomaly-free.
Taking into account the symmetry of the fermion Q, we construct the Lagrangian where Λ represents the UV cutoff of the theory.The expectation values of Φ 1 and Φ 2 yield the fermion mass, m Q = yv 1 v 2 /(2Λ).By integrating out the heavy scalars and employing eqs.(2.4) and (2.6), we derive the effective Lagrangian for axion couplings to fermions 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, the heavy fermion becomes disentangled from the axion, introducing an axion-gluon coupling term where g s is the coupling constant of 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 va fa .However, in this case, the domain wall number for cosmic string solutions is N = 1.This arises because, in the vacuum manifold of U(1) PQ , the minimal winding is achieved by choosing the angle α PQ from 0 to π. α PQ = 0 and α PQ = π are gauge-equivalent, and 0 can return to π through the gauge group U(1) Z (see fig. 1).A similar method of counting domain wall number is observed in the PQWW model [31,32,40], which has three domain walls in an axion string, but va fa = 6.Alternatively, we can construct a model by introducing two sets of quarks, denoted as Q 1 and Q 2 , as discussed by Barr and Seckel [91].These quarks are color-triplets under SU(3) c and are assigned the following chiral charges under U(1) Z × U(1) PQ : Q 1L has charges (1/2, 1/2), Q 1R has charges (−1/2, −1/2), Q 2L has charges (−1/2, −1/2), and Q 2R has charges (1/2, 1/2).This charge assignment ensures that U(1) Z remains anomaly-free.Furthermore, it leads to Q 1 interacting with Φ 1 and Q 2 interacting with Φ 2 though the Lagrangian density, Following a similar procedure, we arrive at the same axion gluon coupling shown in eq.(2.12).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, U(1) B−L is anomaly-free and could serve as a gauge symmetry.Another intriguing possibility is that U(1) Z symmetry can be identified with the U(1) B−L [92].

String solutions
In this section, we study the string solutions of U(1) Z ×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.
For the winding number (1, 0), the classical configurations of the scalar and gauge fields outside the string cores take the form 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 where we take the core size as δ.Since m Z is typically lighter than the scalar field mass m 1 and m 2 , we can choose the core size δ ≃ m −1 Z . 1 The parameter c is determined by minimizing the tension, (2.16) Therefore, the (1, 0) tension outside the core is Next, for the (0, 1) string tension outside the string core, the configurations take the form 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, 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, α Z → α Z + θ.Consequently, (1, 0) and (0, 1) are string and anti-string. 2he (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 µ configuration with c = 1/e cancels the Φ 1 and Φ 2 gradient terms simultaneously.Hence, we conclude that the tension outside the core of (1, 1) string is zero, µ 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, where is the z-component of the "magnetic" field B = ∇ × Z.For a (j, k) string, we choose the Ansatz for the fields: where f 1 (r), f 2 (r), and g(r) are the profile functions and where we have c to minimize the gradient energy (deviation in appendix A).They need to satisfy the boundary conditions when j ̸ = 0 and k ̸ = 0: (2.21) Note that the boundary conditions of the profile function f 1,2 (r) for r → 0 are not necessarily equal to 0 when the corresponding winding number is 0. Substitute the above Ansatz into eq.( 2.19), we find the string tension in terms of f 1 , f 2 , g: 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 Θ(r) is introduced to simplify the profile functions, taking the form as (2.23) 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 , yielding Substituting the above expressions into eq.( 2.22), we get the string tension with the winding number (j, k) as We use + L rc when performing the radial integral.We introduce the Kronecker delta functions, δ j0 , δ 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 section 2.3 with respect to r c , r 1,c , and r 2,c individually.With the two Higgs masses m i = λ i /2v i (i = 1, 2) and the gauge boson mass m Z = e v 2 1 + v 2 2 , the radii are given as Hence, the string tension can be written in terms of the mass of the massive fields (2.28) We define a string global charge q global ≡ j − k.When q 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 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 global | = 1.The tension difference between gauge strings and the two global strings gives rise to the binding energy, represented by the expression For m Z L ≳ 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 µ (1,0) > µ (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 Φ 1 and Φ 2 can carry generic gauge charges q 1 and q 2 , different from +1 that we assign before.The covariant derivative of Φ 1 and Φ 2 is, therefore, Similarly, we deduce the string tension, taking the form where the string global charge q global ≡ jq 2 − kq 1 .

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.(2.20), the equations of motion of Φ 1 , Φ 2 and Z µ leads to the following differential equations: 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 → 1, f 2 → 1, g → 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, To simplify the analysis, we introduce the dimensionless parameter ρ = ev 1 r.The profile functions for (1, 0), (0, 1), and (1, 1) strings are shown in fig. 2.More details of the numerical solutions are given in table 1 of appendix B. For zero winding numbers, we find that the profile functions of the corresponding scalar fields are non-zero at the origin. 3

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 > 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 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 πv 2 1 ln( m 1 m Z ) for the (1, 0) string and πv 2  2 ln( 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 , where we determine x 10 = 1.85 and x 01 = 3.19 by fitting the string tension outside of the Figure 3: String tension and its components vs ρ for winding numbers (1, 0), (0, 1), and (1, 1).The long-dashed lines represent the kinetic energy of complex scalar fields, µ k , while The dashed lines represent the summation of gauge field and potential energy, µ Z + µ U .These three components sum up to be the total tension µ, as shown by solid lines.

Cosmological Implication
The U(1) Z × 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 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 U(1) Z ×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 U(1) Z × 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.

Formation of string network
Two distinct phase transitions occurred during the evolution of the universe, sequentially breaking the U(1) Z × 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 π 1 string, in contrast to those generated after the second phase transition when the symmetry is completely broken.In the second phase transition, the π 1 strings undergo modifications due to the scalar field Φ 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 Φ 1 (x) acquires a non-zero VEV, ⟨Φ 1 (x)⟩ = v 1 / √ 2, while the second scalar field Φ 2 (x) remains at its potential minimum, Φ 2 (x) = 0.As a second-order phase transition, this occurs around the temperature T c,1 ∼ 6λ 1 λ 1 +3e 2 v 1 .The π 1 string with a winding number 1 is therefore formed by the Kibble mechanism [1].The π 1 strings with higher winding numbers, |j| > 1, are rarely formed.Since Φ 2 (x) = 0 remains at its origin, π 1 strings can be regarded as (1, n) string, with the integer n determined by the second phase transition when Φ 2 acquires a VEV.To minimize the energy, the gauge field Z µ compensates for the gradient of the phase of Φ 1 (x), and thus the π 1 string is a 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 π 1 string network enters a scaling regime where a few long π 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 ∼ beginning of the phase transition can be estimated by the Ginzburg length [94], Since this correlation length is much shorter than the separation of the π 1 long strings, H −1 , within the correlation length, Φ 1 (x) can be considered homogeneous.Consequently, we expect the formation of (0, 1) string independently of the pre-existence of the π 1 strings.
In the case of π 1 string, within a distance ∼ ξ 2 , the presence of non-trivial winding number (j = 1) in the Φ 1 fields influences the Φ 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 Φ 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.(2.29) and eq.(2.31), The (1, 0) strings often is lighter than (1, 1) strings, particularly considering v 2 L ∼ O(1) during the second phase transition (see fig. 5).Therefore, we expect the predominant formation of (1, 0) strings.However, it is worth noting that the phase transition may exhibit more complex dynamics, or the adiabatic argument may fail, potentially leading to the simultaneous production of (1, 0), (0, 1) and gauge string.The evolution of (1, 0) and (0, 1) strings is more complicated than the evolution of conventional global string networks.The complexity arises due to the formation of (1, 1) gauge string segments, known as Y-junctions, when the two types of strings encounter and intercommute (see fig. 6).The formation of Y-junctions can be understood as the result of attractive forces between the two strings or µ (1,0) +µ (0,1) −µ (1,1) > 0, as shown in eq.(2.32).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 U(1) strings [95][96][97][98], global 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 U(1) Z × U(1) PQ model, but cosmological simulations for this model is encouraged to validate this conclusion.
In summary, we reveal that (1, 0) and (0, 1) are produced during the second phase transition by considering the adiabatic condition.The (1, 0) and (0, 1) strings form Yjunctions 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 Yjunctions.
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.

QCD axion abundance
In this section, we explore the interesting possibility of the Goldstone modes in the U(1) Z × U(1) PQ being the well-motivated QCD axion.We consider that the symmetry breaking of U(1) Z × 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.( 4.2) by setting L ∼ H −1 (T = Λ 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,63,64,105], 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 ∼ t, leading to the string tension of a long string taking the form 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.(3.2), (3.3).Notably, the first term in eq. ( 4.3) introduces a substantial correction to the string tension, owing to . In the case of (0, 1) strings, their tension closely resembles that of the standard QCD axion strings, as Considering that the string network enters a scaling regime, the number of long strings within one Hubble patch is of the order of O(1).The energy density of the long string within one Hubble volume at time t can thus be expressed as 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 O(1).
To consider the decay of strings into axions, we track the evolution of the density of long strings using the following equation where 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 To achieve the last approximate form, we use the string density evolution function, eq.(4.6), and neglect the time dependence in N i µ i (t).The average energy of axions, 1 ωi (t) , is calculated through the energy spectrum of the emitted axions, dE i dω , 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 ω = 2π/t, i.e., dE dω ∝ δ(ω − 2πt −1 ).In Scenario B, the spectrum takes the form dE dω ∝ 1 ω , with the frequency range ω ∈ [2πt −1 , 2πδ −1 ], where δ represents the UV cutoff.For a (1, 0) string, δ ≃ 1/m 1 ; for a (0, 1) string, δ ≃ 1/m 2 .Although we choose a delta function in the spectrum in Scenario A, the analysis should encompass the case of dE dω ∝ 1 ω q with q > 1, corresponding to an IR-dominant spectrum.Substituting H = 1 2t , we solve eq.(4.7) and find that the number density of axions from strings can be expressed as (4.9) 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 1 ω = t 2π and neglecting the O(1) coefficients, we estimate the axion number density as follows In Scenario B, where 1 ω = t 2π ln t δ −1 , the axion number density is estimated as In the two equations above, the factor of 1 t 3/2 accounts for the dilution due to the expansion of the universe.Following axion production, their number density scales as 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.(4.10) or eq.(4.11) 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 α P Q = α 1 − α 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 The axion mass m a (t) is time-or temperature-dependent.At high temperature (T ≳ 1 GeV), the mass can be approximated as with the dilute instanton gas approximation [106], and at temperatures below 100 MeV Once the axion mass turns on at time t 1 , the strings bound to the wall have an IR cutoff L ∼ m −1 a , comparable to the wall thickness.Hence the energy stored in a (1, 0) or (0, 1) string per unit length can be expressed as 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, The (1, 0) and (0, 1) strings have different transition times due to their distinct string tension µ(t), as shown in eq.(4.15).For walls connected to the light (0, 1) strings, the situation aligns with the standard scenario [47], which approximates t 2 ∼ t 1 .However, when the walls are connected to the heavier (1, 0) strings, where µ (1,0) > µ (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 σ(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 , After time t 3 , the axions produced during the collapse of the wall-string system are boosted.We define an average Lorentz factor γ as the ratio of the energy density to the axion mass at t 3 , denoted as γ ≡ ⟨ω⟩ ma(t 3 ) .The number density of the produced axions is scaled by the volume In the second equality, substituting eq.(4.17), we find that the dependence on time t 3 drops out of the estimate of n DW a (t).According to simulation results in Ref. [47], γ ∼ 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.(4.18) 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.(4.18), allowing us to compare the string decay contribution with the string-wall collapse contribution at t 2 directly.In Scenario A, we compare n str a with n DW a at t 2 , and set the number of (1, 0) long strings is In the last line, we use eq.(4.16) to replace µ(t 2 ).Considering t2 > t 1 and γ ≳ 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 γ = 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, n str a 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 DW a is the number density of axions from domain wall decay bounded by (1, 0) and (0, 1) strings, and ρ vac a = Ω vac a ρ c,0 is the energy density of axions produced by the misalignment mechanism, with ρ c,0 the current critical density of the universe, for f a < 2 × 10 15 GeV, where ⟨θ 2 a ⟩ = π 2 /3 is the average value of the square of the misalignment angle.If one considers axion to explain the 100% of the relic abundance of dark matter observed today, this will require ρ a,0 ρ c,0 h 2 = Ω DM h 2 ∼ 0.12.(4.22) Figure 7 summarizes the three contributions to axion energy density as a function of v 2 for both Scenario A and B. It confirms that the decay from the (1, 0) string with the heavy core is the dominant contribution to dark matter relic abundance for most regions of v 2 .
In the KSVZ model, the axion-photon coupling is linked to the axion mass through the relation g aγγ = 2.0 × 10 −16 C aγ (m a /µeV) GeV −1 , with C aγ = −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 µ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 µ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 ∈ (40, 180) µeV.In contrast, our proposed     [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 powerlaw 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 ≳ v 2 .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.

(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.(2.1), and investigate the axion interactions with classical configurations of Z µ (x), Φ 1 (x) and Φ 2 (x).Subsequently, we perform a gauge transformation, Z µ → Z µ + 1 e ∂ µ α Z .This transformation eliminates the quadratic term of Z µ ∂ µ a and Z µ ∂ µ π z .The Lagrangian density thus takes the form In this equation, we introduce the functions g(ϕ 1 , ϕ 2 ) and f (ϕ 1 , ϕ 2 ) to simplify the notation,

.25)
When considering the vacuum expectation values of ϕ 1 and ϕ 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 (4.23) already reveals that (1, 1) strings interact with axions through the term g(ϕ 1 , ϕ 2 ) e Z µ ∂ µ a, owing to the non-trivial field configurations of Z µ (x), ϕ 1 (x), ϕ 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 ∼ 1/m Z .
To calculate the radiation power of (1, 1) strings to axions, we can employ the Kalb-Ramond field B µν [112][113][114].A duality relationship exists between real massless field, axion, and the two-form antisymmetry gauge field, B µν , given by This relation is satisfied beyond the string core.However, in the presence of ϕ 1 (x), ϕ 2 (x) and Z µ (x), this duality relationship is modified.Based on the field equation of a(x) derived from eq. (4.23), the relation between a and B µν should be replaced with Substituting axion with B αβ or its gauge-invariant field strength H µνλ in the Lagrangian density, we obtain the following expression, where the field strength of B µν is defined as, The last term in each line of eq.(4.29) is included to ensure that the field equations for Z µ and a remain consistent with the original Lagrangian.The field equation for H ναβ is given by 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 µν , the coupling to Z µ term as shown in the last term in eq.(4.29) leads to the radiation of axion from the string.The action for the interaction takes the form Analogous to the derivation of the gravitational wave radiation power [115], the radiation power per solid angle at a frequency ω in a direction k is proportional to the amplitude square of the Fourier transformation of j µν , where |k| = ω for massless axions, and

.34)
We then estimate the radiation power as This estimation considers that the classical field Z µ 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 ∼ f 2 a , and the radiation power of gauge strings to graviton is ∼ Gf 4 a , 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 Ω GW (f ) ∼ O(10 −4 ) Gµ/Γ, with the constant efficiency of gravitational wave emission Γ ∼ 50.

Conclusion
In conclusion, we investigate the string solutions and cosmological implications within the U(1) Z × U(1) PQ model.This model has two hierarchical symmetry-breaking scales given by two complex scalar fields.By assuming that both U(1) symmetries spontaneously broke after the inflation, we find that the U(1) Z × 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 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 µ 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 U(1) Z × 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γγ − 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.configurations.The small ρ behaviors can be found by taking ρ → 0 in the equation of motions section 3.For (1, 0) string, f 1 → 0 and g → 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 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.

Figure 1 :
Figure 1: The cross-section of the vacuum manifold of U(1) Z × U(1) PQ , where the two black dots represent the angle α PQ = 0 and α PQ = π.

z 2 ,Figure 2 :
Figure 2: Profile functions of three string configurations as a dimensionless parameter ρ = ev 1 r.

Figure 6 :
Figure 6: The sketch of a Y-junction.

Figure 7 :
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 × 10 −5 , λ 1 = λ 2 = 1, the Lorentz factor γ = 60 for domain wall decay.

Figure 8 :
Figure8: 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 powerlaw 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 ≳ v 2 .

Table 1 :
Numerical result of string profile functions.