Gravitational wave signals from short-lived topological defects in the MSSM

Supersymmetric theories, including the minimal supersymmetric standard model, usually contain many scalar fields whose potentials are absent in the exact supersymmetric limit and within the renormalizable level. Since their potentials are vulnerable to the finite energy density of the Universe through supergravity effects, these flat directions have nontrivial dynamics in the early Universe. Recently, we have pointed out that a flat direction may have a positive Hubble induced mass term during inflation whereas a negative one after inflation. In this case, the flat direction stays at the origin of the potential during inflation and then obtain a large vacuum expectation value after inflation. After that, when the Hubble parameter decreases down to the mass of the flat direction, it starts to oscillate around the origin of the potential. In this paper, we investigate the dynamics of the flat direction with and without higher dimensional superpotentials and show that topological defects, such as cosmic strings and domain walls, form at the end of inflation and disappear at the beginning of oscillation of the flat direction. We numerically calculate their gravitational signals and find that the observation of gravitational signals would give us information of supersymmetric scale, the reheating temperature of the Universe, and higher dimensional operators.


Introduction
The discovery of the Higgs boson [1,2] makes it clear that the Standard Model has a hierarchy problem between the electroweak scale and the Planck scale. One of the most elegant solutions is a symmetry relating fermions and bosons, called a supersymmetry (SUSY), which results in cancellations of quadratic divergences in quantum collections of scalar field masses including Higgs mass. In addition, the gauge couplings are miraculously unified at the grand unification scale in the minimal SUSY Standard Model (MSSM). This may indicate that the Standard Model gauge interactions are unified as a single gauge interaction at that scale. In addition to these phenomenological successes in particle physics, SUSY theories have plentiful implications in the early Universe. SUSY theories, including the MSSM, usually contain many scalar fields whose potentials are absent in the exact SUSY limit and within the renormalizable level (see e.g., Ref. [3]). Since their potentials are vulnerable to the finite energy density of the Universe through supergravity effects, these flat directions have nontrivial dynamics in the early Universe. One famous example is the Affleck-Dine mechanism, in which the baryon asymmetry is generated through a rotational motion of a flat direction in field-value complex plane [4,5]. In that mechanism, so-called a negative Hubble induced mass term plays an important role for the flat direction to have an initial vacuum expectation value (VEV), which is closely related to the amount of the resulting baryon asymmetry. Flat directions can also be considered as inflaton because their potentials are flat. However, inflation models in supergravity are usually affected by Hubble induced mass terms unless we introduce a shift symmetry [6] or consider a class of D-term inflation [7,8].
In the literature, it is usually assumed that the coefficient of Hubble induced mass terms is constant during and after inflation [5]. However, in Ref. [9], we have pointed out a possibility that the coefficient of Hubble induced mass terms changes at the end of inflation. This is partly because in most of the SUSY inflationary models the energy density of the Universe is dominated by different sources between the era during inflation and the one after inflation. In particular, the coefficient of the Hubble induced mass term can be positive during inflation and negative after inflation. With such a Hubble induced mass term, the flat direction stays at the origin of the potential during inflation and then obtains a large VEV after inflation. 1 Since the flat direction obtains a large VEV after inflation due to the negative Hubble induced mass term, we need to consider higher dimensional operators to stabilize its VEV. In some cases, the higher dimensional operator has U (1) or Z n symmetry. The VEV of the flat direction breaks the symmetry spontaneously and results in formations of topological defects, such as cosmic strings and domain walls, after the end of inflation. Interestingly, a typical width of these topological defects is always of the order of (but smaller than) the Hubble length because the curvature of the potential for the flat direction is of the order of the Hubble parameter. It follows that the energy density of the topological defects never dominate that of the Universe, as shown in this paper. When the Hubble parameter decreases down to the soft mass of the flat direction, which we expect O(1) TeV, the potential for the flat direction is dominated by its soft mass term and it starts to oscillate around the origin of the potential. Since the symmetry is restored at the origin of the potential, the topological defects disappear at that time. Therefore, the topological defects exist only between the end of inflation and the beginning of the oscillation, so that we can not directly observe them at the present epoch. However, we have a chance to observe redshifted gravitational waves (GWs) emitted by the topological defects.
Stochastic GW signals are generated from the dynamics of cosmic strings [10][11][12][13] and domain walls [14][15][16]. In the present case, emission of GWs begins at the end of inflation and ceases when the Hubble parameter decreases down to the soft mass of the flat direction. Since the emitted GW signals have a peak at the wavenumber corresponding to the Hubble length, resulting GW signals have a peak at the wavenumber corresponding to the soft mass of the flat direction. This implies that in principle we can obtain the information of the SUSY scale from the observation of GW signals [9]. In addition, since the energy density of the topological defects is roughly proportional to the square of the VEV of the flat direction, which is related to its higher dimensional operators, its GW signal also has information of higher dimensional operators. Finally, since the GW spectrum is sensitive to what dominates the energy density of the Universe, we could obtain information of the reheating temperature from the observation of GW signals [17,18].
In this paper, we investigate the dynamics of the flat direction which obtains a positive Hubble induced mass term during inflation and obtains a negative one after inflation. We consider three cases for its higher dimensional operators to stabilize its VEV; the one coming from a nonrenormalizable superpotential, the one coming from a Kähler potential with U (1) symmetry, and the one coming from Kähler potentials with Z n symmetry. We numerically follow its dynamics and confirm the above qualitative understanding, i.e., topological defects, such as cosmic strings and domain walls, form after the end of inflation and disappear when the Hubble parameter decreases down to the soft mass of the flat direction. From the simulations, we also evaluate its GW signals and discuss its detectability by future GW detectors.
In the next section, we briefly review flat directions in SUSY theories. In Sec. 3, we explain the potentials for flat directions originating from SUSY breaking and higher dimensional operators. We consider three important cases for higher dimensional operators and show our results of numerical simulations for each case in Secs. 4, 5, and 6. Finally, Sec. 7 is devoted to the conclusion.

Flat directions in the MSSM
A large number of scalar fields are introduced in SUSY theories because SUSY relates fermions and bosons. The potentials for scalar fields are severely restricted by SUSY, which especially results in existence of many flat directions. Flat directions are scalar fields whose potentials are absent within the renormalizable level as long as SUSY is unbroken.
Here, we illustrates how the potentials for scalar fields are absent in SUSY theories by taking a flat direction called the u c d c d c flat direction as an example. Let us focus on a scalar field constructed by right-handed squarks through the following orthogonal matrix: where the lower and upper indices represent flavor and color, respectively. The dots represent other directions, which we are not interested in. Since the inverse matrix is given by the transposed matrix, we obtain In the MSSM, the superpotential is given by 2 within the renormalizable level. Here we have omitted flavor indices for notational simplicity.
It is easy to see that the F -term potential for the field φ is absent such as The D-term potential is also absent like where |D a 3 | 2 and |D 1 | 2 are D-term potentials for SU (3) and U (1) Y , respectively. Therefore, the field φ has a flat potential and is called a flat direction. The above example consists of right handed squarks of u c , d c , d c , and is called the u c d c d c flat direction. It is known that every flat direction is characterized by gauge-invariant monomial in this manner. Table 1 is the list of flat directions in the MSSM [3]. Note that there are many flat directions even in such a simple model.
Although flat directions have no potential within the exact SUSY limit and the renormalizable level, they obtain nonzero potentials through SUSY breaking and nonrenormalizable operators (i.e., underlying higher energy theory). These potentials induce non-trivial dynamics of flat directions. The next section is devoted to discussing this point.

Potentials for flat directions
In this section, we discuss the induced potentials for flat directions through SUSY breaking and nonrenormalizable operators. Flat directions have soft (SUSY breaking) masses of the order of sparticle masses, which are subject to collider experiments and should be larger than O(10 2-3 ) GeV [19,20]. In addition, the finite energy density of the Universe contributes to potentials for flat directions. For instance, scalar fields obtain so-called Hubble induced terms through supergravity effects during inflation because inflation is driven by a finite energy density [5]. This is also the case during the inflaton oscillation dominated era as we explain in Sec. 3.2. Finally, when a flat direction has a large VEV, nonrenormalizable operators become important. 3 In this paper, we consider three important cases for nonrenormalizable operators, which we specify in Sec. 3.3.

Soft terms
Since we have not discovered SUSY particles yet, SUSY particles have to be much heavier than the SM particles. This can be achieved by introducing a SUSY breaking hidden sector, which induces soft mass terms to SUSY particles through some mediation mechanism. The potential for a flat direction can be written like around the weak scale vacuum. Note that not only soft masses but µ-term can also contribute to m φ if the flat direction consists of Higgs field. 4 Here we neglect A terms because they are subdominant and irrelevant in the following discussion. 5 In this paper, we assume that weak scale SUSY is broken at latest earlier than the time at which the Hubble parameter H(t) decreases down to m φ so that the soft mass term of Eq. (3.1) dominates the potential for the flat direction after H(t) m φ .

Hubble induced terms
The potential for a flat direction φ is modified by supergravity effects during and after inflation. The energy density of the Universe is so large that we have to consider its effects 4 In this paper, hereafter, we neglect the flavor indices. This is valid when the Hubble induced terms and/or nonrenormalizable terms breaks flavor symmetry completely. The case with a flavor symmetry is considered in Appendix C. 5 There might be Hubble induced A terms during inflation, which plays important role in the context of the Affleck-Dine baryogenesis. In the present case, the flat direction stays at the origin during inflation, so that such higher dimensional terms are irrelevant. In addition, Hubble induced A terms are irrelevant after the end of inflation because the time-averaged F term F vanishes (see e.g., Ref. [24]). Therefore, in the present scenario, the flat direction can not generate the baryon asymmetry. on the potential for the flat direction. In the supergravity, scalar potentials are written in terms of superpotential, W , and Kähler potential, K, which are functions of scalar fields ψ i specified below. The potential for scalar fields is given by Hereafter, we assume F -term inflation models (see Appendix A for example), in which we can neglect the contribution from D-term potentials. The Kähler potential also determines kinetic terms such as One may assume a minimal Kähler potential of K = ψ i 2 , which results in the canonical kinetic term. However, since we do not have any knowledge about the physics above a cutoff scale, such as the Planck scale, there are expected to be cutoff suppressed terms in the Kähler potential (see Eq. (3.4)).
The relevant fields are a flat direction φ and two fields in the inflaton sector denoted as I and X, so that ψ i are identified as φ, I, and X in Eqs. (3.2) and (3.3). The superfield I contains a field χ which starts to oscillate after inflation and whose energy dominates the energy density of the Universe, that is, χ is identified as inflaton in chaotic inflation models [6] and as water-fall fields in hybrid inflation models [25,26]. The minimal Kähler potential for the field I is given by K min (I) = (I + I * ) 2 /2 for the case of the chaotic inflation model proposed in Ref. [6] and K min (I) = |I| 2 for the case of hybrid inflation models (see Appendix A). The superfield X represents the field whose F -term potential energy drives inflation. In other words, the F term of X satisfies the relation of |∂W inf /∂X| 2 = 3H 2 inf M 2 Pl . In Appendix A, we briefly explain chaotic and hybrid inflation models and identify the fields I and X in those models.
To see how the potential for the flat direction is modified, let us consider the following non-minimal Kähelr potential: where M * is a cutoff scale, and c a (a = X, I) are constants of order one. We expect that M * is less than or equal to the Planck scale M Pl ( 2.4 × 10 18 GeV). We set M * = M Pl in Sec. 4, while we assume M * to be less than M Pl in Secs. 5 and 6. Since the F term of X satisfies the relation of |∂W/∂X| 2 = 3H 2 inf M 2 Pl , the flat direction obtains the following mass, called the Hubble induced mass, during inflation: Let us note that c H inf is of the order of M 2 Pl /M 2 * when c X is of order one. While the coefficient c H inf is assumed to be negative in the context of the Affleck-Dine baryogenesis [4,5], we assume it to be positive in this paper. In this case, the flat direction stays at the origin, i.e., φ = 0, during inflation.
After inflation ends and before reheating completes, the energy density of the Universe is dominated by that of the oscillation of a scalar field (denoted by χ) and the Hubble parameter decreases with time as H(t) ∝ a −3/2 (t), where a(t) is the scale factor. The field χ oscillates with the frequency of its mass m χ . The oscillation amplitude is redshifted as where χ i is the oscillation amplitude just after the oscillation starts around t = t i . Since the oscillation time scale of χ is much smaller than H −1 , that is, H m χ , we can averageχ 2 and χ 2 over the oscillation time scale like The flat direction obtains a Hubble induced mass through higher dimensional kinetic terms as well as the F -term potential of Eq. (3.2) during this oscillation dominated era. From Eq. (3.3), we can see that the relevant term comes from which, together with the one coming from the F -term potentials, results in the Hubble induced mass term of c H H 2 (t) |φ| 2 with during the oscillation dominated era. Let us note that c Hosc is of the order of M 2 Pl /M 2 * when c X and c I are of order one.
As we can see from Eqs. (3.6) and (3.10), the coefficient of the Hubble induced mass term during inflation, c H inf , is generally different from the one during the oscillation dominated era, c Hosc [9]. This is because the field X, whose F term drives inflation, is generally different from the field I, whose oscillation energy dominates the energy density of the Universe after inflation. In Appendix A, by taking the models of chaotic inflation and hybrid inflation in supergravity as examples, we see that the field X is indeed different from the field I. While the nonrenormalizable terms are heavily dependent on the high energy physics beyond the cutoff scale M * , we assume that among many flat directions in SUSY theories (at least) one flat direction has c H inf > 0 and c Hosc < 0. In this case, the flat direction stays at the origin (φ = 0) during inflation, and then, it starts to roll down to a large VEV due to the negative curvature after inflation.
Here we comment on how natural to consider the above flat direction. There is no guiding principle to determine the signs of Hubble induced mass terms unless we specify the underlying high energy model beyond cutoff scale. This is what the well-known Affleck-Dine mechanism is also based on, while it assumes a negative coefficient during inflation. As shown above, the coefficient during inflation and that after inflation do not have the common origin and hence can be different in general. We consider one out of four possible combinations of the signs of the two coefficients. It is natural that one or more flat directions have such a combination because there are a lot of flat directions in the MSSM (see Table 1). This is why our work can find its application in not all but wide class of supersymmetric models.

Higher dimensional terms
Since the flat direction obtains a large VEV due to the negative Hubble induced mass term after inflation, we need to consider its higher dimensional terms to stabilize the VEV of the flat direction. Depending on symmetries that the flat direction possesses, three important cases are considered in this paper.
In the next section, we consider the case that the flat direction has a nonrenormalizable superpotential like where n (≥ 4) is an integer dependent on each flat direction and λ is an positive integer. Let us stress that among possible higher dimensional terms, only the lowest dimensional one is relevant to the dynamics of the flat direction. In the case of the u c d c d c flat direction, for instance, it may come from a term like (see Eq. (2.2)) In this case, the power of superpotential n is 6. Note that n = 3 superpotential is forbidden due to R-parity. When we extend this symmetry to a discrete R-symmetry, n can be 6, 9, 12, . . . for the u c d c d c flat direction. For each flat direction, possible values of n are shown in Table 1 [3]. The F -term potential for the flat direction comes from the superpotential and is given by (see Eq. (3.2)) where we neglect irrelevant higher dimensional terms. Note that the potential of Eq. (3.13) possesses a global U (1) symmetry. This is related to the baryon charge (or the lepton charge or a combination between them), which plays an important role in the context of the Affleck-Dine baryogenesis. Let us emphasize that the symmetry is a global symmetry independent from the local symmetries of the flat direction. The symmetry can be explicitly described as (see Eq. for the case of the u c d c d c flat direction. One can see that this transformation can not be described by the gauge symmetries of SU (3) × U (1) Y , so that the U (1) symmetry is an additional global symmetry. If λ in Eq. (3.11) is sufficiently small or superpotential is absent due to some symmetry, the VEV of the flat direction can be as large as the cutoff scale M * . For example, if we introduce B − L symmetry, it forbids superpotentials for the u c d c d c flat direction like Eq. (3.12). If we assign R-charges as, for example, R(Q) = R(L) = R(e c ) = R(u c ) = R(d c ) = 1 and R(H u ) = R(H d ) = 0, it also forbids superpotentials for the u c d c d c flat direction. 6 In such cases, the potential for the flat direction may be dominated by higher dimensional terms coming from nonrenormalizable Kähler potential [27]. The induced potential is proportional to H 2 like the Hubble induced mass term described above. In Sec. 5, we consider a flat direction with a potential given by where m (≥ 3) is a certain integer. 7 The coefficient a H is given by where we assume a H is positive and of order one. Note that the potential of Eq. (3.15) respects U (1) symmetry.
On the other hand, if we assign R-charges as, for example, R(Q) = R(L) = R(e c ) = R(u c ) = R(d c ) = 0 and R(H u ) = R(H d ) = 2, nonrenormalizable Kähler potentials may not respect U (1) symmetry. In Sec. 6, we consider the same potential with Eq. (3.15) except for an additional U (1) breaking term coming from higher dimensional Kähler potentials: where m = 2, 3, 4, . . . , is an integer. We assume m < 2m − 2 so that the potential for the flat direction is bounded from below. The coefficient b H is given by where we assume b H = O(1). Here we have redefined the field φ such that b H is postitive. This additional term explicitly breaks the U (1) symmetry, but there remains Z m symmetry in that potential.

Summary of this section
Here we summarize the potential for the flat direction. It is given by The R-symmetry has to be broken to obtain nonzero gaugino and higgsino masses and to set the vacuum energy (almost) zero. Its symmetry breaking order parameter is proportional to the SUSY breaking scale, so that R-symmetry breaking terms are suppressed by a factor of m 2 φ /M 2 * . These terms may be comparable to the other terms, such as the soft mass term, around the time when the Hubble parameter decreases down to the soft mass scale. This implies that there is O(1) uncertainty in our result when we omit those R-symmetry breaking terms. On the other hand, it might be possible that R-symmetry breaking terms are suppressed by a factor of m 2 φ /M 2 Pl . In this case, one can easily find that they can be neglected even when the Hubble parameter decreases down to the soft mass scale. After that, higher dimensional terms are irrelevant because the potential for the flat direction is dominated by the soft mass term. Therefore, it is reasonable to consider that R-symmetry may suppress the superpotential for the flat direction. 7 In general, two or more terms can give comparable contributions to the potential for flat direction. While it can change our results by a factor of O(1), in this paper, we consider a nonrenormalizable potential with a single term for simplicity.
where the higher dimensional terms V NR is given by Eqs. (3.13), (3.15), or (3.17). Hereafter, c H inf and c Hosc are collectively denoted as c H . We consider a flat direction which has a positive Hubble induced mass term (c H = c H inf > 0) during inflation and a negative one (c H = c Hosc < 0) after inflation.
In the following sections, we investigate the dynamics of the flat direction and show that topological defects, such as cosmic strings and domain walls, form after the end of inflation. We numerically follow the dynamics of the flat direction and calculate its GW signal. We also discuss the detectability of GW signals and explain that the information of the mass of the flat direction, the reheating temperature of the Universe, and higher dimensional potentials are imprinted on GW signals. In Secs. 4, 5, and 6, we consider the case that the flat direction has a higher dimensional potential of Eqs.  In this section, we consider the case that the flat direction has a nonrenormalizable superpotential of Eq. (3.11). The F -term potential for the flat direction is given by where n ≥ 4.
In the next subsection, we briefly explain the dynamics of the flat direction in the early Universe. We numerically follow its dynamics and calculate its GW signal, and the results are shown in Sec. 4.2. Then we give physical interpretation of the results in Sec. 4.3. Finally, we discuss the detectability of GW signals in Sec. 4.4.

Dynamics of the flat direction
The equation of motion for the flat direction is given bÿ with dots denoting derivatives with respect to t. The potential for the flat direction V (φ) is written as Note that we collectively denote c H inf and c Hosc as c H . As explained in the previous section, we consider a flat direction which has a positive Hubble induced mass term (c H inf > 0) during inflation and a negative one (c Hosc < 0) after inflation. During inflation, the flat direction stays at the origin, φ = 0, due to the positive Hubble induced mass term. Since the coefficient c H becomes negative after the end of inflation, the flat direction rolls down to the stable point and becomes to obtain a large VEV. 8 The VEV is given by

4)
8 One might wonder if the flat direction is affected by thermal effects and is forced to stay at the origin by the effective thermal mass for the flat direction [5,[21][22][23]. However, we can neglect thermal effects on the flat direction by the following discussion. During the inflaton oscillation dominated era, the temperature of the plasma is given by T ∼ (H(t)ΓI M 2 Pl ) 1/4 , where ΓI is the decay rate of the inflaton. The flat direction as long as |c Hosc | H 2 (t) m 2 φ . The potential of Eq. (4.1) has a global U (1) symmetry, which corresponds to the baryon charge (or the lepton charge or a combination of them) in the context of the Affleck-Dine baryogenesis. After the end of inflation, the non-zero VEV of the flat direction breaks the global U (1) symmetry. 9 Since the flat direction stays at the origin during inflation, cosmic string network forms after this phase transition. The energy of cosmic strings per unit length µ is roughly given by depending on dimensionless parameters only logarithmically [28]. By using numerical simulations (see the next subsection), we confirm that cosmic string network reaches a scaling regime well before the string disappears. In this regime, the number of cosmic strings in the Hubble volume is O(1) and the energy density of cosmic strings ρ cs is given by A typical width of cosmic strings is roughly given by the curvature of the potential, which is of the order of the Hubble length in the present case. 10 This means that the width of cosmic strings increases with time. This is one of the outstanding characteristics of cosmic strings in this scenario compared with other cosmic strings considered in the literature.
Since the Hubble induced mass squared decreases with time as ∝ a −3 (t) after the inflation ends, the soft mass term m 2 φ |φ| 2 eventually dominates the potential for the flat direction and the flat direction starts to oscillate around φ = 0. This implies that cosmic strings disappear at the time of |c Hosc | H 2 (t) m 2 φ , which can be rewritten as With m φ = O(1) TeV, cosmic strings disappear much earlier than the bing bang nucleosynthesis (BBN) epoch, so that observational constraints on long-lived topological defects are not applicable to our scenario. The flat direction decays into radiation through gauge interactions soon after cosmic strings disappear. We define the reheating temperature T RH as the temperature when the energy density of the Universe is dominated by that of radiation. The Hubble parameter at that time is given by obtains a thermal mass of the order of the temperature. The ratio of the thermal mass to the Hubble induced mass is thus roughly given by . This is smaller than unity just after inflation if we consider a Planck suppressed decay of inflaton. Thus, the flat direction obtains a large VEV soon after the end of inflation due to the negative Hubble induced mass term. Once the flat direction obtains a large VEV, thermal effects become suppressed and can be neglected because particles in thermal bass obtain large effective masses [22,23]. Therefore, we can neglect thermal effects on the flat direction. 9 It also breaks SM gauge symmetries because flat directions have nonzero SM gauge charges. In the unitary gauge, only gauge fields may obtain a nontrivial configuration. Since the broken gauge field obtains the effective mass of g φ , the gauge field configuration begins to oscillate and becomes trivial (A µ = 0) soon after the phase transition. Therefore, we focus only on the global symmetry of the flat direction. See also the discussion in Appendix C.1.
10 Since a typical width of cosmic strings is of the order of the Hubble radius, the Nambu-Goto approximation, which was used in Ref. [29] for example, is inappropriate to describe these cosmic strings.
where g * is the effective relativistic degrees of freedom for the energy density. Note that in lowscale SUSY models we should require T RH 10 9 GeV to avoid the gravitino problem [30,31]. Even for high-scale SUSY models such as pure gravity mediation, T RH 10 10 GeV is required to avoid an overproduction of LSP (= O(100) GeV) [32][33][34]. In such well-motivated cases, the cosmic strings disappear ( Hereafter we consider such a case.
After they form at the end of inflation and before they disappear at t t decay , cosmic strings emit GWs. We numerically follow the evolution of the flat direction and calculated its GW signal. The results are shown in the next subsection.

Results of numerical simulations
In this subsection, we show results of our numerical simulations for the dynamics of the flat direction and its GW signal. Since we consider the evolution of cosmic strings during the oscillation dominated era (H m χ ), the Hubble parameter decreases with time as By redefining the variables, we rewrite the equation of motion such as Let us emphasize that the parameter λ can be absorbed by the redefinition of φ, so that our numerical simulation can be applied to any cases with an arbitrary value of λ. The variable τ is a dimensionless conformal time and is related with comoving horizon size τ through τ = τ a i H i . Primes denote derivatives with respect to τ . We rewrite important parameters by the conformal time and obtain τ i = 2, (4.14) where we implicitly assume a i ≡ a( τ i ) = 1. The subscript i represents the values at the time of τ i at which we set initial conditions for numerical simulations.
We have performed three-dimensional lattice simulations to solve the classical equation of motion Eq. (4.9) numerically. We use a numerical method similar to the one used in Ref. [16]. They use the 4-th order symplectic integration method to evolve the equation of motion for the flat direction in the expanding background, while we have used the 6-th order symplectic integration method. We ignore the backreaction of metric perturbations on the field evolution, which is negligible for ρ cs H 2 (t)M 2 Pl , that is, for φ M Pl (see Eqs. (4.5) and (4.6)). In the set of simulations, we take a unit of H i ≡ H( τ i ) = 1 and set m φ /H i = 5 × 10 −4 and c Hosc = −15. The final time τ f and the time step ∆ τ is set to be 50 (100) and 0.2 (0.1), for N = 128 3 (256 3 ) grid points, respectively. The comoving box size is L = 100 (200) with N = 128 3 (256 3 ) grid points, which implies that the comoving grid size is ∆x = 100/128(= 200/256) 0.75. Since a typical width of cosmic strings is of the order of the Hubble radius H −1 ( τ ) = ( τ / τ i ) 3 , it is correctly resolved in our numerical simulation. We use the periodicity boundary condition for spatial directions. We require that the horizon scale is smaller than the simulation size (i.e., L/2 ≥ τ ) to avoid the boundary effect on the cosmic string network. The initial condition for simulations is seeded by quantum fluctuations around the origin of the configuration space. To be concrete, we set initial conditions in the following way. First, we draw two samples f k,1 , f k,2 from a normal distribution with variance of for each discretized wavenumber vector k = (2π/L) n, where n is a triplet of integers from −N/2 to N/2 − 1. In the above expression, we set λ such that |φ i | = 10 2 in units of H i . Let us stress that once the cosmic string network follows the scaling evolution, the initial condition (i.e., λ) becomes irrelevant. We generate initial conditions ϕ i , ϕ i by taking a superposition of (discrete) Fourier components given by (4.20) Figure 1 shows one example of cosmic string network evolution obtained by our numerical simulation. The blue region represents the inside of cosmic strings, which we define by |φ(x)| < |φ| /5. One can find that a typical width of cosmic strings increases with time and is roughly proportional to the Hubble length. The black line at the bottom-left corner in Fig. 1 describes a unit of the comoving horizon size (given by τ ) at each time. One can see that each Hubble volume contains O(1) cosmic strings, which means that cosmic string network is in a scaling regime. From our numerical simulation, we find that the spatially averaged magnitude of φ grows with time and eventually reaches of the order of the VEV (Eq. (4.4)) around τ form 10 for the case of c Hosc = −15.
From the numerical simulation for cosmic string dynamics, we can obtain the transversetraceless part of the anisotropic stress by the relation of where the projection operator in Fourier space is given by with k = |k|k and Kronecker's delta δ ij . Given τ RH , we can calculate the energy density of GWs from Eqs. (B.6), (B.7), and (B.14) in Appendix B. We define the energy spectrum of GWs ρ gw such as where ρ tot (τ ) (= 3M 2 Pl H 2 (τ )) is the total energy density of the Universe. Figure 2 shows evolution of GW spectra obtained from our numerical simulations. The GW energy density grows as a typical value of the flat direction increases until it reaches the potential minimum at τ = τ form 10. Then the GW peak energy density decreases with time as Ω gw ∝ τ −12/(n−2) and the peak wavenumber k peak decreases with time as k peak ∝ τ −1 for the case of n ≥ 8. On the other hand, the GW peak energy density decreases with time as Ω gw ∝ τ −2 and its peak wavenumber is constant in time for the case of n = 6. We give physical interpretation of these results in the next subsection.
Taking a non-zero value of m φ into account, we find that the flat direction starts to oscillate around the origin of the potential and cosmic strings disappear eventually around τ = τ decay 40, which is given by Eq. (4.16) with m φ /H i = 5 × 10 −4 and c Hosc = −15. After that, the GW spectrum is freezed out, which means that the GW peak wavenumber becomes constant and the GW spectrum decrease adiabatically as ∝ τ −2 (∝ a −1 ). Note that we have taken m φ = 5 × 10 −4 above just for computational time saving, though it does not change the dynamics of flat direction until τ τ decay . Since H i is roughly given by m χ χ i /M Pl (see Eq. (3.8)), m φ can be much smaller in general. When we take m χ = 10 13 GeV and χ i = M Pl motivated by chaotic inflation model and m φ / |c Hosc | 1/2 = 10 3 GeV, τ decay 4000 is larger than τ form 10 by three orders of magnitude.

Physical interpretations
In this subsection, we discuss how we can understand the results of our numerical simulations shown in Fig. 2.
During the scaling regime, the typical length scale of cosmic string dynamics is the Hubble length. As a result, the emission peak wavenumber of GWs emitted from these cosmic strings is given by the Hubble scale: We can estimate the energy density of GWs in the following way. Its quadrupole moment Q can be roughly given by [35] Q where µH −1 is the total energy of cosmic strings within a Hubble volume. Such cosmic strings emit GWs with the luminosity of where we replace the time derivative with H in the last line. Thus, the produced energy density of GWs can be estimated as where we use Eqs. (4.5) and (4.4) in the second and last equality, respectively. During each Hubble time, cosmic strings emit GWs with the emission peak wavenumber of Eq. (4.25) and with the produced energy density of Eq. (4.29). We have measured ∆Ω gw /∆ log τ from our lattice simulation as shown in Fig. 3. We determine the numerical prefactors for the emission peak wavenumber and produced energy density of GWs from the numerical results such as . (4.31) We find that the numerical prefactors are independent of n and τ (see Fig. 3), so that Eqs. (4.30) and (4.31) can be used for any n and τ . From Fig. 3, we can find that ∆Ω gw /∆ log τ is proportional to a certain power of k for large and small wavenumber modes: where α 2. We also find that these values are almost independent of n and τ from Fig. 3. The GW spectra for large wavenumber modes in Fig. 2 is expected to be given by integrating the GW peak energy density from τ form to τ , such as The scale factor is included in the integrant because the GW energy density decreases with time as ∝ a −1 (τ ) ∝ τ −2 in the oscillation dominated era. We define τ k by the time the GW emission peak wavenumber reaches k, that is, k ∆ peak (τ k ) ≡ k. We implicitly assume τ > τ k in this calculation, which means that Eq. (4.33) is valid for k > k ∆ peak (τ ). The result of Eq. (4.33) at τ = 50 is plotted as the black curves in Fig. 2 and well describes the GW spectra. Note that the GW energy density is roughly proportional to and n ≥ 6. Thus, in principle, we can obtain the value of n from the observation of GWs with high frequencies. Note that for the case of n = 6 and 7, the power-law index of GW spectra is positive (16 − 2n > 0). This implies that the GW spectrum has a peak at τ = τ form . This is because the emitted GW energy density of Eq. (4.31) decreases with time faster than the redshift ∝ a −1 for n = 6 and 7. This implies that the peak energy density emitted at τ < τ form is smaller than the redshifted one emitted at τ = τ form when the spatially averaged magnitude of the flat direction reaches the minimum of the potential. Thus, the GW peak wavenumber is determined by k peak k ∆ peak (τ form ) ∼ τ −1 form . Then the GW peak energy density decreases with ∝ a −1 due to the redshift until reheating completes. The redshift factor is given by a form /a RH = (H RH /H form ) 2/3 . This is the reason that the GW peak wavenumber is constant in time for the case of n = 6 in Fig. 2. The GW spectrum for n = 8 is too flat to determine its peak wavenumber precisely, though the GW spectra are well fitted by the above estimation.
Next, let us consider the GW spectrum for very small wavenumber modes k k peak . For such super-horizon modes, we expect that the leading contribution for the Fourier transformed transverse-traceless part of the anisotropic stress T TT ij is independent of k due to the loss of causality at the large scale [14,36]. As a result, the GW spectrum is proportional to k 2 for the modes smaller than but around k peak , which is in agreement with our numerical simulations in Fig. 3. It is proportional to k for the modes much smaller than k peak , which are out of reach of our numerical simulations. In addition, we have to take into account the effect of the reheating, at which the expansion law of the Universe changes. These effects on super-horizon modes leave their imprints on GW signals at present, which we discuss in the next subsection. We can calculate the GW spectrum for super-horizon modes by substituting a constant T TT ij into Eqs. (B.6) and (B.7) in Appendix B and plot the results at τ = 50 as the gray dotted curves in Fig. 2.
From Fig. 2, we can see that GW spectra drop off for modes k τ −1 form ( k peak ), which enter the horizon before the spatially averaged magnitude of the flat direction reaches VEV |φ| . This is because the produced energy density of GWs is proportional to the fourth power of the field value of the flat direction as shown in Eq. (4.28). The result depends nontrivially on c Hosc because the dynamics of the flat direction around the origin of the potential is mainly determined by the Hubble induced mass term. Let us note that we are less interested in such high frequency modes in light of the detection of GW signals (see Fig. 5).
The GW emission terminates around τ decay , which is defined by |c Hosc | H 2 (t) = m 2 φ (see Eq. (4.16)). This is because the flat direction starts to oscillate around the origin of the potential and cosmic strings disappear at that time. This implies that the GW spectrum has a peak (or at least bends) at the wavenumber corresponding to the soft mass of the flat direction. In fact, the GW peak energy density and frequency at that time is given by for the case of n ≥ 8. Here, we determine the numerical prefactors from the results of our numerical simulations. Although GW spectrum has a peak at k τ −1 form for the case of n = 6 and 7, it bends at the wavenumber of Eq. (4.35) and has the energy density of Eq. (4.34). In either case, the GW spectrum contains information of the soft mass of the flat direction at that wavenumber because of the relation of Eq. (4.35). Hereafter, we denote the wavenumber defined by Eq. (4.35) as the peak wavenumber k peak . The shape of GW spectrum is in fact freezed out around τ = τ decay 40, which is given by Eq. (4.16) with m φ /H i = 5 × 10 −4 and c Hosc = −15. After that, the peak wavenumber does not change and the GW energy density decreases as ∝ a −1 until reheating completes. Figure 4 is a schematic diagram for the resulting GW spectrum emitted from the cosmic strings. We calculate super-horizon modes, i.e., small wavenumber modes by using Eqs. (B.6) and (B.7) with a constant T TT ij . The result is plotted as the blue dashed curve. The GW spectrum for small wavenumber modes near the peak wavenumber is proportional to k 2 , while that for the smaller wavenumber modes is proportional to k. We also take into account the effect of reheating, so that the GW spectrum bends at the wavenumber corresponding to the inverse of the horizon size at the end of reheating [17], We explain this phenomena in the next subsection in detail. We calculate large wavenumber modes by using Eq. (4.33) for the cases of n = 6, 8, 10, and 12, and the results are plotted as the black curves. One can see that the GW spectrum for high frequency modes depends on the value of n. For n ≥ 9, the GW spectrum has a peak at the wavenumber corresponding to the inverse of the horizon scale when cosmic strings disappear, that is, at the time of H(t) m φ / |c Hosc |. For n = 6, 7, and 8, the GW spectrum does not have a peak, but bends at that wavenumber. Note that the magnitude of GW energy density depends on the value of the coupling in the superpotential λ, though we use it to rescale the vertical axis in Fig. 4.  We take into account a nonzero value of reheating temperature, which modifies the GW spectrum for wavenumber smaller than the inverse of the horizon scale at the end of reheating.

GW signals at present
The above calculations and discussions focus on the GW spectrum during the oscillation dominated era. In this subsection, we derive the GW spectrum at the present epoch, taking into account the redshift.
The energy density of GWs at τ decay , Ω gw (τ decay ) is given by our numerical simulations as shown and discussed in the previous two sections. Note that we consider the case that cosmic strings disappear before reheating completes, as we explained in Sec. 4.1. Since the GW energy density Ω gw decreases with time as a −1 during the oscillation dominated era (i.e., the matter dominated era), its present value is given by where t 0 is the present time, Ω r h 2 ( 4.15 × 10 −5 ) is the present energy density of radiation, and g * (g s ) is the effective relativistic degrees of freedom for the energy (entropy) density. In the following calculations, we set g s (t 0 ) = 43/11, g * (t 0 ) = 2 + 21/11(4/11) 1/3 , and g s (t RH ) = g * (t RH ) = 915/4 in the following calculations. Substituting [Ω gw (τ decay )] peak into Eq. (4.37), we obtain where we use Eqs. (4.7), (4.8), and (4.34). The GW peak frequency at present f 0 is given by . We should emphasize that this peak frequency depends only on |c Hosc | −1/2 m φ ( H(τ decay )) and T RH and is independent of n and λ.
Here, let us consider the effect of reheating on the GW spectrum, which we mentioned in the previous section. In the absence of GW source, the metric perturbation is constant for super-horizon modes and decreases with a −1 for sub-horizon modes. Let us define a k as the scale factor at which GW with wavenumber k enters the horizon: k/a k = H(a k ). For H(a) ∝ a −ν , a k scales like a k ∝ k −1/(ν−1) . Noting that Ω gw is proportional to the metric perturbation squared, we obtain Ω gw ∝ k −2/(ν−1) Ω p gw for sub-horizon modes, which have primordial GW spectrum of Ω p gw before entering the horizon. If the primordial GW source is scale invariant, i.e., T T T ij = const. with respect to k, the primordial GW spectrum scales like Ω p gw ∝ k 5 . Since the expansion law of the Universe changes from H ∝ a −3/2 to H ∝ a −2 , we expect Ω gw ∝ k and ∝ k 3 for GWs which enter the horizon before and after reheating, respectively, if the primordial GW source is scale invariant. As a result, the GW spectrum is expected to bend at the wavenumber corresponding to the horizon scale at the end of reheating [9,17,18]. Since we consider the case that cosmic strings disappear before reheating completes, the relevant GW with wavenumber corresponding to the horizon scale at the end of reheating is super-horizon mode when the GW emission terminates (τ τ decay ). Thus, we can use T TT ij = const. in Eqs. (B.6) and (B.7) to estimate the GW spectrum due to the loss of causality for super-horizon scales [14,36]. Then, using j l (x) → 2 l l!x l /(2l + 1)! and n l (x) → −(2l)!/(2 l l!x l+1 ) for x 1 in Eq. (B.8), we obtain Ω gw ∝ k for τ −1 and Ω gw ∝ k 3 for k τ −1 RH as expected. This means that the GW spectrum actually bends at the wavenumber around k τ −1 RH . We denote the frequency corresponding to that wavenumber at present as bend frequency f bend (see Eq. (4.36). It is calculated as Note that this bend frequency is different from the one corresponding to k peak for n = 6 and 7. In these cases, GW spectrum bends at two frequencies: at f = f 0 (k = k peak ) corresponding to the horizon scale at τ = τ decay and at f = f bend (k = k bend ) corresponding to the horizon scale at the end of reheating.
Here, let us explain what kind of information we can obtain through the observation of GW spectrum emitted from cosmic strings. First, we can obtain the value of reheating  Figure 5. GW spectra generated by cosmic strings (dash and dot-dashed curves) and sensitivities of planned interferometric detectors (solid curves). We assume m φ = 10 2 GeV with T RH = 10 7 GeV (blue dashed curve), m φ = 10 2 GeV with T RH = 10 9 GeV (red dashed curve), and m φ = 10 3 GeV with T RH = 10 9 GeV (red dot-dashed curve). We also assume c Hosc = 10 and n = 10. We set the value of λ so that the GW peak energy density is given by 0.01 at τ = τ decay . temperature through detection of the bend frequency of Eq. (4.41). Then we can obtain the soft mass scale of the flat direction through Eq. (4.40). We can also obtain the value of the power n because the slope of the GW spectrum for wavenumber larger than Eq. (4.40) is dependent on n (see Eq. (4.33) or Fig. 4). Finally, we can obtain the value of λ through the measurement of the GW energy density Eq. (4.38). Thus, the observation of the whole GW spectrum gives us the values of |c Hosc | −1/2 m φ , T RH , n, and |c Hosc | −1/4 λ . Figure 5 shows examples of GW spectrum predicted by the present mechanism. The GW peak energy density and peak frequency are given by Eqs. (4.38) and (4.40), which are obtained from our numerical results. While the shape of the GW spectrum for low frequency modes are calculated from Eqs. (B.6) and (B.7) with a constant T TT ij and τ f = τ decay , that for high frequncy modes are calculated from Eq. (4.33). Note that the high frecuency modes corresponding to the horizon scale at the time of cosmic string formation τ form are out of the range of Fig. 5, because the energy scale of inflation is much larger than the mass scale of the flat direction. We take m φ = 10 2 GeV with T RH = 10 7 GeV (blue dashed curve), m φ = 10 2 GeV with T RH = 10 9 GeV (red dashed curve), and m φ = 10 3 GeV with T RH = 10 9 GeV (red dot-dashed curve).
To discuss the detectability of GW signals, we plot single detector sensitivities for LISA [37] and Ultimate DECIGO [38] by using the online sensitivity curve generator in [39] with the parameters in Table 7 of Ref. [40]. Note that in some parameter space we can obtain T RH by using only one single detector for Ultimate DECIGO. One may wonder why our sensitivity curve is different from that in Ref. [18], where they claim that Ultimate DECIGO can reach well below Ω gw ∼ 10 −16 for frequencies of 0.1 − 1 Hz. This is because they assume taking two detectors correlation for a decade to obtain T RH by observation of the inflationary GW background, while we assume a single detector. A single detector is sufficient for our purpose because GW signals from cosmic strings are much stronger than inflationary ones. We also plot cross-correlation sensitivities for Advanced LIGO [41] and ET (ET-B configu-ration) [42], assuming two detectors are co-aligned and coincident. In Fig. 5, we take the signal to noise ratio SNR = 5, the angular efficiency factor F = 2/5, the total observation time T = 1 yr, and the frequency resolution ∆f /f = 0.1. CMB constraints (horizontal lines) are put on the integrated energy density of GWs d log f Ω gw h 2 (t 0 ) [43]. We find that GW signals would be observed by ET or Ultimate DECIGO.
Note that we assume a very small λ in Fig. 5 so that the future GW detector can observe GW signals. However, such a small λ may result in a VEV (Eq. (4.4)) larger than the cutoff scale just after inflation, even though the VEV is smaller than the cutoff scale around the time of H(t) m φ . In this case, the potential for the flat direction should be dominated by higher dimensional superpotential or Kähler potentials to stabilize its VEV at a smaller scale than the cutoff scale for a while. As the Hubble parameter decreases, only the lowest dimensional nonrenormalizable superpotential becomes relevant. Then the GW signals are determined by the discussion in this section. Finally, the potential is dominated by the soft mass term and the GW emission terminates. We should emphasize that the power-law index of GW spectrum for high frequency modes around the peak frequency depends on the value of n as shown in Fig. 4. Thus, for larger wavenumber modes, the GW spectrum bends at the wavenumber corresponding to the horizon scale when the potential becomes dominated by the lowest dimensional nonrenormalizable superpotential. However, such a high frequency modes are beyond the reach of future GW detecters, so that we neglect such effects in Fig. 5.

Case B: W = 0 with U (1) symmetry
In this section, we consider the case that the potential for the flat direction is given by a higher dimensional term of Eq. (3.15), which comes from a nonrenormalizable Kähler potential. It is given by where m ≥ 3. The coefficient a H is given by where we assume a H = O(1).
The following discussion and calculation in this section are essentially the same as the ones given in the previous section. In the next subsection, we briefly explain the dynamics of the flat direction in the early Universe. We numerically follow its dynamics and calculate its GW signal, and the results are shown in Sec. 5.2. Then we give physical interpretation of the results in Sec. 5.3. Finally, we discuss the detectability of GW signals in Sec. 5.4.

Dynamics of the flat direction
The whole potential for the flat direction V (φ) is given by where we assume that c H is positive during inflation (c H inf > 0) and is negative (c Hosc < 0) after inflation. The positive c H inf makes the flat direction stay at the origin during inflation.
Then, after inflation ends, the negative c Hosc makes it obtain a large VEV. The potential takes the minimum at |φ| = |c Hosc | a H (m − 1)  (Eq. (3.16)). The VEV is as large as the cutoff scale M * . Note that the VEV of the flat direction is independent of the Hubble parameter and is constant in time. This is an important difference from the previous section (see Eq. (4.4)).
Since the above potential possesses a global U (1) symmetry, cosmic strings form after the end of inflation. By using numerical simulations (see the next subsection), we confirm that cosmic string network reaches a scaling regime well before cosmic strings disappear. In this regime, the number of cosmic strings in the Hubble volume is O(1) and the energy density of cosmic strings ρ cs is given by Eq. (4.6). The energy density of cosmic strings per unit length µ is given by Eq. (4.5). A typical width of cosmic strings is roughly given by the curvature of the potential, which is of the order of the Hubble length in the present case. This means that the width of cosmic strings increases with time. This is one of the outstanding characteristics of cosmic strings in this scenario compared with other cosmic strings considered in the literature.
When the Hubble induced mass decreases down to the mass of the flat direction m φ , the flat direction starts to oscillate around the origin of the potential and cosmic strings disappear. The time when cosmic strings disappear is given by Eq. (4.7) as in the previous section. For the case around m φ = O(1) TeV, cosmic strings disappear much earlier than the BBN epoch, so that observational constraints on long-lived topological defects is not applicable to our scenario. As we explained in the previous section, we consider the well motivated case that the cosmic strings disappear before reheating completes.
After they form at the end of inflation and before they disappear at t t decay , cosmic strings emit GWs. We numerically follow the evolution of the flat direction and calculated its GW signal. The results are shown in the next subsection.

Results of numerical simulations
In this subsection, we show our results of numerical simulations for the dynamics of the flat direction and its GW signal.
By redefining the variables, we rewrite the equation of motion such as The parameter a H can be absorbed by the redefinition of φ, so that our numerical simulation can be applied to any case with an arbitrary value of a H . We follow the time evolution of the flat direction by solving the equation of motion Eq. (5.5) numerically. Details of the numerical method are explained in the previous section. We confirm that each Hubble volume contains O(1) cosmic strings and a typical width of cosmic strings is roughly proportional to the Hubble radius. From our numerical simulation, we find that the spatially averaged magnitude of φ grows with time and eventually reaches of the order of the VEV (Eq. (5.4)) around τ H i = τ form 10 for the case of c Hosc = −15.
We calculate the transverse-traceless part of the anisotropic stress from the numerical simulation and obtain the energy density of GWs using the equations in Appendix B. The results are plotted in Fig. 6. One can find that the GW peak energy density is almost constant in time and is almost independent of m. We give physical interpretation of these results in the next subsection.
Taking a non-zero value of m φ into account, we find that the flat direction starts to oscillate around the origin of the potential and cosmic strings disappear around τ = τ decay 40, which is given by Eq. (4.16) with m φ /H i = 5 × 10 −4 and c Hosc = −15. After that, the GW peak wavenumber becomes constant (blue dashed line in Fig. 6) and the GW spectrum decrease adiabatically as ∝ τ −2 (∝ a −1 ). As we discussed in Sec. 4.2, the above choice of m φ /H i is convenient for numerical simulations, though it does not change the dynamics of flat direction until τ τ decay . In reality, it can be many orders of magnitude smaller depending on inflation models.

Physical interpretations
In this subsection, we discuss how we can understand the results of our numerical simulations shown in Fig. 6.
Since during the scaling regime the typical length scale of cosmic string dynamics is the Hubble length, the GW emission peak wavenumber is typically given by the Hubble scale like Eq. (4.25). The produced energy density of GWs is estimated as Eq. (4.28) but with the VEV given by Eq. Cosmic strings emit GWs with a constant energy density expressed above. We plot our numerical results of ∆Ω gw /∆ log τ in Fig. 7. It shows that the produced energy density of GWs is constant in time as expected. We can determine the numerical prefactors for the emission peak wavenumber and the produced energy density of GWs from the numerical results like . (5.12) We find that the numerical prefactors are almost independent of m and τ up to by a factor of two. Figure 7 shows that ∆Ω gw /∆ log τ is proportional to a certain power of k for large and small wavenumber modes like Eq. (4.32), where α 2.5 in this case. The GW spectra for  large wavenumber modes (k > k ∆ peak (τ )) in Fig. 6 can be understood by an integration like Eq. (4.33), where we can effectively take n → ∞ in this case. Thus, we obtain where τ k is defined by k peak (τ k ) = k. The result of Eq. (5.15) at τ = 50 is plotted as the black curves in Fig. 6. It is consistent with the large-wavenumber behaviour of GW spectrum found in our numerical simulations, though it is larger than our result by a factor of 2 in the case of m = 4. Note that the GW spectrum is expected to be proportional to k −2 for k ∆ peak (τ form ) k k ∆ peak (τ ). This asymptotic behavior is different from the one obtained in the previous section because the VEV is constant in time in the present case. The gray dotted curves in Fig. 6 show that simulated super-horizon modes are consistent with the spectrum expected from the discussion of causality [14,36] as explained in the previous section.
The blue dashed lines in Fig. 6 show the evolution of GW peak wavenumber in the case of m φ /H i = 5 × 10 −4 . Since the flat direction starts to oscillate around the origin due to its soft mass term, cosmic strings disappear at τ τ decay . The GW peak energy density and frequency at that time is given by [Ω gw (τ decay )] peak 10 1 4 where the numerical prefactors are determined from our numerical results. The shape of GW spectrum is in fact freezed out around τ = τ decay 40, which is given by Eq. (4.16) with m φ /H i = 5 × 10 −4 and c Hosc = −15. After that, the peak wavenumber does not change and the GW energy density decreases as ∝ a −1 until reheating completes.
In summary, the schematic form of resulting GW spectrum is similar to the one described in Fig. 4, except for the power-law index of large wavenumber modes. In the case of this section, the GW spectrum for large wavenumber modes is proportional to k −2 and is independent of m. This is because the potential minimum for the flat direction is constant in time. Thus, in principle, we can distinguish the cases with and without the higher dimensional superpotential through the observation of the GW spectrum with higher frequencies.

GW signals at present
The above calculations and discussions determine the GW spectrum during the oscillation dominated era. In this subsection, we derive the GW spectrum at the present epoch, taking into account the redshift.
The peak energy density of GWs at τ = τ decay , [Ω gw (τ decay )] peak is given by Eq. (5.16) from our numerical simulations. Substituting [Ω gw (τ decay )] peak into Eq. (4.37), we obtain the present value of GW peak energy density such as Note that the peak frequency is independent of m and a H . Since the GW spectrum is sensitive to what dominates the energy density of the Universe, it contains information of reheating temperature as explained in the previous section. In fact, the GW spectrum bends at the frequency given by Eq. (4.41) corresponding to the wavenumber around k τ −1 RH [9,17,18]. Since there are three observables: Eqs. (4.41), (5.18) (see also Eq. (5.4)), and (5.19), in principle, we can obtain three parameters through observation of GW spectrum: |c Hosc | −1/2 m φ , T RH , and M * /M Pl (up to by a factor of order one). Again, let us stress that we can distinguish the case of W = 0 in the previous section from the case of W = 0 in this section. This is because for high frequencies the GW spectrum scales as Ω gw ∝ f (16−2n)/(n−2) ( = f −2 ) for the former case and Ω gw ∝ f −2 for the latter case.  19)). In this case, however, we can not ignore the effect of the backreaction of GW emission when we calculate the spectrum of GWs. In addition, the gradient energy density of the flat direction is comparable to the energy density of the Universe and thus it may spoil the isotropic evolution of the Universe. We can neglect these effects when the cutoff scale is less than the Planck scale. Figure 8 shows examples of GW spectrum predicted by the present mechanism. The GW peak energy density and peak frequency are given by Eqs. (5.18) and (5.19), which are obtained from our numerical results. While the shape of the GW spectrum for low frequency modes are calculated from Eqs. (B.6) and (B.7) with a constant T TT ij and τ f = τ decay , that for large wavenumber modes are calculated from Eq. (5.15). In the figure, we take m φ = 10 2 GeV with T RH = 10 7 GeV (blue dashed curve), m φ = 10 2 GeV with T RH = 10 9 GeV (red dashed curve), and m φ = 10 3 GeV with T RH = 10 9 GeV (red dot-dashed curve). We could obtain rough scales of them from the observation of GW spectrum. To discuss the detectability of GW signals, we plot sensitivities of planned interferometric detectors. The sensitivity curves are the same as in Fig. 5, and details are given in Sec. 4.4. We find that GW signals would be observed by ET or Ultimate DECIGO.
6 Case C: W = 0 without U (1) symmetry In this section, we consider the same potential with the one considered in the previous section, Eq. (5.1), except for an additional U (1) breaking term coming from higher dimensional Kähler potentials like where m ≥ 2. We assume m < 2m − 2 so that the potential for the flat direction is bounded from below. The coefficient b H is given by where we assume b H = O(1). Here we have redefined the phase of the field φ such that b H is postitive. This additional term explicitly breaks U (1) symmetry, but there remains Z m symmetry in that potential.
The following discussion and calculation in this section are essentially the same as the ones explained in the previous two sections, though it is more complicated due to the additional term breaking U (1) symmetry. In the next subsection, we briefly explain the dynamics of the flat direction in the early Universe. We numerically follow the evolution of its dynamics and its GW signal, and the results are shown in Sec. 6.2. Then we give physical interpretation of the results in Sec. 6.3. Finally, we discuss the detectability of GW signals in Sec. 6.4.

Dynamics of the flat direction
Throughout this section, we consider the following potential for the flat direction: where we assume that c H (= c H inf ) > 0 during inflation and c H (= c Hosc ) < 0 after inflation. The positive c H inf makes the flat direction stay at the origin during inflation. After inflation ends, the flat direction obtains a large VEV due to the negative c Hosc . Note that the U (1) breaking term unstabilizes the potential for the flat direction around the cutoff scale. For a small b H , the potential minimum is determined by the Hubble induced mass term as in the previous section (see Eq. (5.4)), while for a small c H inf , it is determined by the U (1) breaking term. Thus, the VEV can be estimated as where p = 0, 1, 2, . . . , m − 1. In the last line, we use Eqs. (3.10), (3.16), and (3.18). One can see that the VEV is as large as the cutoff scale M * . Since the nonzero VEV of the flat direction breaks Z m symmetry, m types of domain walls form after the end of inflation. In general, a typical width of domain walls is given by the inverse of the curvature of the potential at the minimum. In the present case, it is given by the inverse of the mass of the phase direction m θ , which is given by The energy of domain walls per unit area σ can be estimated by These mean that the width of domain walls increases and the energy per unit area decreases with time. These are outstanding characteristics of domain walls in this scenario compared with other domain walls considered in the literature. By using numerical simulations (see the next subsection), we confirm that domain wall system soon reaches the scaling regime, in which the number of domain walls in the Hubble volume is of the order of the domain wall number m . In this regime, the energy density of domain walls is roughly given by We should emphasize that the energy density of domain walls never dominate that of the Universe. This is because the curvature of the potential is given by H(t) and decreases with time as ∝ a −3/2 , which is completely different from the ordinary domain walls [44]. When the Hubble parameter decreases down to the mass of the flat direction m φ , the flat direction starts to oscillate around the origin of the potential and domain walls disappear. Since domain walls disappear well before the BBN epoch, our scenario is free from stringent constraints on long-lived domain walls [44]. As we explained in Sec. 4.1, we consider the well motivated case that the domain walls disappear before reheating completes.
After they form at the end of inflation and before they disappear at t t decay , domain walls emit GWs. In the subsequent subsections, we show our numerical results of domain wall dynamics and spectra of GWs emitted from these domain walls.

Results of numerical simulations
In this subsection, we show our results of numerical simulations for the dynamics of the flat direction and its GW signal.
We redefine the variables in the same way as in the previous subsection (see Eqs. (5.6)-(5.9)) and obtain the equation of motion equivalent to Eq. (5.5) except for the additional U (1) breaking term: The parameter a H can be absorbed by the redefinition of φ and b H , so that our numerical simulation depends on a H and b H only though the combination of b H defined by We follow the time evolution of the flat direction by solving the equation of motion Eq. (6.11) numerically. Details of the numerical method are written in Sec. 4.2. Figure 9 shows one example of domain wall network evolution obtained by our numerical simulation. We set c Hosc = −15, m = 4, and m = 2. The green and cyan regions satisfy −π/8 < θ < π/8 In the blue regions, the VEV of the flat direction satisfies |φ| < | |φ| | /30, where |φ| is the VEV given by Eq. (6.6). In the green and cyan regions, the phase of the flat direction satisfies −π/8 < θ < π/8 and (7/8)π < θ < (9/8)π, respectively. The black line at the bottom-left corner represents the length of the comoving horizon size τ at each time. and (7/8)π < θ < (9/8)π, respectively, where θ is the phase of the flat direction. They represent the inside of domain walls. The blue region satisfies |φ| < |φ| /30 and describes cosmic strings in the limit of b H → 0. The black line at the bottom-left corner in Fig. 9 describes a unit of the comoving horizon size (given by τ ) at each time. One can see that a typical width of domain walls increases with time and is roughly proportional to the Hubble length. One can also see that each Hubble volume contains O(1) domain walls of each type (m = 2 in Fig. 9). From our numerical simulation, we find that the spatially averaged magnitude of φ grows in time and eventually reaches of the order of the VEV (Eq. (6.6)) around τ form 10 in this case.
We calculate the transverse-traceless part of the anisotropic stress from the numerical simulation and obtain the energy density of GWs using the equations in Appendix B. The results are plotted in Figs. 10 and 11. One can find that the GW peak energy density is almost constant in time. We give physical interpretation of these results in the next subsection.
Taking a non-zero value of m φ into account, we find that the flat direction starts to oscillate around the origin of the potential and domain walls disappear around τ τ decay , which is given by Eq. (4.7). After that, the GW peak wavenumber becomes constant and the GW spectrum decreases adiabatically as ∝ τ −2 (∝ a −1 ). In Figs. 10 and 11, the GW spectrum is freezed out around τ = τ decay 40, which is given by Eq. (4.16) with m φ /H i = 5 × 10 −4 and c Hosc = −15. As we discussed in Sec. 4.2, the above choice of m φ /H i is convenient for numerical simulations, though it does not change the dynamics of flat direction until τ τ decay . In reality, it can be many orders of magnitude smaller depending on inflation models.

Physical interpretations
In this subsection, we discuss how we can understand the results of our numerical simulations shown in Figs. 10 and 11.
During the scaling regime, the typical length scale of domain wall dynamics is the Hubble length. As a result, the emission peak wavenumber of GWs emitted from these domain walls is estimated by the Hubble scale like Eq. (4.25). Its quadrupole moment Q can be roughly given by [35] Q ∼ H −2 × m σH −2 , (6.13) where m σH −2 is the total energy of domain walls within a Hubble volume. Such domain walls emit GWs with the luminosity of Eq. (4.27), which is proportional to the square of the quadrupole moment of domain walls. Thus, the produced energy density of GWs can be estimated as (6.14) Domain walls emit GWs with a constant energy density given by Eq. (6.14). Note that cosmic strings should coexist with domain walls because in the limit of b H → 0, the potential for the flat direction restores U (1) symmetry. This implies that GWs are also emitted from cosmic strings. Thus, the total GW energy density is given by where the second term comes from the contribution of cosmic strings and is relevant for small b H . We include a numerical prefactor k br to fit our numerical results. Note that the first term in the right-hand side of this equation is proportional to M 2 * /M 2 Pl (Eq. (6.14)), while the second term is proportional to M 4 * /M 4 Pl (Eq. (5.10)). This results from the fact that the energy density of domain walls is proportional to M * (see Eq. (6.10)), while that of cosmic strings is proportional to M 2 * (see Eqs. (4.5), (4.6) and (5.4)). Since the energy density of domain walls is much larger than that of cosmic strings, domain walls emit GWs more efficiently.
We plot our numerical results of ∆Ω gw /∆ log τ in Figs. 12 and 13. The results are constant in time as expected. We can determine the numerical prefactors for the emission peak wavenumber and the produced peak energy density of GWs from the numerical results like Since the flat direction starts to oscillate around the potential minimum due to its soft mass term, domain walls disappear at τ τ decay . The GW peak energy density and frequency at that time is given by where the numerical prefactors are determined from our numerical results. The shape of GW spectrum is in fact freezed out around τ = τ decay 40, which is given by Eq. (4.16) with m φ /H i = 5 × 10 −4 and c Hosc = −15. After that, the peak wavenumber does not change and the GW energy density decreases as ∝ a −1 until reheating completes.
In summary, the schematic form of resulting GW spectrum is the same as the one described in the previous section, even though domain walls form. It is impossible to distinguish GW signals emitted from cosmic strings and domain walls. However, we can obtain information of SUSY scale and reheating temperature through the observation of GW signals, as we explain in the next subsection.

GW signals at present
The above calculations and discussions determine the GW spectrum during the oscillation dominated era. In this subsection, we derive the GW spectrum at the present epoch, taking into account the redshift. The peak energy density of GWs at τ decay is given by Eq. (6.19) from our numerical simulations. Substituting Ω gw h 2 (τ decay ) peak into Eq. (4.37), we obtain the present value of GW energy density such as  Note that the peak frequency is independent of a H , b H , m, and m . Since the GW spectrum is sensitive to what dominates the energy density of the Universe, it has information of reheating temperature as explained in the previous section. In fact, the GW spectrum bends at the frequency given by Eq. (4.41) corresponding to the wavenumber around k τ −1 RH [9,17,18]. There are three observables: Eqs. (4.41), (6.22), and (6.23), so that we can in principle obtain three parameters through observation of GW spectrum: |c Hosc | −1/2 m φ , T RH , and M * /M Pl (up to by a factor of O(1)). Let us stress that |c Hosc | −1/2 m φ and T RH are not affected by higher dimensional operators and are independent of the sources of GW signals, such as cosmic strings and domain walls. Note that we assume M * M Pl to avoid the effect of the backreaction of GW emission. Figure 14 shows examples of GW spectrum predicted by the present mechanism. The GW peak energy density and peak frequency are given by Eqs. (6.22) and (6.23), which are obtained from our numerical results. While the shape of the GW spectrum for low frequency modes are calculated from Eqs. (B.6) and (B.7) with a constant T TT ij and τ f = τ decay , that for high frequency modes are calculated from Eq. (5.15). In the figure, we take m φ = 10 2 GeV with T RH = 10 7 GeV (blue dashed curve), m φ = 10 2 GeV with T RH = 10 9 GeV (red dashed curve), and m φ = 10 3 GeV with T RH = 10 9 GeV (red dot-dashed curve). To discuss the detectability of GW signals, we plot sensitivities of planned interferometric detectors. The sensitivity curves are the same as in Fig. 5, and details are given in Sec. 4.4. We find that GW signals would be observed by ET or Ultimate DECIGO.

Summary and Conclusions
We have investigated the dynamics of a flat direction in the early Universe, focusing on the case that the coefficient of the Hubble induced mass term is positive during inflation and is negative after inflation (oscillation dominated era). While the flat direction stays at the origin during inflation, it obtains a large VEV after inflation due to the negative Hubble induced mass term. In many cases, flat directions have U (1) or Z n symmetries and thus topological defects, such as cosmic strings and domain walls, form just after the end of inflation. When the Hubble parameter decreases down to the soft mass of the flat direction, it starts to oscillate around the origin of its potential and topological defects disappear. These topological defects disappear well before the BBN epoch and they are free from constraints on long-lived topological defects.
The topological defects emit GWs with an emission peak wavenumber corresponding to the Hubble scale at each time until they disappear when the Hubble parameter decreases down to the soft mass of the flat direction. This means that the resulting GW spectrum has a peak wavenumber corresponding to the soft mass of the flat direction [9]. Therefore, we can obtain the information of SUSY scale through the detection of GW signals. In addition, since the GW spectrum bends at the wavenumber corresponding to the horizon scale at the end of reheating, we can obtain the value of reheating temperature through the observation of the bending frequency. Since the GW energy density depends on the energy density of topological defects, it is related to the VEV of the flat direction. This means that the peak energy density of GW spectrum and the shape of GW spectrum for high frequency modes have information of higher dimensional terms for the flat direction. In this paper, we have calculated the GW spectrum emitted from cosmic strings and domain walls and have shown that the resulting GW spectra actually have the information of SUSY scale, reheating temperature, and higher dimensional terms for the flat direction. We investigate three cases for higher dimensional potentials for the flat direction; the one coming from a nonrenormalizable superpotential, the one coming from a U (1) symmetric Kähler potential, and the one coming from U (1) breaking Kähler potentials. Especially in the first case, we have found that the resulting GW spectrum for high frequnicy modes has information of power-law index for the nonrenormalizable superpotential. We have also found that the resulting GW signals would be measured by future GW observation experiments such as Advanced LIGO, ET, and Ultimate DECIGO.
Here let us comment on the relation between the mass of the flat direction and Q-ball formation [45][46][47][48][49]. Note that m φ is a function of φ through the renormalization. If dm φ (φ)/dφ < 0, the flat direction fragments into long-lived non-topological solitons called Q-balls soon after the flat direction starts to oscillate around φ = 0 at t t decay . This is an obstacle in light of GW detection because it has been shown that GWs are diluted by the existence of Q-ball dominated era [50]. In this paper, therefore, we assume dm φ (φ)/dφ > 0 implicitly. This is realized in a wide class of gravity mediation models [48,51]. In particular, Q-balls never form in SUSY models with a hierarchical spectrum between gauginos and the other sparticles in the MSSM. In models of gauge mediation, since the VEV of the flat direction induces the effective mass of gauge fields, soft masses are effectively suppressed as ∝ φ −2 for g |φ| M s , where M s is a messenger scale [52,53]. Thus, the flat direction need to contain a Higgs field to obtain a sizable value of m φ and satisfy dm φ (φ)/dφ > 0 [54]. If GWs are observed in models of gauge mediation, we can identify m φ as the Higgsino mass parameter µ.
Furthermore, let us comment on the baryon asymmetry (see footnote 5). The flat direction may be kicked by its A term along its phase direction after the inflation. The phase of the flat direction is randomly distributed through the phase transition. Therefore, though the baryon density may be locally generated by the Affleck-Dine mechanism, net baryon asymmetry is absent. The baryon asymmetry is thus never generated through the dynamics of this flat direction. To generate the baryon asymmetry, we may require another flat direction which has a large VEV during inflation and realizes the Affleck-Dine baryogenesis. The leptogenesis is also viable if the gravitino problem can be avoided [32][33][34]55].
In the literature, it is known that cosmic strings emit Nambu-Goldstone modes efficiently, which leads to an effective friction on the dynamics of cosmic strings. As a result, cosmic string loops exist for a long time, which become more efficient source of GW emission than the cosmic string network in the ordinary cosmic strings [56][57][58][59]. On the other hand, in the present case, a typical width of cosmic strings is of the order of the Hubble length. Since the length of cosmic string loops can not be lower than its width, cosmic string loop formation is suppressed. This is why the effect of cosmic string loops can be neglected in the present paper. Note that the effect of cosmic string loops is incorporated in our numerical simulations.
One might wonder if (tachyonic) preheating occurs for the flat direction just after the end of inflation and it affects the resulting GW signals [36,[60][61][62][63][64][65]. In fact, our numerical calculations incoperate the effect of tachyonic preheating for the flat direction. However, let us emphasize that GWs in the observable band are not emitted from preheating which occurs when cosmic strings appear, but emitted from cosmic strings when they disappear. This is because the Hubble scale at the time of the preheating is much larger than the one at the time when topological defects disappear. Since GW signals with such high frequencies are beyond the detectability of future GW observations, we neglect GWs emitted from the preheating. We can also neglect the GW emission from the field I because its typical frequency is far beyond the observable region. This is because the characteristic energy scale of that dynamics is the mass of I, which is of the order of inflation scale. The typical GW frequency emitted from that dynamics is usually much higher than the observable region of 0.1-1 Hz. Furthermore, we consider that the reheating temperature is of the order of or less than 10 7−9 GeV. This implies that inflaton decays into radiation perturbatively via Planck suppressed interactions. In this case, GW emission is not efficient at the time of reheating.

A.2 Hybrid inflation model
In the simplest hybrid inflation model proposed in Ref. [25,26], the superpotential is given by where S,Ψ, and Ψ are superfields, and λ and µ are positive parameters satisfying λ µ. The potential for their scalar components V (S, Ψ,Ψ) is thus given by When all scalar fields are well below the Planck scale (Ψ, Ψ, S 1) but satisfy the relation of S 2µ/ √ λ,Ψ and Ψ are stabilized at the origin of the potential due to large effective masses. In this case, inflation is driven by a false vacuum energy density of µ 4 . Driven by 1-loop induced potential for the field S, it rolls down to the origin of the potential. After the field S reaches the critical value of S cr = 2µ/ √ λ, the water-fall fieldsΨ, Ψ starts to roll down to the potential minimum at Ψ Ψ = µ 2 /λ. Although this simplest model of hybrid inflation in supergravity might predict a spectral index inconsistent with the observation, many successful hybrid inflation models are considered in the literature based on this simplest model. Note that inflation is driven by the F term of the superfield S, while the energy density of the Universe is dominated by that of oscillation energy of water-fall fieldsΨ and Ψ after the end of inflation. In the main part of this paper, the field X is identified with the field S while I collectively represents water-fall fieldsΨ and Ψ.

B Method for the calculation of GWs.
GWs are emitted by cosmic strings [10-13, 15, 71, 72] and/or domain walls [14][15][16] between the two phase transitions, that is, after the end of inflation and before the time of t decay . In this appendix, we explain how to calculate the spectrum of GWs. The calculations are based on Refs. [14,36], which is suitable to our situation compared with the method to calculate GW energy density from localized sources derived in Ref. [68].
The Fourier transformed transverse-traceless part of the spatial metric perturbation h ij obeys the following linearized Einstein equation: where the prime denotes a derivative with respect to conformal time τ defined by dτ = a −1 dt, and T TT ij is the Fourier transformed transverse-traceless part of the anisotropic stress. When the source term T T T ij lasts during the interval of [τ i , τ f ], the solution of Eq. (B.1) for τ f τ τ RH is given by [14] h ij (τ, k) = A ij (k) kτ a j 1 (kτ ) + B ij (k) kτ a n 1 (kτ ) , where j 1 and n 1 are the spherical Bessel and Neumann functions of order one. The coefficients A ij and B ij are given by After reheating completes (τ τ RH ), the solution is described by the spherical Bessel functions of order zero: We match the solution of Eq. (B.2) with that of Eq. (B.5) at τ = τ RH . We then obtain A and B such as where the coefficients are given by The energy density of GWs can be calculated from [35] ρ gw = 1 32πG ḣ ij (t, x)ḣ ij (t, x) , (B.10) where · · · represents an ensemble average for a stochastic background. We replace the ensemble average by an average over a volume V of the comoving box. This is a good approximation for sub-horizon modes. Using the following convention for the Fourier expansion: we rewrite the energy density of GWs as Here we have used h ij kh ij and taken average over the oscillation period since we can observe only sub-horizon modes k τ . Also we have averaged over time in Eq. (B.12) because we are not interested in the resolution of the oscillation of h ij (τ ) with time. We define the GW spectrum as where ρ tot (τ ) (= 3M 2 Pl H 2 (τ )) is the total energy density of the Universe. We can calculate the GW spectrum such as [36] Ω gw (τ ) (B.14) We can calculate T TT ij from lattice simulation for the classical dynamics of a scalar field and then calculate the value of A and B through Eqs. (B.6) and (B.7). Finally, we can obtain the GW spectrum from Eq. (B.14).

C Case with flavor symmetry
In this appendix, we investigate the components of the flat direction in detail. In particular, we consider the case with flavor symmetry, where topological defects may not form when the flat direction obtains a large VEV and phase transition occurs. However, Nambu-Goldstone (NG) modes are randomly distributed over the Hubble scale through the phase transition, so that GWs are emitted from their relaxation dynamics within the Hubble horizon. The resulting GW spectrum is qualitatively the same as the ones obtained in the main part of this paper. This is a realization of so-called self-ordering scalar fields in SUSY theories [69][70][71][72].
In the next subsection, we investigate flavor components of the flat direction in detail and explain what happens in the case with flavor symmetry. Then we briefly review a method suitable to the calculation of GW spectra in this case, which has been derived in the literature in the context of self-ordering scalar fields. Using this method, we calculate the GW spectrum and show the results in Sec. C.3.

C.1 Flavor symmetry
In order to investigate flavor components of flat directions, let us consider the LH u flat direction as an illustration. Since the left-handed lepton doublet L contains three flavors, the LH u flat direction has an ambiguity in terms of the flavor direction, that is, any flavor directions are D-flat directions. This implies that the LH u flat direction can be written as where the coefficient a H is of the order of (M Pl /M * ) 2m−2 and may also break flavor symmetry in general. If we consider, for example, the case of a i 1 j 1 ...i m−n j m−n H = a H δ i 1 j 1 . . . δ i m−n j m−n , 12 the superpotential respects U (3).
Suppose that the Hubble induced mass term is flavor-universal and the higher dimensional superpotential has the U (3) symmetry. As we can see from Eqs. (C.5) and (C.8), the coefficient of the Hubble induced mass term during inflation, c H inf , is generally different from the one during the oscillation dominated era, c Hosc [9]. We consider the case that the coefficient of the Hubble induced mass term is positive during inflation and is negative after inflation. In this case, the flat direction stays at the origin during inflation and has a large VEV after inflation, that is, the phase transition occur just after inflation. At this phase transition, the U (3) symmetry is broken to U (2) symmetry. In this case, topological defects do not form at this phase transition. However, there is another source of GWs. Since NG models have no correlation over distances larger than the Hubble radius due to the causality, they are randomly distributed over the Hubble scale. Therefore, the NG modes obtain gradient energy. The gradient energy turns into the kinetic energy and decreases inside the Hubble radius. As a result, the modes with wavenumber around the Hubble scale (k ∼ aH(t)) dominantly emits GWs [69][70][71][72]. Since the NG modes outside the Hubble radius grow through the nonlinear couplings, the peak frequency of the resultant GW signals correspond to the Hubble scale when the flat direction start to oscillate around the origin at t t decay (Eq. (4.7)). This is another source of GWs which gives a contribution to the GW signal. Note that it is possible that, say, there is O(N )×U (1) symmetry and is broken down to O(N −1) symmetry. In this case, topological defects also form and the resulting GW signal is given by the sum of both contributions. In the following subsections, we calculate the GW energy density and spectrum from the above mechanism.
Hereafter we focus only on the global symmetry of the flat direction for the following reasons. The gauge fields may be randomly distributed at the phase transition. We can take the unitary gauge if the topological defects are not formed. Since, the broken gauge field obtains the effective mass of g φ , the field configuration begins to oscillate and becomes trivial (A µ = 0) soon after the phase transition. For the unbroken gauge symmetry, the corresponding gauge field remains massless. Here we should note that while long wavelength modes of such gauge fields remain constant until they enter the horizon, the NG modes grow through the nonlinear term [72]. Therefore, GW emission from the local symmetry can be neglected.

C.2 Analytical calculation of GW spectra for self-ordering scalar fields
Here we explain a method suitable to the calculation of GW spectra emitted from self-ordering scalar fields. We consider the case with O(N ) symmetry for simplicity. Note that LH u flat direction considered in the previous section has O(6) symmetry rather than U (3) symmetry if gauge and Yukawa couplings are ignored. This is similar to O(4) custodial symmetry in the Standard Model. We expect that the results do not change qualitatively for the case with other symmetries because the discussion in the previous subsection hardly depends on the detail of the manifold where NG bosons reside.

C.3 Results for self-ordering scalar fields
When the Hubble parameter decreases down to the mass of the flat direction m φ , the flat direction starts to oscillate around the origin of the potential. After that, the field value of the flat direction decreases with time due to the expansion of the Universe such as φ ∝ a −3/2 . After the phase transition at the end of inflation and until they start to oscillate around the origin at t t decay , the dynamics of NG modes emit GWs continuously. We obtain the GW spectrum by evaluating Eqs. (C.12), (C.13), and (C.14) numerically. We obtain GW peak energy density and peak wavenumber at t t decay such as Note that we assume M * /M Pl 1 so that we can neglect the backreaction of the flat direction dynamics on the cosmic expansion. Substituting the above formulas into Eqs. (4.37) and (4.39), we obtain the present value of the GW peak energy density such as Ω gw h 2 (t 0 ) 3 × 10 −8 N −1 |c Hosc | −1/2 m φ 10 3 GeV Since a typical length scale of the dynamics of NG modes is the Hubble length, the emission peak wavenumber of GWs is estimated as the Hubble scale. In addition, the GW energy density of Eq. (C.12) is proportional to φ 4 . These observations imply that the resulting GW spectrum emitted from the dynamics of NG modes is qualitatively the same as the ones obtained in Secs. 4, 5, and 6, depending on higher dimensional potentials for the flat direction. In particular, since the GW spectrum is sensitive to what dominates the energy density of the Universe, it has information of reheating temperature as explained in those sections. In fact, the GW spectrum bends at the frequency given by Eq. (4.41) corresponding to the wavenumber around k τ −1 RH [9,17,18]. In summary, the schematic form of resulting GW spectrum is similar to the one described in the main part of this paper, even if there is a large flavor symmetry such as O(N ). Figure 15 shows examples of GW spectrum generated by the self-ordering scalar fields. We assume that φ = M * is constant until t = t decay as in Secs. 5. The GW spectra are calculated from Eqs. (C.12), (C.13), and (C.14). In the figure, we take m φ = 10 2 GeV with T RH = 10 7 GeV (blue dashed curve), m φ = 10 2 GeV with T RH = 10 9 GeV (red dashed curve), and m φ = 10 3 GeV with T RH = 10 9 GeV (red dot-dashed curve). To discuss the detectability of GW signals, we plot sensitivities of planned interferometric detectors. The sensitivity curves are the same as in Fig. 5, and details are given in Sec. 4.4. Detecting these GW signals seems somewhat challenging in the current design of Advanced LIGO, ET, and Ultimate DECIGO unlike those from topological defects in Secs. 4, 5, and 6. However, let us stress that if we can detect these GW signals from somehow, for example, through longer observation time T 1 yr, they provide us invaluable information of |c Hosc | −1/2 m φ , T RH , and M * /M Pl (up to by a factor of O(1)).