Topological defects and nano-Hz gravitational waves in aligned axion models

We study the formation and evolution of topological defects in an aligned axion model with multiple Peccei-Quinn scalars, where the QCD axion is realized by a certain combination of the axions with decay constants much smaller than the conventional Peccei-Quinn breaking scale. When the underlying U(1) symmetries are spontaneously broken, the aligned structure in the axion field space exhibits itself as a complicated string-wall network in the real space. We find that the string-wall network likely survives until the QCD phase transition if the number of the Peccei-Quinn scalars is greater than two. The string-wall system collapses during the QCD phase transition, producing a significant amount of gravitational waves in the nano-Hz range at present. The typical decay constant is constrained to be below O(100) TeV by the pulsar timing observations, and the constraint will be improved by a factor of 2 in the future SKA observations.


Introduction
The Peccei-Quinn (PQ) mechanism is one of the most attractive solutions to the strong CP problem [1,2]. In association with spontaneous breakdown of a global U(1) PQ symmetry, there arises a pseudo Nambu-Goldstone boson called the QCD axion [3,4]. In the conventional scenarios, the axion decay constant F a is of order the U(1) PQ breaking scale. The classical axion window is given by 10 9 GeV F a 10 12 GeV, (1.1) where the lower bound is due to the SN 1987A neutrino burst duration [5][6][7] and the upper bound is due to the cosmological abundance of the axion produced by the misalignment mechanism, barring fine-tuning of the initial misalignment [8][9][10]. The QCD axion is stable in a cosmological time scale, and is a plausible candidate for cold dark matter. See e.g. refs. [11][12][13][14][15][16] for reviews on the QCD axion and the related topics.
Axions are also known to appear in string theory in association with compactification of the extra dimension, and some of them may be so light that they play an important role in cosmology [17][18][19][20]. In particular, a multi-axion inflation model attracted much attention, because multiple axions with sub-Planckian decay constants can conspire to realize a light inflaton with the effective super-Planckian decay constant through the alignment mechanism [21][22][23][24][25][26][27]. A peculiar structure of the charge assignment for the N axions was noted in ref. [28], and a concrete realization was given in refs. [29,30]. If there are many axions with various kinds of shift symmetry breaking terms, they may form an axion landscape [23,24], where eternal inflation and the subsequent quantum tunneling or jump-up to one of the adjacent minima take place repeatedly. In the axion landscape, the slow-roll inflation may be realized via the alignment mechanism, and the inflaton potential is generally JHEP08(2016)044 expected to have small modulations, which may lead to a sizable running of the spectral index [31,32]. The vacuum structure of the axion landscape was studied in refs. [33,34].
Recently, four of the present authors (TH, KSJ, NK, and FT) applied the alignment mechanism to the QCD axion, and studied its cosmological and phenomenological implications [35,36]. In this model, the QCD axion is given by a flat direction composed of multiple axions. Interestingly, even if the decay constant of each axion is around the TeV scale, the QCD axion decay constant can be enhanced up to the classical axion window (1.1). One of the axions or saxions in the model can explain the recently found diphoton excess in LHC [37][38][39][40]. 1 In addition, the required high quality of the PQ symmetry [44,45] can be naturally realized in this model [36].
Cosmological evolution of the QCD axion can be quite involved even in the single axion case. Suppose that the U(1) PQ symmetry is restored in the early Universe. This is the case if the inflation scale or the reheating temperature is sufficiently high. Then cosmic strings are formed when the U(1) PQ is spontaneously broken. Subsequently, the axion potential arises due to the QCD instanton effects, leading to the formation of domain walls bounded by strings. If the domain wall number of the QCD axion is unity, such a string-wall network soon collapses by emitting axions [46,47]. On the other hand, if the domain wall number is greater than unity, the string-wall network is stable, which causes the cosmological disaster. 2 In the aligned axion model with N axions, the associated U(1) N symmetries are likely restored in the early Universe, because the symmetry breaking scale is much smaller than F a . The spontaneous breakdown of the U(1) N symmetries would result in the formation of a complicated string-wall network, which reflects the aligned structure of the multiple axion fields. The main purpose of this paper is to study the structure and evolution of the topological defects in the aligned axion model. As we shall see, in a case with two axions, an isolated string-wall system can be formed, which eventually collapses to a bundle of strings that can be regarded as the QCD axion string. In general, such a bundle of strings glued by walls exists as a solution corresponding to the QCD axion string, and in fact, their tensions are approximately same. However, for N ≥ 3, such a structure is rarely formed in the early Universe, because the solution requires an exponentially large hierarchy in the number of strings of different types. Instead, we expect that the string-wall network survives until a later time, and collapses during the QCD phase transition, emitting gravitational waves. We derive an upper bound on the PQ breaking scale by the pulsar timing experiments [50][51][52][53].
The rest of this paper is organized as follows. In section 2, we review the aligned QCD axion scenario. We discuss the formation and evolution of topological defects in the aligned QCD axion model in section 3. Section 4 is devoted to discussion and conclusions. 1 See refs. [41][42][43] for models of a visible axion as an explanation of the diphoton excess. 2 There is an alternative way to form domain walls in a multi-axion scenario, where the system exhibits a chaotic behavior (axion roulette) during a mass level crossing between two axions [48,49].

JHEP08(2016)044 2 Aligned QCD axion
To enhance the decay constant by the alignment mechanism, one needs multiple axions. Let us introduce N axions, φ 1 , φ 2 , . . . , φ N , each of which respects the shift symmetry, where C i is a real transformation parameter. We assume that N −1 of the shift symmetries are explicitly broken down to their discrete subgroups, giving rise to the potential for (N −1) combinations of the axions. The remaining massless axion is to be identified with the QCD axion, and its anomalous coupling to gluons is generated by including extra PQ quarks as in the Kim-Shifman-Vainshtein-Zakharov axion models [54,55] . The QCD axion is given by a certain combination of the N axions, and it remains massless until the QCD instanton effects are turned on. In this set-up, the QCD axion decay constant depends on how the N − 1 shift symmetries are broken [56], and it can be significantly enhanced if the alignment takes place.
Here we adopt the clockwork axion model [29,30] as a simple and concrete realization of the aligned axion model, and we apply it to the QCD axion following ref. [35]. The following discussion, however, can be straightforwardly applied to more general aligned axion models. In the clockwork axion model, the potential is given by for integer n i , where Λ i and f i are a characteristic energy scale of the shift symmetry breaking and the decay constant of φ i , respectively. Since there are N − 1 shift symmetry breaking terms, there remains a flat direction which is to be identified with the QCD axion. The QCD axion and its decay constant are given by where we have defined n N = 1 for notational convenience. Suppose that all the decay constants f i and the integers n i are comparable to each other, i.e. f i f and n i n. Then, the largest contribution to the QCD axion comes from φ 1 , and the QCD axion decay constant is exponentially enhanced, F a ∼ e N ln n f , for large N . For instance, if f i f = O(1) TeV, one needs an enhancement of order 10 6−9 for F a to be in the classical axion window (1.1), which can be realized for e.g. n = 3 and N = 13 − 19.
One of the virtues of the aligned QCD axion model is that the high quality of the PQ symmetry is naturally explained [35,36]. This is because the actual symmetry breaking scale is much smaller than F a , and so, any Planck-suppressed PQ breaking terms give only JHEP08(2016)044 suppressed contributions to the QCD axion potential. In some case, such Planck-suppressed PQ breaking terms have interesting implications for the QCD axion dynamics [36]. 3 Another cosmological implication is that the PQ symmetry is easily restored in the early Universe. To see this, let us now consider a possible UV completion, where each axion φ i is originally embedded in a phase component of a complex scalar field, Φ i . The model is based on global U(1) N symmetries and the scalar potential has a typical form of In addition, the scalars generically acquire a Hubble-induced mass or thermal mass in the early Universe. Then, if the Hubble parameter or the temperature is higher than m i , the scalars are stabilized at the origin, and the U(1) N symmetries are restored. As the Universe expands, Φ i becomes tachyonic at the origin, and develops a vacuum expectation value, Then, N massless axions φ 1,··· ,N appear as Nambu-Goldstone bosons. At the same time, cosmic strings are produced in association with a non-trivial topological configuration of φ i . Later on, the (N − 1) symmetry breaking terms become important and the (N − 1) axions acquire masses. As an example, we may adopt the following renormalizable potential [30], where i is an order parameter; i 1 implies that the corresponding U(1) symmetry is a relatively good symmetry. This potential generates the axion potential (2.2) with n i = 3. As a result, domain walls appear stretching between the cosmic strings. The strings and walls form a complicated string-wall network. We will study the nature of these topological defects in the next section.
Lastly let us comment on the axion coupling to gluons. In order to couple the QCD axion a to gluons with a coupling suppressed by 1/F a , we introduce the following interaction of Φ N with heavy PQ quarks, taking into account that the φ i -fraction in a quickly decreases with i for f i f . Then, integrating out the PQ quarks, heavy axions and saxions, we are left with the usual QCD axion coupled to gluons. 4

JHEP08(2016)044 3 Topological defects and gravitational waves
Here we study the formation and evolution of cosmic strings and domain walls in the aligned axion model. First we discuss the structure of string-wall network and its correspondence with the ordinary QCD axion string. Then we investigate their evolution based on numerical simulations. Lastly, we evaluate the gravitational waves emitted in the process of the domain wall collisions during the QCD phase transition.

Cosmic strings and domain walls
In the aligned QCD axion model, the actual symmetry breaking scales, f i , are much smaller than F a , and so, the U(1) N symmetries are easily restored in the early Universe. When the U(1) N symmetries are spontaneously broken, the radial component of each complex scalar field develops a nonzero vacuum expectation value, while the phase component is randomly distributed in space. There appear N kinds of cosmic strings corresponding to a non-trivial topological configuration of φ 1,··· ,N . Consider a string of Φ i along the z-axis. Then, one can express the complex scalar field around the string in a cylindrical coordinate system as where θ is the angular coordinate, w i is the winding number and we have neglected the radial component of Φ i . Most of cosmic strings realized in the spontaneous symmetry breaking have w i = ±1, and so, we will focus on this case, and denote the (anti-)string with w i = 1(−1) by S i (S i ) in the following. The energy per unit length, i.e., the sting tension, for an infinitely long global string is given by the sum of the potential energy stored inside the string core and the gradient energy outside the string. The former is roughly µ core,i ∼ f 2 i , and the latter is logarithmically divergent and gives the dominant contribution. For each cosmic string, the tension is estimated as is the typical core radius, and R is the cutoff length corresponding to the distance between strings, which is usually the Hubble radius. Those N kinds of strings individually follow the scaling law until the domain walls are formed.
Later on, the shift symmetry breaking terms such as (2.6) become important, and the N − 1 U(1) symmetries are explicitly broken down to their discrete subgroups, generating discrete minima for the (N − 1) axions. As a result, domain walls appear between strings. There are two kinds of domain walls, W (i,i+1) and W (i,i) , where the former stretches between S i and S i+1 with i = 1, · · · , N − 1, and the latter between the anti-stringS i and string S i . For the potential (2.2), a single wall is attached to S 1 , n i−1 + 1 walls are attached to S i with i = 2, · · · , N − 2, and n N −1 walls are attached to S N . The composition of the walls depends on the initial string configuration, and it also evolves with time as strings and walls annihilate. Specifically, S i with i = 2, · · · , N − 2 has n i−1 + 1 walls whose composition is expressed by

JHEP08(2016)044
with a i = 0, 1 and b i = 0, · · · , n i . Those strings and walls form a complicated string-wall system, whose evolution is numerically studied later in this section.

Cosmic string bundles as QCD axion strings
Here we show that there is a special configuration in which many strings S 1,··· ,N glued by domain walls form an isolated string bundle. The effective tension of the bundle is of the same order of the QCD axion string in a usual axion model without the alignment. Let us begin with the case of N = 2 and adopt the renormalizable potential (2.6). In this case, one possible configuration is such that threeS 1 strings are connected to an S 2 string by the walls W 12 , as schematically illustrated in figure 1a. Once such a configuration is formed, each string gets attracted to each other by the domain wall tension, leading to a single string bundle composed of S 2 + 3S 1 glued by the walls W 12 . Taking account of the fact that the effective winding number w 1 is equal to −3, the effective tension of such a string bundle is where we have used the effective decay constant, F a = 3 2 f 2 1 + f 2 2 (see the relation (2.3)), in the second equality, and we have neglected subdominant contributions from the string cores and the walls. The tension is consistent with that of the QCD axion string with the winding number equal to unity. The above argument can be straightforwardly applied to the case of N ≥ 3. An isolated string bundle is composed of S N + 3S N −1 + 3 2 S N −2 + · · · + 3 N −1 S 1 (S 1 ) for odd (even) N . See figure 1b in the case of N = 3. The effective tension is in general Thus, the string bundle can be regarded as the usual QCD axion string with the decay constant F a . Let us note here that the aligned structure in the axion field space appears in the real space as the string bundle. The string bundle contains exponentially many strings, and as a result, the effective winding number of each string becomes exponentially large. This enhances the string tension by a factor of F 2 a /f 2 . To reinforce the estimation above, here we numerically analyze the structure of an isolated string bundle for N = 2. For this purpose, we have computed a static solution of Φ i that minimizes the action with a potential consisting of eqs. (2.3) and (2.6). The action is computed on a 2-dim lattice with the number of grids N grid = 512 2 . The box size is taken to be about 70 times larger than the typical length scale of strings ( To obtain an isolated string bundle, we have adopted a Dirichlet boundary condition, so that along the boundary eq. (3.1) is satisfied. The minimization of nonlinear actions with such a large N grid is not trivial. In this paper we have adopted the scheduled relaxation Jacobi method [58,59]. Figure 2 shows the configurations of Φ 1 and Φ 2 . Here we have assumed f 1 = f 2 = f and = 0.25. First of all, it can be clearly seen that there are three distinct S 1 strings, which surround theS 2 string at the center. This is consistent with the schematic picture JHEP08(2016)044  in figure 1a. We also note that the contours indicating |Φ 2 | are elongated towards the cores of the S 1 strings. This is because domain walls stretch between S 1 andS 2 strings. For the potential energy in ∆V of eq. (2.6) to be reduced, |Φ 2 | is required to be nonzero along the domain walls. We note that in the right panel in figure 2, the discontinuity in the color scale corresponding to θ 2 = 0 appears wiggling, which results from minimizing the potential energy around domain walls. While it is not as apparent as θ 2 = 0, loci corresponding to θ 2 = 2/3π and θ 2 = 4/3π are also wiggling in a similar way. From our numerical calculation, we can also compute the tension of the string bundle. In figure 3, we plot the energy stored inside the radius from the center of the bundle. Because the gradients of the phases θ i dominate at large radii, the energy approximately logarithmically increases and its asymptotic behavior converges into eq. (3.2).

Formation and evolution of string-wall network
Now the question is if isolated string bundles are produced after the domain walls appear in the early Universe. As we shall argue below, such isolated string bundles are unlikely to be formed if N ≥ 3.
Starting from the random initial condition, cosmic strings are necessarily produced after the spontaneous breakdown of the N U(1) symmetries. The cosmic strings are considered to follow the scaling law in a few Hubble time, and the number of each string per the Hubble horizon is of order unity. At a later time, the shift symmetry breaking terms generate discrete minima for the (N −1) axions, and domain walls are formed between two strings.
Let us take one S N string, and consider a whole object connected to it by various kinds of walls and strings. Assume that such an object contains a finite number of strings. As we start from the S N string, each string contained in this object is identified with either S i or S i with i = 1, · · · , N . Namely, the orientation of the strings is determined by how they are connected to the S N string. The number of strings is monotonically decreasing with time JHEP08(2016)044   such an object will collapse into an isolated string bundle at a sufficiently late time when the whole object is contained in a Hubble horizon. Since there is no essential difference between S N andS N and their number density is comparable, the above condition (3.6) is easily satisfied if such an object with a finite number of strings is formed. However, the object must have a larger bias in the number of the strings of the other types. This can be understood by noting that the object must contain much more S 1 (S 1 ) strings thanS N (S N ) strings for odd (even) N , as it collapses into S N + 3S N −1 + 3 2 S N −2 + · · · + 3 N −1 S 1 (S 1 ) at the end of the day. That is to say, the number of strings in the object must satisfy for 1 ≤ i ≤ N . However, starting with the random initial condition, the number of strings and anti-strings in a part of the whole string-wall network typically satisfy Therefore it is extremely unlikely that the required huge bias in the number of strings is realized. We thus conclude that such an object contains most likely an infinite number of strings for large N , and no isolated objects are produced. 5 Indeed, as we shall see below, we could not find any isolated structure in the numerical simulations with the N = 3 case (see figure 7). On the other hand, in the case of N = 2, the required bias is not huge, and isolated string bundles are formed (see figure 6). Even if isolated string bundles are not formed, the string-wall network will never disappear until the QCD instanton effects are turned on. This can be understood as follows. Note that only N − 1 of the original N U(1) symmetries are explicitly broken, and there remains a massless degree of freedom corresponding to the QCD axion, a. Before the domain wall formation, all possible field values of the N axions are realized in space, and there is no bias in their distribution. However, if the string-wall network disappeared, it would imply that the distribution of the QCD axion is topologically trivial everywhere. This is the case if only a part of its field space is realized, which is inconsistent with the initial condition. Indeed, the QCD axion has a non-trivial topological configuration around the S N string, which would be the core of isolated string bundle. In other words, the S N strings cannot completely disappear in the Universe until the QCD phase transition, irrespective of whether isolated string bundles with S N being the core are formed or a complicated string-wall network remains around it.
In summary, it is highly unlikely that isolated string bundles are formed in the Universe for large N , as long as the initial condition is random. The string-wall network in the aligned axion model likely remains until the QCD phase transition. Figure 4. Two-dimensional lattice simulation for N = 2. Time evolves from left to right. The red (blue) square and triangle points represent S 1 (S 1 ) and S 2 (S 2 ) respectively and the green line represents W 12 .

Numerical simulation of the string-wall network
We have performed two and three-dimensional lattice simulations for the string-wall network in the aligned axion model given by eqs. In the case of N = 2 (figures 4 and 6), S 1 orS 1 becomes the boundary of domain walls W 12 . One can see from the figures that such domain walls tend to shrink and eventually (almost) disappear, leading to the cosmic string bundle, S 2 + 3S 1 orS 2 + 3S 1 , as discussed before. Thus, isolated string bundles are formed in the case of N = 2 and n = 3.
In the case of N = 3 (figures 5 and 7), the domain wall network survives and the isolated structure cannot be seen. In the two-dimensional case, the walls wrap around the lattice simulation box, and both W 12 and W 23 remain. On the other hand, in the threedimensional case, while incomplete string bundles, S 2 + 3S 1 orS 2 + 3S 1 , are formed, there also remain S 3 andS 3 strings. Importantly, their numbers are comparable to each other as expected, which makes it difficult to form isolated string bundles. We could not find any isolated string bundle and the string-wall system remains in the case of N = n = 3.
The string-wall network likely survives until the QCD phase transition and after that it disappears due to the QCD axion potential which serves as an energy bias to break the degeneracy of discrete vacua. During the domain wall collisions, heavy axions are emitted and they will decay into gluons. Depending on the strength of the coupling to gluons and the masses, some of them may be so long-lived that the energetic gluons change the light element abundances. In addition to heavy axions, the QCD axion is produced by the domain walls and strings. Their abundance may be different from the ordinary QCD axion model, but it requires dedicated numerical simulations. The detailed study of those axions produced by the collapse of the string-wall network is left for future work.   JHEP08(2016)044

Gravitational waves from domain wall annihilation
We have seen that the domain walls are long-lived for N ≥ 3 and most probably survive until the QCD axion potential arises. Assuming the scaling regime, the energy density of domain walls evolves as ρ wall ∼ σH. Here σ is the tension of the domain wall, and it is approximately given by where we have set f i = f and i = for simplicity, and m aH is the mass of heavy axions that form the domain walls. Note that those heavy axions are orthogonal to the QCD axion. Since the energy density of domain walls decreases more slowly than radiation or matter, stable domain walls will eventually dominate the Universe and make it unacceptably inhomogeneous, spoiling the standard cosmology. Thus, domain walls must annihilate before they dominate the Universe. In our aligned QCD axion model, the domain walls annihilate after the QCD axion potential arises from the QCD instanton effects at T ∼ 1 GeV if the domain wall number of the QCD axion is equal to unity. This is because the QCD axion potential behaves like a bias V ∼ Λ 4 QCD for the W N −1,N domain walls and they annihilate at σH ∼ Λ 4 QCD [60]. Remaining domain walls W N −2,N −1 , . . . , W 1,2 also annihilate one after another with attached cosmic strings. If σH < Λ 4 QCD is satisfied at T ∼ 1 GeV, the domain walls annihilate soon after the QCD axion potential is turned on. Otherwise, the annihilation temperature is estimated as which implies that the annihilation is delayed significantly for f 100 TeV 1/6 . We have numerically confirmed that the string-wall system collapses after the QCD instanton effects become relevant. Requiring that the domain walls should not dominate the Universe until T = T ann , one obtains a conservative constraint on f , 12) and the annihilation temperature must be T ann 0.1 GeV. 6 In fact, a more stringent constraint comes from the pulsar timing observation, as we shall see below.
In the violent collisions of domain walls, gravitational waves are produced over frequencies ν corresponding to a typical physical length scale [61][62][63][64]. Once the domain-wall network follows the scaling law, a typical curvature radius is of order the Hubble parameter. 6 In the presence of extra PQ breaking terms, each axion is trapped at each potential minima even before the QCD phase transition [36]. In this case, the string-wall network sooner disappears and f can be a larger value.

JHEP08(2016)044
Then, the power spectrum of the gravitational waves has a peak at frequency corresponding to the Hubble parameter at the domain wall annihilation, ν peak (t ann ) H ann . As the Universe expands, the peak frequency is red-shifted, and its present value is given by ν peak,0 1.6 × 10 −7 Hz g * ann 80 which happens to be in the sensitivity range of the pulsar timing experiments. The intensity of the gravitational waves is usually characterized by a dimensionless quantity, Ω gw (t), defined by where ρ gw (t) is the energy density of gravitational waves, ρ c (t) the critical density, and ν is the frequency. The peak value of Ω gw at the annihilation is estimated as where˜ gw 0.7 ± 0.4 [64] is an efficiency parameter of the gravitational wave emission, and A parametrizes the energy density of domain walls as ρ wall = Aσ/t in the radiationdominated era. In the case of a simple Z 2 potential, A 0.8 ± 0.1, and it is considered to increase in proportion to the number of the potential minima [46,47]. In our aligned axion model, the vacuum structure is much more complicated than the case of Z 2 , and we expect A = O(1 − 10). In fact, we have obtained A ∼ N from two dimensional lattice simulations.

JHEP08(2016)044
As mentioned before, it may take a longer time for the string-wall system to annihilate completely in the aligned axion model, compared to the usual case of a single PQ scalar. As a result, the gravitational power spectrum may be modified, and, in particular, the peak frequency may be slightly lowered. Taking account of this uncertainty, the upper bound on f may be tightened by a factor of several. Detailed lattice numerical simulations for the produced gravitaional waves are warranted. The future observation by SKA will reach Ω gw h 2 ∼ 10 −13 [65] and the upper bound on f will be improved by a factor of 2. It implies that the gravitational waves produced in the aligned axion model with a decay constant of O(10 − 100) TeV will be probed by the pulsar timing experiment in future.

Discussion and conclusions
Let us discuss the case in which the U(1) N symmetries are only partially restored by some of the scalar fields which develop VEVs after inflation. Suppose that Φ N develops a nonzero VEV during inflation. In this case, a complicated string-wall network does not remain for a long time. This is because, for the cosmic string to remain in the end, the phase of Φ N should rotate from 0 to 2π about the string, but if Φ N has a non-zero VEV already, such configuration does not appear. It seems that only the cosmic strings associated with the scalars which develop non-zero VEVs after inflation will appear for the moment, but they are attached to domain walls and so, they will disappear as soon as the domain walls appear. As for the isocurvature fluctuations, only some fraction of the field space is populated by this partial spontaneous symmetry breaking. Therefore, after the string-wall network disappears, some amount of isocurvature perturbations may be left. But it is difficult to know what actually happens unless the evolution of the string-wall network is numerically studied. In any case, complete symmetry restoration seems more plausible, in which case no isocurvature perturbations are generated.
In this paper we have studied the formation and evolution of the complicated stringwall network in the aligned axion model. Since the actual PQ symmetry breaking scale is much smaller than the classical axion window, the symmetry restoration is more plausible in this scenario. Such axion model has several virtues: the high quality of the PQ symmetry is naturally explained, and no isocurvature perturbations are generated if the symmetry is restored during or after inflation. We have shown that there exists a solution of the isolated string-wall system, which can be identified as the QCD axion string as their string tensions are approximately equal to each other. However, the formation probability of such isolated string bundles is extremely suppressed for N ≥ 3, because it requires a highly biased population of each type of strings (S i ) over anti-strings (S i ), which cannot be realized if the strings follow the scaling law. Thus, the string-wall network is considered to be extremely long-lived, and most probably survives until the QCD axion potential appears. During the violent domain wall collisions, a significant amount of gravitational waves are emitted. We have derived an upper bound on the axion decay constant, f i O(100) TeV, by the pulsar timing experiments, and the bound will be improved by a factor of 2 in the future observations by SKA.