$Z_3$ meta-stable states in PNJL model

We study the $Z_3$ meta-stable states in the Polyakov loop Nambu-Jona-Lasinio (PNJL) model. These states exist for temperatures above $T_m \sim $ 362 MeV and can decay via bubble nucleation. We numerically solve the bounce equation to compute the nucleation rate. Our results suggest that, in the context of heavy-ion collisions, the likely scenario for the decay of the meta-stable states is via spinodal decomposition.


I. INTRODUCTION
In pure SU (N ) gauge theories, the energy density increases sharply across critical temperature (T c ). It is believed that this is due to deconfinement of the constituents (gluons) of low energy excitations of the theory. This transition from confined phase to deconfined state, known as confinement-deconfinement (CD) transition, has been extensively studied in the literature [1][2][3][4][5]. The CD phase transition is found to be second order for N = 2 [6][7][8][9][10][11] and, first order for N ≥ 3 [12,13]. The Polyakov loop, which transforms as a Z N spin, plays the role of order parameter. It is real valued for N = 2 and complex for N > 2. Above the critical temperature, in deconfined phase, it acquires a non-zero expectation value, spontaneously breaking the Z N symmetry. This leads to N degenerate vacua. This non-trivial nature of the deconfined state allows for the existence of topological defects such as, domain walls for N = 2, and domain walls connected by strings for N > 2 [14][15][16].
In a realistic theory such as quantum chromo dynamics (QCD), there are fermions (quarks) in the fundamental representation. The presence of these fermions lead to explicit breaking of the Z N symmetry. The strength of the explicit symmetry breaking depends on the quark masses as well as the number of quark flavors [17][18][19][20][21]. It affects the nature of the CD transition [17,19] as well as the transition temperature. For large explicit symmetry breaking, the CD transition turns into a cross-over while the transition temperature tends to decrease. Furthermore, there are no N degenerate vacua in the deconfined phase. Out of the previous N vacua, all but one becomes the ground state. With explicit symmetry breaking, the topological defects can still exist, but far above T c and most of them are time dependent (non-static) [22].
The explicit breaking of the Z N symmetry due to matter fields has been studied by calculating the the partition function or the effective potential of the Polyakov loop, when the gauge and matter field fluctuations are small [23][24][25][26]. These perturbative calculations are reliable for high temperatures (T >> T c ), i.e in the early Universe when the gauge coupling is expected to be small. Calculations which include fluctuations up to second order (one loop) show the presence of meta-stable states [20,23,25,27]. These states have been studied extensively in the context of cosmology. They are found to be long lived, and can leave observable imprints while decaying via a first order phase transition in the early Universe [27,28]. However, in the early Universe the number of effective quark flavors is larger than 3, in which case, the free energies of the meta-stable states, at one loop, are positive and hence lead to negative pressure and entropy [20]. This problem doesn't arise in QCD near the critical temperature as the number of flavors is effectively ≤ 3.
The study of Z 3 meta-stable states for small temperatures, in particular near T c is important as they may play an important role in the the evolution of quark-gluon plasma (QGP) in heavy-ion collision experiments. Near T c perturbative calculations are expected to break down due to large coupling and fluctuations. There are very few studies of Z N symmetry using non-perturbative lattice QCD simulations. Lattice QCD results for 2 flavors show that out of the previously 3 degenerate vacua only one remains stable, while the rest become meta-stable states [21]. The two meta-stable states are degenerate, related via Z 2 symmetry. Further, the meta-stability depends on the temperature, with the meta-stable states becoming unstable below 750 MeV [21]. In general, non-perturbative lattice simulations are essential for a quantitative estimate of the explicit symmetry breaking; however, the mean field approaches provide a qualitative understanding. In a recent study of Z 3 symmetry in the Polyakov loop quark meson (PQM) model [29], it is found that the meta-stable states exist above 310 MeV.
In heavy-ion collisions, the initial conditions are far from equilibrium. The system quickly thermalizes in less than a f m time. In such a scenario, it is possible that the whole or part of the system can get trapped in one of the meta-stable states [28]. Also, if the system somehow thermalizes to a state of super hot hadron gas, which is a possibility at high baryon density, it will decay through bubble nucleation and some of the bubbles will have meta-stable cores.
In the present work, the Z 3 meta-stable states are studied in the PNJL model. In this model they exist above T m ∼ 362 MeV. If such a state is present, it can either become unstable (when temperature drops below T m ) or decay through nucleation of bubbles which grow in real time converting the meta-stable state to stable state. To compute the nucleation rate, the Euler-Lagrange equation for the bubble/bounce solution [30][31][32][33] is numerically integrated. The action as well as other properties of the bounce solution are found to depend strongly on the temperature. This study finds that the likely scenario for the evolution/decay of a meta-stable state in heavy ion collision is spinodal decomposition. In such an evolution large angular oscillations of the Polyakov loop are expected.
The paper is organized as follows. In section II, Z 3 symmetry in pure SU (3) gauge theory is discussed. We briefly go through the explicit breaking of Z 3 symmetry in the PNJL model and compute the thermodynamic properties of the meta-stable states in section III.
In section IV, we present the calculation of the bounce solution. In section V, we discuss the evolution of the meta-stable states in heavy-ion collisions and present our conclusions in section VI.

II.
Z 3 SYMMETRY IN PURE GAUGE THEORY.
In path integral formulation, gauge fields which are periodic in the temporal direction only contribute to the partition function, i.e where β = 1 T . This boundary condition allows for the gauge transformations, U ( x, τ ), to be periodic up to a factor z ∈ Z N , such as Though the partition function is invariant under the above gauge transformation, the Polyakov loop transforms as a Z N spin. The Polyakov loop is defined as Here 'P' denotes path order, 'g' is the gauge coupling and A 0 = A a is related to the free energy FQ Q (r) of a static (infinitely heavy) quark-anti-quark pair at infinite separation [34].
In the following, we briefly describe the Z 3 symmetry in the effective potential for the Polyakov loop which describes the CD transition in pure SU(3) gauge theory [35]. We consider the following Landau-Ginsburg effective potential for the complex field Φ [35][36][37], Different forms of the effective potential in terms of the field Φ have been proposed [36,[38][39][40][41][42]. The Z 3 symmetry and the first order nature of the CD transition require a cubic term in the effective potential. The factor T 4 takes care of the dimension of the effective potential [35]. In the mean-field approximation the minimum (minima), Φ th , of the effective potential gives the thermal average of the Polyakov loop, L(T ). Φ th = L(T ). In this approximation the pressure P is given by, The above effective potential with the following form of b 2 (T ) with a 0 = 6.75, a 1 = -1.95, a 2 = 2.625, b 3 = 0.75, b 4 = 7.5 [38,43] and T 0 = 270 MeV [13,[44][45][46] reasonably reproduces the lattice QCD results around the transition point. Across the critical temperature T 0 , the Polyakov loop expectation value jumps discontinuously. For T > T 0 , there are three degenerate minima in the effective potential, which can be seen in Fig. 1(a), the contour plot of the effective potential in the complex Φ = (Φ 1 , Φ 2 ) plane at T = 1.5T 0 . We also plot the variation of the potential along a circle going through the three minima, i.e Φ = |Φ th |e iθ in Fig. 1(b). The three degenerate vacua are situated at θ = 0, 2π 3 and 4π 3 , related by Z 3 rotation.
When the dynamical quarks are included, the Z 3 symmetry is not exact any more. While the pure gauge part of the action is Z 3 symmetric, the quark part of the action is not invariant Here D ν is the covariant derivative, This term takes into account the interaction between the gauge and quark fields. m f and µ f are quark mass and chemical potential of quark flavor f respectively. G s is the four quark contact interaction strength. The thermodynamic potential in the mean field approximation for the above theory with two quark (u, d) flavors is given by [38,59,60], Here In the calculations G s = 5.02GeV −2 and m 0 = 5 MeV. For T 0 = 190 MeV in the effective potential U , the PNJL results are consistent with the lattice results [38,60].
For the temperature dependence of the quark condensates and expectation value of the Polyakov loop we minimize the thermodynamic potential Ω(σ, m, T 0 , T ) by numerically solving the following set of equations, Φ 1 and Φ 2 are real and imaginary parts of Φ. The numerical program requires initial trial values of Φ and σ. It updates the trial values so that the thermodynamic potential decreases. The process stops once a minimum is reached within a certain numerical accuracy.
This method can not find all the minima at once. For each minima the numerical procedure is repeated, by suitable choices of initial conditions.
One can show that, in the ground state Φ is real valued, i.e Φ 2 = 0 (θ = 0). Hence, for the ground state, we solve the above equations with initial values of |Φ| > 0, and θ = 0. We took the zero temperature value of σ as its initial value. For the meta-stable states, the Z 3 rotated values of the ground state Φ as initial value works well. Since the Z 3 symmetry is explicitly broken, Z 3 rotated Φ do not solve the equations. Though the meta-stable states are found to be close to the initial values. The value of σ in the meta-stable is found to differ from the stable one.
With above numerical procedure, in the temperature range of T c (176) MeV -T m (362) MeV, only one solution is found with Φ 1 > 0 and Φ 2 = 0. Z 3 rotated values of this solution as initial condition does not result in any new solution. As the temperature is increased, for T > 361.7 MeV, two local minima appear around θ = 2π 3 and θ = − 2π 3 . The thermodynamic potentials for these two states are found to be same, but higher than that for the ground state for which θ = 0. The values of Φ for the meta-stable and stable states are no more related by Z 3 rotation.
In Fig. 2, the thermodynamic potential versus θ has been plotted for a given value of σ and |Φ| for T =380 MeV and 570 MeV. In Fig. 2(a) the temperature is close to T m when the effective potential develops two saddle points. Comparing Fig. 2(a) and Fig. 2(b) one can clearly see signs that the barrier between meta-stable and stable states increases with temperature. These features can also be seen in Fig. 3, the contour plot of the thermodynamic potential in the Φ 1 -Φ 2 plane.    the radial direction given by r 2 = |x| 2 + τ 2 in the Euclidean space. It has been shown that the problem is equivalent to calculating the classical evolution of a particle in the Euclidean space in presence of the inverted potential −V (φ), where the particle rolls down from the stable vacuum and bounces up to the meta-stable one. The decay rate then can be written as the summation of all such "bounces" [30]. At high temperatures, owing to the periodicity of the field theory in the "time" direction, the field configurations will have O(3) symmetry on a time slice [32]. For a single scalar field theory with a meta-stable state, the bounce can then be calculated by solving the equation of motion, with the boundary conditions φ → φ m as r → ∞, where φ m is the value of the field in the meta-stable state. For r ∼ 0 the field is expected to be close to the stable state. If r were the time variable, equation (12) would be the equation motion of a particle in an inverted potential with a damping term. The required boundary conditions are equivalent to the trajectory of a particle starting from the maximum of the inverted potential, rolling down and climbing up to the local maximum which corresponds to the meta-stable point.
As the particle approaches the local maximum its velocity approaches zero. The critical bubble nucleation probability rate per unit volume at finite temperature is proportional to exp(−S/T ), where S is the action of the bounce solution.
A. The bounce The bounce Eq. 12 is non-linear in φ, which makes it difficult to solve analytically. Only in the thin-wall approximation, when the stable and meta-stable are almost degenerate, the bounce can be calculated analytically. Such an approximation will not be valid in the present case as the difference in the thermodynamic potential between the stable and meta-stable states increases with T and dominates over the barrier height. Hence, numerical integration is the only way to find the bounce/bubble profile. The numerical integration is straight forward when there is a symmetry, for example in U (1) theory, where only the radial mode of the field appears in the bounce equation. The phase is taken to be uniform, for minimum action bubble profile.
In the PNJL model there is no such symmetry, the real and imaginary parts of Polyakov loop field and the sigma field are expected to have non-trivial profiles. Since evolving all the three fields simultaneously proved extremely difficult, we kept sigma field constant throughout the trajectory, that is, σ is independent of r. Later we will consider sample profiles for σ to estimate the corrections to the action. In the present case, the thermodynamic potential Ω(Φ 1 , Φ 2 , σ) replaces V (φ) in Eq. 12. The equations to be solved simultaneously are given by, The boundary conditions are Φ i → Φ m i as r → ∞. Φ m i , i = 1, 2, are the values of Φ 1 , Φ 2 in the meta-stable state. For the numerical integration r is discretized as r → r n = nδ, where δ is the lattice spacing. δ must be small compared to the length scale of typical variations of Φ. The integration starts from r = r 0 ∼ 0. Two types of discretizations of Eq. 13 are considered. In the first approach, the values of Φ 1 , Φ 2 at r = r 0 and r = r 0 + δ are used to generate the trajectory. In the second approach the two equations are rewritten as four first order equations. In this case, the values of Φ 1 , Φ 2 as well as their derivatives at r 0 determine the trajectory. It has been checked that both these methods of integration give same results.
In Fig. 5 we show the plot of the inverted potential in the Polyakov loop plane. There is a ridge which connects the stable and meta-stable states. The height of the ridge initially drops from the stable peak but eventually rises to the meta-stable peak. The bounce profile must start near the stable state and approach the meta-stable state. This can happen only for a unique choice of initial conditions, i.e. position (Φ 1 , Φ 2 ) and velocity ( dΦ 1 dr , dΦ 2 dr ). For wrong choices, the trajectory will either fall into the local minimum at the center or cross the ridge and fall off to infinity. Hence, the initial conditions must be tuned which is achieved by standard bisection method. We have checked that the results do not change for smaller δ. i.e., 364.8 MeV to 475 MeV. Fig. 6 shows the temperature dependence of the bubble profiles.
The values of Φ 1 and Φ 2 approach asymptotically to their corresponding meta-stable values.
For temperatures just above T m the barrier between the stable and meta-stable states is small, while Ω s − Ω ms is large. Starting the field close to Φ s i at r 0 ∼ 0 will always lead to overshooting. Hence the initial values of the fields (at the center of the bubble) must be farther away from the stable point. Since the field starts already on a higher slope for small r, damping dominates the profile giving a broad "wall" profile for the bounce. For higher temperatures, the initial point is closer to the stable point. The force term is small, so is acceleration. The field gets to spend more time near the stable state. Therefore the core radius of the bubble increases as we go towards higher temperatures. The bubble "wall" is slimmer because when the particle eventually starts rolling the damping is small.   7 shows the radii of these bubbles as a function of temperature. We define the radius of the bubble as the radial distance from the center to the point where the field drops half way to the meta-stable value. Here we notice that the radii for the two different fields are not the same. The radius of the Φ 1 profile is slightly higher than that of Φ 2 . This is because the curvature of the potential along Φ 1 and Φ 2 , or in other words, their mass scales, are different. Fig. 8 (left), shows a plot of the θ profile of the bubble profile at T = 475 MeV. We also plot the magnitude of the Polyakov loop versus radius in Fig. 8 Here α is a constant given by 2N/g 2 [36], where N is the number of colors and g is the gauge coupling constant. For g/4π = 0.3, α = 1.6. Fig. 9(a) shows the plot of the bubble action in units of temperature vs T /T c . Let us recall here that σ was kept constant at the meta-stable value in the bubble. For an estimate of the change in the action, an approximate σ profile is computed by minimizing Ω with respect to σ for a given Φ 1 , Φ 2 profile. With the σ profile the bubble action increases slightly. The increment is below 20% for 365 MeV, and decreses to 10% at 420 MeV . We also checked with σ profiles similar to both Φ 1 and Φ 2 profiles. In both cases the action increased only by few percent. The bubble nucleation probability per unit time per unit volume, or in other words, the decay rate of the meta-stable state is given by [32] Γ = T 4 S 2πT One can see from Fig. 9(b) that this value is insignificant for temperatures higher than 2T c . It is expected that the decay rate will decrease if we take into account effects of σ profile.
We do a quick calculation of the decay rate for the case of heavy-ion collisions, assuming the thermalization time to be t i = 0.38 fm and the initial temperature of the order of T i = 750MeV. We consider the QGP to be cylindrical with radius 8 fm and length 3 fm.
The number of bubbles nucleated within this volume when the system cools down to a temperature T is estimated as follows. The bubble action S as a function of temperature is obtained by fitting our data points. We used a longitudinally boost invariant 2 + 1D hydrodynamic simulation with Galuber optical initial conditions following [68], to fit the temperature evolution. The number of bubbles nucleated in volume V during the time t when the temperature drops to T is given by We find that N (T m ) = N (t m ) ∼ 0.005, where t m is the time when the temperature drops to T m . The value of N (t) rapidly decreases for higher temperatures. If we assume a larger equilibration time, that is, a smaller initial temperature, the value of N (t) decreases further. The effect of the σ field profile will only enhance the suppression in the nucleation rate. Thus the bubble nucleation probability is very small in heavy ion collisions. This leads to an interesting scenario, known as the spinodal decomposition. When the meta-stable states become unstable, the field will roll down to the stable state resulting in large angular fluctuations. This will have interesting consequences in the dynamics of heavy-ion collisions including flow and also may possibly lead to coherent emission of particles.

VI. CONCLUSIONS
We have studied the Z 3 meta-stable states in PNJL model. The meta-stable states exist at and above the temperature T m ∼ 362 MeV which is well above the deconfinement transition temperature T c = 176 MeV. We discuss the probability of the decay of these meta-stable states by calculating the stable bubble nucleation probability in the meta-stable regions using bounce solution. Our calculations suggest that, the probability of these states decaying by tunnelling into stable states is very small in the case of heavy-ion collisions.
This leads to large angular oscillations of the Polyakov loop field, which can have interesting consequences in the dynamics of flow and coherent emission of particles.

ACKNOWLEDGMENTS
We thank A. P. Balachandran, Shreyansh S. Dave and Ajit Srivastava for important comments and suggestions.