Cosmic R-string, R-tube and Vacuum Instability

We show that a cosmic string associated with spontaneous $U(1)_R$ symmetry breaking gives a constraint for supersymmetric model building. In some models, the string can be viewed as a tube-like domain wall with a winding number interpolating a false vacuum and a true vacuum. Such string causes inhomogeneous decay of the false vacuum to the true vacuum via rapid expansion of the radius of the tube and hence its formation would be inconsistent with the present Universe. However, we demonstrate that there exist metastable solutions which do not expand rapidly. Furthermore, when the true vacua are degenerate, the structure inside the tube becomes involved. As an example, we show a"bamboo"-like solution, which suggests a possibility observing an information of true vacua from outside of the tube through the shape and the tension of the tube.

In [2,3,9], by exploiting the Nelson-Seiberg theorem [1], a connection between metastability and R-symmetry was demonstrated in the context of generalized Wess-Zumino models with generic superpotential. From more phenomenological viewpoint, the U(1) R symmetry must be broken explicitly or spontaneously to generate Majorana gaugino masses. Gaugino masses are not induced by SUSY breaking without the U(1) R symmetry breaking.
Through the cosmological phase transition, there may appear solitonic objects such as domain walls, cosmic strings and monopoles [25] through the Kibble-Zurek mechanism [26,27]. When a global U(1) symmetry is spontaneously broken, a global string appears [28]. Thus, when the U(1) R symmetry is broken spontaneously in SUSY models, there would appear a global string, which we refer as an R-string. 1 The R-string would be stable in those models in which the U(1) R breaking vacuum is a global minimum, and that would lead to several cosmologically interesting aspects [31].

R-string
Consider the following simplest spontaneous R-symmetry breaking model. Superpotential is linear in a chiral superfield X which will be a trigger for SUSY breaking, To stabilize the pseudo-moduli X in this SUSY breaking vacuum, we introduce the following non-canonical Kähler potential by hand, (2.1) Thus, the potential of this theory is given by This model can be viewed as a low energy effective theory of one of the O'Raifeartaigh models studied in [12]: When the pseudo-moduli space is stable everywhere along messenger directions, by integrating out the messengers, one obtains non-trivial corrections to the Kähler potential. Expanding the Kähler potential up to O(|X| 6 ), one can reproduce a theory similar to (2.2). On the other hand, when the pseudo-moduli space has a tachyonic direction at a point in the space, which is phenomenologically interesting situation in gauge mediation models [3,20], the existence of messengers is crucial and two-field model is required. This is the main topic in the next section.
When µ 2 X > 0 and λ X > 0, the X field develops its vacuum expectation value and the Rsymmetry is broken (X has the R charge 2). The minimum of the potential V (X) is obtained at |X| = X min ≡ 2µ 2 X /λ X . Note that this vacuum is the global minimum of the potential V (X). Since the global U(1) R symmetry is spontaneously broken, the R-string would be formed.
Let us introduce dimensionless variables as Then the effective Lagrangian is given by (2.4) Here a positive definite metric at the minimum (T = 1) requires λ > 1/4. We use the dimensionless cylindrical coordinate (ρ, θ,z) for constructing a straight R-string along the z-axis. We make the following standard Ansatz, where f (0) = 0 and f (ρ) → 1 at ρ → ∞. We numerically solve the equation of motion for a minimal winding solution (n = 1). The solution is shown in Fig. 1.
For later convenience, let us estimate a size of R-string, R, by using the following simple approximation f (ρ) = ρ R n for ρ ≤ R and f (ρ) = 1 for ρ > R, (2.6) where the power of ρ is determined by requiring smoothness of the configuration at ρ = 0.
The total energy of this configuration, E, per the string length ∆z is estimated as where an R independent constant c n (λ) which should be numerically determined is introduced, and an IR-cutoff Λ is also introduced in order to regularize a well-known logarithmic divergence of a global vortex. The above energy takes the minimum at R ≈ R string , which is the transverse size of the R-string. Here m T is the dimensionless mass of T in the vacuum T = 1. For instance for λ = 1/2 and n = 1, we obtain m T = 1 and R string = 2 √ 3.
To check the approximation (2.6) by comparing with numerical calculations, we introduce another definition of a transverse size of the R-string as Here K T gives a finite contribution from the kinetic term along the ρ direction into the energy and is useful to define the transverse size. Note that this quantity does not include the cut-off dependence. We compute this both analytically with the linear approximation and numerically, see Fig.2. For instance we observe R T = 1.95 for λ = 1/2 by a numerical calculation. As can be seen in Fig. 2, the linear approximation nicely reproduces the numerical results (we need to pay attention to errors of 10% − 30% in the linear approximation).

Tube solution
Here, in order to illustrate the tube solution, we study a non-supersymmetric bosonic theory as a toy model. Let us study the model with the following scalar potential, This potential has two degenerate vacua, that is, |X| = 0 and v. At the former vacuum, the global U(1) symmetry is unbroken, while the U(1) symmetry is broken at the latter vacuum.
Then, a global string would be formed. Again, let us rescale the fields and coordinates as then the Lagrangian becomes We make the Ansatz for the minimally winding string, where f (0) = 0 and f (r) → 1 at r → ∞. The solution is again obtained numerically which is shown in Fig. 3. As can be seen in Fig. 3, the string has a substructure that is a hole inside the string. Thus, we refer it as the tube. It is the asymmetric phase (X = 0) outside the tube while it is the symmetric (X = 0) phase inside it.
This tubelike string solution can be regarded as a tube of a domain wall. Indeed, there also exists a domain wall in this model. For instance, a solution interpolating the two vacua T = 1 at x 1 = ∞ and T = 0 at x 1 = −∞ is given by 14) with a dimensionless tension T wall = 1/2. Thus, assuming that the field configuration of the tube along the radial direction is well described by this solution, the total energy E per the tube length ∆z of the tubelike solution with a radius ρ = R can be estimated by as long as the "thickness" of the wall is much smaller than the radius R. Minimizing this, we get the transverse size of the tube solution Note that the stabilization mechanism of the tube solution is different from that of the Rstring (without a hole) where the kinetic energy and the potential energy are balanced. As a result, the transverse sizes have different dependences on the winding number n as R ∝ n for the R-strings and R ∝ n 2 for the tubes, respectively.
We can define a transverse size R T of this tube-like string similar to Eq.(2.9) with K T = (f ′ (ρ)) 2 and observe R T = 2.06 for the minimal winding tube by a numerical calculation.
Ratios R T /2n 2 with higher winding solutions are listed in Fig.4. It suggests that the above approximation works well.

Metastable R-tube
In section 2.1, we studied the single field SUSY breaking model as a toy model of one of the O'Raifeartaigh models discussed in [12] in which classical pseudo-moduli space is stable everywhere. In this section, we move on to phenomenologically interesting situation where pseudo-moduli space has a tachyonic direction, in which large gaugino masses are generated by gauge mediation [20,38]. In such models, the R-symmetry breaking vacuum is metastable, thus the R-string solution can be a tube-like domain wall with winding number as showed in section 2.2.
Here, we study the illustrating supersymmetric model with two superfields, X and φ. These superfields have R-charges, R[X] = 2 and R[φ] = 0, and the superpotential is given by In addition, we consider the following effective Kähler metric, This model has the Z 2 symmetry, under which X and φ are Z 2 even and odd, respectively. This model has the discrete SUSY vacua, and the SUSY breaking vacuum, where the U(1) R symmetry is also broken. The former is the true vacuum, while the latter is the metastable vacuum whose vacuum energy is V = µ 4 V(1) = µ 4 1 − 1 4λ . For later convenience, let us introduce dimensionless variables by then the Lagrangian is of the form This Lagrangian is characterized by two dimensionless parameters λ and ǫ. For instance, in the SUSY breaking vacuum (T, s) = (1, 0), dimensionless masses for T and s are wriiten by respectively, and that is, existence of the SUSY breaking vacuum requires If one is interested in vacuum selection, a simple criterion is a ratio of tachyonic masses at the origin (T, s) = (0, 0) where we have Since in the early universe field values are assumed to be around the origin, if tachyonic mass of T is larger, one may expect that supersymmetry breaking model is preferable 2 . An inequality is required for selecting the SUSY breaking vacuum. Now we are ready to construct the R-tube in this two-scalar model. To this end, we make the Ansatz Similar to the model in section 2.2, it is the symmetric phase inside the tube while it is the asymmetric phase outside the tube. A sharp contrast among two models is that the outside is the true vacuum in section 2.2 and is the false vacuum in this section. One may guess that stable solution does not exist since the core of the tube has lower energy than its outside and hence larger radius would be favored energetically, which causes the "roll-over" problem.
However, because the tension of the wall acts as a centripetal force for the R-tube, we will find that there exist metastable tube-like field configurations. In order to see a typical R-tube numerical solution in this model, here we show an example in Fig. 5 with ǫ = 1, λ = 0.27. Note that the profile function of the winding field T , whose mass is very small, has a very long tail compared to that of the solution in section 2.2. On the other hand, the unwinding field s, whose dimensionless mass is of order 1, converges exponentially. In the next subsection, we will investigate stability of the R-tube by varying those parameters.

Instability of R-string and Broken Z 2 Symmetry
If we set s = 0 to keep the Z 2 symmetry, the model discussed in this section reduces to just the model discussed in section 2.1 except for an overall factor. Therefore the R-string solution (s = 0) without a hole inside is also a solution in this model. However, such an R-string would be almost always unstable and transforms into an R-tube with non-zero s inside. Since non-vanishing s means the broken Z 2 symmetry, we observe below that the Z 2 symmetry inside the R-tube in this model is almost always broken. Let us consider an infinitesimal fluctuation s(≪ 1) around the R-string solution discussed in section 2.1 and study whether a direction along s is tachyonic or not. A linearized equation for s is given with an eigenvalue q 2 as For instance we observe tachyonic masses of s numerically for many sets of parameters {λ, ǫ} as shown in Fig.6. We therefore make a conjecture that for almost all the winding number n and the almost whole parameter region of {λ, ǫ} satisfying the inequalities (3.8). This conjecture means that the R-string with s = 0 is always unstable and a stable R-tube solution, if it exists, must have the following property s| ρ=0 = 0. (3.14) In this paper we will assume this conjecture holds and will not consider the constraints from the stability of R-string configuration.

Rough Estimation for R-tube
If the domain wall consisting the R-tube is sufficiently thin and resides at ρ = R, its total energy E per the length ∆z can be estimated as as has been done in Sec. 2.2. See Fig.7. Note that the total energy has divergence terms with an IR-cutoff Λ proportional to Λ 2 and log Λ. The former is the energy density V(1) = 1 − 1 4λ of the SUSY breaking vacuum, and the latter is that for the well-known global string tension.
Here, T wall is a (dimensionless) tension of the domain wall. This potential has a local minimum (maximum) at ρ = R tube (R max ) with , (3.16) if the dimensionless tension T wall is sufficiently large as T wall > 2n. (3.17) This is, therefore, a necessary and sufficient condition for existence of the R-tube as long as the approximation Eq. (3.15) is valid. Such configurations that satisfy the inequality R ≥ R max can not avoid to spread out toward the infinite of the space. Note that comparing a thickness L wall of the domain wall with R tube , if L wall ≪ R tube holds, the above estimation (3.15) works well and we will observe the SUSY vacuum inside the R-tube, namely Using an approximation discussed in the next subsection, we can show the lower limit of the Therefore R tube can not be very small. If R tube is comparable with L wall , a configuration of R-tube approaches one of the R-string, but s| ρ=0 keeps non-vanishing even there in almost all the cases as we discussed.

Linear Approximation for the Domain Wall
In order to estimate the transverse size of the R-tube and its stability following the discussion in the previous subsection, we need the data {T wall , L wall }. We here evaluate them assuming that the domain wall in the R-tube can be well approximated by a flat domain wall interpolating the SUSY vacuum at x = −∞ and the SUSY breaking vacuum at x = ∞. Let us consider this configuration in the following. Note that, however, there is an ambiguity for definition of T wall and profile functions for the domain wall since the flat domain wall itself is unstable. It is natural to set a relation between the total energy of the system E wall and the tension T wall of the domain wall sitting at x = x with IR-cutoff Λ ± as which gives a force (pressure) from the SUSY vacuum to the domain wall Moreover, it is natural to require for the relation, is hold near the domain wall solution. When V(1) = 0 holds, the above relation can be derived from Derrick's theorem [28]. Then, we define the tension T wall and a position x of the wall in terms of only kinetic terms K without using a potential V as (3.23) The equation (3.22) is enough to estimate data {T wall , L wall } as the following. Let us approximate the profile functions (T, s) for the domain wall by piecewise-linear functions as and (T, s) = (1, 0) for x ≥ L wall and (T, s) = (0, 1) for x < 0. By inserting this approximation to Eq.(3.22) we find that the l.h.s (r.h.s) is proportional to L −1 wall (L wall ). Note that the tension of the domain wall can be expressed as Minimizing it in terms of L wall , we get where In the purple region,R-tube with n = 1 is stable but others n ≥ 2 are unstable. In the light purple region, R-tubes with n = 1, 2 is stable but others n ≥ 3 are unstable. The region below the red line represents (3.10). The SUSY breaking vacuum is unstable in the yellow region (see (3.8)).
Here we took Λ − < 0 and Λ + > L wall . Especially we find inequality (3.28) We have used this for deriving the inequality in Eq. (3.19).
Finally, using the result T wall in the linear approximation, we show the stability condition

Numerical calculation for stability
In the previous subsection, exploiting linear approximation, we found the stability condition of the R-tube with winding number n. Here, we try to check the parameter dependence of the stability by numerical calculation. We adopt a kind of relaxation method to find a configuration of the R-tube. See Appendix A for details. Since we are interested in a parameter region close to borders of the stability of two winding numbers, so we have to treat relatively unstable configuration, which require careful analysis. Because of this complexity, we focus on a couple of examples for the numerical analysis.
As a first example, we take a parameter ǫ = 1, λ = 27/100 where according to the linear approximation, winding number n = 1 is stable but n = 2 is unstable (see Fig 8). Following the relaxation method, we take an appropriate initial function and finite relaxation time τ , then we calculate minimum energy configurations. As we show in Fig 9 and Fig 10, energy convergences of the configurations have a clear difference in two cases. Here we removed a contribution E vev of the vacuum energy density from the total energy E and calculated the following dimensionless energy and we take the IR-cutoff of the energy as Λ = 50. The configuration with n = 2 is monotonically loosing the energy and in a sufficiently late relaxation time τ , the energy decreases as a linear function, which clearly suggests instability of this configuration. On the other hand, as for the configuration with n = 1, the energy seems to converge to a constant value. This sharp difference nicely matches with the result of the liner approximation. However, it is worth noting that our numerical analysis is done with a finite precision which is appropriately chosen by reasonable calculation time. Thus, beyond our calculation precision, there may exist an unstable mode which may yield slight energy loss. Thus, as long as we use a kind of relaxation method with a fixed initial condition, it may be, in general, hard to conclude complete stability of the configuration. However, even if the small instability exists, the life-time of R-tube can be longer than the decay time of the R-tube originating from an explicit U(1) R breaking effect. In many phenomenological models, the global R-symmetry is already broken by adding gravity due to the constant term in superpotential. Thus, at a point of the early universe, R-tubes disappear by generating axion domain walls [31,[39][40][41][42][43]. Therefore, as long as the stability is long enough compared with its lifetime, we can treat the R-tube as a stable solution.
As a second example, we choose ǫ = 1/50 and λ = 6/10. Small ǫ is favorable in model building, partially because the longevity of the false vacuum. So from phenomenological point of view, it is important to study a configuration of R-tube with small winding number in this parameter region. Again, using the relaxation method, we numerically calculate the energy convergence of two cases, n = 1 and n = 2. As shown in Fig 11, the total energy of R-tube with n = 2 converges to a constant value. Thus, within our calculation accuracy, the configuration looks stable. Also, the confituration with n = 1 is similarly stable. With these numerical results and linear approximation shown in the previous subsection, it may be plausible that in small ǫ region R-tubes are relatively stable and the roll-over process does not occur.  Figure 11: Energy against the relaxation time τ . Energy convergence of the configuration with winding number n = 2 for ǫ = 1/50, λ = 6/10.

Effective potential for light mode
As emphasized above, there may exist a very light mode which may cause instability of a configuration. Although treatment of such light modes in the relaxation method is not an easy task, but we would like to propose a method to uncover the existence of a light mode.
The most interesting mode is a fluctuation of the size of the R-tube. Generally speaking a zero mode (moduli) around a solution is frozen in the relaxation method, and a light mode of which dependence in the total energy is quite small, seems to be almost frozen even if it exists. To detect such light mode and search the true stable solution, we need to take a lot of different initial conditions for the relaxation method. To be concrete, we show initial conditions for the fields T, s, 3 With various values for ρ 0 which roughly indicates a transverse size of R-tube, we calculate minimum energy configurations with finite relaxation time. For any value of ρ 0 , the energy converges like Figs. 9 and 11 for those values of n, ǫ and λ. However, final configurations can have small differences of the energy and the tube-size. To represent the size of the tube, it would be useful to introduce the following definition similar to (2.9), (3.31) Here we defined two sizes of the R-tube, R T and R s . A reason for introducing two sizes can be seen in a discrepancy between the linear approximation in section 3.4 and numerical results shown below. As has mentioned in section 2.1, these quantities do not include the cut-off dependence and are well-defined candidates for the size of the tube.
As an example, we take λ = 27/100, ǫ = 1. Varying the initial position ρ 0 , we calculate the minimum energy configuration with finite relaxation time. First of all, we show a correspondence between the initial condition ρ 0 and the size R s in Fig 12. R s is evaluated with a converged configuration. Since it is one-to-one correspondence, varying the initial condition ρ 0 represents varying the size of the R-tube. Now let us show the low energy effective potential for the fluctuation mode of the size. We plot the energy (3.29) at the relaxation time τ = 50, which we will denote as T tube ≡ 2πE(τ =  Figure 12: One to one correspondence between the initial condition ρ 0 and the size R s for ǫ = 1, λ = 27/100. T tube ≡ 2πE(τ = 50). Gradients in the energy with respects to τ at τ = 50(the right panel).
Finite gradients indicate instability of configurations and noises for small R s imply the limit bound of precision of calculations.

50)
, with respect to the position of R s in Fig.13. Surprisingly, we observe a monotonically decreasing potential in terms of R s for the model with λ = 27/100, ǫ = 1 which we explained. A large plateau with tiny gradient in Fig.13 is consistent with Fig.9 where R s can be almost regarded as a massless mode. Fig.13 implies that the R-tube in this case is unstable and will expand to the infinity. This clearly suggests that the border line of instability of minimum winding R-tube shown in Fig 8 does not match with the numerical analysis above. We find that for small R s the two sizes behave differently as shown in Fig.14 and this fact tells us why the estimation (3.15) with a single size shown in section 3.4 does not work for small R s . It would be interesting to study two-scale linear approximation for better approximation.
Here, we proposed a way to analyze the low energy effective theory for a very light mode  Figure 14: Relation between the two transverse sizes R T , R s with λ = 27/100, ǫ = 1. R T is decreasing for small R s whereas R T for large R s is proportional to R s . Therefore R T is not a good quantity for parametrizing the potential.
corresponding to the fluctuation of the tube-size and checked the instability of the mode. However, this light mode is significantly affected by various corrections such as thermal effects, supergravity effects and quantum corrections. We would study these corrections elsewhere.

Bamboo solution: Tube junction
In the previous section, we showed that an R-string can form a tube-like domain wall with winding number and inside of the wall is in the SUSY preserving vacuum. It would be wonderful if we could extract any evidences of the existence of the SUSY vacuum from outside of the tube 4 . Toward this goal, in this section we demonstrate that quantum number in the SUSY vacua significantly affects the shape and the tension of the string. Concretely, we construct a junction of the R-tube. To the best of our knowledge, this kind of soliton has not been known so far. In order to demonstrate an explicit solution, let us again take the two-scalar model in Eq. (3.6). As shown in Eq. (3.11), reflecting the Z 2 symmetry of the model, there are two different R-tubes 5 . The one has γ = +1 and the other has γ = −1. The skin of the R-tube, namely the profile of T field, is independent of the choice of γ, so that one can naturally imagine that the two R-tubes with different γ can be smoothly connected. The This R-bamboo solution may be created when two R-tubes collide. If the two tubes are different kind, the domain wall must be created at the junction of the two tubes. At the same time, the anti R-bamboo may be created. This is very similar phenomenon to monopole and anti-monopole creation associated with the non-Abelian string reconnection. This is interesting issue but is beyond the scope of this paper, so we leave it as a future work. Finally, it is worthy to note that stability of bamboo configuration is not guaranteed by our numerical approach. As emphasized in the previous section, our analysis is a kind of relaxation method. Actually, if there is no domain wall inside the R-tube, the configuration seems to have a small instability shown in Fig 13. However, in the bamboo configuration, the existence of the domain wall increases the stability because the energy of the domain wall is proportional to R 2 . Thus, if the number of the domain walls is large, the bamboo

Conclusion
In this paper, we have demonstrated a fascinating role of a cosmic R-string/R-tube by using two toy models with spontaneous R-symmetry breaking. The first example shown in section 2.1 is a single field model. The model can be regarded as a toy model of one of O'Raifeartaigh models studied in [12] in which pseudo-moduli space is stable everywhere. In this model, string-like defect generated by the Kibble-Zurek mechanism is stable and very close to the known global strings. Winding number dependence of the size of the string is linear in n. On the other hand, in the two-field model shown in section 3, a string-like object is a tube-like domain wall interpolating a false vacuum and a true SUSY vacuum. One specific feature of two-field model is existence of a tachyonic direction at R-symmetry restoring point X = 0. Because of this, the core of the string is not stable under fluctuation toward the tachyonic direction. Thus, inside the core in which R-symmetry is restored, can be filled out by the true vacuum through the tachyonic direction. Naively, one may think that such an R-tube is unstable. However, as we shown in the main text numerically, there exist metastable Rtube solutions for certain parameters. An interesting property of such R-tube is the winding number dependence of the size. By the linear approximation, we estimated the dependence and found that its dependence is n 2 rather than n. Therefore, there is a tendency that tube configuration with larger winding number is more unstable. Numerically, we checked the higher winding instability at a sample parameter ǫ = 1 and λ = 27/100: We showed that the configuration with winding number n = 2 is unstable and the total energy decreases monotonically.
If an unstable R-tube is created by the Kibble-Zurek mechanism, it rapidly expands and our universe will be completely filled by the true vacuum. This process gives constraints for models building. However, it is worthy to emphasize that such roll-over process can be protected by D-term contribution or thermal potential. As is demonstrated in [21,45,46], when a D-term contribution cannot be negligible, it can lift the tachyonic direction and stabilize the pseudo-moduli space. In such models, the roll-over process does not occur. Also, when the amplitude of (tachyonic) messenger mass at the origin is sufficiently smaller than that of R-symmetry breaking field, the vacuum selection is successfully realized by exploiting the thermal potential or Hubble induced mass. If such thermal potential keeps lifting the tachyonic direction until the time of the R-string decay [31,[39][40][41][42][43], the roll-over process can be successfully circumvented. However, even if such an early stage scenario is assumed, there are sever cosmological constraints on R-axion density as studied in [31].
Finally we comment on our numerical analysis done in section 3. Since we adopted a kind of relaxation method with a fixed initial condition to find energetically minimum configurations, it is not easy to conclude the stability of the configuration sharply. There may exist a very light mode which does not change the energy significantly. To uncover the existence of the very light mode, we proposed a method for studying the low energy effective potential for such light mode. By changing the initial conditions and observing converged configurations, we can estimate the effective potential. Using this method, at the parameter ǫ = 1 λ = 27/100, we find emergence of a light model and find a large plateau in the effective potential. It would be important to apply the method to a wide range of parameter space and study borders of instabilities of configurations with various winding numbers numerically. This is beyond the our scope, so we will leave it as a future work.

A Relaxation method
In this section we review a kind of relaxation methods applied to constructing numerical solutions in this paper. Let us consider the following general Lagrangian for scalars φ α in where g αβ is the metric for the field space. Our goal is to find an numerical static solution φ α (t, x i ) = φ α sol (x i ) for this system. The shooting method is a good strategy for a system with a single scalar fields and a single spatial coordinate x 1 , but only in that case. For systems with multi fields in higher dimensions, the shooting method does not work very well and we need the relaxation method explained bellow. Let us introduce a 'relaxation time' τ instead of the real time t and suppose that φ α depend on τ as φ α (τ, x i ). Then τ dependence of φ α is defined by, with Neumann condition at the boundary of the region Σ ⊂ R d The added term in the r.h.s of Eq.(A.2) i the so-called friction term. Actually, due to this term we can show that the ordinary total energy E with integral region Σ is no longer constant but a monotonic decreasing function E = E(τ ) of τ as If the energy has the global minimum, therefore, we can obtain, at least, one static solution φ α sol by using the relaxation method. If there exists multiple local minima of the energy, an initial condition of φ α chooses one of them.
In actual numerical computation we need a cutoff of the relaxation time at τ = τ fin and we regard φ α (τ fin , x i ) as a solution. A relation between τ fin and precision of a solution φ α = φ α (τ fin , x i ) can be discussed in the following. If deviations of φ α from a solution φ α sol are efficiently small, they can be expanded as f n (x i )a n (τ ), (A. 6) with profile functions f n (x i ) for massive modes of mass m n around the configuration φ α sol . Eq.(A.2) gives development of coefficients a n (τ ) as a n (τ ) = a 0 n e −m 2 n τ , (A.7) and behavior of the energy is controlled by the lightest mass m 1 as Here a constant A ∈ R >0 depends on an initial condition we took and is assumed to be the same order as E sol . To get precision 10 −p , therefore, we have to take a time τ fin for this relaxation method as where m 1 can be roughly guessed by a typical mass scale of the system. With large τ , we sometimes observe a random behavior of E(τ ) which is a signal of the limit bound of machine precision. See Fig.17 for an example. If a configuration of φ α is accidentally near to a saddle point, we also observe an exponential decay of E(τ ), but after that it collapses like a waterfall as with a tachyonic massm 2 = −|m 2 |. Therefore an exponential behavior of the energy do not always guarantee that a stable solution is obtained. Taking multi initial conditions of φ α can avoid this technical error.  Figure 17: Typical behavior of an energy in the relaxation method. Here we get accuracy of 10 −7.5 for a solution. Calculations after τ ≈ 32 turns out to be meaningless.