Metastable vacuum decay and $\theta$ dependence in gauge theory. Deformed QCD as a toy model

We study a number of different ingredients related to the $\theta$ dependence, metastable excited vacuum states and other related subjects using a simplified version of QCD, the so-called"deformed QCD". This model is a weakly coupled gauge theory, which however preserves all the relevant essential elements allowing us to study hard and nontrivial features which are known to be present in real strongly coupled QCD. Our main focus in this work is to test the ideas related to the metastable vacuum states (which are known to be present in strongly coupled QCD in large $N$ limit) in a theoretically controllable manner using the"deformed QCD"as a toy model. We explicitly show how the metastable states emerge in the system, why their life time is large, and why these metastable states must be present in the system for the self-consistency of the entire picture of the QCD vacuum. We also speculate on possible relevance of the metastable vacuum states in explanation of the violation of local $\cal{P}$ and $\cal{CP}$ symmetries in heavy ion collisions.


I. INTRODUCTION AND MOTIVATION
A study of the the QCD vacuum state in the strong coupling regime is the prerogative of numerical Monte Carlo lattice computations. However, a number of very deep and fundamental questions about the QCD vacuum structure can be addressed and, more importantly, answered using some simplified versions of QCD. In the present paper, we study a set of questions related to metastable vacuum states and their decay to the true vacuum state using the so-called "deformed QCD" toy model wherein we can work analytically. This model describes a weakly coupled gauge theory, which however preserves many essential elements expected for true QCD, such as confinement, degenerate topological sectors, proper θ dependence, etc. This allows us to study difficult and nontrivial features, particularly related to vacuum structure, in an analytically tractable manner.
The fact that some high energy metastable vacuum states must be present in a gauge theory system in the large N limit has been known for quite some time [1]. A similar conclusion also follows from the holographic description of QCD as originally discussed in [2]. The fundamental observation on the emergence of these excited vacuum states was made in a course of studies related to the resolution of the U (1) A problem in QCD in the large N limit [3][4][5]. In the present work we do not introduce quarks (which play an essential role in the formulation of the U (1) A problem) into the system, but rather, study pure gluodynamics, and the metastable vacuum states which occur there. Nevertheless, the key object relevant for the resolution of the U (1) A problem, the so-called topological susceptibility χ, still emerges in our discussions in pure gluodynamics because it plays an important role in understanding of the spectrum of the ground state and multiple metastable states. Indeed, the topological susceptibility is defined as χ(θ) = ∂ 2 Evac(θ) ∂θ 2 . Therefore, the information about the ground (or in general metastable) states E vac (θ) is related to the θ behaviour of the system formulated in terms of the topological susceptibility χ(θ).
When some deep questions are studied in a simplified version of a theory, there is always a risk that some effects which emerge in the simplified version of the theory could be just artifacts of the approximation, rather than genuine consequences of the original underlying theory. Our study using the "deformed QCD" as a toy model is not free from this potential difficulty with misinterpretation of artifacts as inherent features underlying QCD. Nevertheless, there are few strong arguments suggesting that we indeed study some intrinsic features of the system rather than some artificial effects. The first argument is discussed in the original paper on "deformed QCD" [6] where it has been claimed that this model describes a smooth interpolation between strongly coupled QCD and the weakly coupled "deformed QCD" without any phase transition. In addition, there are a few more arguments based on our previous experience with the "deformed QCD" model, see below, which also strongly suggest that we indeed study some intrinsic features of QCD rather than some artifact of the deformation.
Our arguments are based on the computation [7] of the contact term in the "deformed QCD", see also [8] with some related discussions. The key point is that this contact term with a positive sign (in the Euclidean formulation) in the topological susceptibility χ is required for the resolution of the U (1) A problem [3][4][5]. At the same time, any physical propagating degrees of freedom must contribute with a negative sign, see [7] with details. In [3] this positive contact term has been simply postulated while in [4,5] an unphysical Veneziano ghost was introduced into the system to saturate this term with the "wrong" sign in the topological susceptibility. This entire, very non-trivial picture, has been successfully confirmed by numerical lattice computations. More importantly for the present studies, this picture has been supported by analytical computations in "deformed QCD" in which all the nontrivial crucial elements for the reso-lution of the U (1) A problem indeed emerge in analytical analysis.
Indeed, the non-dispersive contact term in topological susceptibility can be explicitly computed in this model and is given by [7] where q(x) is the topological density operator. It has the required "wrong sign" as this contribution is not related to any physical propagating degrees of freedom, but is rather related to the topological structure of the theory, and has a δ(x) function structure as it should. In this model χ is saturated by fractionally charged weakly interacting monopoles describing the tunnelling transitions between topologically distinct, but physically equivalent topological winding sectors. Furthermore, the δ(x) function in (1) should be understood as total divergence related to the infrared (IR) physics, rather than to ultraviolet (UV) behaviour as explained in [7] The singular behaviour of the contact term has been confirmed by the lattice computations where it has been found that the singular behaviour at x → 0 is an inherent IR feature of the underlying QCD rather than some lattice size effect [9][10][11][12].
In addition, one can explicitly see how the Veneziano ghost postulated in [4,5] is explicitly expressed in terms of auxiliary topological fields which saturate the contact term (1) in this model [13]. In other words, the η field in this model generates its mass (which is precisely the formulation of the U (1) A problem) as a result of a mixture of the Goldstone field with the topological auxiliary field governed by a Chern-Simons like action, see [13] for the details.
All these features related to the θ dependence which are known to be present in strongly coupled regime also emerge in the weakly coupled "deformed QCD" toy model. Therefore, we interpret such behaviour as a strong argument supporting our assumption that the "deformed QCD" model properly describes, at least qualitatively, the features related to the θ dependence and vacuum structure of QCD, including the presence of metastable states which is main subject of the present work.
The specific computations we perform related to the metastable vacuum states have never been performed using numerical lattice (or any other) methods. Therefore, we do not have the same level of luxury present in our previous studies of the contact term [7] in which our results were supported by numerous lattice computations. Nevertheless, as the specific questions about the metastable states are closely related to much more generic studies of the θ dependence in the system, as reviewed above, we are still confident that our results presented below, based on the "deformed QCD" model, are inherent qualitative properties of QCD rather than some artificial effects which may occur due to the deformation.
Our presentation is organized as follows. We start in section II by reviewing a simplified ("deformed") version of QCD which, on one hand, is a weakly coupled gauge theory wherein computations can be performed in theoretically controllable manner. On other hand, this deformation preserves all the elements relevant to our study such as confinement, degeneracy of topological sectors, nontrivial θ dependence, presence of non-dispersive contribution to topological susceptibility, and other crucial aspects pertinent to the study of the metastable states. In section II B we explicitly demonstrate the presence of metastable states in this model. In section III we review the general strategy to compute a decay of metastable vacuum states to the true vacuum in the path integral formulation. Finally, in section IV we present our numerical analysis on the life time of the metastable states as a function of a "semi-classicality" which is a parameter determining the region of validity of our semiclassical computations. We conclude in section V with speculations on possible consequences and manifestations of our results for physics of heavy ion collisions where a metastable state might be formed as a result of collision, and the system, which is order the size of a nuclei, might be locked in this state for sufficiently long period of time ∼ 10 fm/c.

II. DEFORMED QCD
Here we overview the "center-stabilized" deformed Yang-Mills developed in [6]. In the deformed theory an extra "deformation" term is put into the Lagrangian in order to prevent the center symmetry breaking that characterizes the QCD phase transition between "confined" hadronic matter and "deconfined" quark-gluon plasma, thereby explicitly preventing that transition. Basically the extra term describes a potential for the order parameter The basics of this model are reviewed in section II A, while in section II B we classify the metastable states which is inherent element of the system.

A. The Model
We start with pure Yang-Mills (gluodynamics) with gauge group SU (N ) on the manifold R 3 × S 1 with the standard action and add to it a deformation action, built out of the Wilson loop (Polyakov loop) wrapping the compact dimension The parameter L here is the length of the compactified dimension which is assumed to be small. The coefficients of the polynomial P [Ω(x)] can be suitably chosen such that the deformation potential (4) forces unbroken symmetry at any compactification scales. At small compactification L the gauge coupling is small so that the semiclassical computations are under complete theoretical control [6]. As described in [6], the proper infrared description of the theory is a dilute gas of N types of monopoles, characterized by their magnetic charges, which are proportional to the simple roots and affine root α a ∈ ∆ aff of the Lie algebra for the gauge group U (1) N . For a fundamental monopole with magnetic charge α a ∈ ∆ aff (the affine root system), the topological charge is given by and the Yang-Mills action is given by The θ-parameter in the Yang-Mills action can be included in conventional way, withF µν ≡ µνρσ F ρσ . The system of interacting monopoles, including the θ parameter, can be represented in the dual sine-Gordon form as follows [6,7], where ζ is magnetic monopole fugacity which can be explicitly computed in this model using the conventional semiclassical approximation. The θ parameter enters the effective Lagrangian (9) as θ/N which is the direct consequence of the fractional topological charges of the monopoles (6). Nevertheless, the theory is still 2π periodic. This 2π periodicity of the theory is restored not due to the 2π periodicity of Lagrangian (9). Rather, it is restored as a result of summation over all branches of the theory when the levels cross at θ = π(mod 2π) and one branch replaces another and becomes the lowest energy state as discussed in [7]. Finally, the dimensional parameter which governs the dynamics of the problem is the Debye correlation length of the monopole's gas, The average number of monopoles in a "Debye volume" is given by The last inequality holds since the monopole fugacity is exponentially suppressed, ζ ∼ e −1/g 2 , and in fact we can view (11) as a constraint on the region validity where semiclassical approximation is justified. This parameter N is therefore one measure of "semi-classicality". For our studies in what follows it is convenient to express the action in terms of dimensionless variables. We rescale x as follows x = x /m σ such that x becomes a dimensionless coordinate. All distances now are measured in units of m −1 σ . With this rescaling the potential term is explicitly proportional to the parameter of semiclassicality N . The coefficient on the kinetic term in the above action (9) also esquires the same factor N such that the action (9) assumes a very nice form: with σ N +1 identified with σ 1 . In formula (12) we used x rather than x to simplify notations. The Lagrangian entering the action (12) is then dimensionless with a prefactor N , such that the rest of the action depends on N (in the number of fields), but no longer depends on g, the gauge coupling, or L, the compactification scale. This is the form of the action we use in our calculations. Note that the large N limit in the so called "double scaling" limit had been discussed previously in [6]. It has been argued that in this limit we can consider any finite N but cannot consider the formal limit N → ∞ since the parameter space shrinks to a point. This will not matter for the present work as we always carry out the computations for large, but finite N .

B. Metastable Vacuum States
Here we concentrate on the Euclidean potential density for the σ fields at θ = 0, where again σ N +1 is identified with σ 1 . To simplify notations we skip a large common factor N in our discussions which follow. We restore this factor in our final formula. Also, we have added a constant (N ) so that the potential is positive semi-definite. The lowest energy state, denoted by σ (−) , is the state with all σ fields sitting at the same value (σ n = σ n+1 ) and has zero energy. This is clearly the true ground state of the system, but there are also potentially some higher energy metastable states.
For an extremal state we must have ∂U ∂σ n = 0 (14) for all n, which gives immediately A necessary condition for a higher energy minimum of the potential is thus that the σ fields are evenly spaced around the unit circle or (up to a total rotation), where m is an integer. A sufficient condition is then again for all n. This gives us which using (16) gives So, we get a constraint on m in the form of (19), and also on N . From (19) it is quite obvious that metastable states always exist for sufficiently large N , which is is definitely consistent with old and very generic arguments [1]. In our simplified version of the theory one can explicitly see how these metastable states emerge in the system, and how they are classified in terms of the scalar magnetic potential fields σ(x). One should also remark here that a non-trivial solution with m = 0 in (19) does not exist 1 in this simplified model for the lowest N = 2, 3, 4. Therefore, in our study we always assume N ≥ 5.
Looking back at the potential (13), the lowest energy of the possibilities are given by m = ±1, so that the lowest energy metastable states, denoted by σ (+) , are given by (again up to a constant rotation) To understand the physical meaning of the solutions describing the nontrivial metastable vacuum states, we recall that the operator e iαa·σ(x) is the creation operator 1 N = 4 deserves a special consideration as at m = ±1 the second derivative (17) vanishes. It may imply a presence of the massless particles in the spectrum for these excited vacuum states. It may also correspond to a saddle point in configuration space. We shall not elaborate on this matter in present work.
for a monopole of type a at point x, as it was explicitly demonstrated in [7], Therefore, the vacuum expectation value M a (x) describes the magnetization of a given metastable ground state classified by the parameter m. As one can see from (16), the corresponding vacuum expectation value M a (x) always assumes the element from the centre of the SU (N ) group. Specifically, for the lowest metastable vacuum states given by (20), the magnetization is given by The fact that the confinement in this model is due to the condensation of fractionally charged monopoles has been known since the original paper [6]. Now we understand the structure of the excited metastable states also; mainly, these metastable vacuum states can be also thought of as a condensate of the monopoles. However, the condensates of different monopole types, n from (20), are now shifted by a phase such that the corresponding magnetization receives a non-trivial phase (22). A different, but equivalent way to classify all these new metastable vacuum states is to compute the expectation values for the topological density operator for those states. By definition The imaginary i in this expression should not confuse the readers as we work in the Euclidean space-time. In Minkowski space-time this expectation value is obviously a real number. A similar phenomenon is known to occur in the exactly solvable two dimensional Schwinger model wherein the expectation value for the electric field in the Euclidean space-time has an i, see [14] for discussions within present context. The expectation value (23) is the order parameter of a given metastable state.
The crucial point we want to make here is that a metastable vacuum state with m = 0 in general violates P and CP invariance since the topological density operator itself is not invariant under these symmetries. Precisely this observation inspires our suggestion, to be discussed in conclusion, that such metastable states could be the major source of the local P and CP violation observed in heavy ion collisions at the Relativistic Heavy Ion Collider (RHIC) at Brookhaven, and the Large Hadron Collider (LHC). Now we come back to our discussions of the lowest metastable states (20). Putting the metastable configuration back into the potential (13) we find that the energy density separation between the true ground state and the lowest metastable states (20) is given by 2 The choice of sign in (20) is irrelevant for the purposes of calculating the vacuum decay since the two states m = ±1 are degenerate in terms of energy and have the same energy splitting with respect to the ground state. These states, however, are physically distinct as the expectation value of the gauge invariant operator (23) has opposite signs for these two metastable vacuum states. This implies that all P and CP effects will have the opposite signs for these two states, while the probability to form these two metastable states is identical, as is the decay rate. So, while our fundamental Lagrangian is invariant under these symmetries, the metastable vacuum states, if formed, may spontaneously break that symmetry.

III. METASTABLE VACUUM DECAY
In this section we briefly review the general theory and framework for calculating metastable vacuum decay rates in Quantum Field Theory. For a more thorough discussion see [15] [16]. The process for the decay of a metastable vacuum state to the true vacuum state is analogous to a bubble nucleation process in statistical physics. Considering a fluid phase around the vaporization point, thermal fluctuations will cause bubbles of vapor to form. If the system is heated beyond the vaporization point, the vapor phase becomes the true ground state for the system. Then, the energy gained by the bulk of a bubble transitioning to the the vapor phase goes like a volume while the energy cost for forming a surface (basically a domain wall) goes like an area. Thus, there is some critical size such that smaller bubbles represent a net cost in energy and will collapse while larger bubbles represent a net gain in energy. Once a bubble forms which is larger than the critical size it will grow to consume the entire volume and transition the whole of the sample to the vapor phase. To understand the lifetime of such a 'superheated' liquid state, the important calculation is, therefore, the rate of nucleation of critical bubbles per unit time per unit volume (Γ/V ). Similarly, we aim to calculate this decay rate for our system with from the metastable state σ (+) to the ground state σ (−) , though through quantum rather than thermal fluctuations. Consider a general system with a ground state field configuration, φ (−) , and metastable field configuration, 2 One should comment here that the vacuum energy of the ground state E (±) ∼ N in this model scales as N in contrast with conventional N 2 scaling in strongly coupled QCD. However, the ratio /E (±) ∼ N −2 shows the same scaling as in strongly coupled QCD.
Energy density φ-field configuration Qualitative diagram for potential A very qualitative picture for the potential of a general system with a global ground state, φ (−) , and a higher energy metastable state, φ (+) , with an energy splitting between the two given by . φ (+) , with an energy density difference between the two given by . Qualitatively the potential for the system should be understood as something like FIG. 1. Classically, a system in the configuration φ (−) is stable, but quantum mechanically the system is rendered unstable through barrier penetration (tunneling).
The semiclassical expression for the tunneling rate per unit volume is given by [15] where the Euclidean action, S E , is the action upon analytically continuing to imaginary time and is given by We have explicitly left in (25) to emphasize the semiclassical expansion. The action, (26), in the exponent of (25) is evaluated in the field configuration called the "Euclidean bounce" which we have denoted φ b . The Euclidean bounce is a finite action configuration which solves the classical equations of motion and interpolates, in Euclidean time, from the metastable state to a configuration "near" the ground state and back. Making reference to the potential depicted in FIG. 1, continuing to Euclidean time essentially describes a system with the sign of the potential flipped. As such, the bounce describes a path starting at the now local maximum φ (+) at t → −∞ rolling down into the valley and up to the classical tuning point near the higher peak φ (−) , then reversing and traveling back to φ (+) at t → +∞. In order for the action to be finite the bounce must also tend to φ (+) as the spacial coordinates go to infinity in any direction.
Additionally, we have glossed over one technicality by representing φ as a single dimension while in fact it is not. In principle there is a classical turning surface, call it Σ, rather than a single point and so there may be many paths from the peak at φ (+) to the surface Σ. The resolution however is straightforward. Each such path contributes as (25) and so the path of minimal action is the dominant path. For details see [17]. Furthermore, the minimal action path is spherically symmetric so that the action can be written in terms on a single radial dimension, Thus, the bounce we should find is the minimal action configuration which solves the equation of motion, subject to the boundary conditions φ → φ (+) as ρ → ∞ and dφ/dρ → 0 as t → 0.
In the limit of small separation energy the bounce approaches the peak φ (−) more closely and spends longer in the region around the peak, so that the bounce configuration resembles a bubble with the interior at φ (−) , the exterior at φ (+) , and a domain wall surface interpolating between the two 3 . If the bubble is very large, corresponding to very small , then the curvature at the interpolating surface is small and the surface appears flat. Alternately, simply note that the second term in (28) goes like 1/ρ, so if the fields only change appreciably around a thin surface at large ρ, the second term can be neglected and the equation of motion reduces further to the much simpler form Therefore, if the separation energy, , between the two states is small, we need only solve for the one dimensional soliton interpolating between φ (+) and φ (−) which solves (29). This is called the thin-wall approximation, and is the framework in which we will work in the next section. In the deformed model discussed in the previous section, the separation ∼ 1/N , so that the thin-wall approximation coincides with the large N approximation. For the thin wall approximation the full action action reduces to where S 1 is the one dimensional action across the domain wall given by What remains then is to determine the size of the bubble, R. The stipulation that the bounce configuration describes a classical path implies that it extremizes the action (30). Thus, by variation, which yields R = 3S 1 / . Notice again the similarity to a bubble nucleation problem. This extremal action with respect to the bubble size is in fact a maximum, and as such the action increases with R for smaller size and decreases with R for larger. Hence, the bounce configuration which saturates the decay rate is essentially a bubble of critical size as discussed when making this analogy to bubble nucleation. We now have everything required to calculate the exponent for the vacuum decay (25) assuming we can solve for a classical path associated with the one dimensional action (31) interpolating between the two states φ (+) and φ (−) . We have not discussed the coefficient A, and indeed it is a much more complicated problem related to the functional determinant of the full differential operator, δ 2 S/δσ 2 . This calculation is beyond the scope of this work for the deformed model as we are primarily interested in the dependence on N and we expect the determinant factor to vary much more slowly than the exponential. The only strong dependence we expect for the determinantal factor is through a factor of S E /2π for each zero mode [15].

IV. COMPUTATIONS
We must now translate the results of the previous section for a general four dimensional system to the particular three dimensional effective theory for deformed gauge theory discussed in section II A. As discussed, we shall look for a one dimensional classical path of minimal action interpolating between σ (−) and σ (+) described by the 1d action, setting θ = 0, The 3d action in the thin wall approximation, analogous to (30) is then given by where the last step is again computed by using variational analysis to get R = 2S 1 / .
The condition for the validity of the thin wall approximation is essentially that the interior of the bubble is very near the true ground state σ (−) so that it is nearly stable and stays near σ (−) for large ρ. We want Rµ 1, where µ 2 = ∂ 2 U/∂σ 2 n (σ (−) ) is the curvature of the potential at the ground state; here µ = √ 2. Thus, we need where, again, So, we should solve the equations of motion with σ N +1 identified with σ 1 , derived from the action (33) subject to the boundary conditions and The ϕ in (39) is a relative rotation angle between the two boundaries, since each of the two states are only defined up to a rotation. The angle is determined by demanding a minimal action interpolation. That is, we should minimize the action with respect to the interpolating field configuration and also with respect to this angle ϕ. The final solution thus obtained will then be defined only up to an arbitrary total rotation which will be important later. Additionally, we expect the solution to be a soliton (instanton-like) in the sense that it should be well contained with only exponential tails away from the center so that the interpolation occurs in an exponentially small region. The characteristic size of this region is given by m −1 σ in the original model, or just 1 in dimensionless notations used here. Then, we can calculate the vacuum decay rate as where we have put in the part of the coefficient that we can calculate related to the zero modes in the system. There are six zero modes: three translations, two spacial rotations, and the one global σ-rotation discussed earlier.
We now discuss in section IV A the numerical technique employed to solve this problem, and in section IV B the results of these numerical calculations.

A. Numerical Technique
Numerical techniques can be broadly classified as Iterative or Explicit Finite Difference (EFD). EFD methods involve approximating the derivatives as an N × N matrix and then inverting the latter only once to find the solution. They work well for simple linear problems where there is a unique non-trivial solution to the equations of motion. In iterative methods, we first define an initial guess for the solution and then relax that solution until the residual is below a certain user defined precision. This is a better technique for non-linear systems where several non-trivial solutions may exist, and is therefore the method we employ here. The equations of motion (37) are a coupled version of the sine-Gordon equation and are quite non-linear.
The sine-Gordon equation for a single field, u = sin(u), has a soliton solution given by interpolating between 0 and 2π, which seems like good starting point. As such, we choose a similar form for the initial guess at the solution for the coupled equations, and hence define our initial guess to be of the form This initial guess has two important properties; it satisfies the boundary conditions (38) and (39), and tails off toward those boundaries as decaying exponentials for x → ±∞, which is the type of behaviour expected, as discussed previously. The equations of motion (37) are on an infinite domain and must be truncated to be solved numerically. We want to truncate the domain to a region beyond which changes in the σ i (r) are numerically insignificant. Given that the tails of the σ n (r) (and we expect the final solution) are decaying exponentials, choosing the domain [− 16,16] means that the boundary values are within ∼ 10 −7 of their final values and is suitable for our purpose.
In order to promote numerical stability particularly around the boundary values we employ Chebyshev spectral methods for integrals and derivatives, as described in [19] [20], using an unevenly spaced Chebyshev grid given by where N p is the number of grid points. Notice that the Chebyshev grid defines a domain [−1, 1] and so we scale x → x 16 in order to express functions with the chosen domain on the spectral grid.
From now on, we will use the following notation: σ i n denotes the i th grid point of the n th field, where N is the total number of fields while N p is the number of grid points. The differentiation matrix (N p × N p ) is given by [19] (page 570) as where, p j = 2 j = 0 or N 1 otherwise.
Any higher derivative is then given by repeated multiplication by D. This differentiation matrix (44) is basically just the linear operator describing interpolating a function on the grid points by an N th p order polynomial and differentiating that polynomial. Since it uses knowledge of the entire function rather than just the few nearby points like a finite difference, the accuracy of the derivative is generally much better than any small order finite difference. Furthermore, using a grid spaced in this way provides much more numerical stability, counteracting the Runge phenomenon. [20] The algorithm we employ to minimize the action with respect to the field configuration is essentially to treat the action as a potential over the configuration space formed by each σ field at each grid point, then to take steps in the negative gradient direction. Essentially, iterating the expression where δ is a chosen step size, which we start as δ = 1. At each iteration, we enforce the boundary conditions and check if the action (33) applied to the new configuration is in fact smaller than the old configuration. If so, we move to the new configuration. If not, the step was too large and we have overstepped the section of the potential with a downward slope, so we go back to the old configuration and reduce the step size by δ → δ/2 and iterate this procedure until we find a good step or reduce the step size below our desired precision. We then reset δ = 1 and continue until we cannot find a good step within our desired precision. Once reaching a position from which no step reduces the action, we are within the defined precision of a minimum of the action. In order to probe more of the configuration space we took a Monte-Carlolike approach wherein we adjusted our initial guess by adding some Gaussian noise in an envelope, (1 − |x|) 2 . We chose an envelope of this form because we expect that the solution has the sort of exponential tails of our initial guess (42), while we do not know the form of the core of the domain wall. Thus, it is sensible to probe more of the configuration space related to the specific details of the core.

B. Results
The first issue we address is the question about the favoured angle, ϕ, between the two boundaries, (38) and  to arbitrary N we set The particular fit parameters are given for completeness but should not be regarded as terribly important. We are after all working in a toy model, and as such should expect a good qualitative picture but not take the details too seriously. It is however interesting that the final form for the decay rate is given as, putting the parameter N back in, with both a and b positive. Thus, the decay rate does drop off exponentially in N and our other semiclassical parameter N , and indeed faster than any perturbation term would describe as previously conjectured. It is a semiclassical calculation, but the behaviour is fundamentally non-perturbative, and it is only parametrically justified when N 1.

V. CONCLUSION
Our conclusion can be separated on two different parts: 1). solid theoretical results within the "deformed QCD" model, and 2). some speculations related to strongly coupled QCD realized in nature.
We start with the first part of the conclusion in which our basic result is as follows. We have demonstrated that the "deformed QCD" model shows (once again) that some qualitative features expected to occur in strongly coupled regime in large N limit as argued in [1] do emerge in the simplified version of the theory as well. We demonstrated the existence of metastable vacuum states with energy density higher than the ground state by ∼ 1/N , and have shown that the life time of the metastable states is exponentially suppressed in this model with respect to semi-classicality parameter N . Suppression increases even further with increasing number of colours N for a fixed N , and it is given by (48).
In this simplified system one can explicitly see these metastable states, how they are classified, and the microscopic dynamics which govern the corresponding physics. We shall not repeat this analysis, instead referring to section II B. However, the important remark we would like to make here which is relevant for the speculative portion of our conclusion is that the P and CP invariance is generally violated in these metastable vacuum states.
We believe that this feature of spontaneous breaking of the P and CP invariance in metastable states is quite a generic feature which is shared by strongly coupled pure gauge theories (for sufficiently large N ).
Therefore, we now speculate that precisely this spontaneous symmetry breaking effect is responsible for the asymmetries observed at the RHIC (Relativistic Heavy Ion Collider) and the LHC (Large Hadron Collider). To be more specific, the violation of local P and CP symmetries has been the subject of intense studies for the last couple of years as a result of very interesting ongoing experiments at RHIC [21,22] and, more recently, at the LHC [23][24][25][26]. The main idea for explaining the observed asymmetries is to assume that an effective θ( x, t) ind = 0 is induced on a sufficiently large scale as a result of collision [27]; see [28] for a recent review and introduction to the subject with a large number of references to original papers. The induced θ( x, t) ind = 0 obviously violates local P and CP symmetries on the same large scales where θ( x, t) ind = 0 is correlated. It may then generate a number of P and CP violating effects, such as Charge/Chiral Separation (CSE) and Chiral Magnetic (CME) Effects. One of the critical questions for the applications of the CME to heavy ion collisions is a correlation length of the induced θ( x, t) ind = 0. Why are these P odd domains large?
We would like to speculate that the crucial element in understanding of this key question might be related to the metastable states which are the subject of the present work. To be more specific, we suggest that the system at high temperature might be locked in one of these metastable states during the cooling stage. It then implies a spontaneous symmetry breaking of P and CP symmetries as (23) suggests. If this happens one should obviously expect a number of P and CP effects to occur coherently in the entire system characterized by a large scale of order the size of nuclei of order L Λ −1 QCD . In other words, we identify θ( x, t) ind = 0 from [27] with the effective theta parameter 2π/N which enters (23) and which manifests a spontaneous violation of the P and CP symmetries in the system.
The presence of such long range order (which itself is a consequence of a spontaneous selecting of a metastable vacuum state in the entire system) may explain why CME is operational in this system and how the asymmetry can be coherently accumulated from entire system. This identification would justify the effective Lagrangian approach advocated in [27,29] wherein θ( x, t) ind is treated as slow background field with correlation length much larger than any conventional QCD fluctuations, L Λ −1 QCD . It is important to emphasize that the P and CP symmetries are good symmetries of the fundamental QCD. The asymmetries can only be observed in heavy ion collisions in event by event analysis when the system might be locked for sufficiently long period of time τ ∼ L/c Λ −1 QCD in a metastable state in one collision with one specific sign for the topological density (23). Because the metastable states with opposite signs for the topological density operator (23) have the same energy, which state is chosen for a particular event is random and evenly distributed. Thus, it is clear that if one averages over large number of events, the asymmetry will be washed out as the probability to form these metastable states is identical and the lifetime for the two is the same as we mentioned in section II B. However, in the event by event studies the asymmetry will be evident in the system. Apparently, this is precisely what has been observed, see the recent review paper [28] for the details.