A review of gravitational waves from cosmic domain walls

In this contribution, we discuss the cosmological scenario where unstable domain walls are formed in the early universe and their late-time annihilation produces a significant amount of gravitational waves. After describing cosmological constraints on long-lived domain walls, we estimate the typical amplitude and frequency of gravitational waves observed today. We also review possible extensions of the standard model of particle physics that predict the formation of unstable domain walls and can be probed by observation of relic gravitational waves. It is shown that recent results of pulser timing arrays and direct detection experiments partially exclude the relevant parameter space, and that a much wider parameter space can be covered by the next generation of gravitational wave observatories.


Introduction
The progress of direct observations [1,2] of gravitational waves (GWs) is bringing about drastic developments in astrophysics and cosmology. We expect to obtain a lot of important information about physics at very high energies from direct observations of GWs due to the fact that they interact very weakly with matter and hence preserve almost all the features characterizing astrophysical or cosmological events [3,4]. Experimental sensitivities for the direct detection of GWs have been improved substantially during the past decades, and many new GW observatories are now planned to be built in the world. In this context, it is worth investigating various possible sources of GWs and clarifying to what extent we can extract the information about new physics from future observations. So far, various cosmological sources of relic GWs are discussed in the literature, such as the primordial amplification of vacuum fluctuations [5][6][7][8][9][10], cosmological phase transitions [11,12], cosmic strings [13][14][15], and preheating after inflation [16][17][18][19][20]. Furthermore, the recent starting of the GW astronomy can provide a distinctive way of testing General Relativity and other theories of gravity [21].
In this article, we consider domain walls as possible cosmological sources of GWs. Domain walls are sheet-like topological defects, which might be created in the early universe when a discrete symmetry is spontaneously broken [22]. Since discrete symmetries are ubiquitous in high energy physics beyond the Standard Model (SM), many new physics models predict the formation of domain walls in the early universe. By considering their cosmological evolution, it is possible to deduce several constraints on such models even if their energy scales are much higher than that probed in the laboratory experiments.
In general, the formation of domain walls is regarded as a problem in cosmology [23], since their energy density soon dominates the total energy density of the universe, which conflicts with the present observational results. However, we can consider the possibility that domain walls are unstable and collapse before they overclose the universe [24][25][26]. Their unstability might be guaranteed if the discrete symmetry is only approximate and explicitly broken by a small parameter in the theory. In such a scenario, a significant amount of GWs can be produced during the process of collisions and annihilations of domain walls, and they may remain as a stochastic GW background in the present universe. Observations of such relic GWs will enable us to trace the events in the very early universe and provide a new way of investigating physics at very high energies.
The purpose of this article is to review the physics of cosmic domain walls and to evaluate the detectability of GWs produced by them in present and future observations. In particular, we summarize the results of recent theoretical developments including various particle physics motivations so far proposed in the literature and the methods to estimate the GW signatures based on field theoretic lattice simulations.
We note that the production mechanism of GWs discussed in this article is different from that discussed in the context of phase transitions [e.g. Refs. [11,12,[27][28][29][30]]. In the latter case, a strong first order phase transition is assumed, and GWs are produced due to the collision of bubbles and subsequent turbulences. On the other hand, in the former case it is not necessary to assume the first order phase transition, and GWs are produced due to the late-time motion and interaction of domain walls. The key ingredient of this scenario is the existence of quasi-degenerate vacua after the phase transition.
This review is organized as follows. The theoretical basics of domain walls and cosmological constraints on them are described in Sec. 2. The semi-analytical approach to estimate the production of GWs from domain walls is discussed in Sec. 3. Section 4 deals with particle physics models that predict the production of a significant amount of GWs from domain walls. In Sec. 5, we compare the GW signatures with sensitivities of present and future experiments. Finally, we conclude in Sec. 6.
2 Domain walls and cosmology 2.1 Field theory As an illustrative example, let us consider the following toy model of a real scalar field φ, Note that the potential V (φ) has two degenerate minima at φ = ±v. In this theory there is a discrete Z 2 symmetry, under which the field transforms as φ → −φ. This discrete symmetry is spontaneously broken when the scalar field acquires a vacuum expectation value (VEV), φ = ±v.
The scalar field takes one of the two discrete values (+v and −v) after the spontaneous symmetry breaking, which means that two different domains can appear. Domain walls are produced around the boundary of these two domains. Consider a static planar domain wall configuration lying perpendicular to the z-axis in the Minkowski space, φ = φ(z). Solving the field equation We see that φ(z) approaches ±v as z → ±∞, and it rapidly changes around z = 0. The width of the domain wall δ can be estimated as a typical length scale of the spatial variation of φ(z), The energy-momentum tensor for the static solution φ = φ(z) is given by Integrating T 00 over the direction perpendicular to the wall, we obtain its surface energy density, Note that the integration of the spatial components results in the same quantity: dzT 11 = dzT 22 = −σ. Therefore, σ is also referred to as the tension of domain walls.
Another interesting example is the following model of a real scalar field a, where the field a is defined within a finite domain [0, 2πv], and N is a positive integer. This kind of potential naturally arises in the context of axion models, which will be discussed in Sec. 4.2. In this model, there is a discrete Z N symmetry under which the field transforms as a/v → a/v + 2πk/N with k = 0, 1, . . . , N − 1. This symmetry is spontaneously broken once the field a settles down to one of N degenerate minima of the potential, and domain walls are formed around the boundary of these vacua. If we consider a planar wall orthogonal to z-axis a = a(z), the solution of the classical field equation reads This configuration interpolates between two vacua, a/v = 2πk/N at z → −∞ and a/v = 2π(k + 1)/N at z → +∞. From the above solution, we can estimate the thickness of the wall as The tension of domain walls reads (2.10) As we have seen in above two examples, the properties of domain walls can be characterized by two model-dependent quantities, the tension σ and the thickness δ. In general, the wall thickness is roughly given by Compton wavelength of the field which causes the spontaneous breaking of discrete symmetry, while the tension is estimated in terms of the height of the potential energy V 0 separating the degenerate minima, In the following, we do not specify the magnitude of σ and provide some model-independent arguments about the evolution of domain walls. We will come back to the mode-dependent issues in Sec. 4.

Cosmological evolution
Domain walls exist if different vacua are populated in the universe. Whether such a distribution of vacua appears or not depends on the cosmological initial conditions. In particular, it is widely believed that the universe underwent a period of exponentially rapid expansion, called inflation. Depending on the conditions at the inflationary epoch, the formation of domain walls must be seriously taken into account.
Suppose that the toy model scalar field φ introduced in the previous subsection stayed at a certain vacuum before the inflationary period. In such a setup, we naively expect that domain walls do not exist in the present universe, since a domain on which φ takes an uniform value exponentially glows during inflation and the size of such a domain is much larger than the present horizon size. However, such a naive expectation is not necessarily true. During inflation, the field φ acquires vacuum fluctuations of order δφ ∼ H inf /2π if its effective mass m φ [m 2 φ = 2λv 2 in the model given by Eq. (2.2)] is smaller than H inf [31][32][33], where H inf is the Hubble parameter during inflation. Once such a condition is satisfied, the φ field easily jumps into other domains within one Hubble time, and as a result many different domains can exist after inflation, which leads to the formation of domain walls. Furthermore, even if H inf is sufficiently small such that the φ field never acquires large fluctuations during inflation, it can thermalize and have thermal fluctuations due to the reheating after inflation. If this is the case, the discrete symmetry is thermally restored when the maximum temperature after inflation T max becomes larger than m φ . After that, domain walls are formed when the universe cools below some critical temperature. Therefore, we expect that the formation of domain walls can happen if either the Hubble parameter during inflation H inf or the maximum temperature after inflation T max is sufficiently larger than the mass m φ of the field φ. 1 After the formation of domain walls, their dynamics can be described by two kinds of forces. One is the tension force, which is given by where R wall represents a typical curvature radius of the walls. The other is the friction force, which appears if there is a interaction of the field composing the core of domain walls and particles in thermal bath. It can be estimated as [36] p F ∼ ∆p · n ∼ vT 4 , (2.13) where ∆p ∼ T v is a typical momentum transfer due to the collision with a particle, n ∼ T 3 is the number density of particles, and v is the velocity of domain walls. These two forces are balanced, p T ∼ p F , and from this condition we obtain where we have used T 4 ∼ M 2 Pl /t 2 assuming the radiation dominated background, and M Pl ≃ 2.435 × 10 18 GeV is the reduced Planck mass. Let t r denote the time at which domain walls become relativistic. Equations (2.14) and (2.15) imply that their curvature radius becomes comparable to the horizon size at that time, R wall ∼ M 2 Pl /σ ∼ t r . We also see that their energy density ρ wall dominates over the total energy density of the universe ρ c at that time, In other words, they remain non-relativistic as long as their energy density is subdominant. The friction force is exponentially damped when the temperature of the background radiations becomes less than the mass of particles that interact with domain walls. Therefore, we expect that the effect of the friction force becomes negligible at sufficiently late times if domain walls only interact with massive particle states. On the other hand, if they interact with lighter particles such as those in the SM, we must carefully evaluate the effect of the friction force in order to describe their late time evolution. In what follows, we focus on the case where the friction force becomes negligible at sufficiently early times.
Once the friction force becomes irrelevant, the dynamics of domain walls is dominated by the tension force, which stretches them up to the horizon size. Many numerical studies [37][38][39][40][41][42][43] confirmed that the evolution of domain walls in this regime can be described by the scaling solution, in which their energy density evolves according to the simple scaling law ρ wall ∝ t −1 , and their typical size is given by the Hubble radius ∼ t. 2 An analytic method to calculate the evolution of domain walls was also proposed in Refs. [44,45], which again showed the existence of the scaling solution.
It will be convenient to parameterize the energy density of domain walls in the scaling regime as  [47,48]. It was shown that the value of A for the case with N = 2 agrees with Eq. (2.18), and that it increases proportionally with N . The energy density of domain walls in the scaling regime ρ wall ∝ t −1 decays slower than that of cold matters ρ matter ∝ R −3 (t) and radiations ρ rad ∝ R −4 (t), where R(t) is the scale factor of the universe. Therefore, they gradually dominate the energy density of the universe. From the condition ρ c (t) = ρ wall (t), we estimate the time at which the wall domination occurs, (2.19) Here we assumed that the total energy density of the universe is dominated by radiations before t = t dom , i.e., ρ c = 3M 2 Pl /4t 2 . Once domain walls dominate the energy density of the universe, the subsequent evolution of the universe is drastically altered. The equation of state for an isotropic gas of non-relativistic domain walls is given by w = −2/3 [36], which implies that the scale factor in the wall dominated universe evolves as (2.20) Such a rapid expansion is incompatible with standard cosmology. Even if the energy density of domain walls is subdominant at the present time, they may cause another problem. Since their typical curvature radius is comparable to the Hubble radius, they introduce large scale density fluctuations, whose magnitude is estimated as where G is Newton's gravitational constant, and we used t 0 ∼ H −1 0 with t 0 and H 0 being the present cosmic time and the Hubble constant, respectively. The observation of the cosmic microwave background radiation implies δρ/ρ O(10 −5 ), from which we obtain the following condition This constraint was first discussed in Ref. [23], and it is referred to as the Zel'dovich-Kobzarev-Okun bound. We see that domain walls with a tension as large as σ > O(MeV 3 ) must not exist in the universe at the present time.

Biased domain walls
One possible solution to the domain wall problem is to introduce an energy bias in the potential, which lifts the degenerate minima [24][25][26]. 3 Let us consider the model for the real scalar field φ discussed in Sec. 2.1. Here we artificially introduce the following term in addition to Eq. (2.2), where ǫ is a dimensionless constant. The modified potential is shown in Figure 1. This potential has minima at φ = ±v, but there is an energy difference between them, Because of the existence of this energy difference, domain walls become unstable and eventually collapse. Note that the additional term (2.23) explicitly breaks the discrete Z 2 symmetry. Therefore, this solution works if the discrete symmetry is not exact, but holds only approximately. We note that domain walls cannot be created from the beginning if the energy deference between two vacua V bias is sufficiently large [25]. In order to clarify the condition to have domain walls in the presence of the energy bias, let us consider the probabilities p + and p − in which the scalar field ends up in the plus vacuum (φ = +v) and the minus vacuum (φ = −v), respectively, after the phase transition. The ratio between these two probabilities is given by where ∆F = V bias · ξ 3 is the difference of the free energy between two vacua, ξ is the correlation length at the epoch of the phase transition, and we estimate T as the Ginzburg temperature, T ≃ V 0 · ξ 3 , with V 0 being the height of the potential barrier between two minima. The above equation implies that the spatial distribution of two vacua after the phase transition becomes asymmetric if V bias = 0. According to the prediction of percolation theory, the critical value above which an infinite cluster of the minus vacuum appears in the space is given by p c = 0.311, if the system is treated as a three dimensional cubic lattice [50]. Requiring that a large cluster of the false vacuum appears in the space (p − > p c ), we obtain In other words, large scale domain walls are expected to be formed as long as the above condition is satisfied. Even if V bias is sufficiently small such that domain walls are created at the phase transition, the false vacuum region tends to shrink due to the existence of the energy difference: There is a volume pressure force acting on the walls, whose magnitude is estimated as p V ∼ V bias . The collapse of domain walls happens when this pressure force becomes greater than the tension force p T ∼ σ/R wall . If we assume that domain walls have reached the scaling regime beforehand, their typical curvature radius is given by R wall ≃ t/A [see Eq. (2.17)]. Hence, from the condition that two forces become comparable, p V ∼ p T , we can estimate their annihilation time: where C ann is a coefficient of O(1). If the annihilation occurs in the radiation dominated era, the temperature at t = t ann is given by where g * (T ) is the relativistic degrees of freedom for the radiation energy density at a given temperature T . The value of C ann (or C d in Ref. [48]) can be determined from numerical simulations. It typically takes the value of C ann ≃ 2-5, depending on N for the model with the Z N symmetry [Eq. (2.7)]. It also depends on the choice of the criterion to determine the decay time of domain walls in the simulations. For more details, see Ref. [48]. Note that the lifetime t ann is inversely proportional to V bias . If the energy bias is sufficiently small, domain walls live for a long time. Requiring that their collapse occurs before they overclose the universe t ann < t dom [see Eq. (2.19)], we obtain the lower bound 4 on the magnitude of the energy bias, V bias > 4C ann A 2 σ 2 /3M 2 Pl , or (2.29) By using Eq. (2.28), this condition can be rewritten in terms of the annihilation temperature, Even if domain walls are annihilated before they overclose the universe, their decay products may behave as dangerous relics, which places additional constraints on the magnitude of the energy bias. In particular, if domain walls decay into the SM degrees of freedom, the decay products can destroy light elements created at the epoch of BBN, which conflicts with the standard cosmological scenario. The ratio between the energy density of domain walls and the entropy density around that time is estimated as where g * s (T ) is the relativistic degrees of freedom for the entropy density at the temperature T corresponding to the cosmic time t. According to the constraints on energy injection at the epoch of BBN [53,54], we must require that the lifetime should be shorter than t ann 0.01sec, if we assume that a significant fraction of the energy density of domain walls is converted into energetic particles. This condition leads to another lower bound on the magnitude of the energy bias, (2.32) We note that this constraint is derived under the assumption that the decay products strongly interact with SM particles, and hence it depends on details of underlying particle physics models. If some stable relics are produced from long-lived domain walls, they would behave as dark matter and contribute to the energy density of the present universe. In this case, the tension of domain walls and the magnitude of the energy bias are further constrained from the observed dark matter abundance. Such a constraint is particularly relevant to axion models, which will be discussed in Sec. 4.2.

Estimation of gravitational waves from domain walls
Domain walls having a tension larger than the bound (2.22) must not exist at the present time, but there is a possibility that they are annihilated before they overclose the universe due to the existence of the energy bias. It is expected that such collapsing domain walls produce GWs, which are potentially observable today.
The production of GWs from cosmic domain walls was discussed by several authors [24,55,56], while the first quantitative study aiming at comparing the GW signatures and the sensitivities of experiments was carried out in Ref. [57]. In Ref. [57], the energy density of the relic GWs was estimated by solving the evolution of collapsing domain walls numerically and specifying some ansatzes for initial field configurations. A more improved estimation was performed in Refs. [58,59], where the production and the evolution of domain walls in the expanding universe were investigated based on the field theoretic lattice simulations, and the spectrum of GWs was computed by applying the method introduced in Ref. [60]. The results of Refs. [58,59] were updated in Ref. [46] by correcting some error in the numerical code and improving the dynamical range of the simulations.
Let us roughly estimate the energy density of GWs produced by domain walls. Here we assume that the energy density of domain walls obeys the scaling law (2.17), and that the typical time scale of the gravitational radiation is given by the Hubble time ∼ t. 5 According to the quadrupole formula, the power of the gravitational radiation is given by P is the quadrupole moment of domain walls, and M wall ∼ σAt 2 is their mass energy. Therefore, the energy density of GWs ρ gw ∼ P t/t 3 reads From this estimate we expect that the energy density of GWs produced by domain walls is proportional to the square of their tension σ 2 and remains almost constant. Strictly speaking, the above quadrupole formula cannot be directly applied to domain walls, since it is only valid in the far-field regime, while the domain wall network is an spatially extended medium. In order to check the validity of this estimate, an alternative formalism must be employed. The formalism to compute the production of GWs from dynamical scalar fields is descried in Ref. [60]. In this approach, the spectrum of GWs can be numerically computed by using Green functions with transverse-traceless parts of the stress-energy tensor of the scalar field. In this way, the features shown in Eq. (3.1) can be checked by performing detailed numerical simulations.
In Refs. [46,58,59], numerical simulations of domain walls were performed with the aim of computing the spectrum of GWs produced by them. In the numerical studies, the evolution of the real scalar field φ in the simple toy model with Z 2 symmetry [Eqs. (2.1) and (2.2)] was investigated by solving the field equation in the Friedmann-Robertson-Walker background, where the potential V (φ) is given by Eq. (2.2). The simulations were executed in the 3D cubic lattice with the periodic boundary condition. The radiation dominated background [R(t) ∝ t 1/2 ] was assumed in Refs. [46,58], while the evolution in the matter dominated background [R(t) ∝ t 2/3 ] was also investigated in Ref. [59]. From the configuration of the scalar field in the numerical simulations, one can estimate the spectrum of GWs produced by them. The spectrum can be computed by using Green functions with transverse-traceless parts of the stress-energy tensor of the scalar field [60]. The results obtained in Ref. [46] are shown in Figure 2, where S k represents the spectrum of GWs per unit logarithmic frequency interval, with V being the volume of the comoving simulation box. The spectrum has a peak at the scale corresponding to the Hubble radius. We note that the horizontal axes in Figure 2 represents the comoving wavenumber, and hence the location of the peak k peak shifts according to k peak /R(t) ∼ H(t). Since the smallest scale of the structure of domain walls is given by their core width δ [see Eq. (2.4)], the spectrum falls off at a large wavenumber corresponding to that scale, k/R(t) ∼ δ −1 . Furthermore, S k increases as ∼ k 3 for k < k peak , and decreases as ∼ k −1 for k > k peak . The behavior S k ∝ k 3 at small k can be deduced from causality [46,62]. The estimate in Eq. (3.1) can be checked by computing the following quantitỹ where the subscript "peak" means that the quantity is evaluated at the peak of S k . The results of numerical simulations clearly show that the value ofǫ gw remains almost constant after domain walls enter into the scaling regime, and it is estimated as [46] ǫ gw ≃ 0.7 ± 0.4, (3.5) where the error corresponds to the statistical uncertainty. Furthermore, the value ofǫ gw hardly depends on the choice of the value of the parameter λ, which determines the tension σ [see Eq. (2.6)]. These facts are consistent with the expectation that the amplitude of GWs is given by Eq. (3.1) during the scaling regime. Let us estimate the peak amplitude of GWs produced by long-lived domain walls. The spectrum of GWs at the cosmic time t is characterized by the following quantity [3,4]: where f = k/2πR(t) is the frequency corresponding to the comoving wavenumber k. From Eq. (3.4), we have the peak amplitude at the annihilation time of domain walls, Here, we assume that the production of GWs is suddenly terminated at t = t ann , 6 and that it happens during the radiation dominated era. Then, the peak amplitude of GWs at the present time t 0 is given by  (3.9) Note that the large GW amplitude is predicted if the tension σ is large and the annihilation temperature T ann is low, which corresponds to the case where domain walls lived for a long time. We can also estimate the peak frequency of GWs in terms of the Hubble parameter at the annihilation time of domain walls, Hz g * (T ann ) 10 The high annihilation temperature T ann results in the high peak frequency. Note that there is a cutoff frequency corresponding to the width of domain walls, which is much higher than the peak frequency. The results of numerical simulations imply that the spectrum of GWs behave as Ω gw ∝ f −1 for the intermediate frequency range f peak < f < f δ . So far we have assumed that the annihilation of domain walls happens during the radiation dominated era. If it happens before reheating, the above estimates are modified accordingly. Let us assume that the energy density of the universe is dominated by that of the inflaton, which behaves as non-relativistic matter, when domain walls are annihilated. Since the Hubble parameter evolves as H 2 ∝ R −3 (t) at that stage, instead of Eq. (3.8) we have where T reh is the reheating temperature, and H reh and H ann are the Hubble parameters at T = T reh and T = T ann , respectively. In the case of the perturbative decay of the inflaton, the Hubble parameter at this stage is given by [63] (3.13) Using this relation, we obtain for T ann > T reh . (3.14) Furthermore, the peak frequency reads f peak ≃ 8.9 × 10 −2 Hz g * s (T reh ) 100 −1/3 g * (T ann ) 100 T reh 10 4 GeV The spectrum of GWs for the Z N symmetric model [Eq. (2.7)] was also analyzed in Ref. [47]. Similar to the above model with a real scalar field, the spectrum has a peak at the scale corresponding to the Hubble radius, and the peak amplitude agrees with the estimate based on Eq. (3.1). However, the shape of the spectrum at f > f peak differs from that in the real scalar field model, and it slightly changes according to the value of N . This N dependence might be caused by the fact that many configurations whose sizes are smaller than the Hubble radius are produced in the model with large N , which results in the enhancement of the amplitude of GWs at high frequencies.

Particle physics models
In the previous section, we have shown that the amplitude of GWs is determined by two parameters, the tension of the domain wall σ and the temperature at the domain wall annihilation T ann . The latter is related to the energy bias V bias for quasi-degenerate vacua [see Eq. (2.28)]. The values of σ and V bias depend on underlying particle physics models, and hence the prediction for the peak amplitude and its frequency differs according to the details of the models. In this section, we briefly review various particle physics models proposed in the literature that predict the formation of unstable domain walls and the production of GWs from them.

Standard Model Higgs field
Intriguingly, there is a possibility that the dynamics of the SM Higgs field induces the formation of unstable domain walls, which can produce a significant amount of GWs [64]. According to the recent analysis of the effective potential of the Higgs field based on the measured values of the Higgs boson mass and the top quark mass, our electroweak vacuum is likely to be metastable in the framework of the SM [65,66]. Indeed, the solution of the renormalization group equation for the Higgs self coupling implies that it becomes negative at some energy scale Λ, which is much higher than the electroweak scale. It is probable that there exists some new physics around that scale, which lifts the Higgs potential. If this is the case, the effective potential of the Higgs field can have two minima, which leads to the formation of Higgs domain walls in the early universe.
The effect of new physics can be modeled by introducing a higher dimensional operator in the Higgs potential, where ϕ represents the SM Higgs field value, and λ(ϕ) is the field-dependent Higgs self coupling obtained by solving the renormalization group equation and treating the Higgs field value as the renormalization scale, λ(µ) = λ(ϕ). Here we ignore the quadratic term which leads to the electroweak vacuum, since its effect on the dynamics of the Higgs field is negligible at high energies. The instability scale Λ is sensitive to the top quark mass, and it can take a value from 10 10 GeV to the Planck scale within the error of the measured top quark mass. Up to the detailed values of Λ and the self coupling, the Higgs potential can have two minima, one of which is our electroweak vacuum, and the other one is at a higher energy scale ϕ = ϕ f determined by minimizing Eq. (4.1). Furthermore, one can consider the possibility that the high-scale minimum ϕ f is just a local minimum, and that two minima are quasi-degenerate, i.e., the energy difference V bias between two minima is much smaller than the height of the potential energy V 0 separating them (see Figure 1). Let us assume that the Higgs potential has quasi-degenerate minima as described above. If the inflationary scale is sufficiently high, the Higgs field acquires large quantum fluctuations during inflation. After inflation, the Higgs field takes different values in different patches of the universe, which results in the formation of domain walls. These domain walls are annihilated when the effect of the energy bias V bias becomes relevant, and subsequently the electroweak minimum dominates the universe. The tension of the domain wall can be estimated as [64] σ ∼ ϕ 2 where the width of the domain wall δ ∼ ϕ f /V 1/2 0 can be fixed by minimizing the tension. The precise values of ϕ f , V 0 , and V bias should be obtained by solving detailed renormalization group equations, and they depend on the values of various parameters in the SM such as the top quark mass, the strong gauge coupling, and the Higgs boson mass. In Ref. [64], the magnitude of the energy bias V bias was treated as a free parameter as it can be adjusted by tuning the value of Λ, and it was shown that there exists a parameter region in which a significant amount of GWs is produced by long-lived domain walls. The typical peak frequency reads f peak ∼ 10 −3 -10 2 Hz, which is relevant to future direct detection experiments. The peak frequency cannot be lower than this range, since a smaller value of ϕ f is required, which cannot be realized in this framework.
The production of GWs from Higgs domain walls was also investigated in Ref. [67]. Contrary to the above discussions, it was concluded that the amplitude of GWs is too small to be observed in the planned detectors. However, it should be noted that the scenario considered in Ref. [67] is different from the above scenario in the sense that the high-scale minimum is located at a superplanckian value and that two minima are non-degenerate. In such a case, domain walls do not enter into the scaling regime and collapse soon after the formation, leading to the small amplitude of relic GWs.

Axion models
The axion [68,69] appears in the extensions of the SM with the Peccei-Quinn (PQ) mechanism [70,71], which has been proposed as a solution to the strong CP problem of quantum chromodynamics (QCD). It arises as a pseudo Nambu-Goldstone boson when a hypothetical global U(1) symmetry (called the PQ symmetry) is spontaneously broken. Its interaction with other particles is suppressed by a large decay constant F ∼ O(10 9 -10 11 ) GeV, and hence it is regarded as one of the best motivated candidates of cold dark matter [72][73][74]. Furthermore, string theory suggests the existence of many axion-like particles (ALPs) [75,76]. For more comprehensive reviews, see Refs. [77][78][79][80].
The crucial feature of the axion models is that they predict the formation of domain walls if the PQ symmetry is restored and broken after inflation [81]. The global U(1) PQ symmetry is explicitly broken to its subgroup Z N due to topological charge fluctuations in the QCD vacuum [82,83], and the effective potential for the axion field a at low energies is described by that in Eq. (2.7), where the axion mass is given by m ∼ F π m π /F ∼ 6 µeV (10 12 GeV/F ), F π ≃ 92 MeV is the pion decay constant, m π ≃ 135 MeV is the pion mass, and F = v/N is the axion decay constant. In the early universe, first the line-like objects, called global strings, are formed due to the spontaneous breaking of the U(1) PQ symmetry when the temperature of the universe becomes T ∼ v. Subsequently, domain walls are formed around the epoch of QCD phase transition. At that time strings are attached by N domain walls, and the hybrid networks of strings and domain walls, called stringwall systems are formed. 7 The tension of axionic domain walls is given by Eq. (2.10), where the approximation implies that there would be some finite corrections in the zero-temperature effective potential [83,86], which we ignore for simplicity. The evolution of string-wall systems differs according to the number of degenerate minima N . If N = 1, the systems collapse soon after the formation due to the tension of domain walls [87], and the present energy density of GWs produced from them is too small to observe. On the other hand, they are stable if N > 1, and we need to introduce explicit symmetry breaking terms in order to guarantee that they are annihilated before they overclose the universe [56,81]. For instance, Planck-suppressed higher dimensional operators can induce sufficiently small energy bias between N degenerate minima. It should be noted that such Planck-suppressed operators induce a large CP-violating effect which spoils the original PQ solution to the strong CP problem [88][89][90][91][92][93], and that the dimension of those operators must be sufficiently high in order to avoid the experimental limit on the CP violation [94]. These kinds of higher dimensional operators naturally arise if we assume that the PQ symmetry is an accidental symmetry of an exact discrete symmetry [95][96][97]. In Ref. [98], it was argued that the long-lived domain walls in the axion models with N > 1 can produce a significant amount of GWs with the peak frequency f peak ∼ 10 −11 Hz. However, it turned out that such a parameter region is excluded since the abundance of cold axions produced by long-lived domain walls exceeds the observed cold dark matter abundance [47,48]. There still 7 In some exceptional cases, domain walls may be formed even if strings do not exist. For instance, if the initial value of the axion field is tuned to the location which is very close to the top of the cosine potential (2.7) and its fluctuations are sufficiently large, domain walls without strings can be formed around the time of QCD phase transition. Domain walls without strings can also be formed due to the level crossing between the axion and an ALP [84,85], if there exists an ALP whose mass is comparable to the axion mass around the epoch of QCD phase transition. remains some parameter region which avoids all observational constraints, and in such a region the predicted amplitude of GWs is very small, Ω gw h 2 10 −20 [47].
A similar argument can be applied to the models with ALPs, but in such models the ALP mass is not necessarily related to its decay constant. Therefore, one can treat them as two independent parameters in low energy phenomenology. In particular, if there exist some couplings between ALPs and SM particles and the ALP mass is sufficiently large, it is possible to avoid the dark matter overclosure bound, since ALPs produced by long-lived domain walls can decay into radiations. In Ref. [99], it was pointed out that domain walls in the ALP models can produce baryon asymmetry of the universe as well as GWs. A very high peak frequency f peak ∼ O(100) kHz is predicted in this scenario, since the temperature at the domain wall annihilation must be high, T ann 10 11 GeV, in order to generate sufficiently large baryon asymmetry.
The formation of domain walls and the production of GWs are also predicted in the context of the aligned axion models [100][101][102], which have been built explicitly by applying the clockwork mechanism discussed in Ref. [103]. In the aligned axion models, the axion a is described in terms of the flat direction of plural axion-like fields φ i (i = 1, . . . , N ax ), where N ax is the total number of the axion-like fields. A large decay constant F for the axion can be realized even though the actual decay constants F i for N ax axion-like fields are much smaller than F . Since the symmetry breaking scales F i are much smaller than the usual PQ scale F ∼ O(10 9 -10 11 ) GeV, the CP violating effects from Planck suppressed operators remain small, which naturally explains the high quality of the PQ symmetry.
As an explicit ultraviolet completion of the aligned axion model, one can consider a model based on N ax complex scalar fields associated with N ax global U(1) symmetries. It is assumed that (N ax − 1) U(1) symmetries are explicitly broken due to the operators proportional to some small

Supersymmetric models
The rich structure of supersymmetric theories gives rise to various possibilities of the formation of domain walls in the early universe. In Refs. [104,105], the formation of domain walls and the production of GWs in the context of the spontaneous breaking of discrete R symmetries were discussed. Let us assume that there exist some hidden SU(N c ) gauge interactions in addition to those in the SM. In the supersymmetric extensions of gauge theories, there exist fermionic partners of gauge bosons, called gauginos. Such gauginos may settle down in a condensate in the early universe due to the corresponding strong gauge forces. Topological charge fluctuations associated with such gauge forces break the global U(1) R symmetry of the theory down to its discrete subgroup Z 2Nc , and this Z 2Nc symmetry is spontaneously broken further down to the Z 2 subgroup due to the gaugino condensation (see e.g. Ref. [106]). It is known that the effective potential has N c degenerate vacua after the gaugino condensation, and domain walls with the tension are formed around that time [107,108], where Λ c represents the scale at which the gauge interactions become strong.
In order to avoid the cosmological domain wall problem, it is necessary to introduce a term that induces the energy bias between degenerate vacua. Such a energy bias is obtained if we assume that there exists a constant term w 0 in the superpotential, which explicitly breaks the discrete Z 2Nc symmetry. Note that the magnitude of the constant term should be w 0 ∼ m 3/2 M 2 Pl in order to cancel the positive contribution to the cosmological constant associated with supersymmetry breaking effects, where m 3/2 is the mass of gravitinos. This constant term results in the energy bias in the effective potential, , we see that the annihilation of domain walls occurs when the Hubble parameter becomes comparable to the gravitino mass, This result implies that there is a possibility to probe the gravitino mass from the observation of GWs [104]. For instance, if the domain wall annihilation occurs during the radiation dominated era, the peak frequency is given by f peak ∼ 10 3 Hz (m 3/2 /1 TeV) 1/2 . The formation of domain walls is also predicted in the next-to-minimal supersymmetric SM (NMSSM). The NMSSM is a possible extension of the minimal supersymmetric SM (MSSM), in which an additional gauge singlet superfield is introduced in order to provide a solution to the µ-problem [109] of the MSSM. Here, µ is a dimensionful parameter appearing in the superpotential of the MSSM, µH u H d with H u and H d being two Higgs doublet superfields, and its magnitude should be of the order of the soft supersymmetry breaking scale rather than the natural cutoff scale such as the Planck scale. In the NMSSM, a discrete Z 3 symmetry is imposed in order to forbid all dimensionful quantities in the superpotential, and the µ term with the appropriate magnitude is induced due to the dynamics of the scalar component S of the singlet superfield, see Refs. [110,111] for reviews.
The Z 3 symmetry is spontaneously broken when the S field acquires expectation values, and domain walls are formed around that time. The tension of domain walls depends on various parameters including the singlet-Higgs couplings and soft supersymmetry breaking parameters, but typically σ 1/3 ∼ O(TeV) if the singlet-Higgs couplings are relatively large. On the other hand, in the decoupling limit where the singlet-Higgs couplings become much smaller than unity, the tension can take much larger values [112], where κ is a dimensionless coupling appearing in the superpotential, W ⊃ (1/3)κS 3 , S ∼ m soft /κ is the VEV of the singlet scalar, and m soft represents the soft supersymmetry breaking mass scale. Due to the large tension, the amplitude of GWs produced from domain walls can be enhanced accordingly.
It was argued that the domain wall problem in the NMSSM cannot be solved just by introducing Planck-suppressed operators which provide the energy bias among different vacua [113], since such interactions radiatively induce a large tadpole operator that destabilizes the VEV of the singlet field. One possible solution is to impose various additional symmetries to arrange the form of Planck-suppressed interactions such that the scalar potential only contains a small bias term, V bias ∼ ζm 3 soft S + h.c., where ζ is a loop suppression factor [114]. Another possibility is to assume that the Z 3 symmetry is anomalous for QCD or some hidden strong gauge interactions [115]. In this case, the energy bias is given by V bias ∼ Λ 4 s , where Λ s is the scale at which the corresponding gauge interactions become strong. In any case domain walls must be annihilated before the epoch of BBN, since their decay products would destroy light elements [see Eq. (2.32)]. This fact implies that the peak frequency of GWs produced from domain walls should be higher than f peak ∼ 10 −9 Hz. Assuming that the domain wall annihilation occurs just before BBN, one can constrain the parameters of the NMSSM in the decoupling limit from pulsar timing observations [112].
In Ref. [116], the production of GWs from domain walls formed after thermal inflation was discussed. Thermal inflation [117,118] is introduced in order to suppress the abundance of harmful light long-lived scalar fields, called moduli, appearing in the context of string cosmology [119,120]. In this model, a short period of inflation is driven by the potential energy of a scalar field called flaton, which is trapped at the origin of the scalar potential because of thermal effects. A discrete Z n symmetry is imposed in order to guarantee the flatness of the flaton potential, where n is an integer satisfying n ≥ 4. After thermal inflation, the thermal effects become irrelevant and the flaton acquires non-zero VEVs, which break Z n symmetry and lead to the formation of domain walls. The domain walls can be annihilated if we introduce an additional term that explicitly breaks Z n symmetry [121,122] or if we assume that the Z n symmetry is anomalous for QCD [55,116]. It was shown that such domain walls can produce a significant amount of GWs if their lifetime is sufficiently long, and the peak frequency typically lies in the range of 10 −6 -10 −3 Hz [116].

Implications for present and future observations
In this section, we discuss implications of GWs from domain walls for ongoing and planned experimental searches. Here we marginalize the model-dependences described in the previous section, and treat σ and T ann as free parameters to clarify the parameter region that is relevant to present and future observations.
The leading ground-based interferometer is Advanced LIGO [123], whose first observing run placed limits on the amplitude of the stochastic GW background Ω gw < 1.7 × 10 −7 with 95% confidences for 20-86 Hz by assuming a flat GW spectrum [124], and these limits are about 33 times tighter than the previous limits set by Initial LIGO and Virgo [125]. In addition to them, ground-based interferometer KAGRA [126] will soon start to run in Japan. In Europe, the more advanced ground-based observatory, Einstein Telescope (ET) [127], is planned with the aim of achieving further improvement in sensitivity. Space-borne interferometers such as eLISA [128,129] and DECIGO [130] are planned to be launched in the future, and they will enable us to explore lower frequency ranges, which cannot be probed in the ground-based experiments. Much lower frequencies ∼ 10 −9 -10 −8 Hz are probed by using the pulsar timing array (PTA). Recently, European Pulsar Timing Array (EPTA) set a limit on the amplitude of a flat stochastic GW background Ω gw h 2 < 1.2 × 10 −9 at a reference frequency of f = 1yr −1 [131], and it is about one order of magnitude tighter than that obtained in previous PTA searches [132][133][134]. The sensitivity will be improved in future PTA projects such as SKA [135] and FAST [136].
Sensitivities of various ongoing and planned experiments are summarized in Figure 3. For the sensitivities of Advanced LIGO, we plot the constraint on arbitrary power spectra ("Advanced LIGO O1") and the design sensitivity with the assumption of two years observations in Advanced LIGO and Virgo ("Advanced LIGO design") reported in Ref. [124]. These lines imply that GW signals above (below) them correspond to SNR ≥ 2 (SNR ≤ 2). For other interferometers, we assume one year cross-correlation searches and plot the lines with SNR = 2. The sensitivity curve for ET is produced by using a fitting function in [137]. The sensitivity of eLISA depends on the detailed detector configurations, and here we assume the C1 configuration, whose parameters are specified in Refs. [128,129]. For instrumental noises of DECIGO and Ultimate DECIGO, we used the parameters specified in Ref. [138]. It should be noted that GWs produced from white-dwarf (WD) binaries may lead to a significant confusion noise, which decreases sensitivities at lower frequencies of f 0.1 Hz. In Figure 3, we adopt the fitting formula for the WD confusion noise specified in Refs. [139,140] in addition to instrumental noises of DECIGO and Ultimate DECIGO. The sensitivities of EPTA and SKA are taken from [141,142]. Advanced LIGO design (blue), and ET (purple). The sensitivity curves for DECIGO and Ultimate DECIGO contain both the instrumental noise and the WD confusion noise. Light colored regions represent typical spectra of GWs from domain walls for σ 1/3 = 10 5 GeV and T ann = 0.1 GeV (light red), σ 1/3 = 10 9 GeV and T ann = 10 4 GeV (light green), and σ 1/3 = 10 11 GeV and T ann = 10 8 GeV (light blue).
In Figure 3, we also plot the GW signatures from cosmic domain walls for three choices of parameters. In these plots, we used Eqs. (3.9) and (3.10) to estimate the peak amplitude and frequency. The spectra are extrapolated based on the frequency dependences implied by the results of numerical simulations, Ω gw ∝ f 3 for f < f peak and Ω gw ∝ f −1 for f > f peak . We see that sufficiently large GW signatures are predicted according to the values of σ and T ann . Following Ref. [61], in Figure 4 we specify the parameter region of T ann and σ 1/3 relevant to observations. The colored regions in Figure 4 correspond to the parameter values for which the peak amplitude of GWs from domain walls [Eq. (3.9)] exceeds the sensitivity curves plotted in Figure 3. In Figure 4, we also plot the parameter region denoted by "Wall domination", which corresponds to the potential uncertainties since the energy density of domain walls dominates the total energy density of the universe [see Eq. Up to the value of V 0 , this condition gives an upper limit on T ann . Here we take V 0 = σ 4/3 as a typical estimate of the height of the potential barrier. 9 The corresponding parameter region is denoted by "No domain walls" and shown in Figure 4.  Figure 3 and in Figure 4, we have assumed that the annihilation of domain walls occurs during the radiation dominated era. If it occurs before reheating, the amplitude and frequency of GWs are modified according to the value of the reheating temperature [see Eqs. (3.14) and (3.15)].

Both in
From Figure 4, we see that the significantly large GW amplitude is predicted if the energy density of domain walls is close to dominate the total energy density of the universe. It should be noted that recent observational results by EPTA and Advanced LIGO already exclude some parameter spaces. Future observations with improved sensitivities are expected to probe much wider ranges of the parameter space, and from such observations we will obtain richer information about high energy physics beyond the SM.

Conclusion
Various well-motivated particle physics models predict the formation of unstable domain walls in the early universe, and it is possible to probe such models by observing GWs produced by them. The signatures of GWs can be characterized by two quantities, the tension of domain walls σ and the temperature at the annihilation of them T ann . Values of these parameters depend on the details of models, and they range over many orders of magnitude. Accordingly, future broadband observations of GWs including PTA, ground-based, and space-borne interferometers will allow us to explore new physics at various energy scales, some of which cannot be reached in the conventional laboratory experiments. Assuming that the annihilation of domain walls occurs during the radiation dominated era, we have shown that the ranges of 10 −2 GeV T ann 10 9 GeV and 10 3 GeV σ 1/3 10 12 GeV can be covered by future experiments.
So far we have estimated the spectrum of GWs from domain walls based on a naive extrapolation of the results obtained in the field theoretic lattice simulations. However, it is difficult to estimate the spectrum of GWs accurately over broad frequency ranges due to the limitation of the dynamical ranges of the simulations. In order to resolve this difficulty, it will be necessary to develop some alternative method to compute the spectrum of GWs analytically. Such an approach would enable us to estimate the signatures of GWs more quantitatively, and it can be used to distinguish the signal of domain walls from that of other sources. Furthermore, the results of numerical simulations of domain walls associated with Z N symmetric models imply that the shape of the spectrum slightly differs depending on the value of N [47]. The N -dependent feature in the spectrum of GWs might be regarded as an additional information to distinguish different models, and it deserves further investigation.
The search for cosmological GWs would have a great impact on high energy physics and cosmology. The collapse of domain walls will provide a possible way to interpret the results of forthcoming GW experiments. Even if there is no evidence of a signal, such information can be used to constrain various particle physics models beyond the SM.