New type of string solutions with long range forces

We explore the formation and the evolution of the string network in the Abelian Higgs model with two complex scalar fields. A special feature of the model is that it possesses a global U(1) symmetry in addition to the U(1) gauge symmetry. Both symmetries are spontaneously broken by the vacuum expectation values of the two complex scalar fields. As we will show the dynamics of the string network is quite rich compared with that in the ordinary Abelian Higgs model with a single complex scalar field. In particular, we find a new type of string solutions in addition to the conventional Abrikosov-Nielsen-Olesen (local) string solution. We call this the uncompensated string. An isolated uncompensated string has a logarithmic divergent string tension as in the case of the global strings, although it is accompanied by a non-trivial gauge field configuration. We also perform classical lattice simulations in the 2 + 1 dimensional spacetime, which confirms the formation of the uncompensated strings at the phase transition. We also find that most of the uncompensated strings evolve into the local strings at later time when the gauge charge of the scalar field with a smaller vacuum expectation value is larger than that of the scalar field with a larger vacuum expectation value.


Introduction
Cosmic strings (also known as vortex solutions) are one-dimensional topological defects which appear in various contexts of particle physics and condensed matter physics. The simplest theoretical framework to describe the string (vortex) formation is the Abelian Higgs model, where the U(1) gauge symmetry, U(1) local , is spontaneously broken by a vacuum expectation value (VEV) of a complex scalar field. The string appearing in this model is called the Abrikosov-Nielsen-Olesen (ANO) string [1]. At the phase transition of the U(1) local breaking, the strings are formed and make up a web-like structure, socalled the string network (see e.g. the textbook [2]). Numerical simulations have widely investigated the properties of the interactions and the evolution of the string network (see e.g. ).
In this paper, we discuss the formation and the evolution of the string network in the Abelian Higgs model with two complex scalar fields. In particular, we consider a model with an additional global U(1) symmetry, U(1) global . Such an additional global symmetry naturally arises when the ratio of the absolute U(1) local charges of the two complex scalar fields is larger than 3 (or less than 1/3). In this case, there is no gauge-invariant interaction term which breaks U(1) global at the renormalizable level, and hence, it emerges as an accidental symmetry. This class of models has been considered, for instance, in the context JHEP02(2020)058 of the axion models where the Peccei-Quinn symmetry [29][30][31][32] is protected by (abelian) gauge symmetries [33][34][35][36][37][38] from the quantum gravity effects [39][40][41][42][43][44]. This type of models is also discussed in condensed matter physics to understand the physics of high-temperature superconductivity [45]. 1 As we will see, the dynamics of the string network is much richer compared with that in the ordinary Abelian Higgs model with a single complex scalar field. In particular, we observe the existence of a new type of string solutions in addition to the ANO (local) string. We call the new type of the string solutions, the "uncompensated" strings. Around an uncompensated string, the covariant derivatives of the complex scalar fields are not canceled by the gauge field configuration at the large distance from the core of the string. Accordingly, an isolated uncompensated string has a logarithmically divergent string tension as in the case of the global strings. In the early universe, those strings can be formed at the phase transition, where the divergence is cutoff by a typical distance between the strings of the order of the Hubble length. Our classical lattice simulations (in the 2 + 1 spacetime dimensions) confirm the formation of the uncompensated strings at the phase transition.
As a characteristic feature of this system, the strings with winding numbers larger than 1 appear ubiquitously. This feature should be contrasted with the string formation in the Abelian Higgs model with a single complex scalar field in which it is quite rare to have a cosmic string with a large winding number in the cosmological context. We also find that the uncompensated strings have long-range repulsive and attractive forces. With the long-range forces, we find that the uncompensated strings tend to combine into the local ANO strings in the evolution of the string network when the gauge charge of the scalar field with a smaller VEV is larger than that of the scalar field with a larger VEV.
The organization of the paper is as follows. In section 2, we discuss the static string solutions analytically. In section 3, we show the results of the classical lattice simulations of the formation and the evolution of the string network in the 2+1 dimensional spacetime. In section 4, we derive the long-range force between the uncompensated string and the global string in an analytical way. The final section is devoted to our conclusions and discussions.

String solutions in Abelian Higgs model
In this section, we discuss the string solutions in the Abelian Higgs model with two charged scalar fields. We assume that both the scalar fields obtain non-vanishing VEVs. A specific feature of the model is that it possesses the U(1) global symmetry in addition to the U(1) local symmetry. We here discuss the static field configurations in an analytical way with the expansion of the universe being neglected.

Model
The action of the model is given by 2

JHEP02(2020)058
with the Lagrangian densities Here, φ n (n = 1, 2) denote the two complex scalar fields and is the field strength of the gauge field A µ of the U(1) local symmetry. The covariant derivative is defined by where q n is the charge of φ n , and e is the gauge coupling constant. In what follows, we normalize q 1 and q 2 so that they are relatively prime integers without loosing generality. The scalar potential is taken to be where λ 1,2 (> 0) and κ are real valued dimensionless coupling constants, and η 1,2 are real valued constants with a mass dimension. This potential possesses two U(1) symmetries under the phase rotations of the two complex scalar fields. One of the linear combination of these symmetries is identified as the U(1) local gauge symmetry, while the other symmetry is a global U(1) global symmetry. The global U(1) global symmetry naturally appears as an accidental symmetry in the renormalizable theory when |q 1 | + |q 2 | is larger than 4. 3 We consider the case with λ 1 λ 2 > 4κ 2 . In this case, both the U(1) symmetries are spontaneously broken by the VEVs [53], φ n = η n , (n = 1, 2) . (2.7) The Goldstone modes in this system are decomposed into the gauge-invariant Goldstone boson and the would-be Goldstone boson [36,37]. To see this, let us define the phase component fieldsã n (n = 1, 2) by where the decay constants are given by f n = √ 2η n . The domains of the phase component fields are givenã respectively. The U(1) local gauge symmetry is realized by the shifts ofã n , where α is a local parameter of the U(1) local transformation.
JHEP02(2020)058 The kinetic terms ofã n are given by (2.11) In the final expression, we use the Goldstone bosons defined by (2.12) The field b is the would-be Goldstone boson which is absorbed by the U(1) local gauge boson in the unitary gauge. The mass of the gauge boson is given by, The field a is the Nambu-Goldstone boson of the global U(1) global symmetry, which is invariant under the U(1) local symmetry (see eq. (2.10)).
In figure 1, we show the gauge orbits of U(1) local for (q 1 , q 2 ) = (1, 4) as the blue lines. As seen in the figure, the shift of b from 0 to 2π corresponds to the shift of a 1 from 0 to 2π q 1 and the shift of a 2 from 0 to 2π q 2 . The direction perpendicular to the blue lines corresponds to the gauge-invariant Goldstone boson, a. The domain of a is given by the interval between the blue lines, (2.14) JHEP02(2020)058

Field equations of the static string solutions
The Euler-Lagrange equations of φ n and the gauge field are given by, Hereafter, we take the temporal gauge, A 0 = 0, which reduces the field equations tö and yields the constraint equation, The dot denotes the time-derivative. 4 As an ansatz for the static string solutions, we assume

21)
A θ (r) = 1 e ξ(r) , A r = A z = 0 , (2.22) in the cylindrical coordinate where (r, θ, z) is the radial distance, the azimuth angle, and the height, respectively. The integers n 1,2 denote the winding numbers of the strings consisting of φ 1 and φ 2 , respectively. Under the ansatz, the field equations in eqs. (2.17) and (2.18) are reduced to

JHEP02(2020)058
Here, we rescale the radial coordinate r to R by R = r/r 0 where The dash in the field equations denotes the derivative with respect to R. We also define for n = 1, 2.

Static string solutions
To obtain the static string solutions, we consider the following boundary conditions, for R → 0, and for R → ∞. Those boundary conditions are the same in the case of the static string solution of the Abelian Higgs model with one complex scalar field. In the region of R → ∞, h 1 (r) and h 2 (r) get close to unity, the field equation of ξ in eq. (2.25) becomes for R → ∞. This field equation has an asymptotic solution, and hence, the gauge field becomes The asymptotic behaviors of the covariant derivatives of φ n are then given by As the covariant derivatives contribute to the string tension via the string tension diverges at r → ∞ unless both D θ φ 1 and D θ φ 2 vanish. Eqs. (2.33) and (2.34) show that this is possible only when the winding numbers satisfy, The string solution which satisfies eq. (2.36) has a finite string tension. The string solution which does not satisfy eq. (2.36), on the other hand, has a logarithmically divergent tension as in the case of the global string. In the following, we call the former strings the compensated (local) strings, while the latter the uncompensated strings. It should be emphasized that the uncompensated strings are new-type of the string solutions which are absent in the conventional Abelian Higgs model with a single complex scalar field.
To show examples of the string profiles in our model, we solve the equations (2.23)-(2.25) with the boundary conditions given in eqs. (2.28)(2.29) and (2.31). We impose the asymptotic boundary conditions at R = 80. In figure 2, we show the static string solutions for various values of (n 1 , n 2 ) in the case of (q 1 , q 2 ) = (1, 4). We take other parameters as e = 1/ √ 2, η 1 /η 2 = 4, λ 1,2 = 1, and κ = 0.  uncompensated strings. The figure shows that the string core of φ 2 becomes smaller for a larger n 2 . The figure also shows the gauge field configuration converges to the asymptotic value in eq. (2.32). We see that the gauge field configuration clings to the string solution of φ 1 for a small R and it converges to the asymptotic value at a large R.
The solution for n 2 = 0 requires an explanation. In this case, the condition, h 2 (r) = 0, at R → 0 is not required and we assume the Neumann boundary condition for h 2 . Even in this case, h 2 has a non-trivial configuration at around the string solution with n 1 = 0 (see eq. (2.24)).
In figure 3, we show how each type of the string solutions winds in the domain map of (ã 1 ,ã 2 ) in the case of (q 1 = 1, q 2 ) = (1, 4). The figure shows that the configuration of the compensated string, i.e. (n 1 , n 2 ) = (1, 4), coincides with the orbit of the gauge transformation (see figure 1). The uncompensated strings, on the other hand, wind also in the direction of the gauge-invariant Goldstone boson. In fact, the string configuration with (n 1 , n 2 ) winds the Goldstone direction as which vanishes only for the compensated strings (see eq. (2.12)).

Preparation for simulations
Here, let us summarize the setup of our numerical simulations. To take into account the cosmic expansion in the radiation dominated universe, we use the conformal time, τ , with which the metric is given by, where a(τ ) is the scale factor. In the expanding universe, the field equations are modified tö where H denotes the conformal Hubble parameter and the dot denotes the derivative with respect to the conformal time.
We also add the thermal mass term to the scalar potential The thermal mass term stabilizes the symmetry enhancement point, φ 1,2 = 0, and hence, the U(1) local × U(1) global symmetry is restored at the high temperature, T η 1,2 . In the followings, we discuss the formation/evolution of the strings for, as reference values. We also discuss the parameter dependencies at the end of this section. At the initial time, we impose the random values for φ i (t in , x) = δφ(x) so that the distributions of the scalar fluctuations are given as the Planck distribution with the temperature T = √ 3η 1 . Then we can determineȦ i (t in , x) by solving eq. (2.19) with the Fourier transformation. At this stage, we can freely choose A i (t in , x), so we simply assume A i (t in , x) = 0. Throughout the simulations, we assume the radiation-dominant universe. We performed the simulations in a comoving box with periodic boundaries, and 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. The string configuration at the spot 1 has the winding numbers, (n 1 , n 2 ) = (1, 4). This configuration corresponds to the compensated (local) string (see eq. (2.36)). As discussed in section 2, the non-trivial configuration of the gauge field cancels the covariant derivatives of the both complex scalars. The string configurations at the spots 2 and 3 are, on the JHEP02(2020)058 In the present setup, we take η 1 /η 2 = 4, and hence, the phase transition of the U(1) local ×U(1) global breaking takes place in stages. At the first stage, the U(1) local breaking takes place at T ∼ η 1 . At the second stage, the U(1) global breaking takes place at T ∼ η 2 . At the first stage, the string configurations with n 1 = 0 (mostly n 1 = 1) are formed in the similar manner to the original Abelian Higgs case. Then, around the string configurations with n 1 = 1, φ 2 also takes non-trivial configurations with various n 2 at the second stage. These configurations become either the compensated and the uncompensated strings. At the second stage, it is also possible that φ 2 takes non-trivial configurations with n 2 = 0 in the region where φ 1 is trivial, i.e., φ 1 = η 1 . These configurations are the global strings, in which φ 1 = η 1 is barely affected by the non-trivial configurations of φ 2 due to the hierarchy η 1 η 2 .

Evolution of string network
After the formation of the string configurations, the string network exhibits complicated evolution. In figure 5, we show the time evolution of the numbers of the strings with various n 2 in the computational domain. The results are obtained by averaging over 10 realizations, and the error bars show the statistical variance of the number of strings arising from the randomness of the initial conditions. To identify each string and compute its winding number for φ 1 and φ 2 , we use the algorithm developed by ref. [25], and see also ref. [27]. Throughout this paper, we regard two strings separated within 20η −1 1 as one string. We define N c as the number of compensated strings which have (n 1 , n 2 ) = (1, 4), N u as that of uncompensated strings with n 2 q 1 − n 1 q 2 = 0 with n 1 = 0, N g,1 as that of global strings JHEP02(2020)058 with (n 1 , n 2 ) = (0, 1), and N g,>1 as that of global strings with n 1 = 0 and n 2 > 1. The total number of string, N tot , is given as N tot = N c + N u + N g,1 + N g,>1 . Note that in our simulations we observed no strings with n 1 > 1. We also show the evolution of the number fractions of the compensated strings, f c = N c /N tot (red line), the uncompensated strings, f c = N u /N tot (green line), and the global strings, f g = (N g,1 + N g,>1 )/N tot (blue line) in the middle panels. The right panels show the number ratio of the compensated strings to all the strings, (3.6) and the number ratio of the strings with |n 2 q 1 − n 1 q 2 | > 1, Here, R dw includes contributions both from the uncompensated and the global strings with |n 2 q 1 − n 1 q 2 | > 1. This ratio is particularly important when the gauge-invariant Goldstone boson a plays the role of the QCD axion. In fact, the so-called the axion domain wall number around the cosmic strings are given by (a multiple of) |n 2 q 1 − n 1 q 2 |, and the axion domain wall problem can be avoidable only for |n 2 q 1 − n 1 q 2 | = 0, 1, if we assume that the phase transition of the U(1) symmetries take place after the end of inflation (see the appendix A). The figure shows that the number fraction of the compensated local strings increases in time, while those of the uncompensated and the global strings decrease. This means that the uncompensated strings and the global strings tend to combine together to form the compensated strings. After a while, the number of compensated strings tends to be converged, since the cosmic expansion separates the strings from each other and makes the compensation difficult. In section 4, we discuss how the long-range force between the uncompensated and the global strings appears.
The simulation result in figure 5 also shows that the fraction of the compensated strings converges to a value smaller than 1, and hence, there remain uncompensated/global strings at later times. We also found that R dw → 0 in this simulation, which indicates that all the remaining strings in the network have the winding number either 0 or 1 in the gaugeinvariant Goldstone direction. This result has important implications for the axion model building (see the appendix. A).
In the simulation of figure 5, we do not include the thermal mass term in eq. (3.4). Accordingly, the string formations of φ 1 and φ 2 take place simultaneously. In a realistic situation, however, the string formation of φ 2 takes place much later than that of φ 1 . In figure 6, we show the same figures in figure 5 with the effects of the thermal mass in eq. (3.4) taken into account. In this setup, the formation of the second string is delayed. Hence we performed the simulations with a larger box whose initial size is given by 80H −1 in and the grid size being 4096 2 . Even in this case, the winding number of most of strings is n 2 = 4 and thus they are compensated. So we can conclude that the string network at late time is not sensitive to the presence of the thermal effect. From now on, we do not take the thermal mass into account for the following simulations.

Parameter dependence
So far, we have fixed the model parameters to the values in eq. (3.5). Here, we briefly discuss parameter dependencies of the behavior of the string network.
To see the dependence on the charge ratio, we have performed numerical simulations for (q 1 , q 2 ) = (1, 1) and (1,2), while the other parameters are fixed. In these cases, we have confirmed that the string networks show similar behavior to the case of (q 1 , q 2 ) = (1, 4), and the number ratio of the compensated strings increases in time as shown in figure 7. In fact, the compensated strings are dominated in the late time and thus most of strings take n 2 = q 2 .
We have also performed numerical simulations for (q 1 , q 2 ) = (4, 1) with the other parameters fixed. In this case, we found that the string configurations have either (n 1 , n 2 ) = (1, 0) or (n 1 , n 2 ) = (0, 1) even at late times, and they do not combine into the compensated strings unlike the case of (q 1 , q 2 ) = (1, 4) as shown in the bottom panels of figure 7. This can be understood by the fact that the long-range force between the uncompensated and the global strings are proportional to the covariant derivative of φ 2 around the uncompensated strings (see section 4). In fact, D θ φ 2 around the string configurations of (n 1 , n 2 ) = (1, 0) is smaller for (q 1 , q 2 ) = (4, 1) than that for (q 1 , q 2 ) = (1, 4). As a result, we find that the string network consists of the strings with (n 1 , n 2 ) = (1, 0) and (n 1 , n 2 ) = (0, 1) for (q 1 , q 2 ) = (4, 1) even at a later time.
Altogether, our lattice simulations in the 2 + 1 dimensional spacetime suggest that the string network is dominated by the compensated strings R c 1 for q 1 < q 2 and η 1 > η 2 , while it is dominated by the uncompensated one for q 1 > q 2 and η 1 > η 2 . We also find that R dw → 0 for q 1 < q 2 and η 1 > η 2 , We will discuss the behavior of the string network in the 3 + 1 dimensional spacetime in a separate paper, where the R c and R dw can be drastically different from those in the present simulations, since the reconnection rates of the strings can be drastically different in the 3 + 1 dimensional spacetime.
We have also performed the simulation for κ = 0.3. In this case, the mass term of φ 2 becomes positive at around φ 1 = 0 which stabilizes the symmetry enhancement point, φ 2 = 0 at the center of string core of φ 1 . This feature is expected to make the correlation between φ 1 and φ 2 stronger. We have numerically confirmed that the number of φ 2 strings formed around the φ 1 string tends to be larger for κ > 0 than the case with κ = 0, although we do not discuss details of those tendencies in this paper (see also [54]).

Interaction between strings
In this section, we discuss a long-range force between a compensated/uncompensated string with (n 1 , n 2 ) = (n 1 = 0, n 2A ) and a global string with (n 1 , n 2 ) = (0, n 2B ). We call the former to be the string A and the latter to be the string B. As our simulation was performed in 2 + 1 dimensional spacetime, we also consider the string configurations in the 2-dimensional space. In the 2-dimensional space, the spatial coordinates of the centers of the string A and B are given by r A and r B , respectively.

JHEP02(2020)058
In our discussion, we assume η 1 η 2 . 8 In this limit, the configurations of φ 1 and the gauge field are hardly affected by the configuration of φ 2 , and hence, we can treat them as the fixed background fields. Besides, we may neglect the φ 2 contribution to the Euler-Lagrange equation of the gauge field in eq. (2.25) as c 1 ( 1) c 2 . With such a gauge field configuration, D θ φ 1 vanishes at the large distance from the string A as in the case of the ANO string.
When the string A and B are well separated, the configuration of φ 2 can be approximated by, where φ 2A and φ 2B are the isolated string configurations of the string A and B. The longrange force between the string A and B can be read off from the difference of the string energy per unit length [15], where the energy functional is given by In this expression, we have taken r A = 0, and (r, θ) are the radius and the azimuthal angle around the string A. In this coordinate system, A θ is the gauge field configuration around the string. When the string A and B are well separated, ∆E is approximately given by Besides, by remembering that at the large distance, the contributions from the first term in eq. (4.4) is subdominant compared to those of the second term. As a result, we find Here, θ B is the azimuth angle around the string B, and its derivative with respect to θ is given by (4.8)

JHEP02(2020)058
Therefore, the energy difference is given where F (x) is a decreasing function of x.
The long-range force between the uncompensated and the global strings is given by the derivative of ∆E with respect to |r A − r B |. 9 Thus, we find the long-range force is attractive for while it is repulsive for With this force, the uncompensated string with n 1 = 1 and n 2 < 4 and the global string with n 1 = 0 and n 2 = 1 combine together rather quickly for (q 1 , q 2 ) = (1, 4). 10 The expression in eq. (4.9) also explains the weakness of the attractive force between the string with (n 1 , n 2 ) = (1, 0) and (n 1 , n 2 ) = (0, 1) for (q 1 , q 2 ) = (4, 1).

Conclusions and discussions
In this paper, we studied the formation and the evolution of the string network in the Abelian Higgs model with two complex scalar fields. As a special feature of the model, the model possesses a global U(1) symmetry in addition to the U(1) gauge symmetry. Both symmetries are spontaneously broken by the VEVs of the two complex scalar fields. In this model, the dynamics of the string network is quite rich compared with that in the ordinary Abelian Higgs model with a single complex scalar field. In particular, we found the existence of a new type of string solutions, the uncompensated strings, in addition to the ANO local string. The isolated uncompensated string has a logarithmic divergent energy per unit length as in the case of the global strings, although it is accompanied by non-trivial gauge field configuration.
We also performed the classical lattice simulations in the 2 + 1 dimensional spacetime, which confirmed the formation of the uncompensated strings at the phase transition. We also found that most of the uncompensated strings evolve into the compensated strings at later time when the gauge charge of the scalar field with the smaller VEV is larger than that of the scalar field with a larger VEV. Such a behavior can be understood by long range forces between the uncompensated string and the global string.
Finally, we comment on the implications of the results in the present paper to the axion models where U(1) global is identified with the Peccei-Quinn symmetry. As we have 9 To avoid the logarithmic divergence of F (x), we need to regularize the integration in eq. (4.7) at a large radius. The long-range force is, on the other hand, free from the divergence. 10 The long-range attractive force presents even for n2A = 0 where h2 has a non-trivial configuration due to the non-trivial A θ configuration.
JHEP02(2020)058 mentioned in subsection 3.3, we found that all the strings in our simulation for (q 1 , q 2 ) = (1, 4) have the effective winding number either of 0 or ±1 in the direction of the gaugeinvariant Goldstone boson at late time. When we apply this model to the axion model, at most one domain wall is attached to the cosmic string below the QCD scale [55]). Such string-wall network does not cause the infamous domain wall problem. Thus, the axion model with q 1 = 1 and q 2 ≥ 4 and η 1 η 2 might be free from the domain wall problem even when the U(1) local and U(1) global takes place after inflation. The result of the present paper is, however, based on the simulation in the 2 + 1 dimensional spacetime, and hence, not conclusive. We will perform the classical lattice simulation in the 3 + 1 dimensional spacetime in a separate paper.

A Domain wall problems in QCD axion models
Here, we briefly summarize the axion domain wall problem when we apply the accidental global symmetry, U(1) global , to the Peccei-Quinn symmetry [29][30][31][32] (see also [36][37][38]). The Peccei-Quinn mechanism is one of the most successful solutions of the Strong CP problem, where the effective θ-angle of QCD is canceled by the VEV of the pseudo-Nambu-Goldstone boson, the axion ( [55] for review).
To see how the axion appears in the present setup, let us consider the KSVZ axion model [56,57] and introduce N 1 and N 2 of the vector-like multiplets of the fundamental representation of the SU(3) c gauge group of QCD which couple to φ 1 and φ * 2 through, Here, (Q i ,Q i ) collectively represents the N i vector-like multiplets (i = 1, 2). We suppress the Yukawa coupling constants for notational simplicity. Below the U(1) local and U(1) global breaking scales, the vector-multiplets become massive. 11 Then, by integrating out the vector-like multiplets, the couplings of the Goldstone bosons in eq. (2.12) to QCD is given by, Here, g s is the QCD coupling constant, G andG are the gauge field strength of QCD and its dual. We suppress the color and the Lorentz indices. To satisfy the anomaly free condition of the U(1) local gauge symmetry, we require, N 1 q 1 − N 2 q 2 = 0 , (A.3) 11 We assume that η1,2 are much higher than the QCD scale.

JHEP02(2020)058
and hence, the effective Goldstone-QCD interactions are reduced to Here, we have introduced a parameter N m by N m = N 1 /q 2 = N 2 /q 1 (N m ∈ Z >0 ). Therefore, the gauge-invariant Goldstone can be identified with the QCD axion. Now, let us discuss how the axion evolves in the cosmological history. Let us assume that the U(1) symmetry breaking takes place after inflation. Then, the axion locally settles down to a point at the bottom of the scalar potential in eq. (2.6) after the phase transition. The axion, however, winds around the cosmic string with an effective winding number, n eff = n 1 q 2 − n 2 q 1 , (A. 5) in the domain of the axion [0, 2πF a ) as given in eq. (2.14). Below the QCD scale, the axion coupling to QCD in eq. (A.4) results in a periodic axion potential with a period, F a /N m . Therefore, the axion potential gives rise the energy contrasts around the cosmic strings and the cosmic strings are attached by N dw = N m × |n eff | domain walls. For |N dw | > 1, the string-wall network is stable and immediately dominates the energy density of the universe, which conflicts with the Standard Cosmology. 12 Around the compensated strings, i.e. n eff = 0, on the other hand, no walls are formed, and hence, it does not cause cosmological problems. Besides, the string-wall network with N dw = 1 also does not cause cosmological problems, since the network is not stable and decays immediately after the QCD phase transition [58,59]. The latter possibility requires N m = 1 and |n eff | = 1. Therefore, the domain wall problems in the axion model can be avoided if the cosmic strings evolves into either R c → 1 or R dw → 0 at a later time. 13 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.