$Z_3$ meta-stable states in PNJL model

TL;DR

This paper investigates Z3 meta-stable states in the PNJL model, analyzing their existence above a certain temperature, their decay via bubble nucleation, and potential implications for heavy-ion collisions.

Contribution

It provides a numerical analysis of the nucleation rates of Z3 meta-stable states in the PNJL model and discusses their decay mechanisms in high-temperature conditions.

Findings

Meta-stable states exist above T_m ~ 194 MeV.

Decay occurs via bubble nucleation, with computed nucleation rates.

Decay in heavy-ion collisions likely via spinodal decomposition.

Abstract

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 ~ 194 MeV and can decay via bubble nucleation. We numerically solve the bounce equation to compute the nucleation rate. We speculate that, in the context of heavy-ion collisions, the likely scenario for the decay of the meta-stable states is via spinodal decomposition.

Equations

$$A_{\mu }(\vec{x},0)=A_{\mu }(\vec{x},\beta ),$$

$$U(\vec{x},0)=zU(\vec{x},\beta ).$$

$$L(\vec{x})=\frac{1}{3}\mathrm{Tr}\left(\mathcal{P}\exp \left[ig{\int }_0^{\beta }d\tau A_0(\vec{x},\tau )\right]\right)$$

$$L(T)=⟨\frac{1}{V}\int L(\vec{x})d^{3}x⟩,$$

$$|L(T){|}^2=\exp \left[-\beta F_{\bar{Q}Q}(r=\infty )\right].$$

$$U(\bar{\Phi },\Phi ,T)=b_4{T}^4\left[-\frac{b_2(T)}{4}(|\Phi {|}^2+|\bar{\Phi }{|}^2)-\frac{b_3}{6}({\Phi }^3+{\bar{\Phi }}^3)+\frac{1}{16}(|\Phi {|}^2+|\bar{\Phi }{|}^2{)}^2\right].$$

$$\mathrm{P}=-U({\bar{\Phi }}_{th},{\Phi }_{th},T).$$

$$b_2(T)=\left(1-1.11T_0/T\right){\left(1+0.265T_0/T\right)}^2{\left(1+0.3T_0/T\right)}^3-0.487$$

$${\mathcal{L}}_{{}_{PNJL}}=\sum _f({\bar{\Psi }}_f(i{\gamma }_{\nu }{D}^{\nu }-m_f){\Psi }_f+G_s[({\bar{\Psi }}_f{\tau }_a{\Psi }_f{)}^2+({\bar{\Psi }}_{f}i{\gamma }_5{\tau }_a{\Psi }_f{)}^2])+U(\bar{\Phi },\Phi ,T)$$

$$\begin{array}{r}\Omega =-\sum _{f=u,d}{\int }_0^{\infty }\frac{d^{3}p}{(2\pi {)}^3}(2T\ln [1+3\Phi e^{-\beta (E_f-{\mu }_f)}+3\bar{\Phi }{e}^{-2\beta (E_f-{\mu }_f)}+e^{-3\beta (E_f-{\mu }_f)}]\\ +2T\ln [1+3\bar{\Phi }{e}^{-\beta (E_f+{\mu }_f)}+3\Phi e^{-2\beta (E_f+{\mu }_f)}+e^{-3\beta (E_f+{\mu }_f)}])\\ -6\sum _{f=u,d}\int \frac{d^{3}p}{(2\pi {)}^3}{E}_f\theta (\Lambda -|\vec{p}|)+\sum _{f=u,d}{G}_s{\sigma }_f^2+U(\bar{\Phi },\Phi ,T).\end{array}$$

$$\frac{\partial \Omega }{\partial {\Phi }_1}=0,\frac{\partial \Omega }{\partial {\Phi }_2}=0,\frac{\partial \Omega }{\partial \sigma }=0.$$

$$\frac{d^2\varphi }{d{r}^2}+\frac{2}{r}\frac{d\varphi }{dr}=\frac{\partial V}{\partial \varphi }$$

$$\frac{d^2{\Phi }_1}{d{r}^2}+\frac{2}{r}\frac{d{\Phi }_1}{dr}$$

$$\frac{d^2{\Phi }_2}{d{r}^2}+\frac{2}{r}\frac{d{\Phi }_2}{dr}$$

$$S=\int{4\pi r^{2}dr\Big{[}\frac{1}{2}\alpha T^{2}\left\{\left(\frac{d\Phi_{1}}{dr}\right)^{2}+\left(\frac{d\Phi_{2}}{dr}\right)^{2}\right\}+\frac{1}{2}G_{s}^{2}\left(\frac{d\sigma}{dr}\right)^{2}+\Omega(\Phi_{1},\Phi_{2},\sigma)\Big{]}}$$

$$\Gamma=T^{4}\Big{(}\frac{S}{2\pi T}\Big{)}^{3/2}\exp(-S/T)$$

$$F_b=-\frac{4\pi }{3}\rho R_d^3+4\pi \delta R_d^2,R_d=\frac{2\delta }{\rho }.$$

$$N(\tau )=V{\int }_{{\tau }_i}^{\tau }\Gamma (t)dt.$$

Full text

I Introduction

In pure SU(N) gauge theories, energy density increases sharply across the 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 to deconfined state of gluons, known as confinement-deconfinement (CD) transition, has been extensively studied in the literature Kuti:1980gh ; Celik:1983wz ; Susskind:1979up ; Engels:1980ty ; Karsch:2001cy . The CD phase transition is found to be second order for $N=2$  Damgaard:1987wh ; Engels:1988ph ; Christensen:1990qs ; Damgaard:1994np ; Engels:1994xj ; Engels:1998nv and, first order for $N\geq 3$ Boyd:1995zg ; Boyd:1996bx . The Polyakov loop, which transforms as a $Z_{N}$ spin, plays the role of an 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$ Balachandran:2001qn ; Layek:2005fn ; Gupta:2010pp .

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 Green:1983sd ; Oevers:1997yf ; Karsch:2000zv ; Belyaev:1991np ; Deka:2010bc . It affects the nature of the CD transition Green:1983sd ; Karsch:2000zv 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) Gupta:2011ag .

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  Gross:1980br ; Weiss:1980rj ; Weiss:1981ev ; Guo:2018scp . These perturbative calculations are reliable for high temperatures ($T>>T_{c}$), 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 Gross:1980br ; Weiss:1981ev ; Belyaev:1991np ; Dixit:1991et . These states have been studied extensively in the context of cosmology. In the early Universe, they are found to be long lived, and can leave observable imprints while decaying nucleation of bubbles of true vacuum as in a first order transition Dixit:1991et ; Ignatius:1991nk . However, 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 Belyaev:1991np . This problem doesn’t arise in QCD near the critical temperature as the number of flavors is effectively $\leq 3$.

The study of $Z_{3}$ meta-stable states for small temperatures, in particular near $T_{c}$ is important as they may affect 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 gauge coupling constant 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 other two become meta-stable states Deka:2010bc . 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 2$T_{c}$  Deka:2010bc . 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 Mishra:2016ipq , 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 $fm$ 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 Ignatius:1991nk . 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 at zero baryon chemical potential. In this model, they exist above $T_{m}\sim 194$ MeV. If such a state exists, 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 Coleman:1977py ; Callan:1977pt ; Linde:1981zj ; Gleiser:1991rf 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. This will lead to large oscillations of the Polyakov loop. We mention here that in heavy-ion collisions the baryon chemical potential is small but non-zero. At finite chemical potential the thermodynamic potential has an imaginary part. There are several papers that discuss how to include the effect of non-zero $\mu$ Dumitru:2005ng ; Rossner:2007ik . With $\mu$, the contribution of the fermions to the free energy will increase. Since the fermions break the $Z_{N}$ symmetry, we expect that finite $\mu$ will lead to more explicit breaking. Following Roessner et.al., for small $\mu$, we calculated the thermodynamic potential to the leading order, i.e, keeping only the real terms and found that $T_{m}$ increases slightly with $\mu$ Mishra:2016ipq . With increase in $T_{m}$ there is lesser time available for the nucleation of bubbles which enhances the likelihood of meta-stable states becoming unstable.

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

[TABLE]

where $\beta={1\over T}$. This boundary condition allows for the gauge transformations, $U(\vec{x},\tau)$, to be periodic up to a factor $z\in Z_{N}$, such as

[TABLE]

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

[TABLE]

Here ‘${\cal P}$’ denotes path order, ‘g’ is the gauge coupling and $A_{0}=A_{0}^{a}{\tau^{a}\over 2}$ is the temporal gauge field. Here $\tau_{a}$ are the Pauli matrices with ‘$a$’ denoting the color indices. Under a $Z_{3}$ gauge transformation, Eq. (2), the Polyakov loop, transforms as $L(\vec{x})\rightarrow zL(\vec{x})$. The thermal as well as volume average of the Polyakov loop $L(\vec{x})$,

[TABLE]

is related to the free energy $F_{\bar{Q}Q}(r)$ of a static (infinitely heavy) quark-anti-quark pair at infinite separation Svetitsky:1985ye .

[TABLE]

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 Pisarski:2000eq . We consider the following Landau-Ginsburg effective potential for a complex field $\Phi$ Pisarski:2000eq ; Dumitru:2000in ; Dumitru:2001vc .

[TABLE]

Different forms of the effective potential in terms of the field $\Phi$ have been proposed  Dumitru:2000in ; Ratti:2005jh ; Ratti:2006wg ; Roessner:2006xn ; Bhattacharyya:2010jd ; Ghosh:2007wy . Across the critical temperature $T_{0}$, the Polyakov loop expectation value jumps discontinuously. 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 Pisarski:2000eq . In the mean-field approximation the minimum (minima), $\Phi_{th}$, of the effective potential gives the thermal average of the Polyakov loop, $L(T)$. $\Phi_{th}=L(T)$. In this approximation the pressure $P$ is given by,

[TABLE]

The above effective potential with the following form of $b_{2}(T)$

[TABLE]

and the coefficients $b_{3}=2.0$ and $b_{4}=0.6016$ reproduces the pressure of the pure gauge theory computed from non-perturbative lattice method(s). Following Dumitru:2000in , for QCD, $b_{4}$ is rescaled as $b_{4}=0.6061\times 47.5/16$. Here the coefficients are chosen such that the expectation value of the order parameter $\Phi$ approaches unity for $T\rightarrow\infty$. Hence, the fields and the coefficients in the above potential are rescaled as $\Phi\rightarrow\Phi/x$, $b_{2}(T)\rightarrow b_{2}(T)/x^{2}$, $b_{3}\rightarrow b_{3}/x$, $b_{4}\rightarrow b_{4}/x^{4}$, where $x=b_{3}/2+{1\over 2}\sqrt{b_{3}^{2}+4b_{2}(T=\infty)}$ for temperature $T\rightarrow\infty$. Note that we have used a polynomial form of the Polyakov loop potential which is different from the more commonly used form as in Ratti:2005jh . The $T_{m}$ for the latter case is higher than the initial temperature of QGP in heavy ion collisions which is not in agreement with lattice study Deka:2010bc . Also it is not clear that these models will be valid at such high temperatures.

For $T>T_{0}$, there are three degenerate minima in the effective potential, which can be seen in Fig. 1, the contour plot of the effective potential in the complex $\Phi=(\Phi_{1},\Phi_{2})$ plane at $T=1.3T_{0}$. We also plot the variation of the potential along a circle going through the three minima, i.e $\Phi=|\Phi_{th}|e^{i\theta}$ in Fig. 1. The three degenerate vacua are situated at $\theta=0$, ${{2\pi}\over 3}$ and ${{4\pi}\over 3}$, related by $Z_{3}$ rotation.

When dynamical quarks are included, the $Z_{3}$ symmetry is broken. While the pure gauge part of the action is $Z_{3}$ symmetric, the quark part of the action is not invariant under the $Z_{3}$ gauge transformations. This is because the gauge transformed quark fields are no longer anti periodic along the temporal direction. The non-trivial $Z_{3}$ gauge transformations can act only on the gauge fields. The situation is similar to the presence of an external (explicit breaking) field in spin systems. For example, in the presence of an external field the Ising model Hamiltonian has both $Z_{2}$ symmetric and broken terms. As the magnetization still describes the Ising transition, the field $\Phi$ too describes the CD transition Engels:1993tp . In the following section, we discuss the PNJL model which provides a prescription to include the effect of quarks on the $Z_{3}$ symmetry and the Polyakov loop effective potential.

III Meta-stable states in PNJL model

The PNJL model is an extension of the Nambu-Jona-Lasinio (NJL) model. The NJL model is a phenomenological model formulated on the basis of the chiral symmetry of QCD and describes the dynamics of low energy excitations as well as the chiral transition Nambu:1961tp ; Nambu:1961fr ; Klevansky:1992qe ; Vogl:1991qt ; Scavenius:2000qd ; Hatsuda:1994pi ; Buballa:2003qv ; Nebauer:2001rb . Since there are no gauge fields in this model, it can not describe the CD transition. The PNJL model attempts to include the gauge fields by adding the effective potential $U(\bar{\Phi},\Phi,T)$ to the NJL Lagrangian Ratti:2005jh ; Mukherjee:2006hq ; Roessner:2006xn ; Costa:2008gr ; Ghosh:2008et ; Deb:2011en . Further, in the fermion part of NJL model, covariant derivative substitutes the standard one. The PNJL Lagrangian is given by,

[TABLE]

Here $D_{\nu}$ is the covariant derivative, $D_{\nu}=\partial_{\nu}-i(gA_{\nu}+{\delta_{0}^{\nu}}\mu_{f})$, $A_{\nu}=A_{\nu}^{a}{\tau^{a}\over 2}$. Here subscript $f$ refers to the $u$, $d$ quark flavors. This term takes into account the interaction between the gauge and quark fields. $m_{f}$ and $\mu_{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 Ghosh:2007wy ; Ratti:2005jh ,

[TABLE]

Here $\sigma_{f}=\langle{\bar{\Psi}}\Psi\rangle_{f}$ is the quark condensate. $\mu_{u}$ and ${\mu_{d}}$ are $u$ and $d$ quark chemical potentials respectively. For the present calculations chemical potentials are set to $\mu_{u}=\mu_{d}=0$. The masses of $u$ and $d$ quarks are taken to be degenerate, i.e $m_{u}=m_{d}=m_{0}$. $E_{u,d}=\sqrt{p^{2}+\Sigma_{u,d}^{2}}$ are the single particle energies with $\Sigma_{u,d}=m_{0}-G_{s}\sigma_{u,d}-G_{s}\sigma_{d,u}$ where $G_{s}=10.08$ GeV*-2* and $m_{0}=5$ MeV. For $T_{0}=190$ MeV in the effective potential ${\it U}$, the thermodynamic potential given by Eq. 10 results in qualitatively similar thermodynamic behaviour as that in Ghosh:2007wy ; Ratti:2005jh .

For the temperature dependence of the quark condensates and expectation value of the Polyakov loop, we minimize the thermodynamic potential $\Omega(\sigma,m,T_{0},T)$ by numerically solving the following set of equations,

[TABLE]

$\Phi_{1}$ and $\Phi_{2}$ are real and imaginary parts of $\Phi$. The numerical program requires initial trial values of $\Phi$ and $\sigma$. It evolves the trial values such 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 $\Phi$ is real valued, i.e $\Phi_{2}=0$ ($\theta=0$). Hence, for the ground state, we solve the above equations with initial values of $|\Phi|>0$, and $\theta=0$. We took the zero temperature value of $\sigma$ as its initial value. For the meta-stable states, the $Z_{3}$ rotated values of the ground state $\Phi$ as initial value works well. Since the $Z_{3}$ symmetry is explicitly broken, $Z_{3}$ rotated $\Phi$ does not solve the equations. However, the meta-stable states are found to be close to $Z_{3}$ rotations of the stable state. The value of $\sigma$ in the meta-stable state is found to differ from the value in the stable state.

With above numerical procedure, in the temperature range of $T_{c}(170)$ MeV - $T_{m}(194)$ MeV, only one solution is found with $\Phi_{1}>0$ and $\Phi_{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>194$ MeV, two local minima appear around $\theta={2\pi\over 3}$ and $\theta=-{2\pi\over 3}$. The thermodynamic potentials for these two states are found to be same, but higher than that for the ground state for which $\theta=0$. The values of $\Phi$ for the meta-stable and stable states are no more related by $Z_{3}$ rotation.

In Fig. 2, the thermodynamic potential versus $\theta$ has been plotted by fixing $|\Phi|$ in the stable state ($|\Phi|=0.79763$ for $T=$199.5 MeV and $|\Phi|=0.92962$ for 247 MeV) and minimising thermodynamic potential with respect to $\sigma$. In Fig. 2 the temperature is close to $T_{m}$ when the effective potential develops two saddle points. Comparing Fig. 2 and Fig. 2 one can clearly see signs that the barrier between meta-stable and stable states increases with temperature. In Fig. 3, the contour plots of the thermodynamic potential in the $\Phi_{1}$-$\Phi_{2}$ plane are shown. At each point $\Phi_{1}$, $\Phi_{2}$ in the contour plot, we have fixed $\sigma$ at the value which minimises the thermodynamic potential. For temperatures above $T_{m}$ there are two meta-stable states and one stable state. We denote the meta-stable state for which $\Phi_{2}<0$ by M1 and the other by M2. The stable state is denoted by SS. Subscripts $s$ and $ms$ on variables denote their values in the stable state and the meta-stable states respectively.

Fig. 4 shows the difference in the thermodynamic potential of the meta-stable to the stable state ($\Omega_{ms}-\Omega_{s}$) vs $T$. This difference increases with temperature, which suggests enhancement in the explicit breaking for larger temperatures. In Fig. 5 we show the Polyakov loop in the M1, M2 and SS states for small values of $\mu$. We find that with increase in $\mu$ the absolute value of the Polyakov loop decreases in the meta-stable states while it increases in the stable state. $\mu$ contribution effectively tilts the thermodynamic potential towards the stable state. We find that the phase of the M2 state shifts toward higher value. The phase of M1 decreases as it is the complex conjugate of M2. We also find that the barrier height between (M1, M2) and SS decreases while $T_{m}(\mu)$ increases slightly.

IV Bounce solution for the decay of meta-stable states

In PNJL model, with 2 quark flavors, meta-stable states exist above $T_{m}\sim 194$ MeV. Even though the value of $T_{m}$ in this model is too small compared to the lattice result  Deka:2010bc , we believe that it will give qualitative results for the effect of meta-stable states in heavy-ion collisions. We mention here that, with other Polyakov loop potentials, higher $T_{m}$ can be achieved by tuning the integration cut-off in Eq. 10 though the results do not qualitatively differ from the present case.

As mentioned before if a system is in the meta-stable state, it will eventually decay to the stable state. A meta-stable state can either become unstable if the temperature drops below $T_{m}$ or decay via nucleation of bubbles like in a first order phase transition. At finite temperature, there will be fluctuations in the form of bubbles, with stable states in their core. The free energy of a bubble consists of two components, the volume component and the surface component. The volume component comes from the free energy difference between the stable and the meta-stable states. The surface component comes from the fact that, the fields ($\Phi,\sigma$) have to interpolate between stable values at the center to meta-stable values outside. For a critical bubble, these two components balance and a small fluctuation can make it grow or collapse. Thus the critical bubble and its nucleation rate play an important role in a first order phase transition. For decay of the state M1(M2) the fields $\Phi$ and $\sigma$ will have values corresponding to the SS inside the bubble. Both these fields vary smoothly across the bubble wall and approach the values corresponding to M1(M2). Given that the free energy does not depend on the sign of $\Phi_{2}$, the bubbles interpolating M1 and SS will have same action as the other interpolating M2 and SS. The critical bubble is obtained from the bounce solution which is a saddle point of the Euclidean action. We must mention here that the bubble nucleation picture here is not related to any phase transition but to the fact that the theory allows existence of meta-stable states above a certain temperature and they can tunnel into the stable state.

The decay rate of the false vacuum (meta-stable state) can be calculated in the semi classical approximation where the dominant contribution comes from the configurations with the least action Coleman:1977py ; Callan:1977pt , i.e bounce solutions. It is shown that such configurations, at zero temperatures have $O(4)$ symmetry, reducing the problem to one degree of freedom along the radial direction given by $r^{2}=|\bf x|^{2}+\tau^{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(\phi)$, 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” Coleman:1977py . 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 Linde:1981zj . For a single scalar field theory with a meta-stable state, the bounce can then be calculated by solving the equation of motion,

[TABLE]

with the boundary conditions $\phi\rightarrow\phi^{m}$ as $r\rightarrow\infty$, where $\phi^{m}$ is the value of the field in the meta-stable state. For $r\sim 0$ the field is expected to be close to the stable state. If $r$ were the time variable, Eq. 12 would be the equation of 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 (which corresponds to the stable state), 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.

IV.1 The bounce

The bounce Eq. 12 is non-linear in $\phi$, 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 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, $\sigma$ is independent of $r$. Later we will consider sample profiles for $\sigma$ to estimate the corrections to the action. We also calculate the lower bound of the action. In the present case, the thermodynamic potential $\Omega(\Phi_{1},\Phi_{2},\sigma)$ replaces $V(\phi)$ in Eq. 12. The equations to be solved simultaneously are given by,

[TABLE]

The boundary conditions are $\Phi_{i}\rightarrow\Phi_{i}^{m}$ as $r\rightarrow\infty$. $\Phi_{i}^{m}$, $i=1,2$, are the values of $\Phi_{1},\Phi_{2}$ in the meta-stable state. For the numerical integration $r$ is discretized as $r\rightarrow r_{n}=n\delta$, where $\delta$ is the lattice spacing. $\delta$ must be small compared to the length scale of typical variations of $\Phi$. The integration starts from $r=r_{0}\sim 0$. Two types of discretizations of Eq. 13 are considered. In the first approach, the values of $\Phi_{1},\Phi_{2}$ at $r=r_{0}$ and $r=r_{0}+\delta$ 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 $\Phi_{1},\Phi_{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. A few other approaches of calculating bounce solution for multiple field cases are discussed in Konstandin:2006nd ; Wainwright:2011kj ; Piscopo:2019txs ; Masoumi:2016wot ; Athron:2019nbd ; Sato:2019axv ; Profumo:2010kp ; Chigusa:2019wxb .

In Fig. 6 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 $(\Phi_{1},\Phi_{2})$ and velocity $(\frac{d\Phi_{1}}{dr},\frac{d\Phi_{2}}{dr})$. For wrong choices, the trajectory will fall off to infinity either through the center along $\theta=\pi$ or by crossing the ridge to $\theta\sim-\pi/3$. Hence, the initial conditions must be tuned which is achieved by standard bisection method. The basic idea is that with the given initial conditions, position ($\Phi_{1}(r_{0}),\Phi_{2}(r_{0}))$ and velocity (${d\Phi_{1}\over dr}|_{r_{0}},{d\Phi_{2}\over dr}|_{r_{0}}$), if a trajectory undershoots (overshoots) the meta-stable point (peak), we start with an initial choice (position) closer to (farther from) the “global” minimum (global peak). Further the bisection method is used to find the direction of initial velocity such that the undershoot trajectory makes a $180^{o}$ turn or the overshoot trajectory passes through the meta-stable state peak. We have checked that the results do not change for smaller $\delta$.

IV.2 The bubble

The bubble profiles are computed for temperatures in the range $T=1.05T_{0}~{}-~{}1.2T_{0}$, i.e., 199.5 MeV to 228 MeV. Fig. 7 shows the temperature dependence of the bubble profiles. These bubble profiles represent the decay of M1 to SS. The values of $\Phi_{1}$ and $\Phi_{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 compared to $\Omega_{s}-\Omega_{ms}$. Starting with initial values of the field close to $\Phi_{i}^{s}$ at $r_{0}\sim 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 thinner because the particle/field spends larger time near the stable maximum and when it eventually starts rolling, the damping is small.

Fig. 8 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 $\Phi_{2}$ profile is slightly higher than that of $\Phi_{1}$. This is because the curvature of the potential along $\Phi_{1}$ and $\Phi_{2}$, or in other words, their mass scales, are different.

Fig. 9 (left), shows a plot of the $\theta$ profile of the bubble at $T=218.5$ MeV. We also plot the magnitude of the Polyakov loop versus radius in Fig. 9 (right).

V Evolution of meta-stable state in heavy-ion collisions

The decay rate of these meta-stable states depends on the bubble action, which is given by,

[TABLE]

Here $\alpha$ is a constant given by $2N/g^{2}$ Dumitru:2000in , where $N$ is the number of colors and $g$ is the gauge coupling constant. For $g/4\pi=0.3$, $\alpha=1.6$. Fig. 10 shows the plot of the bubble action in units of temperature vs $T$. Let us recall here that $\sigma$ was kept constant at the meta-stable value in the bubble. For an estimate of the change in the action, an approximate $\sigma$ profile is computed by minimizing $\Omega$ with respect to $\sigma$ for a given $\Phi_{1},\Phi_{2}$ profile. With the $\sigma$ profile the bubble action decreases slightly. The decrement is below $10\%$ for all the temperatures calculated. We also checked with $\sigma$ profiles scaled like both $\Phi_{1}$ and $\Phi_{2}$ profiles, interpolating between $\sigma_{s}$-$\sigma_{ms}$. In the case where $\sigma$ profile was scaled like $\Phi_{1}$ the action was minimum (less than $20\%$ decrement).

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 Linde:1981zj

[TABLE]

One can see from Fig. 10 that this value is as small as $10^{-6}$/fm for the smallest temperature above $T_{m}$ (199.5 MeV) and grows insignificant for higher temperatures. Though it is difficult to solve for the bubble in the full case by including equation corresponding to $\sigma$ in Eq. 13, we can compute the lower bound of the full action. The corresponding nucleation rate will then be the upper bound. Note that $\Gamma$ has a peak in $S/T$. However, for the range of $S/T$ in our calculation, $\Gamma$ is a monotonically decreasing function of $S/T$. To compute the upper bound on the nucleation rate we fix the $\sigma$ field at the SS value ($\sigma_{s}$). To see how $\sigma(r)=\sigma_{s}$ leads to this bound, we write the free energy ($F_{b}$) of a critical bubble of radius $R_{d}$ as

[TABLE]

where $\rho$ is the free energy difference between the stable and meta-stable state, $\delta$ is the free energy cost (surface tension) as the fields vary smoothly between the two states. The position of the meta-stable state when the $\sigma$ is fixed at $\sigma_{s}$ deviates from the full case such that $|\Phi_{s}-\Phi_{ms}|$ decreases. This effectively leads to decrease in $\delta$. Also with $\sigma=\sigma_{s}$ the meta-stable states always have higher $\Omega_{ms}$ compared to the full case. This leads to increase in $\rho$. Hence the free energy for $\sigma=\sigma_{s}$ is lower compared to the full case. Note that the barrier height also plays a role in determining $\delta$ which is the reason the bubble action grows with temperature. However for a given temperature the barrier would slightly decrease as $\sigma$ changes from $\sigma_{ms}$ to $\sigma_{s}$.

We do a quick calculation of the decay rate for the case of heavy-ion collisions, assuming the thermalization time to be $\tau_{i}=0.38$ fm and the initial temperature of the order of $T_{i}\sim 550$ MeV. We consider the QGP to be cylindrical (the midrapidity region) 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 Glauber optical initial conditions following Kolb:2000sd , to fit the temperature evolution. The number of bubbles nucleated in volume $V$ during the time $\tau$ when the temperature drops to $T$ is given by

[TABLE]

We find that $N(T_{m})=N(\tau_{m})$, where $\tau_{m}$ is the time at which temperature is $T_{m}$, is vanishingly small. Here we have used the profile with sigma scaled as $\Phi_{1}$. If any of the other profiles are used, this number only decreases. The value of $N(\tau)$ rapidly decreases for higher temperatures. If we assume a larger equilibration time, that is, a smaller initial temperature, the value of $N(\tau)$ decreases further. We have also computed the free energy of the bubble (action) and the nucleation rate by fixing $\sigma=\sigma_{s}$. We find that the bubble action is smaller by an average factor of $\sim 0.45$ compared to the case when $\sigma$ is fixed at $\sigma_{ms}$. In this calculation the nucleation rate increases by a factor of hundred, though still remaining negligibly small. We have considered the effects of $\mu$ upto $\sim 100$ MeV. For $\mu=100$ MeV $T_{m}$ increases by $\sim 2$ MeV. Since the barrier height decreases with $\mu$, the bubble action $S(T,\mu)<S(T,\mu=0)$. However the system spends lesser time above $T_{m}(\mu)$, hence the results will not change qualitatively. Hence, we do not expect any bubble nucleation in heavy ion collisions. This leads to an interesting scenario, known as the spinodal decomposition. When the meta-stable states become unstable below $T_{m}$, the field will roll down to the minimum resulting in large angular fluctuations. Fig. 11 shows the evolution of the Polyakov loop and chiral condensate at the centre of the quark gluon plasma after the temperature falls below $T_{m}$ at zero chemical potential. Since the meta-stable states become unstable, the fields roll down and oscillate around the minimum. These oscillations will have interesting consequences in the dynamics of heavy-ion collisions including flow, jet energy loss Lin:2013efa and also may possibly lead to coherent emission of particles. As discussed above, at a small finite chemical potential, $T_{m}(\mu)$ increases, and spinodal decomposition is expected to occur earlier.

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}\sim 194$ MeV. For small values of $\mu$ we found that $T_{m}$ increases slightly, i.e by $\sim 2$ MeV at $\mu=100$ MeV. The barrier height between the meta-stable and stable states decreases. We have discussed 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. The bubble action measured in the units of temperature increases roughly linear in temperature. For small $\mu$, relevant for heavy-ion collisions, the bubble action decreases slightly though the system spends lesser time above $T_{m}$. Our results suggest that, the probability of these states decaying by tunnelling into stable states is very small in the case of heavy-ion collisions. Ultimately the meta-stable state will become unstable and the fields will start rolling towards the minimum. This will lead to large oscillations of the Polyakov loop field, which may have interesting consequences to the dynamics of flow, jet energy loss, and also may lead to coherent emission of particles.

Acknowledgements.

We thank A. P. Balachandran, Shreyansh S. Dave, Ajit Srivastava, Rajarshi Ray and Ranjita Mohapatra for important comments and suggestions.

REFERENCES