{'fraction_non_alphanumeric': 0.04210701045758867, 'fraction_numerical': 0.015617219056050463, 'mean_word_length': 4.055454545454546, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Axion cosmology with long-lived domain walls

We investigate the cosmological constraints on axion models where the domain wall number is greater than one. In these models, multiple domain walls attached to strings are formed, and they survive for a long time. Their annihilation occurs due to the effects of explicit symmetry breaking term which might be raised by Planck-scale physics. We perform three-dimensional lattice simulations and compute the spectra of axions and gravitational waves produced by long-lived domain walls. Using the numerical results, we estimated relic density of axions and gravitational waves. We find that the existence of long-lived domain walls leads to the overproduction of cold dark matter axions, while the density of gravitational waves is too small to observe at the present time. Combining the results with other observational constraints, we find that the whole parameter region of models are excluded unless an unacceptable fine-tuning exists.


Introduction
Axion [1] is a hypothetical particle which arises as a consequence of the Peccei-Quinn (PQ) mechanism [2], the most attractive solution to the strong CP problem of quantum chromodynamics (QCD) [3,4]. It arises as a Goldstone boson when a global U (1) PQ symmetry (so called PQ symmetry) is spontaneously broken at some high-energy scale denoted as F a . The phenomenological searches indicate that axions are "invisible" due to the weak coupling with matter which is suppressed by 1/F a . This invisibleness is realized by two kinds of models.
One is the KSVZ model [5], which is realized by introducing additional quarks to the standard model. Another is the DFSZ model [6], which is realized without introducing additional fermions by assuming two Higgs doublets. These invisible axions lead to rich implications to astrophysics and cosmology, such that they are good candidates of dark matter filled in the universe [7]. The remarkable feature of the axion models is that they predict the formation of topological defects in the early universe [see [8] for reviews]. When the PQ symmetry is spontaneously broken, vortex-like defects, called strings, are formed. These strings are attached by surface-like defects, called domain walls, when the temperature of the universe becomes comparable to the QCD scale. The structure of the domain walls is determined by an integer number N DW which is referred as the "domain wall number". The value of N DW is related to the U (1) PQ -SU (3) c -SU (3) c anomaly which is given by the formula [4,9,10] where Q PQ (q i ) is the U (1) PQ charge for quark species q i , i=L(R) represents the sum over the left(right)-handed fermions, and T a are the generators of SU (3) c normalized such that TrT a T b = Iδ ab where I = 1/2 for fundamental representation of SU (3) c . Here we choose the unit Q PQ = 1 for minimum non-vanishing magnitude of PQ quantum numbers. In the KSVZ models we obtain N DW = 1, while in the DFSZ models it is predicted that N DW = 2N g where N g is the number of generations. As we discuss in Sec. 2, the cosmological scenario is different between models with N DW = 1 and N DW > 1. Recently, we studied the cosmological implication of the models with N DW = 1 through the careful numerical simulations [11]. In this paper, we turn our research into the models with N DW > 1.
For the case N DW > 1, it is well known that domain walls are stable and they eventually overclose the universe, which conflicts with the standard cosmology [9,12]. So far, several solutions to this problem are proposed. One is to introduce a small symmetry breaking parameter, called the bias [13], in the Lagrangian, which makes the walls unstable and leads to their annihilation at late times. Another solution is to embed the discrete subgroup of U (1) PQ into another continuous group [14]. In such models, the degenerate vacua are connected by symmetry transformations, leading to the same cosmological property with the case for N DW = 1. It is also possible to avoid the domain wall problem by assuming that inflation has occurred after the PQ symmetry breaking. In this case, the inflationary scale is constrained by the observation of isocurvature fluctuations in the cosmic microwave background (CMB) [15]. Among these possibilities, here we focus on the first solution, the biased domain walls.
It was pointed out that the bias term, which causes the annihilation of domain walls, can be found in effects of gravity [16]. If we regard the low energy field theory as an effective theory induced by a Planck-scale physics, we expect that the global symmetry is violated by higher dimensional operators suppressed by the Planck mass M P . This small symmetry violating operators cause the decay of domain walls. However, in the case of PQ symmetry, the situation is more complicated. It was argued that these Planck-scale induced terms easily violate the PQ symmetry and cause the CP violation [17][18][19]. In particular, the operators with dimension smaller than d = 10 are forbidden [18] by a requirement of CP conservation. If it is true, the biased domain walls should be long-lived, and disappear at late time due to the effect of highly suppressed operators.
In our previous numerical study [20], we found that the late-time annihilation of domain walls actually occurs. We also pointed out that the gravitational waves produced by domain walls may be detectable in the future experimental searches. However, the previous results strongly depended on the production rate of axions and gravitational waves, which was treated as free parameters. In this work, extending the previous study [20], we aim to determine the amount of radiated axions and gravitational waves precisely. This requires the estimation of the spectra of radiated particles from domain walls, which can be investigated by performing field-theoretic numerical simulation of scalar fields.
In this paper, we study the formation and evolution of string-wall networks in the models with N DW > 1 and calculate the spectra of axions and gravitational waves radiated from them. This work might be regarded as the first study in which the particle radiations from long-lived domain wall networks are estimated quantitatively. Although we focus on the axionic models, part of our findings might be applicable to other particle physics models which predict the existence of long-lived domain walls, such as domain walls formed by discrete Rsymmetry breaking [21], thermal inflationary models [22], and Z 3 symmetry breaking in NMSSM [23].
The organization of this paper is as follows. We introduce the model and review the property of topological defects in Sec. 2. In Sec. 3, we show the results of numerical simulations. Using the numerical results, we calculate the relic density of axions and gravitational waves in Sec. 4. We show that the overabundance of axions produced by domain walls gives severe constraints on the model parameters. We also review other observational constraints. Finally, we make a summary in Sec. 5. The analysis method that we used in the numerical studies are described in Appendices A, B, C, and D. We make a brief comment on the determination of the surface mass density of domain walls in Appendix E.
In this paper, we work in spatially flat Friedmann-Robertson-Walker (FRW) universe with a metric where R(t) is the scale factor of the universe. We denote the cosmic time as t and the conformal time as τ , where dτ = dt/R(t). A dot represents a derivative with respect to the cosmic time, while a prime represents a derivative with respect to the conformal time, i.e. = ∂/∂t, and ′ = ∂/∂τ .

Model of domain walls bounded by strings
We consider the model of a complex scalar field Φ with the Lagrangian density given by where the scalar potential is given by Here, θ is the angular part of the scalar field (i.e. Φ = |Φ|e iθ ). The first term in Eq. (2.2) is the usual mexican-hat potential, which causes spontaneous breaking of U (1) PQ symmetry. In the early universe, the U (1) PQ symmetry is broken when the temperature of the universe falls below the energy scale ∼ η. At this time, the scalar field gets a vacuum expectation value | Φ | 2 = η 2 , and cosmic strings are formed. Then, the angular part of the scalar field can be identified as axion field a, Φ = ηe ia/η . The second term in Eq. (2.2) may arise due to the non-perturbative effect of QCD. Because of the presence of this term, the U (1) PQ symmetry is explicitly broken into its subgroup Z N DW . This term becomes non-negligible when the mass of the axion m becomes greater than the friction term which comes from the cosmic expansion, m H, where H is the Hubble parameter. Then, the axion field begins to roll down to one of N DW degenerate minima represented by At this stage, domain walls are formed, separating the N DW vacua [9]. These domain walls are attached to cosmic strings around their boundaries [8].
We note that there are two different energy scales in this model. One is the scale of the PQ symmetry breaking, given by the axion decay constant Phenomenologically, it takes a value F a ≃ 10 9 -10 12 GeV (see e.g. [3]). Another is the QCD scale, Λ QCD ≃ O(100)MeV, which induces the second term of the effective potential (2.2). The mass of the axion is given by m ∼ O(0.1)Λ 2 QCD /F a . Since there is a large hierarchy between these two scales F a ≫ Λ QCD , the height of the potential barrier of the first term ∼ η 4 in Eq. (2.2) is much larger than that of the second term ∼ Λ 4 QCD . Because of the existence of this large hierarchy, the formation of strings occurs much earlier than the formation of domain walls. Once strings are formed, they evolve into the "scaling" regime, in which the population of strings in the Hubble volume tends to remain in O(1). In order to keep this scaling property, the strings lose their energy by emitting closed loops, and these loops decay by radiating axions [24][25][26][27][28][29]. Subsequently, domain walls are formed when the Hubble parameter becomes comparable to the mass of the axion H ∼ m. Shortly after that, the condition is satisfied at the time t w , where σ wall = 9.23mF 2 a is the surface mass density of domain walls [see Eq. (E.7)], µ string (t) ≃ πη 2 ln(t/δ s ) is the mass energy of the strings per unit length, and δ s ≃ 1/ √ λη is the width of strings. After the time t w , the dynamics is dominated by the tension of domain walls.
The fate of the string-wall networks is different between the case with N DW = 1 and the case with N DW > 1. If N DW = 1, domain walls bounded by strings are unstable, and they disappear immediately after the formation [30]. On the other hand, if N DW > 1, the stringwall networks are stable, and they eventually dominate the energy density of the universe. If it occurs at sufficiently late time in the universe, it may conflict with standard cosmology [12].

Domain wall problem and its solution
Now, let us look closer at the model with N DW > 1. After the time t w given by Eq. (2.5), domain walls are straightened by their tension force up to the horizon scale. In the past numerical studies [31], it was found that these networks of domain walls bounded by strings evolve into the scaling regime, which is similar to the property of the networks of strings.
In this regime the typical length scales of defects, such as the wall curvature radius and the distance of two neighboring walls, are given by the Hubble radius H −1 = 2t. Then, the energy density of domain walls behaves as Since the dilution of the energy density ρ wall is slower than that of matters ∼ 1/R 3 and radiations ∼ 1/R 4 , domain walls eventually dominate the energy density of the universe. This occurs at the time where G is the Newton's gravitational constant. In other words, the wall domination occurs when the temperature of the universe becomes less than T ∼ O(10)keV [20].
In order to avoid this overclosure problem, Sikivie [9] phenomenologically introduced a term δV = −Ξη 3 (Φe −iδ + h.c.) (2.8) in the potential (2.2), where Ξ is a dimensionless parameter which is assumed to be much less than unity 1 . This term, called a bias [13], explicitly breaks the discrete symmetry and lifts degenerate vacua represented by Eq. (2.3). If this kind of term exists in the potential, it affects as a pressure on walls, which eventually annihilate them. We can naively estimate the lifetime of the networks t dec by the ratio between the surface mass density of domain walls σ wall and the difference of the potential energy between two neighboring vacua [20], Requiring that the decay of walls occurs before the wall domination (i.e. t dec < t WD ), we obtain a lower bound on Ξ (2.10) In our previous study [20], we found that string-wall networks actually disappears in the time scale given by Eq. (2.9) and pointed out that significant amounts of gravitational waves might be produced by these defect networks. We also found that the cosmological abundance of cold axions produced by string-wall networks gives another constraint on the model parameters. In this work, we are going to give a concrete analysis of gravitational waves and axion radiations from these string-wall networks by numerically computing the spectrum of these radiations.

Numerical simulations
We performed field-theoretic lattice simulations with the model of PQ scalar field Φ whose potential given by Eq. (A.1). We use the three dimensional homogeneous grid with the size N 3 = 512 3 . In the numerical studies, we normalize the dimensionful quantities in the unit of η. For example, x → xη, Φ → Φ/η, etc. We choose the comoving size of the simulation box as L = 80 in this unit. Then, the spatial resolution is estimated as ∆x = L/N ≃ 0.16. We give the time evolution in terms of the conformal time τ rather than the cosmic time t. We set the initial time as τ i = 2.0 and the final time as τ f = 40.0. The dynamical range of the simulation is therefore τ f /τ i = 20. We normalize the scale factor at the initial time R(τ i ) = 1. The other parameters are summarized in table 1. The technical details on the analysis method are described in Appendices A, B, C, and D. With the parameters shown in table 1, the physical scale of the Hubble radius H −1 , the width of the wall δ w = m −1 , and the width of the string δ s = λ −1/2 divided by the lattice spacing ∆x phys = R(t)∆x = R(t)L/N are given by At the end of the simulation τ f = 40, these ratios become H −1 /∆x phys ≃ 256 < N , δ w /∆x phys ≃ 3.2, and δ s /∆x phys ≃ 1.01. Therefore, these length scales are marginally resolvable at the final time of the simulation. We naively guess that domain walls formed in the model with N DW = 2 have a common property with that formed in the model of real scalar field with spontaneous breaking of a discrete Z 2 symmetry. Based on this observation, we also consider the real scalar field theory with the Lagrangian density given by where φ is the real scalar field. In Appendix E, we show that this model mimics the axionic model with N DW = 2, in the sense that the surface mass density and the width of domain walls are same. It should be emphasized that we do not perform the simulation including the bias term (2.8). This means that we do not simulate the annihilation process of topological defects. In order to simulate the annihilation process, much larger dynamical range is required, since each of N DW vacua is not homogeneously percolated in the simulation with the small dynamical range [20]. In our previous study [20], we could improve dynamical range by performing two-dimensional lattice simulations where the size of data is given by N 2 . In the present work, our aim is to investigate the spectrum of axions and gravitational waves produced by defect networks, while it is impossible to calculate the gravitational waves properly in the two-dimensional field theoretic simulations. Therefore, we just calculate the spectra in the onset of the scaling regime, and extrapolate the results for the case in which string-wall networks survive for a long time. Then we use the results of two dimensional simulations of the annihilation process of the networks [20] when we estimate the decay time of them [see Eq. (4.16)].

Evolution of string-wall networks
Now let us show the results of numerical simulations. We show the visualization of one realization of the simulation in Fig. 1. Using the identification scheme described in Appendix B, we confirmed that N DW domain walls attached to strings are formed at late time of the simulations. The structure of defect networks becomes more complicated when we increase the number N DW . Figure 2 shows the spatial distribution of the phase of the scalar field θ in the model with N DW = 3 at selected time slices. Since we give the initial configuration of the scalar field as Gaussian random amplitude (see Appendix A.3), the scalar field randomly oscillates around the minima |Φ| 2 = η 2 of the potential, given by the first term of Eq. (2.2). Strings are formed when this initial random configuration becomes relaxed. Subsequently, the value of the phase of the scalar field θ is separated into N DW = 3 domains. Domain walls are located around the boundaries of these domains. We note that the initial stage of simulations does not correctly follow the history of the universe, since domain walls are formed immediately after the formation of strings. As we noted in Sec. 2.1, in the realistic case, strings are formed much earlier than the formation of domain walls. This is due to the hierarchy between two different scales, Λ QCD /F a ∼ 10 −11 . On the other hand, in our choice of the parameters in numerical simulations, this ratio becomes as small as Λ QCD /F a ∼ √ m ∼ O(0.1) (recall that we used the unit η = 1). Because of the limitation of the dynamical range of the numerical simulations, we cannot choose this ratio as an arbitrary small value. Therefore, we concentrate on the evolution at late times where the string-wall networks evolve along to the scaling regime, and ignore the effect of the initial relaxation regime. We will discuss this point again in Secs. 3.2 and 3.3.
We define the length parameter ξ of strings 4) and the area parameter A of domain walls These parameters take a constant value when the defect networks enter the scaling regime (i.e. ρ string ∼ 1/t 2 and ρ wall ∼ 1/t). Henceforth we refer to these parameters as the "scaling parameters". We plot the time evolution of these parameters in Fig. 3. From this figure, we see that the defect networks are in scaling regime in the time scale τ 15. The length parameter ξ does not strongly depend on N DW , and its value is consistent with ξ ≃ 0.5-1.3 obtained in previous studies [11,28,29,32]. On the other hand, the value of the area parameter A increases for large N DW . This fact agrees with the intuitive argument that  In these figures, we used the same data that are used to visualize the result with N DW = 3 in figure 1. The value of θ varies from −π (blue) to π (red). At late times, the value of θ is separated into three domains represented by blue, red, and green region. Domain walls are located around the boundary of these three regions, θ = π/3, θ = −π/3, and θ = ±π. Strings, which are represented by white lines, pass through the point where three regions meet each other.
simulation box is unchanged. The right panel of Fig. 3 also implies that the area of domain walls in the real scalar field theory given by Eqs. (3.2) and (3.3) evolves similarly to the axionic model with N DW = 2. This means that the behavior of the networks for the case with N DW = 2 is similar to that produced by the model with Z 2 symmetry, except that they contain strings on the surface of walls. We note that the bump-like behavior of A in the real scalar field model at the initial stage of the simulation should be regarded as a numerical fake. We identify the position of domain walls as the point on which scalar field changes its sign in the model with the real scalar field, which is different from the identification scheme used for the model with the complex scalar field (see Appendix B.2). These different identification schemes cause the different behavior of A at the initial stage, where domain walls do not relax into the scaling regime.

Production of axions
From data of the scalar field obtained by numerical simulations, we compute the power spectrum P (k) of axions radiated by string-wall networks, which is defined by whereȧ(t, k) is the Fourier component of the time derivative of the axion field, and . . . represents an ensemble average. Here, we mask the position where topological defects exist, and extract the pure radiation component in the axion field produced by the defect networks. See Appendix C for details on the analysis method. Figure 4 shows the results for P (k, t f ) for various values of N DW evaluated at the final time of the simulation. We find that the spectrum has a peak around the momentum scale determined by the mass of the axion, k/R(t f ) ∼ m, or k ∼ R(t f )m ∼ 2. This implies that most of the axions are produced by the self-interaction of domain walls whose width is given by the inverse of the mass of the axion ∼ m −1 . However, this peak does not seem to fall off, but makes a tail at higher momenta. This tail-like feature was also found in another numerical study [11] in which we investigated the spectrum of axions produced by the annihilation of N DW = 1 domain walls bounded by strings. The appearance of this high-frequency tail can be understood in terms of the conjecture made in Refs. [26,27,33,34]. The strings whose radius smaller than the horizon size quickly decay due to the tension of domain walls in the time scale comparable to (or less than) the Hubble time. This fast process does not change much the energy spectrum of the axion field composing the string configuration, which has a form dE/dk ∼ 1/k with a high-momentum cutoff of order δ −1 s (the width of strings) and a low-momentum cutoff of order t −1 (the Hubble radius) [26]. This 1/k behavior makes a tail in the spectrum of radiated axions, as we see in Fig. 4. This argument cannot be applied for the spectrum of axions produced by global strings [28,29], which has a sharp peak at the horizon scale since in the absence of domain walls the strings lose their energy in the time scale much longer than the Hubble time. The high-momentum cutoff which is given by the width of the string ∼ δ −1 s is not clearly shown in Fig. 4 because of the poor resolution at the final time of the simulation.  Note that the amplitude of the high-frequency tail becomes small if we increase the domain wall number N DW . In the case of large N DW , the decay of loops of strings with small curvature radius might be "frustrated", because of the tension of domain walls acting on the strings from various directions. Then, the radiation of hard axions with spectrum ∼ 1/k is suppressed, compared with the case of small N DW .
We emphasize that our conclusion will be different from that made in [34] in support of [26,27]. In Ref. [34] it was claimed that the spectrum of axions radiated by domain walls behaves as ∼ 1/k in the whole frequency domain, and hence the mean energy of radiated axions is proportional to ln( √ λη/m) ∼ 60, which suppresses the relic abundance of cold dark matter axions. Our result, however, shows that the spectrum does not behave like ∼ 1/k except for the high-frequency modes, and most axions are mildly relativistic, making a peak at the low frequency. Hence there is no correction like ln( √ λη/m) ∼ 60, and the relic abundance of axions turns out to be large compared with the result of [34]. This fact leads to severe constraints on the axion models, as we will see in Sec. 4.
The above observation can be confirmed by estimating the energy of radiated axions. Roughly speaking, the product P (k)k represents the fraction of the energy density of axions with comoving wavenumber k against the total energy of axions. Hence the ratio of the energy between hard axions which contribute to the tail of the power spectrum and the energy of axions which contribute to the peak of the power spectrum can be estimated as On the other hand, the ratio between the energy of domain walls and the energy of strings at the final time of the simulations becomes E string /E wall ∼ [t f µ string ]/[t 2 f σ wall ] ∼ 10 −1 . Therefore, strings have 10% of the energy of domain walls at the final time, which contribute to the hard component ∼ 1/k of the power spectrum of radiated axions. Since the fraction E string /E wall decreases with time, we expect that the amplitude of the high-frequency tail would be negligible compared with the height of the peak of the spectrum if we follow the evolution for much larger dynamical ranges.
Since the population of the axions radiated by string-wall networks is dominated by lowmomentum modes which has the frequency comparable to the mass of the axion m, we expect that the mean energy of radiated axions is determined by the parameter m. Here, using the result of P (k, t 2 ), where t 2 is chosen at the late time of the simulation [see Eq. (C.12)], we compute the mean comoving momentum of radiated axions and its ratio to the axion mass The results are shown in table 2. We see that the value of the mean momentum is not strongly depend on the number N DW . However, the uncertainty becomes large for the case with large N DW . This is because a large portion of the simulation box is masked when we evaluate the power spectrum, which causes a large systematic errors when we evaluate the spectrum in the large scale 2 , as shown in Fig. 4. The claim is that the parameter defined by Eq. (3.8) does not depend on the time t 2 . To confirm it, we calculate the parametersk and ǫ a at different time steps, while fixing the reference time τ 1 = 14 in Eq. (C.12). The results are shown in table 3. We find that the value ofk changes proportionally to R(t 2 ). On the other hand, the ratio ǫ a between the physical momentumk/R and the axion mass m remains constant. Therefore, we use the parameter ǫ a as a constant of O(1) when we perform the analytic investigations in Sec. 4. We close this subsection by showing Fig. 5, the time evolution of the total energy density of axions and topological defects in the simulation box. From this figure, we see that the domain wall networks evolve as ρ wall ∼ 1/t, or ρ wall ∼ 1/τ 2 in conformal time. This behavior is expected from the property of the scaling solution, given by Eq. (2.6). We also find that the energy density of axions radiated by string-wall networks behaves similarly to that of topological defects at the late time. We will confirm this behavior in Sec. 4.1.

Production of gravitational waves
Next, we investigate the production of gravitational waves from string-wall networks. The spectrum of gravitational waves is represented by the quantity where ρ c (t) is the critical density of the universe at the time t, and ρ gw is the energy density of gravitational waves. We compute Ω gw using the data of the scalar field Φ obtained from the numerical simulation. The method of the calculation is summarized in Appendix D. Figure 6 shows the results for Ω gw (k, t) for various values of N DW . The slope of the spectrum changes at two characteristic scales. One is the Hubble radius (k ∼ 2πR(τ f )H(τ f ) ∼ 0.1) and another is the width of the domain wall (k ∼ R(τ f )m ∼ 2). This result agrees with the general argument made by two of the present authors (M. K. and K. S.) in [35]. However, the amplitude of Ω gw at high momenta is enhanced for the case with large N DW . This can be interpreted as follows. Since the number of domain walls is large in the model with large N DW , domain walls are "frustrated", and there are many small scale configurations whose correlation length is shorter than the Hubble radius. Therefore, the small scale anisotropy becomes large in the model with large domain wall number N DW . We also plot the spectrum of gravitational waves produced by domain walls in the real scalar field theory given by Eqs. (3.2) and (3.3). Although the evolution of domain wall networks shown in Fig. 3 is similar between the model with N DW = 2 and the model with real scalar field, the form of the spectrum of gravitational waves is slightly different between these two models at high momenta. This might be caused by the different form of the effective potential which determines the small scale structure of the defects.
As we discussed in previous subsections, at the early times (τ 15) topological defects do not relax into the scaling regime. However, the nonlinear dynamics of scalar fields during this initial stage produces a "burst" of gravitational waves, which may contaminate the spectrum of gravitational waves produced at late times where the defects networks evolve along to the scaling solution. We subtract this contribution with the same manner that we used for the subtraction of axions produced at the initial stage [see Eq. (D. 19)]. Figure 7 shows the time evolution of the spectrum of gravitational waves, where we compare the effect of this subtraction procedure. We see that the form of the spectrum is different at early times  . Time evolution of energy densities for various values of N DW . ρ defects represents the energy density of the string-wall networks which is evaluated by summing up the data of scalar field placed near the defects whose position is identified by using the algorithm described in Appendix B. ρ a,kin , ρ a,grad , and ρ a,pot are the kinetic, gradient, and potential energy density of radiated axions, respectively. ρ a,tot is the sum of ρ a,kin , ρ a,grad , and ρ a,pot . These energy densities of axions are evaluated by subtracting the data of scalar field which contribute to ρ defects .
if we make a subtraction. In particular, there is a bump at the scale given by the Hubble radius k/R ∼ H, and a plateau which falls off at the scale given by k/R ∼ m. This form is different from what we found in our previous studies [35,36], where a nearly "flat" spectrum extends between two characteristic scales. It seems that this difference in the form of spectra is caused by the contamination of radiations produced at the initial relaxation regime. We also note that the dynamical range of the past numerical studies (τ f /τ i ≃ 12) [35,36] is shorter than that of the present studies (τ f /τ i = 20). Therefore, the bump-like feature is less apparent in the past numerical studies with shorter dynamical ranges. The form of the spectrum at the final time of the simulation is similar between the result with subtraction and that without subtraction, which indicates that contaminations of the initial stage are diluted due to the cosmic expansion. The spectrum shown in Fig. 6 is normalized in terms of the peak amplitude in order to give the relative form of the spectrum between different theoretical parameter N DW . On the other hand, the total amplitude of gravitational waves is determined by the following arguments. Since the spectrum of gravitational waves has a peak at the momentum which is given by the Hubble scale, we regard that the most of the gravitational waves are generated at the length scale ∼ t. Suppose that the radiation of gravitational waves occurs with the time scale comparable to the Hubble time. The power of gravitational waves radiated from domain walls is estimated by using the quadrupole formula [37] where Q ij ∼ M DW t 2 is the quadrupole moment, and M DW ∼ ρ wall t 3 ∼ σ wall At 2 is the mass energy of domain walls, where we used Eq. (3.5). The energy of gravitational waves is estimated as We see that the energy density of gravitational waves does not depend on the time t Now, we define the efficiency ǫ gw of gravitational waves where ρ (sim) gw is the total energy density of gravitational waves computed from numerical simulations. We plot the time evolution of the parameter ǫ gw in Fig. 8. Although the value of ǫ gw does not remain in constant exactly, it seems to converge into a universal value of O(5-10).
We note that the time dependence of ǫ gw is strongly affected by the choice of the time τ g at which we subtract the contribution of gravitational waves produced in the initial relaxation regime [see Eq. (D. 19)]. This is shown in Fig. 9, where we plot the time evolution of ǫ gw for various values of τ g . If we choose τ g at an early time of the simulation, ρ (sim) gw contains the contribution of gravitational waves produced during the initial relaxation regime. Hence the behavior of ρ (sim) gw might be different from what we expect from Eq. (3.13), where we assumed that the energy density of domain walls is given by the scaling formula (3.5). This additional contribution is not negligible even at the late time of the simulations because the dynamical range of the numerical simulation is short. This constancy of ǫ gw should be tested in future numerical studies with improved dynamical ranges. Since the time variation of ǫ gw is small for large value of τ g , as shown in the plot of τ g = 25 in Fig. 9, we assume that ǫ gw takes a constant value when string-wall networks enter into the scaling regime, and use the value ǫ gw ≃ 5 when we estimate the abundance of gravitational waves in the next section.

Observational constraints
In this section, we give constraints on the parameters of axion models with long-lived domain walls, using the results of numerical studies described in the previous section. First, we calculate the relic abundance of cold axions and gravitational waves produced by string-wall networks in Sec. 4.1. Then we briefly review some experimental constraints in Secs. 4.2 and 4.3. We combine these constraints, and discuss their implications on the models with N DW > 1 in Sec. 4.4. We find that the overabundance of the cold axions produced by domain walls gives a stringent bound on the model parameters, which excludes the most part of models, except for tiny loopholes. Finally we comment on some exceptions in Sec. 4.5.

Axion cold dark matter abundance
Here we assume that the energy of strings becomes negligible compared with that of domain walls, and axions and gravitational waves are only produced by domain walls. In this approximation, the evolution of energy densities for domain walls, axions, and gravitational waves is described by the following coupled equations  where A and ǫ gw are constants. Substituting these expressions to the above equations, we obtain From Eqs. (4.4), (4.6), and (4.7), we find Note that the sign of the right hand side of this equation becomes negative for t (Gσ wall ) −1 . This is not problematic, since domain walls have to collapse before the time t ∼ (Gσ wall ) −1 otherwise they overclose the universe [see Eq. (2.7)]. Let us define the energy and number of axions per comoving box where n a (t) is the number density of axions at the time t. Integrating Eqs. (4.2) and (4.8), we find where t r is the time when domain walls begin to radiate axions, which might be chosen as the time of the QCD phase transition. In the regime t ≫ t r , it reduces to For sufficiently early times t ≪ (Gσ wall ) −1 , the first term of Eq. (4.12) dominates over the second term, which leads to the behavior ρ a ∝ 1/t. This agrees with what we observed in numerical simulations, shown in Fig. 5. The ratio between E a (t) and N a (t) can be determined by the result of numerical simulations where ǫ a is defined by Eq. (3.8). Suppose that domain walls disappear at the time t dec and the radiation of axions and gravitational waves is terminated at that time. In this case, the energy densities of axions and gravitational waves at the present time t 0 are given by where the subscript "dec" in the left hand side of Eq. (4.14) represents that the axions produced by the decay of domain walls. The annihilation of domain walls occurs in the time scale estimated by Eq. (2.9) (4.16) The numerical coefficient can be fixed to α ≃ 18, according to our previous study [20]. Note that the relation ) where h is the reduced Hubble parameter today: H 0 = 100hkm·sec −1 Mpc −1 . We can fix the numerical values of A, ǫ a and ǫ gw from the result of numerical simulations. Then, the relic abundance of axions and gravitational waves is determined by three theoretical parameters: N DW , Ξ, and F a . The second term in the square bracket in Eq. (4.19) represents the effect of gravitational radiation, which is negligible at early times but becomes relevant when domain walls survive for a long time. In particular, when domain wall survives enough time such that they overclose the universe, t dec ∼ (Gσ wall ) −1 , the contribution of the second term in the square bracket in Eq. (4.19) becomes O(1). In such a case, the gravitational field around walls cannot be regarded as a small perturbation [38], and we must use the full general relativistic treatment, which is out of the scope of this paper. Hence, the expression (4.19) is valid only for the time t ≪ t WD ∼ (Gσ wall ) −1 . This condition is automatically satisfied as long as we consider the parameter region where domain walls do not overclose the universe. We note that there are other contributions to the axion cold dark matter abundance, in addition to Ω a,dec . One is the coherent oscillation of the homogeneous axion field (zero mode) [7] and another is the production of axions from global strings [24]. Their contributions are estimated as [11] Ω a,0 (t 0 )h 2 ≃ 1.2 × F a 10 12 GeV where the subscripts "0" and "string" represent that they are the contributions of zero modes and global strings, respectively 3 . In the above expressions, we fixed the QCD scale such that Λ QCD = 400MeV. The total abundance of cold dark matter axions is given by the sum of Eqs.  We require that Ω a h 2 should not exceed the observed value of the abundance of cold dark matter Ω CDM h 2 = 0.11 [39]. This gives an upper bound on F a and a lower bound on Ξ.

Neutron electric dipole moment
Next we consider the effect of the Ξ term on the degree of strong CP-violation, following Refs. [17,18]. The QCD Lagrangian contains a term where G aµν is the gluon field strength,G a µν is the dual of it, andθ is a dimensionless parameter. This term contributes to the neutron electric dipole moment (NEDM) with the amplitude d n = 4.5 × 10 −15θ ecm [3]. The recent experimental bound on the NEDM gives |d n | < 2.9 × 10 −26 ecm [40], which requires θ < 0.7 × 10 −11 . (4.25) The original idea of Peccei and Quinn is to setθ = 0 by introducing a symmetry [2]. However, if the additional term (2.8) exists in the effective potential, it shifts the CP-conserving minimum with a magnitude controlled by the parameters Ξ and δ.
Substituting the parametrization Φ = ηe ia/η into the effective potential given by Eqs. (2.2) and (2.8), and expanding for small a/η, we obtain the effective potential for a where m phys is the effective mass of the axion Requiring that this should not exceed the experimental bound (4.25), we obtain an upper bound for Ξ.

Astrophysical bounds
It is known that various astrophysical observations give a lower bound on F a . Here we just summarize the constraints obtained so far. See [41,42] for more comprehensive reviews. First of all, the emission of axions in the sun gives a constraint on the axion-photon coupling, g aγγ ≡ (α em /2πF a )c aγγ , where α em is the fine-structure constant, and c aγγ is a numerical factor of O(1). The helioseismological consideration on the anomalous solar energy loss due to the axion emission gives a bound g aγγ < 10 −9 GeV −1 [43]. There are also direct experimental searches. The Tokyo Axion Helioscope [44] gives a bound g aγγ < (5.6-13.4) × 10 −10 GeV −1 for 0.84eV < m < 1.00eV. Phase I of the CERN Axion Solar Telescope (CAST) [45] gives g aγγ < 8.8 × 10 −11 GeV −1 for m 0.02eV. It is improved in Phase II [46] as g aγγ < 2.2 × 10 −10 GeV −1 for m 0.4eV. The sensitivity is expected to be improved up to g aγγ < few × 10 −12 GeV −1 in the next generation helioscopes [47].
The observation of globular clusters gives another bound. The helium-burning lifetimes of horizontal branch (HB) stars give a bound for axion-photon coupling g aγγ 0.6 × 10 −10 GeV −1 [48]. Furthermore, the delay of helium ignition in red-giant branch (RGB) stars due to the axion cooling gives g aee < 2.5 × 10 −13 [49], where g aee ≡ C e m e /F a is the axionelectron coupling, m e is the mass of the electron, C e = cos 2 β/3, and tan β is the ratio of two Higgs vacuum expectation values.
The axion-electron coupling is also constrained by the observation of white-dwarfs. The cooling time of white-dwarfs due to the axion emission gives a bound g aee < 4 × 10 −13 [50]. Recently, it is reported that the fitting of the luminosity function of white-dwarfs and the pulsation period of a ZZ Ceti star is improved due to the axion cooling, which implies the axion-electron coupling g aee ≃(0.6-1.7)×10 −13 [51] or g aee ≃(0.8-2.8)×10 −13 [52].
Finally, the axion-nucleon coupling is constrained from the energy loss rate of the supernova 1987A [41]. This gives the strongest astrophysical bound on the axion decay constant compared with other bounds described above, F a > 4 × 10 8 GeV. (4.29) Hence in the following we take it as a lower bound on F a when we combine it with other observational constraints. Note that the estimation of the axion emission rate suffers from various numerical uncertainties which may modify the bound by a factor of O(1) [41].

4.4
Implication for models with N DW > 1 In Fig. 10, we plot the observational constraints described in the previous subsections. Based on the estimation given by Eq. (4.23), we find that the overabundance of cold axions gives a stringent lower bound on the bias parameter, Ξ 10 −50 -10 −52 . This lower bound is much stronger than that obtained from the overclosure of domain walls (2.10). Combined with the NEDM constraint, it completely excludes the parameter region if we assume that the phase δ of the bias term (2.8) is O(1). This bound might be weakened if we assume the extremely small value of δ, but this assumption spoils the genius of the original PQ solution to the strong CP problem. We must fine-tune δ in order to solve the fine-tuning problem ofθ. We find that the amplitude of gravitational waves produced by domain walls is weak Ω gw h 2 10 −20 in the parameter space considered above. This is irrelevant to any observations, since even the ultimate phase of DECIGO [53], which is a space-borne interferometer planned to launch in the future, will have a sensitivity Ω gw h 2 ∼ 10 −18 [54]. The weakness of the signal of gravitational waves is due to the fact that domain walls annihilate in early epoch in order to avoid the overclosure of axions produced from them, and the dilution factor of cosmic expansion becomes large.
The difference of this result from our previous study [20] arises due to the different estimation on the mean momentum of radiated axions. In the present analysis described in Sec. 3.2, it turns out that the radiated axions are mildly relativistic with the momentum comparable to their mass, k/R ≃ m, or ǫ a ≃ 1. However, in the previous study [20], we assumed that the radiated axions have a mean energy ω ≃ γm where γ ≃ 60, according to the result of the past numerical study 4 [34]. This leads to an underestimation of the number of axions produced by domain walls by a factor ǫ 2 a + 1/γ ∼ 10 −2 . In this case, the bound which comes from the overclosure of cold axions is weakened. The large factor of γ originated in the assumption made in [34], that the radiated axions have a hard spectrum dE/dk ∼ 1/k which leads to a logarithmic correction γ ∼ ln( √ λη/m) ∼ 60. In Sec. 3.2 we have shown that this contribution is negligible compared with the contribution of the peak located at the axion mass. We note that there is a parameter region around F a ∼ 10 6 GeV, called the "hadronic axion window" [56][57][58], that is not excluded by observations. This occurs due to the ambiguity in light quark masses which leads to the cancellation in the axion-photon coupling g aγγ for KSVZ models. In such a case, the astrophysical bounds on g aγγ do not significantly constrain the value of F a . Here we do not consider this possibility, as it strongly depends on the quantum number assignment of the PQ-charged fermions [58]. Also, we note that in this case the axion becomes a candidate of hot dark matter [57] which can be constrained by observation of the large-scale structure [59].
From the results shown in Fig. 10, we conclude that the axion models with N DW > 1 are excluded if the PQ symmetry is broken after inflation. The exception is the case in which the value of δ is suppressed. As we see in Fig. 10, the degree of tuning in δ is < 1.1×10 −2 in order to avoid all observational constraints. Furthermore, by using Eqs. (4.19) and (4.28), one can show that the degree of tuning does not strongly depend on the value of N DW . Hence, about 1% fine tuning is required regardless of the value of N DW .

Scenario with extremely small δ
Although there is no theoretical reason to consider the extremely small value of δ, it is possible, by accident or design of some fundamental physics, to yield δ whose value is small enough to avoid the NEDM bound. Let us comment on such a situation. In this case, the value of Ξ can take larger value than that we considered so far. In particular, if Ξ becomes as large as the axion mass is not determined by QCD instanton effect, but by the correction term m * [see Eq. (4.27) and also Ref. [18]]. Then, the cosmological history is significantly modified. When the condition (4.30) is satisfied, the magnitude of the bias term (2.8) exceeds that of the QCD potential [the second term of Eq. (2.2)]. Hence, "domain walls" with the tension σ wall ≃ 8m * F 2 a are formed when the Hubble parameter becomes H ≃ m * . The subsequent history is categorized into two possibilities, analogously to the scenarios described in Sec. 2. One possibility is that strings are get attached by domain walls at the time t ∼ m * , and they quickly disappear due to the tension of walls. Let us call this scenario case I. On the other hand, if the structure of the term (2.8) is the form ∝ Φ n , where n is an integer, n domain walls are attached to strings, and they survive for a long time. Let us call this scenario case II. These string-wall networks annihilate due to the effect of a "bias" term which arises from the QCD instanton effect, δV ∼ Λ 4 QCD ∼ m 2 F 2 a . This occurs at the time t dec, * ∼ σ wall /Λ 4 QCD ∼ 8m * /m 2 . Requiring that domain walls should disappear before the epoch of big bang nucleosynthesis (BBN), t dec, * < t BBN ∼ 1sec, we obtain (4.31) In case I, the abundance of axions is estimated in the usual way, as a sum of the coherent oscillation, strings, and decay of domain walls bounded by strings. The number of axions is fixed at t ∼ t * ≡ m −1 * . At that time, the energy density of axions is given by ρ a (t * ) ∼ 10 × m 2 * F 2 a , where the factor 10 arises from the fact that the abundance of axions becomes larger by an order of magnitude if we include the contribution of strings and domain walls [11,29]. The present energy density becomes ρ a (t 0 ) ∼ 10(R(t * )/R(t 0 )) 3 m 2 * F 2 a , which leads where the subscript "nt" indicates that they are produced non-thermally. Axions are also produced from thermal bath of the primordial plasma. The difference from the usual scenario is that the axion mass is given by m * so that their relic energy density can be large. The number of axions is fixed when their interactions with the thermal bath with the process mediated by gluons freeze out. Their number at the present time is given by n th a (t 0 ) ≃ 7.5cm −3 [60], which leads Ω th a h 2 ≃ 3.2 × 10 −7 × N DW F a 10 10 GeV This bound is more severe than that obtained from the lifetime of domain walls, given by Eq. (4.31). However, if Ξ is much larger than the values considered above, axions become heavy and no longer stable. Their dominant decay channel is a decay into two photons. The lifetime is given by where Γ γ is the decay late of the process a → γγ, and we fixed the numerical coefficient c aγγ = 1 in the axion-photon coupling for simplicity. Since the axion mass is not given by m, we must treat m * and g aγγ as different parameters. Cosmological bounds on such models are investigated in [61]. The heavy axions which decay at early times lead to various effects on the cosmological observations, such as the distortion in the spectrum of CMB, and the disagreement of the theoretical prediction with the observed light element abundance. Each of the observations gives a severe constraints on the coupling g aγγ . In order to avoid these constraints, the lifetime of axions must be short enough to decay earlier than the onset of the BBN, t γ 10 −2 sec. This requirement leads a bound (4. 38) In summary, if the value of Ξ is sufficiently large, the mass of axions is determined by the term proportional to Ξ, and this parameter is again constrained by various observational considerations. In order to avoid the observational bounds, axions should be light enough to avoid the overclosure of the universe, whose condition is given by Eq. (4.34) or (4.36), or heavy enough to decay before the BBN epoch, whose condition is given by Eq. (4.38).
- 25 -In this paper, we have investigated the axion cosmology with N DW > 1, where domain walls survive for a long time in the early universe. We performed three-dimensional lattice simulations to follow the evolution of string-wall networks in the scaling regime, and computed the spectra of axions and gravitational waves produced from them. The radiated axion is mildly relativistic with the momentum comparable to the axion mass, which is similar to the result of the model with N DW = 1 [11]. The ratio of the mean momentum of axions to their mass, given by Eq. (3.8), has a value ǫ a ≃ 1. The spectrum of gravitational waves has a peak at the momentum comparable to the Hubble scale k/R ∼ H, with a plateau extending up to the momentum comparable to the inverse of the width of domain walls k/R ∼ m. The magnitude of gravitational waves is given by Eq. (3.13) with the numerical coefficient ǫ gw ≃ 5.
Using the numerical results, we calculated relic densities of axions and gravitational waves produced by long-lived string-wall networks. We found that the significant amounts of axions are produced from topological defects, and the requirement that they should not overclose the universe gives a stringent constraint on the model parameters. Combined with bounds given by astrophysical observations and the experimental result on the NEDM, we showed that the most part of parameter regions are excluded unless the value of the phase of the bias term δ is extremely fine-tuned. In such parameter regions, the amplitude of gravitational waves is too weak to observe even in future gravitational wave detectors with high sensitivities.
Cosmological domain wall problem leads to serious constraints on the axion models. The options that we can choose to avoid this problem are listed as follows: 1. The axion is realized such that N DW = 1. This can be obtained by introducing a heavy quark 5 [5]. Domain walls disappear shortly after their formation, and there is no cosmological problem. In this case, axion mass is constrained into m ≃ 10 −3 -10 −2 eV [11].
2. The PQ symmetry is broken during inflation. In this case, the population of topological defects is wiped away beyond the horizon scale, and domain walls are not problematic. However, the observation of the isocurvature fluctuations in anisotropies of CMB gives constraints on the model parameters [15].
3. The existence of Ξ term annihilates domain walls. This possibility is allowed only if the phase δ of the bias term is extremely suppressed. However, the value of Ξ is still constrained by cosmology, such as Eqs. (4.34), (4.36), and (4.38).
If the PQ symmetry is a global symmetry arising accidentally in the low energy effective theory, it is easily violated by higher dimensional operators which is suppressed by powers of the Planck mass M P [19]. These operators might be identified as the Ξ term, and their magnitude must be suppressed in order not to shift the value ofθ. The existence of longlived domain walls makes a problem much worse, since the value of Ξ cannot be arbitrary small. Hence we conclude that the axion models with N DW > 1 seems to be excluded from cosmological considerations, except for the cases 2 and 3 enumerated above. However, in the cases 2 and 3, some tunings on other theoretical parameters such as the inflationary scale (for case 2) or the value of δ (for case 3) are required. Furthermore, even in the case 1, there still remains a question of how Ξ is suppressed, which should be addressed from some fundamental theories.

A Numerical schemes
In this appendix, we describe the analysis method which we used in the numerical studies. Our numerical code is the combined version of that we used in past numerical studies [11,20,29,62].

A.1 Field equations
The equation of motion for Φ is obtained by varying the Lagrangian density given by Eq. (2.1) with the potential The second term of Eq. (A.1) is different from that of Eq. (2.2). The original potential (2.2) is not well defined at Φ = 0, which causes instabilities in the numerical studies. We avoid this singularity by multiplying a factor |Φ| in front of the cosine term. The difference between Eq. (2.2) and Eq. (A.1) is not important in the bulk region on which |Φ| = η, and we observe that the quantitative behavior such as time evolution of the scaling parameter of topological defects is not so much affected by this modification.
In the numerical studies, we follow the evolution of two real scalar fields φ 1 and φ 2 which are defined by Φ = φ 1 + iφ 2 . We normalize all dimensionful quantities in the unit of η, i.e. η = 1. With this normalization, the equations of motion for two scalar fields in the FRW background becomë If we use the conformal time τ (dτ = dt/R(t)), the above equations become where f ≡ RΦ = f 1 +if 2 , and a prime denotes a derivative with respect to the conformal time τ . Here, we assume the radiation dominated background. This model has three parameters: λ, N DW and m. In the numerical studies, we fix the values of λ and m, and vary the value of N DW .

A.2 Discretization of differential equations
To execute the time evolution which follows from Eqs. (A.4) and (A.5), we used the symplectic integration scheme developed in [63]. In this numerical scheme the (implicit) Hamiltonian of the system is exactly conserved and there is no secular growth of the conserved quantity. In practice, we compute the following sequence of mappings, from r = 1 to r = 4, in each time step. The values of coefficients c r and d r are given by [63] It turns out that this integration scheme has fourth-order accuracy in time.
The spatial derivative and the Laplacian of the scalar fields are computed by using the finite difference method. For example, the Laplacian is given by where f i,j,k = f (τ, x i , y j , z k ), and ∆x = L/N is the grid spacing. This formula has fourthorder accuracy in ∆x. We imposed the periodic boundary condition on the boundaries of the simulation box.

A.3 Initial conditions
We generate initial conditions in the similar way with Refs. [20,35]. We treat φ 1 and φ 2 as two independent real scalar fields, and assume that each of them has massless quantum fluctuations at the initial time. With this assumption, we generate initial conditions as Gaussian random amplitudes in momentum space. Then we Fourier transform them into the configuration space to give the spatial distribution of the fields. We emphasize that this choice of initial conditions does not take account of the correct circumstances in the QCD epoch. However, we expect that the result in the late time of the simulations are qualitatively unchanged if we used the different initial conditions, since the evolution of defect networks in the scaling regime is not so much affected by the initial field configurations. See appendix of Ref. [20] for more discussions.

B Identification of topological defects
Using the lattice algorithm described in Appendix A, we can solve the time evolution of two real scalar fields φ 1 (t, x) and φ 2 (t, x) in the comoving simulation box. From these data of scalar fields, we estimate various physical quantities, such as the length of strings, the area of domain walls, and the spectra of axions and gravitational waves radiated by defect networks.
In order to calculate such physical quantities, we have to identify the position of topological defects in the simulation box. In this appendix, we summarize our identification method of topological defects, and describe how to calculate scaling parameters of strings ξ and domain walls A defined in Sec. 3.1.

B.1 Identification of strings
We used the method developed in [29] to identify the position of strings. Consider a quadrate which consists of four neighboring grids in the simulation box (see points A, B, C, and D in Fig. 11). At each vertex point, we estimate the phase of the complex scalar field θ from the data of two real scalar fields φ 1 and φ 2 . Then, we obtain a mapping from the quadrate in the real space to the phase distribution in the field space, as shown in Fig. 11. Let us denote the minimum range of the phase which contains all the images of four vertices as ∆θ. The string exists within the quadrate if ∆θ > π is satisfied, as long as the value of the phase changes continuously around the core of strings. We can also determine the position at which a sting penetrates the quadrate. At the boundary of the quadrate, we determine the positions of points on which φ 1 = 0 or φ 2 = 0 is satisfied, by using linear interpolation of the values of φ 1 and φ 2 . If the string correctly penetrate the quadrate, there are two points with φ 1 = 0 and two points with φ 2 = 0 on the boundary of the quadrate. Then we determine the position of the string as the point at which the line connecting two points with φ 1 = 0 intersects with the line connecting two points with φ 2 = 0 (see left-hand panel of Fig. 11).
This simple criterion breaks down when the quadrate is penetrated by more than two strings. However, we observed that such region is at most 1% of the whole simulation box when the system relaxes into the scaling regime. Therefore, this scheme enables us to determine the position of strings with at least 99% accuracy.

B.2 Identification of domain walls
The identification of domain walls is simpler than that of strings. We separate the region of the phase of Φ, 0 ≤ θ < 2π, into N DW domains. For example, if N DW = 3, we obtain three domains with 0 ≤ θ < π/3 and 5π/3 ≤ θ < 2π (vac. 1), π/3 ≤ θ < π (vac. 2), and π ≤ θ < 5π/3 (vac. 3), as shown in Fig. 11. At each grid point, we compute the phase of the scalar field and assign the number of the vacuum domain (i.e. vac. 1, vac. 2, or vac. 3) which contains that point. Let us call the neighboring grid points that differ by a unit lattice spacing ∆x = L/N as the "link". We identify that the domain wall intersects the link if the number of vacuum is different between the two ends of the link.

B.3 Calculation of scaling parameters
The scaling parameters of strings and domain walls are given by where ρ string is the energy density of strings, ρ wall is the energy density of domain walls, µ string is the mass energy of strings per unit length, and σ wall is the surface mass density of domain walls. It might be difficult to determine the values of µ string and σ wall directly from the data of scalar fields since they are quantities obtained by integrating ρ string and ρ wall over the coordinates which are perpendicular to the direction along which the defects extend. Instead of estimating µ string and σ wall , we calculate the length l of strings and area A of domain walls in the comoving coordinates, such that where V = L 3 is the comoving volume of the simulation box. Then, the scaling parameters are given by The length of strings l can be computed by connecting the loci of strings identified by the method described in the previous subsection. The computation of the area of domain walls A is not so straightforward, since it depends on the way how we define the segment of A in each lattice. Here, we use the simple algorithm introduced by Press, Ryden, and Spergel [64]. Define the quantity δ ± which takes the value 1 if the number of vacuum is different between the two ends of the link, and the value 0 otherwise. We sum up δ ± for each of the links with the weighting factor where ∆A = (∆x) 2 is the area of one grid surface, and θ ,x , etc. is a derivative of θ(x) with respect to x. The weighting factor given by the function of ∇θ gives the average number of links per area segment. This method takes account of the orientation of the surface of walls in the lattice. See [64] for details.

C Calculation of power spectrum of radiated axions
The calculation of power spectrum of axions radiated from topological defects is executed by using the masking analysis method developed in Refs. [11,29]. In this appendix, we describe the main line of the calculation scheme. Readers may also refer to [11,29] for complete description on the analysis method. The power spectrum P (k, t) of axions is defined by where . . . represents an ensemble average.ȧ(t, k) is the Fourier component of the time derivative of the axion fieldȧ (t, k) = d 3 xe ik·xȧ (t, x), where the value ofȧ(t, x) is obtained from the simulated data of Φ andΦ If a(t, x) is a free field, it can be shown that P (k, t) is proportional to the energy spectrum of axions [11]. The data ofȧ obtained by numerical simulations contain the contaminations from the core of strings or domain walls. We mask the contribution of axion field near strings or domain walls by introducing a window function Then, we compute the power spectrumP (k) of the masked axion fielḋ Unfortunately, the masked spectrumP (k) is not equivalent to the true spectrum of free axions. This discrepancy can be resolved by using the pseudo-power spectrum estimator (PPSE) Here, M (k, k ′ ) is the window weight matrix defined by where W (k) is the Fourier transform of the window function W (x), Ω k is a unit vector representing the direction of k, and dΩ k = d cos θdφ. It turns out that the ensemble average of P PPSE (k) is equivalent to the spectrum of radiated axions [29]. Naively, M (k, k ′ ) in Eq. (C.8) can be evaluated by averaging |W (k − k ′ )| 2 /V 2 over discrete k and k ′ in the lattice. This however requires O(N 2 ) arithmetics, which is practically prohibited when the grid number N is large as is in our analysis. Instead, we approximate M (k, k ′ ) by partially taking the continuous limit, which we here describe briefly. With a transformation P = k − k ′ and R = k+k ′ 2 , Eq. (C.8) can be rewritten as (C.9) In the continuous limit, the integration over R can be performed analytically.
where P = |P|, R = |R| and µ ≡ R·P RP . We thus obtain Evaluation of Eq. (C.11) in a lattice has a good scaling with only O(N ) arithmetics. On the other hand, the approximation can cause a systematic error in estimation of P PPSE (k), which shall be discussed shortly later.
The power spectrum P PPSE (k, t) contains the contributions of axions produced before the time t. This means that it contains the contributions of axions radiated at initial epoch where the defects do not relax into the scaling regime. Since we are interested in the spectrum of axions radiated from topological defects which evolve along the scaling solution, we must subtract the contribution of axions produced at initial relaxation regime. This can be achieved by evaluating the power spectra for selected time steps t 1 and t 2 > t 1 , and subtracting the spectrum obtained at t 2 from that obtained at t 1 [11] ∆P (k, where (ω(k, t 2 )/ω(k, t 1 ))(R(t 1 )/R(t 2 )) 3 is a dilution factor which comes from the effect of the cosmic expansion, and ω(k, t) = m 2 + k 2 /R(t) 2 (C.13) is the energy of axions with momentum k/R(t). Here, t 1 should be chosen as the time of the onset of the scaling regime. Throughout this paper, we use τ 1 = 14 as the subtraction time, at which the system begins to follow the scaling regime. Figure 12 shows the effect of this subtraction procedure. We see that the position of the peak is falsified, and the amplitude of the spectrum is overestimated if we do not perform the subtraction. Figure 12 also shows that the spectrum is overestimated at higher momentum scale when we do not mask the data near the position of topological defects.  The red plot shows the spectrum evaluated with masking the contamination of the core of topological defects and subtracting the contribution of axions produced at initial stage (τ ≤ 14). The green plot shows the spectrum obtained by the same model without masking the contamination of defects. The blue plot shows the spectrum evaluated without masking and subtraction.
We note that the form of the power spectrum is affected by the choice of the region where the masking takes place. In Eq. (C.4), we mask the field data over grid points within a distance d width from the position of the core of defects which we identified by using the methods described in Appendix B. Here we choose the critical distance as d width = n width δ s = n width λ −1/2 for strings and d width = n width δ w = n width m −1 for walls, where δ s and δ w are widths of strings and domain walls, respectively, and n width is an arbitrary number. Figure 13 shows the power spectrum for various values of n width . We see that the power spectrum in large k is overestimated when we use the smaller value of n width . Since the spectra agree within the error bars for the results with n width ≥ 2, we use n width = 2 for numerical analyses presented in the text. Finally, we comment on the estimation of errors of the power spectrum. Error is estimated by the following procedure based on the Monte Carlo sampling. After we compute P PPSE (k), we generateȧ(k) as Gaussian random field with a variance P PPSE (k)/k 2 in the Fourier space. Then we execute the inverse Fourier transformation ofȧ(k) to obtain the data ofȧ(x) in the real space. Let us denote this spurious field data asȧ MC (x). We apply the same masking procedure to the fieldȧ MC (x) and obtain the power spectrum P PPSE,MC (k). We iterate this procedure by n MC times, and compute the mean power spectrum and the covariance matrix of P PPSE,MC (k) where P i PPSE,MC (k) is the power spectrum obtained at i-th sampling step. Although PPSE should be unbiased as formally shown in Ref. [29], we observed thatP PPSE,MC (k) does not converge into P PPSE (k) obtained by the field-theoretic simulation. This bias (P PPSE,MC (k) − P PPSE (k)) arises mainly from use of the approximation M (k, k ′ ) in Eq. (C.11). In our analysis, we regard the bias as the systematic error. Then the covariance matrixC(k, k ′ ), which accounts for both the statistical and systematic errors, can be given bỹ C(k, k ′ ) = C(k, k ′ ) + (P PPSE,MC (k) − P PPSE (k))(P PPSE,MC (k ′ ) − P PPSE (k ′ )) ). (C. 16) In addition to these errors in estimation of P PPSE (k) in each realization, P PPSE (k) also varies from one realization to another, variance from which should be also included as a statistical error. Thus energy spectrum averaged over realizations and its covariance matrix can be finally given bŷ where N r is the number of realizations, and the subscript (j) indicates a realization. We execute the numerical simulations with N r = 20 different realizations of initial field configurations. The results we present in the text is obtained fromP PPSE (k) andĈ(k, k ′ ).

D Calculation of gravitational waves
We use the method introduced in [65] to calculate the spectrum of gravitational waves. The evolution of gravitational waves is described by spatial metric perturbations h ij around the FRW background ds 2 = −dt 2 + R 2 (t)(δ ij + h ij )dx i dx j . In the Fourier space, the equations of motion for metric perturbations reduce to the form [35] ∂ 2 ∂y 2 + 1 χ ij = S ij k 2 , (D. 4) where we defined the variable y ≡ kτ and the source S ij (k, τ ) ≡ 16πGR(τ )T TT ij (k, τ ). (D.5) The transverse traceless (TT) part of the stress-energy tensor is computed by applying the projection operator T TT ij (k, τ ) = Λ ij,kl (k)T ij (k, τ ) = Λ ij,kl (k){∂ k Φ * ∂ l Φ}(k, τ ), (D.6) Λ ij,kl (k) = P ik (k)P jl (k) − 1 2 P ij (k)P kl (k), (D.7) P ij (k) = δ ij −k ikj , (D.8) wherek = k/|k|, and {∂ k Φ * ∂ l Φ}(k, τ ) is the Fourier transform of ∂ k Φ * (x, τ )∂ l Φ(x, τ ). There are two homogeneous solutions of Eq. (D.4), cos y and sin y. The special solution is given by the time integral of the source term convoluted with two homogeneous solutions ij (k, τ ), (D.10) and C (2) ij (k, τ ) = y y i dy ′ cos kτ S ij (k, τ ′ ) k 2 ≡ 16πG k 2C (2) ij (k, τ ). (D.11) The energy density of the gravitational waves is given by where a prime denotes a derivative with respect to conformal time τ . We replace the ensemble average in Eq. (D.12) by an average over a volume V of the comoving simulation box, Substituting the solution given by Eq. (D.9), and ignoring the terms with higher order in RH, we obtain (D.14) where we averaged over a period of the oscillation of χ ij (k, τ ) with time. The spectrum of gravitational waves is given by and ρ c (τ ) is the critical density of the universe at conformal time τ ρ c (τ ) = 3H 2 8πG .
(D. 18) In radiation dominated background, we expect that Ω gw (k, τ ) ∝ S k (τ ), since R 4 H 2 remains constant. Therefore, we calculate S k (τ ) in the numerical studies. Note that Ω gw (k, τ ) given by Eq. (D. 16) does not correspond to the spectrum observed today. The spectrum of gravitational waves observed today is obtained by multiplying Ω gw (k, τ ) with a damping factor which arises from the following matter-dominated epoch.
As we discussed in Sec. 3.3, we subtract the contribution of gravitational waves produced at the initial epoch where the defects do not relax into the scaling regime. This subtraction can be executed as follows. Since the energy density of gravitational waves shifts as ∝ R −4 , the quantity S k (τ ) ∝ R(τ ) 4 dρ gw (τ )/d ln k is proportional to the number of gravitons with the momentum k per comoving volume. Therefore, the difference of the spectral function ∆S k (τ ) ≡ S k (τ ) − S k (τ g ) (D. 19) gives the spectrum of gravitational waves produced during the time interval [τ g , τ ]. Here τ g is a reference time which we treat as a free parameter in the numerical simulations.

E Surface mass density of domain walls
In this appendix, following Ref. [8], we shortly comment on the surface mass density of domain walls in the models used in numerical simulations. First, consider the model of scalar field φ with the potential given by where λ r and η r are constants. This model has a discrete Z 2 symmetry in which the scalar field transforms as φ → −φ. In the flat Minkowski background, the classical field equation of this model has a static solution φ(z) = η r tanh λ r 2 η r z , which describes a domain wall lying perpendicular to the z-axis. By integrating the energy density of the field over the direction perpendicular to the surface of the wall, we obtain the surface mass density Next, let us consider the axionic model with the Lagrangian given by Eq. (2.1). If we restrict ourselves to low energy configurations at the QCD scale, we can put Φ = ηe iθ , which gives the effective potential for θ L θ = η 2 2 (∂ µ θ) 2 + m 2 η 2 N 2
We note that the surface mass density of axionic domain wall, given by Eq. (E.6) should be modified if we include the structure of the neutral pion field which varies inside the wall. In this case, the surface mass density is given by [66] σ wall = 4.32f π m π η/N DW ≃ 9.23mη 2 /N 2 DW , (E.7) where f π is the pion decay constant, and m π is the mass of the pion. In the second equality, we used m u /m d ≃ 0.48 [3], where m u and m d are the mass of up quark and down quark, respectively.