Cosmic string in Abelian-Higgs model with enhanced symmetry — Implication to the axion domain-wall problem

In our previous work, we found new types of the cosmic string solutions in the Abelian-Higgs model with an enhanced U(1) global symmetry. We dubbed those solutions as the compensated/uncompensated strings. The compensated string is similar to the conventional cosmic string in the Abrikosov-Nielsen-Olesen (ANO) string, around which only the would-be Nambu-Goldstone (NG) boson winds. Around the uncompensated string, on the other hand, the physical NG boson also winds, where the physical NG boson is associated with the spontaneous breaking of the enhanced symmetry. Our previous simulation in the 2+1 dimensional spacetime confirmed that both the compensated/uncompensated strings are formed at the phase transition of the symmetry breaking. Non-trivial winding of the physical NG boson around the strings potentially causes the so-called axion domain- wall problem when the model is applied to the axion model. In this paper, we perform simulation in the 3+1 dimensional spacetime to discuss the fate of the uncompensated strings. We observe that the evolution of the string-network is highly complicated in the 3+1 dimensional simulation compared with that seen in the previous simulation. Despite such complications, we find that the number of the uncompensated strings which could cause can be highly suppressed at late times. Our observation suggests that the present setup can be applied to the axion model without suffering from the axion domain-wall problem.


Introduction
Recently, we discussed new types of string (vortex) solutions in the Abelian-Higgs model with two complex scalar fields [1]. As a peculiar feature, the model has an accidental U(1) global symmetry, U(1) global , enhanced by the hierarchical charge assignment of the U(1) gauge symmetry, U(1) local . By the spontaneous breaking of U(1) local and U(1) global , one of the Nambu-Goldstone (NG) boson is absorbed by the Higgs mechanism, while the other one appears as a physical NG boson. This class of setup has been discussed in the context of the QCD axion models where a global Peccei-Quinn symmetry [2][3][4][5] is identified with U(1) global [6][7][8][9][10]. In those applications, the gauge symmetry protects the PQ symmetry from the explicit breaking caused by the quantum gravity effects [11][12][13][14][15][16].
In [1], we call the new solutions, the "compensated" and the "uncompensated" strings. The compensated string is more like the local string [17], in which only the would-be NG boson winds non-trivially. The uncompensated string has, on the other hand, properties in between the local string and the global string solutions [17]. Around the uncompensated strings, not only the would-be NG boson but also the physical NG boson wind. 1 We also found that the uncompensated strings have long-range repulsive and attractive forces. 2 1 Similar to the compensated/uncompensated string, the composite string solutions in the de Sitter space have been discussed in [18]. 2 This type of the string solutions has been discussed in [19], where the evolution of the string network is considered in the presence of the explicit breaking of the global U(1) symmetry. In this paper, we discuss how the network evolves in the absence of the explicit breaking, and hence, no domain walls are attached until we apply the model to the axion model. We thank Guy Moore for letting us know this previous work.

JHEP09(2020)054
In our previous work [1], we also performed classical lattice simulation of the timeevolution in the spacetime of the 2+1 dimension (2+1 D). (The earlier numerical simulation of the vortex formation in a similar setup has been performed in [20], which exhibits consistent features with our simulations in [1].) The simulation confirmed the ubiquitous formation of the uncompensated strings at the phase transition. As a remarkable feature, the compensated/uncompensated strings have multiple winding numbers. This should be contrasted with the string formation in the conventional Abelian-Higgs model, where the strings with multiple winding numbers are hardly formed at the phase transition. We also found a tendency that the uncompensated strings evolve into the local strings by the long-range forces mentioned above. Indeed, at a later time, we observed that most of the uncompensated strings end up being the local string in the 2+1 D simulation.
In this paper, we explore the formation and the evolution of the string network by performing 3+1 D classical lattice simulation (see the pioneer works using lattice simulations for the Abelian-Higgs strings and the global strings, e.g., [21][22][23][24][25][26][27][28][29][30][31][32]). As we will see, the collision and the reconnection of the uncompensated strings are much more complicated than those seen in the 2+1 D simulation. 3 In particular, we observe the formation of "bridges" between two strings, which is a peculiar feature of the present model. The bridge formation makes the cosmological evolution of the string network highly complicated compared with the conventional string networks. The formation of the bridges also makes the combination process of the uncompensated strings into the compensated strings less efficient in the 3+1 D than that in the 2+1 D. The uncompensated strings potentially cause the axion domain-wall problem when this model is applied to the QCD axion. Thus, the 3+1 D simulation is imperative for the axion domain-wall problem where the PQ symmetry is identified with U(1) global . 4 The organization of the paper is as follows. In section 2, we review the setup of the model and the results in our previous work. In section 3, we simulate the collision of two strings and show the results. In section 4, we describe the setup and analysis to study the evolution of the string networks. In section 4.2, we perform simulation of the time-evolution of the string network. In section 5, we discuss the implications to the axion domain-wall problem. The final section is devoted to the summary.

Review of new string solutions in Abelian-Higgs model
In this section, we summarize the model setup.
We also review the compensated/uncompensated string solutions, as well as the results of the simulation in the 2+1 D spacetime [1].

Model
The action of the model is given by, 5 In the 2+1 D simulation, the strings appear as "particles" in the 2 D plane, and hence, the reconnection just looks like the binding of two particles. 4 See [33] for the string-wall network in the QCD axion model appearing in the non-Abelian gauge symmetry [34]. 5 The metric (−, +, +, +) is employed in this paper.

JHEP09(2020)054
Here, φ n (n = 1, 2) are the two complex scalar fields with U(1) local charge q n , A µ is the U(1) local gauge field, and e denotes the gauge coupling constant of U(1) local . We take q 1 and q 2 to be relatively prime numbers without loss of generality. The scalar potential V (φ 1 , φ 2 ) is given by where λ 1,2 and κ are real dimensionless constants, and η 1,2 are real constants with mass dimension one. The model possesses an U(1) global symmetry in addition to the U(1) local symmetry. Such U(1) global symmetry naturally appears as an accidental symmetry in the renormalizable theory when |q 1 | + |q 2 | is larger than 4. We assume λ 1 λ 2 > 4κ 2 so that both the φ 1,2 obtain non-vanishing vacuum expectation values (VEV's), The VEV's in eq. (2.4) break both the U(1) local and U(1) global symmetries spontaneously. Accordingly, there appear two NG modes, where one corresponds to the would-be NG boson (b), and another one is the gauge-invariant NG modes (a). We can extract these two modes from the phase components of the complex scalar fields, Here,ã 1,2 are pseudo-scalar fields, and f n ≡ √ 2η n (n = 1, 2) are the decay constants of them. The gauge-invariant and the would-be NG modes, a and b, are given by [8,9], We can see the component a is indeed invariant under the U(1) local gauge transformation, where α denotes the U(1) local gauge transformation parameter. The domain of the gauge invariant axion is given by, (2.8)

Uncompensated string
At the phase transition, we expect string formation due to the non-trivial homotopy group of the vacuum manifold, π 1 (U(1) local × U(1) global ) = Z × Z. In fact, we can obtain the JHEP09(2020)054 static string solutions under the ansatz, Here, (r, θ, z) are the radial distance, the azimuth angle, and the height in the cylindrical coordinate, respectively. The integers n 1 and n 2 denote the winding numbers of the strings consisting of φ 1 and φ 2 , respectively. The boundary conditions are for r → 0, and for r → ∞. See [1] for the numerical configurations of the static string solutions satisfying these boundary conditions. The asymptotic behaviors of the covariant derivatives of the complex scalar fields are given by, When they are not vanishing at r → ∞, those terms lead to logarithmic divergences of the string tension at r → ∞ via The string tension becomes finite only when D θ φ 1 and D θ φ 2 vanish, i.e.
In [1], we call the string solutions which do not satisfy eq. (2.17) but with n 1 > 0, the uncompensated strings. The string solutions satisfying eq. (2.17) are the local strings which we call the compensated strings. As in the case of the global string, the uncompensated string has a diverging string tension. Unlike the global string, however, it also has non-vanishing U(1) local gauge flux at the string core. Therefore, the uncompensated strings have properties in between the global string and the local string. The divergence of the string tension of the uncompensated string is cutoff by a typical distance between the strings. Thus, it can be formed at the phase transition in the early universe, although it would have an infinite tension if it is completely isolated.

JHEP09(2020)054
Around an uncompensated string with (n 1 , n 2 ), not only the would-be NG boson but also the physical NG winds non-trivially As we will discuss in section 5, N dw corresponds to the domain-wall number in the application to the axion models. Obviously, the compensated string has the vanishing domain-wall number, N dw = 0.

2+1 D simulation results
In [1], we studied the formation and the evolution of the strings in the radiation dominated universe by the classical lattice simulation in the 2+1 D spacetime. We took the parameters, as reference values. Due to η 2 /η 1 1, φ 1 obtains the VEV earlier than φ 2 in the evolution. We observed the formation of the compensated, the uncompensated, and the global strings after U(1) symmetry breaking. Here, the global string corresponds to the strings with n 1 = 0, which appears when φ 2 obtains the expectation value. As the gauge boson has become massive by that time, the global string does not have gauge flux at its core. In our simulation, we found that the winding number of the formed global strings are only ±1. Accordingly, the formed global strings have the domain-wall number N dw = 1. 6 The simulation also showed that most of the uncompensated strings evolve into the compensated strings by absorbing or releasing the global strings. Such behavior is understood by the long-range forces between uncompensated and global strings, as discussed in [1]. Thus, the 2+1 D simulation suggested that the string network evolves into the one with only the compensated strings (N dw = 0) and the global strings (N dw = ±1).
In our 2+1 D simulation, we found that the strings with |N dw | > 1 disappear at late times, which implies that the model might be free from the axion domain-wall problem. However, this is not conclusive, because strings reconnect in the 3+1 D spacetime, which cannot be simulated in the 2+1 D simulation. In the following, we perform the 3+1 D simulation.

Colliding two strings
In this section, we simulate the collision between a compensated string and an uncompensated string with various winding numbers as a demonstration. We take the parameters 7 JHEP09(2020)054

Simulation setup
The Euler-Lagrange equations of the scalar fields and the gauge field are given by, For the colliding string simulations, we take the Lorenz gauge, ∂ µ A µ = 0. We prepare the initial configurations for the scalar fields and the gauge field, realizing the compensated and the uncompensated strings with n 1 = 1 and arbitrary n 2 . Following our previous paper [1] (see its section IIB), we first solve the one-dimensional boundaryvalue problem for φ n and A µ , and obtain the axially symmetric static solution of the strings along the z-axis. We put two of the resultant string solutions on the x-axis with a separation d in the three-dimensional box, such as (d/2, 0, z) and (−d/2, 0, z). We rotate them by ±α/2 with respect to the x-axis so that the relative angle between them becomes α. We also perform the Lorentz-boost for each string with velocity (v 1 , v 2 ) = (−v, +v) respectively so that they collide with each other at the coordinate origin of the box. Taking d much larger than the width of the strings, i.e., d O η −1 1,2 , we approximate the initial configuration with the two Lorentz-boosted strings by φ n (t 0 , µ (x) (p = 1, 2) represent each field configuration in a Lorentz-boosted string at (d/2, 0, z) and at (−d/2, 0, z), respectively. In the following simulations, we take d = 20η −1 1 , α = 0.2π and v = 0.7, except for figure 3 in which we set α = 0.06π.
In our simulations, we impose the boundary conditions on each side of the 3 D spatial box so that the scalar and gauge fields on them are always given by the superposition of the Lorentz-boosted strings constructed above. Hence, we must terminate the simulation when any effects from the impact point of the collision come to the boundaries.

Result
Case 0: (n 1 , n 2 ) = (1, 2). Before discussing the collision with the uncompensated string collision, let us see the collision between two compensated strings. The snapshots of the collision are shown in figure 1. 8 The upper (lower) figures show the string configurations made by φ 1 (φ 2 ). We observe that the two compensated strings simply reconnect and go away from each other as in the conventional Abelian-Higgs model. Case 1: (n 1 , n 2 ) = (1, 1). In figure 2, we show the snapshots of the collision between a compensated string and an uncompensated string with (n 1 , n 2 ) = (1, 1) at the late time. As shown in the left panel, the φ 1 -strings 9 reconnect and move away as in the previous case. On the other hand, as shown in the right panel, the φ 2 -strings do not simply go away. JHEP09(2020)054 After the collision, the n 2 = 1 φ 2 -string partially reconnects with n 2 = 2 φ 2 -string and form a combined system connected by a n 2 = 1 φ 2 -string.
A schematic picture of this final state is shown in figure 3. One red rod and one green rod are the φ 1 -string with n 1 = 1 and the φ 2 -string with n 2 = 1, respectively. The direction of the arrow denotes the direction of the magnetic flux of the strings. We call the string segment connecting two strings the "string bridge". The string bridge is formed by the reconnection of the bunch of the strings with different n 2 .
When the collision angle, α, is small, we find that even the second reconnection takes place (figure 4). By the tension of the string bridge, two bunches of the strings do not leave away, and they are attracted to each other after the first reconnection. Once they collide and reconnect again (t = 87.5 in the figure), they divide into two bunches of the strings, which are the same string bunches before the first collision. Then, they eventu-JHEP09(2020)054 Figure 3. A schematic picture of the string system after the collision. The red rod is a φ 1 -string with n 1 = 1. The green rod is a φ 2 -string with n 2 = 1. The arrows on the strings denote the magnetic flux in the strings. We call the string segment connecting two strings the string bridge. The string bridge is formed by the reconnection of the bunch of the strings. The strings are cut for drawing figures.
ally move away from each other as if they pass through each other. See figure 5 for the schematic diagram.
Case 2: (n 1 , n 2 ) = (1, 0). We also simulate the case with (n 1 , n 2 ) = (1, 0) uncompensated string. As shown in figure 6, at the collision, the bunch of the φ 2 -strings reconnect and the string bridges consisting of two φ 2 -strings are formed. 10 After the formation of the bridges, the two bridges are separated and then keep a distance between them. This phenomenon is understood by the repulsive force between the two φ 2 -strings. See figure 7 for the schematic diagram.
The colliding simulations have shown the formation of the string bridges, implying that the complicated string network emerges in the cosmological setup. To see this and the fate of the uncompensated strings in the network at the late time, we perform the 3+1 D classical lattice simulations in the radiation dominated universe in the next section.

Cosmological evolution of string network
We explore the formation and the evolution of cosmic strings by the 3+1 D classical lattice simulations. In this section, we take the parameters, as reference values. These parameters are the same as those in the previous work of the 2+1 D simulation [1] and also see eq. (2.20).

Setup
We perform simulations in the radiation dominated universe whose spacetime metric is given as   where τ is the conformal time and a(τ ) is the scale factor. In the simulations, we use the temporal gauge, A 0 = 0, and hence, the field equations are given bÿ   Here, H = aH denotes the conformal Hubble parameter, i, j, k = 1, 2, 3, and the dot denotes the derivative with respect to the conformal time.

3)
The initial conditions and the spatial periodic boundary conditions are the same as in our previous study [1]. As the initial conditions, we impose the random values for each grid point, φ(t 0 , x) = δφ(x), which obey the Planck distribution with the temperature T = √ 3η 1 . The time derivatives,Ȧ i (t 0 , x), are determined by solving the constraint equation, and as A i (t 0 , x) can be freely chosen, we set A i (t 0 , x) = 0. The time-evolution is calculated by the Leap-Frog method, and the spatial derivatives are approximated by the second-order finite-difference. The simulation parameters are summarized in table 1.

Result
We show the 3 D spatial snapshots of the φ 1 -and the φ 2 -strings at the simulation end in figure 8. The left (right) panel shows the configuration of the φ 1(2) -strings. 11 As we assume η 1 = 4η 2 , the width of φ 1 -string is roughly 1/4 of that of φ 2 -string and thus the φ 1 -strings look thinner. Note that the physical size of the lattice spacing reaches to the width of φ 1 -string at the simulation end, given as ∆x ∼ 1.2η −1 1 . We observe the larger number of strings in the φ 2 sector than in the φ 1 sector, indicating that there are many global strings with n 1 = 0 and n 2 = 0. See, e.g., a φ 2 -string connecting from the bottom surface to the left one near the coordinate origin in the right panel; we cannot observe the accompanying string in the φ 1 sector in the left panel. This is again due to our arrangement of the size of the VEVs, η 1 η 2 . In the right panel of figure 9, we show the string core of the φ 2 -string. The higher winding strings with n 2 > 1 are represented as a bundle of strings with n 2 = 1, thus looking thick. The higher winding strings are expected to be accompanied with φ 1 -strings. To confirm this behavior, we count the number of φ 2 -strings along with a φ 1 -string, as shown in the left panel of figure 9. Hereafter, n 1 denotes its absolute value. In the simulation, we find all the φ 1 -strings have n 1 = 1. In this plot, a segment of a φ 1 -string coloured in red indicates that there are four φ 2 -strings along it distancing within 4η −1 1 . 12 Along the φ 1 -string, we find that the sign of the winding numbers of the φ 2 -strings are the same with that of the φ 1 -string. Thus, hereafter, we also take n 2 as its absolute value. The segments colored in red are "compensated" since they satisfy (n 1 q 2 − n 2 q 1 ) = 0. On the other hand, there are three and five φ 2 -strings along the blue and green segments, respectively. The segments colored in blue or green are "uncompensated" with n 1 q 2 − n 2 q 1 = ±1, 11 See also http://numerus.sakura.ne.jp/research/open/NewString3D/ for supplemental material. 12 For the choice of the distance, see the appendix A.  respectively. At the point where the color is changed, e.g. from red to blue, it is expected that there are either the string bridges as discussed in section 3 or thin ambient global strings with no accompanying φ 1 -strings. See figure 10 for the schematic picture of an example string faction.
Let us see the time-evolution of the string segments of various kinds. From the left panel of figure 9, we compute the lengths of strings with n 1 = 1 as L ng , which measures the total length of the non-global strings. We also compute L n 2 , the length of the φ 1string segment which is accompanied by the φ 2 -strings with n 2 . In figure 11, we show the time-evolution of the fractions, The results are averaged over ten realizations and the error bars indicate the standard deviation. In this plot, the red curve shows the time-evolution of the length of the compensated strings, we call it R c in the following. At late times, we find that about 60% of the total non-global strings settle into the compensated strings. This result is rather surprising in viewing the complicated reconnection processes in the 3+1 D spacetime discussed in the previous section. In the 2+1 D simulation [1], the fraction R c settles to R c ∼ 0.8 at the end of the simulation. The reduction of R c from R c ∼ 0.8 (in 2+1 D) to R c ∼ 0.6 (in 3+1 D) shows the effects of the complicated reconnection processes which make the combination process less efficient.
In addition to the compensated/uncompensated strings, there are many φ 2 -strings which are not accompanied by the φ 1 -strings. Those are the global strings, i.e., n 1 = 0. In the simulations, we find no global strings with n 2 > 1. To see the fraction of the global strings, we compute the length of the global string, L g , and define the fractions, Here, the total string length is defined by In figure 12, we show the time-evolution of the fractions. The figure shows that the stringnetwork is dominated by the global strings although the fraction decreases at late times. The latter fact indicates that some of the global strings are absorbed by the uncompensated strings to form the compensated strings, and some of them are reconnected with each other to form loops. 13

Physical Nambu-Goldstone winding
Let us discuss the time-evolution of the winding numbers of the gauge-invariant NG mode, N dw , defined in eq. (2.19). As we will discuss in the next section, the strings with |N dw | > 1 potentially cause the axion domain-wall problem if they remain abundantly at late times.
In the 2+1 D simulation in [1], we found that the number of the strings with |N dw | > 1 vanishes at late times. If this is the case, the axion model based on the present set up can 13 The reconnection and the loop formation of the global strings with n2 = 1 take place in the usual way, regardless of the existence of φ1-strings. Hence, the global strings can disappear faster than the non-global strings. be free from the axion domain-wall problem (see discussion in section 5). In the 3+1 D, however, the string combination processes are less efficient, and hence, it is imperative to estimate how large fractions of the strings have |N dw | > 1 at late times in the 3+1 D simulation.
To compare with the results of the 2+1 D simulation, we compute the ratio, in 3+1 D, where L others is the total length of the strings with |N dw | > 1. 14 This quantity is analogous to R dw in [1]. In figure 13, we show the time-evolution of R dw and R c . The figure shows that the R dw is suppressed over the simulation time, but it does not vanish at late times. The non-vanishing R dw is due to the complexity of the string network evolution in the 3+1 D space-time.
In the figure, we also show the evolution of the rescaled fraction, just for convenience. As L tot is dominated by the contribution of the global strings, R others is more suitable to see the fraction of the strings which have non-trivial NG winding around the non-global strings. A caveat here is that R dw and R others do not count the uncompensated strings with |N dw | > 1 precisely. In the 2+1 D simulation, the strings are point-like objects in the 2 D space. In this case, R dw can be defined as the number fraction of the separate strings with |N dw | > 1. In the 3+1 D simulation, on the other hand, the strings are connected as a network, and hence, the number of the φ 2 -strings, or the winding number n 2 , changes along a φ 1 -string. This is the reason why we define the various fractions based on the length of the segments in eqs. (4.5), (4.6), (4.8), and (4.9). Besides, the winding numbers of the segments are determined by counting the number of φ 2 -strings with n 2 = 1 along a φ 1 -string within 4η −1 1 . Thus, for example, R  Figure 13. The time-evolution of R c , R dw , and R others , for (q 1 , q 2 ; η 2 /η 1 ) = (1, 4; 0.25). n 2 = n 2 , when the φ 2 -string around the φ 1 -string is loosely bounded or when there are some φ 2 -strings in the vicinity accidentally. Therefore, R dw should be taken as an upper limit of the fraction of the segment with N dw > 1, since some portion of R 3,4,5 which contribute to |N dw | ≤ 1 are counted as R dw . 15 Finally, let us show the parameter dependence. In figure 14, we show the time-evolution of R c , R dw and R others for η 2 /η 1 = 0.75, 0.50, 0.25 (from left to right). We observe the increase of R c for a smaller η 2 /η 1 . This shows that more uncompensated strings are bounded into the compensated strings with N dw = 0 thanks to more abundant global strings at the phase transition of φ 2 for η 2 η 1 . On the other hand, we do not observe significant tendencies on η 2 /η 1 for R dw and R others compared with the increase of R c . These results indicate that R dw and R others could have sizable contributions from the contamination discussed above.
We also show dependence on the charge ratio. figure 15 shows the time-evolution of R c , R dw and R others for q 2 = 2, 3, 4 (from left to right). We observe that R c decreases for a larger q 2 /q 1 . This is expected because the uncompensated strings require more φ 2 -strings to evolve into the compensated strings for a larger q 2 /q 1 .
Finally, let us comment that the case with the opposite charge ratio, q 2 /q 1 = 1/4, while 15 As R3 + R4 + R5 R dw , the relative uncertainty of R dw caused by the contamination is larger than that of R3 + R4 + R5 (see also the appendix A). Rc, Rdw, Rothers tη1 Rc Rdw Rothers Figure 15. The time-evolution of R c,dw,others for various q 2 . We take q 2 = 2, 3, 4 from left to right. The rightmost panel is the same with figure 13. keeping the ratio of the VEVs, η 2 /η 1 = 0.25. In this case, we find that the strings with (n 1 , n 2 ) = (1, 0) and (0, 1) remain at the end of the simulations, although we do not show the results. Similar results are also obtained in the 2+1 D simulation. This behavior can be understood that we need four φ 1 -strings and one φ 2 -string to form the compensated string. This is difficult as φ 1 -strings are formed at T ∼ η 1 , while the φ 2 -strings are formed at much lower temperature, T ∼ η 2 . As the global strings, formed at T ∼ η 2 , have N dw = 4, we expect R dw is close to 1, for a model with η 1 > η 2 and q 1 > q 2 .

Implication to axion domain-wall problem
As we have seen in the previous section, R dw becomes small at late times for q 1 = 1, q 2 > 1, and η 1 η 2 . Besides, we found that the global strings with n 2 = 1 is the most abundant at late times. In this section, we discuss implications of these results on the axion domain-wall problem.
First, let us consider the axion model in which the Peccei-Quinn (PQ) symmetry is identified to the U(1) global symmetry [8][9][10] ( [35] for review). As in the KSVZ axion model [36,37], we introduce N 1 and N 2 of the QCD colored vector-like multiplets, (Q if i ,Q if i ) (i = 1, 2, f i = 1-N i ), and couple them to φ 1 and φ * 2 via, Here, we suppress the Yukawa coupling constants and the flavor indices. The anomaly-free condition of the U(1) local gauge symmetry requires, which can be solved by 16 The integer N m is a free parameter of the axion model. 16 We define q1 and q2 being relatively prime numbers.

JHEP09(2020)054
Once the U(1) symmetries are broken spontaneously, the vector-multiplets become massive, and they leave the coupling of the gauge-invariant NG boson in eq. (2.6) with the QCD gluon, This shows that the gauge-invariant NG boson a plays the role of the QCD axion.
As we have discussed in the previous section, the string network for q 1 = 1, q 2 > 1, and η 1 η 2 is dominated by the global strings (with only n 2 = 1). The network also has a small contribution from the uncompensated segments with |N dw | > 1. Around these strings, the axion (i.e. the gauge-invariant NG boson) winds N dw times in its domain defined in eq. (2.8).
When the temperature of the universe decreases below the QCD scale, the axion obtains a periodic potential by the effect of eq. (5.4). The period of the potential is F a /N m , and hence, it has a discrete Z Nm shift symmetry. Due to the axion winding around the strings, the axion potential causes the energy contrasts around the strings. By remembering the periodicity of the axion potential, the energy contrasts have the periodicity of N m × |N dw | which results in N m × |N dw | domain-walls attached to the strings.
For N m > 1, the string-wall network is stable since the domain-walls formed in the above process are associated with the spontaneous breaking of the Z Nm symmetry. In this case, all the strings are attached by at least N m of the domain walls, and they are pulled into multiple directions. Once stable domain-walls are formed, they immediately dominate the energy density of the universe and lead to the overclosure. To avoid such a problem, we consider a model with N m = 1 in the following.
For N m = 1, most strings are attached by only one domain-wall, since the string network is dominated by the global strings which have |N dw | = 1 in the present setup. Since a string which are attached by only one domain-wall is pulled into one direction, it can annihilate with another string with opposite winding which is attached to the other side of the domain-wall [38,39]. Thus, the string network with |N dw | = 1 does not cause the domain-wall problem.
Even for N m = 1, however, the string network of the present model potentially results in stable domain-walls since it has segments with |N dw | > 1. In fact, if the multiple domain-walls connect string segments with |N dw | > 1, the strings are pulled into multiple directions and form a stable string-wall network. As we have seen, however, R dw 1 and the most strings are the global string with N dw = 1 for η 2 /η 1 1, q 1 = 1, and q 2 > 1. Thus, the stable string-wall networks made by the string segments with |N dw | > 1 are surrounded by many global strings with |N dw | = 1. Therefore, the stable string-wall networks are immediately shredded when the global strings pass through the domain-wall, and disappear before they dominate the energy density of the universe. 17 These arguments strongly suggest that the axion model based on the model with η 2 /η 1 1, q 1 = 1, q 2 > 1, and N m = 1 is free from the axion domain-wall problem.

JHEP09(2020)054
Let us finally comment that the absence of the domain-wall problem cannot simply result from the fact that the model has no discrete symmetry for N m = 1. For example, if we take q 1 /q 2 = 1/4, we expect R dw is close to 1, and the string network is dominated by the global strings with N dw = 4. Thus, in this case, the model suffers from the domain-wall problem even in the absence of the discrete symmetry.

Summary
In this paper, we studied the formation and evolution of the string network in the model with one local U(1) and one global U(1) symmetries by performing the classical lattice simulations in the 3+1 D spacetime. As shown in section 3, the reconnection processes of the present model are highly complicated. In particular, the uncompensated strings are not simply reconnected, but yielding bridges which connect the two bunch of the strings (figure 3). We also studied the time-evolution of the string network in the 3+1 D spacetime. Despite its complexity, we observed that more than 50% of the uncompensated strings settle into the compensated strings even in 3+1 D ( figure 11 and figure 12).
As a result of the network simulations, we found that the string-network is dominated by the global string (figure 12) at late times. Besides, we also found that the fraction of the string segments with |N dw | > 1 is suppressed, i.e. R dw 1 for η 2 /η 1 1, q 1 = 1 and q 2 > 1 (figures 13 and 15). These results strongly suggest that the PQ axion model based on the present setup can be free from the axion domain-wall problem when we take η 2 /η 1 1 with q 1 = 1, q 2 > 1, and N m = 1. 18 This type of the axion model naturally explains the origin of the PQ symmetry from the gauge symmetry [6][7][8][9][10]34]. The absence of the domain-wall problem makes this type of the axion more attractive.
There are several open questions. For the evolution of the string network, it is important to see whether the network follows the scaling law or not. Due to the formation of the string bridges, the scaling law may not be observed even for the time duration simulated in this paper. Another question is the production of the gravitational waves. Their spectrum may have a feature arising from the complexity of the network associated with the combination process of the uncompensated strings, and therefore there is a possibility that we can distinguish them from the gravitational waves from the other kind of strings like Abelian-Higgs strings. In the present work, we do not discuss these issues further and leave them for future studies.

A Computation of winding number
In section 4, we computed the winding number of φ 2 -string along a φ 1 -string as shown in figure 9. The winding number is given as the number of φ 2 -strings with n 2 = 1 lying along a φ 1 -string within a separation s. In the main text, we fixed s = 4η −1 1 . In figure 16, we show the dependence of the computed winding numbers on the parameter s. The bottom-left panel in this figure shows the case with s = 4η −1 1 and is the same as figure 11. The top-left, top-right, and bottom-right figures show the cases with s = 16η −1 1 , 6η −1 1 , and 2η −1 1 , respectively. By increasing s from 4η −1 1 to 16η −1 1 , we find that R others increases and R 3,4,5 decrease. This trend can be explained by the overcounting of global strings around the φ 1 -strings for a larger s. On the other hand, for s = 2η −1 1 , R 4 decreases compared to the case of s = 4η −1 1 . By remembering that four φ 2 -strings within a distance of O η −1 1 around a φ 1 -string should combine into a compensated string, the decrease of R 4 for s = 2η −1 1 indicates the underestimate of the φ 2 -strings in a compensated string due to too small choice of s. In our analysis, we adopt s = 4η −1 1 as an optimal choice for R 4 .
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.