Bubble dynamics in the Polyakov quark-meson model

1. Introduction

At sufficiently high temperatures and densities, a theory of strong interaction, quantum chromodynamics (QCD), supports an idea that the hadronic matter at low energy undergoes a phase transition into a deconfined and chiral-symmetry restored quark matter denoted as quark-gluon plasma (QGP). In experiment, how to explore and illuminate the basic properties of QGP is a current topic related to the physics of the heavy-ion collisions (HICs) in laboratories [1,2,3] as well as to the astrophysics of compact stars [4,5,6]. On the theoretical side, although QCD is applicable for determining the properties of strongly interaction matter at a very high energy, the understanding of the nature of the hadron quark phase transition remains a challenge in actual calculations, in particular when finite chemical potentials are concerned. Since an ab initio approach, Monte Carlo techniques on the lattice, is severely hindered by the notorious fermion sign problem [7]. Therefore, adopting effective models provides an alternative therapy to capture the nontrivial features of QCD, such as two salient features intimately connected with the nonperturbative structure of the QCD, chiral symmetry and the confinement. Among these effective models are the Nambu-Jona-Lasinio (NJL) model [8,9], the linear sigma model (LSM) [10] and their extensions formulated with effects of the deconfinement, the Polyakov Nambu-Jona-Lasinio model (PNJL) [11,12,13] and the Polyakov Quark Meson Model (PQM) [14,15,16] for quark matter. While others are the relativistic mean field (RMF) model and its extrapolations incorporated with chiral symmetry, the chiral nucleon meson model [17,18,19,20], the parity doublet model [21,22,23,24,25,26] and the relativistic mean field including chiral potential and confinement effects [27,28,29] for nuclear matter.

Since the first direct detection of gravitational waves by LIGO Collaboration [30], a new era of gravitational wave astronomy offer a new way to study cosmological first-order phase transition [31,32,33]. The subject has gathered fresh interest due to the fact that much of the present interest in primordial first-order phase transition can be explained by detecting a stochastic gravitational wave background in the current and near-future gravitational wave detectors [34,35,36]. Once the signal of its spectrum is to be observed, it would dramatically increase our capability to probe the earliest Universe cosmology. Meanwhile, with the observation of gravitational waves from the binary neutron star merger event GW170817 [37,38], astrophysics has now entered the multimessenger era. Combined with other observations, such as the radius, the mass and the tidal deformability of neutron star, we have so far the most stringent constraint on the equation of states (EOS) of dense matter, which play a key role in the dynamical evolution and structure of neutron star [39,40,41,42,43]. So that any efforts to study strong interaction as well as the possible first-order hadron quark phase transition are continuously deserved to be made in the field of astrophysics and cosmology.

For QCD phase transition, most models in description of strong interactions predict a first-order phase transition from a chiral-symmetry broken hadron phase to a chiral symmetric quark phase at high density and low temperature. Whereas, there is a smooth crossover conversion at low chemical potential but high temperature. At the endpoint of the first-order coexistence line there should exist a critical endpoint (CEP) where the first-order phase transition is to be of the second-order. In laboratories, to explore and locate the CEP through ultra-relativistic heavy-ion collision are the main goals of the Beam Energy Scan (BES) programs at Relativistic Heavy-Ion Collider (RHIC) and other near future experimental facilities, the empirical confirmation of the existence of the CEP will be a milestone of studying QCD phase diagram and strong interactions [44,45]. On the theoretical side, the existence of the CEP guarantees the occurrence of the first-order phase transition, and vice versa. Since much attentions have been drawn to the crossover phase conversion of QCD over past decades, in the present work, most of our studies are focused on the topic of a first-order hadron quark phase conversion at high density and low temperature.

In a first-order phase transition, a system changes from a metastable high energy phase to a stable lower energy phase through the nucleation of bubbles, which is separated from a surrounding metastable phase by an interface. The appearance of a bubble is to be considered as a natural consequence of the thermal and quantum fluctuations of the thermodynamical system sufficiently close to a first-order phase transition. According to these fluctuations, bubbles are expected to create. Among them, some will grow up while others will shrink and disappear, depending on its energy budget with regard to the homogeneous preexisting metastable phase. Only the bubbles larger than a critical size will ultimately survive and then have a decisive role in the process of bubble nucleation. In the end, these critical bubbles will expand, collide and eventually coalesce to complete a fist-order phase transition. The modern theory of bubble nucleation was pioneered by the works of Langer in the context of classical statistical mechanics in the late 1960s [46,47], and applied to field theory by Kobzarev et al [48]. These initial results were put onto firm theoretical ground in seminal papers by Callan and Coleman for relativistic quantum field theory at zero temperature [49,50,51]. Shortly afterwards, the studies were generalized to nonzero temperature by Affect [52] and Linde [53,54]. The remarkable goal of a nucleation theory is to determine the bubble nucleation rate per unit time per unit volume based on the semi-classical saddle point approximation. Within this method, the key point is to find the bounce solution, which is an unstable instanton field configuration that extremizes the action and interpolate between the false vacuum and true vacuum on the other side of the barrier. Once the bounce is found, one can compute the crucial issues intimately related to the dynamics of a first-order phase transition, such as the critical radius, the surface tension and the saddle point action.

Traditionally, in description of a first-order hadron quark phase transition, the equation of states can be technically obtained by matching the pressures in terms of quark and hadron matters via the Maxwell construction or the Gibbs condition. The reason is that for each aspect of matters what we needed is the global minimum of the effective potential, which gives the pressure of the system. However, in order to calculate the bounce solution, the knowledge of the entire potential is necessarily requested. In other words, a unified shape of the effective potential that contains hadron and quark phases is necessary for the existence of exact bounce solution. For a realistic nuclear matter, the bubble nucleation and the surface tension have been calculated and obtained in Ref. [19], where the chiral nucleon-meson (CNM) model with the nucleon degree of freedom is employed [17,18]. The CNM model incorporated with chiral symmetry has been fit to reproduce the nuclear properties at saturation in symmetric matter. Recently, it has been used to study the possible mixed phase at the chiral phase transition in neutron stars [20] and a first-order nuclear liquid-gas phase transition in nuclear matter [19]. Due to the lack of quark degree of freedom, the model does not include the effect of confinement and the chirally restored phase can only be treated in a very rough approximation of dense quark matter. On the contrary, with a similar chiral potential, a quark meson model can be successfully used to study the chiral phase transition at finite temperature and density for two flavors [55,56] and three flavors [57,58]. In particular, within the QM model, the order parameters for chiral phase transition are well defined by quark condensates. Besides phase diagram and thermodynamics, the dynamics of a first-order hadron quark phase transition are also investigated in the QM model in Refs. [59,60,61,62,63]. Moreover, the model can be extended by coupling light chiral quarks to the Polyakov loop fields in order to study the chiral and deconfinement phase transitions simultaneously and confront with the Lattice data directly [14,15,16]. In the thin-wall approximation, the surface tension and phase diagram have been obtained in the three flavors PQM model in Refs. [64,65].

Aside from a few works in Refs. [59,63], most of studies in the QM model or the PQM model are applied in the thin-wall approximation, and the sea-quark contribution due to the presence of the fermion vacuum fluctuation in the effective potential are normally ignored. Because the applicability of the thin-wall approximation is limited and deteriorates when two minima of effective potential become separated. In particular, there exists a weak first-order phase transition in phase diagram at high density and low temperature, where the main mechanism of a phase transition is changed to a spinodal decomposition [66] as temperatures are close to the spinodal lines. Furthermore, the sea quark contribution in the pressure plays a significant role in the study of the hadron quark phase transition for the sake that it will soften the first-order phase transition and reduce the surface tension of hadron quark interface dramatically. Therefore, by including the fermion vacuum fluctuation and the decofinement, we study the bubble nucleation features related to the dynamics of first-order chiral phase transition in the PQM model from nuclear matter to quark matter or vice versa.

The paper is organized as follows. In the following section we briefly describe the Polyakov quark-meson model and present the effective potential in the mean field approximation. Sect. III describes phase diagram and effective potentials at two different directions of order parameter fields. In Sec. IV, we give a brief review of the theory of homogeneous thermal nucleation and discuss methods to search the local minima of the effective potential. Our results and discussions are presented in Sec. V, whereas we gives our summary in the last section.

2. Formulation of the Model

The Lagrangian density of the quark meson model (QM) which has an explicit realization of chiral $S{U}_L\left(2\right)\otimes S{U}_R\left(2\right)$ symmetry is given by [10,55,56]

$$\mathcal{L}=\overline{\psi }\left[i{\gamma }^{\mu }{\partial }_{\mu }-g(\sigma +i{\gamma }_5\overrightarrow{\tau }·\overrightarrow{\pi })\right]\psi +\frac{1}{2}\left({\partial }_{\mu }\sigma {\partial }^{\mu }\sigma +{\partial }_{\mu }\overrightarrow{\pi }·{\partial }^{\mu }\overrightarrow{\pi }\right)-U(\sigma ,\overrightarrow{\pi }),$$

(1)

where the spinor field $\psi$ is a flavor doublet

$$\psi =\left(\genfrac{}{}{0pt}{}{u}{d}\right).$$

(2)

The pseudoscalar field $\overrightarrow{\pi }=({\pi }_1,{\pi }_2,{\pi }_3)$ is an isotriplet pion fields and $\sigma$ is the spin-0 isosinglet field. g is a flavor-blind Yukawa interaction of strength that couples the quarks with four mesons.

The pure mesonic potential, which exhibits both spontaneously and explicitly broken chiral symmetry, is given by the expression

$$U(\sigma ,\overrightarrow{\pi })=\frac{\lambda }{4}{\left({\sigma }^2+{\overrightarrow{\pi }}^2-{\vartheta }^2\right)}^2-H\sigma .$$

(3)

Here $\lambda$ is a positive coupling constant for the mesonic fields, the constant quantity $\vartheta$ is usually called the vacuum expectation value of scalar field. The last term in above equation is due to the finite current-quark masses and is fixed by the partially conserved axial-vector current relation which gives $H=f_{\pi }{m}_{\pi }^2$.

Besides the chiral symmetry, it is well-known that the basic quantum chromodynamics (QCD) Lagrangian has another salient feature, the confinement. The confinement is connected with the center symmetry of the color gauge group and the center symmetry is believed to get spontaneously broken in the high-temperature/high-density regime of quark phase. In order to incorporate the physical aspect of the confinement-deconfinement phase transition in the QM model, the Polyakov loop operator that is related to the free energy of a static color charge is usually introduced and applied. In such an extension, the Polyakov loop operator is introduced as a path ordered Wilson loop of the gauge field in temporal direction,

$$\begin{array}{c}\hfill L\left(\overrightarrow{x}\right)=\mathcal{P}\exp \left[i{\int }_0^{\beta }d\tau A_4(\overrightarrow{x},\tau )\right],\end{array}$$

(4)

where $\mathcal{P}$ denotes path ordering, $\beta =1/T$ being the inverse of temperature and $A_4=i{A}^0$. Accordingly, the color traced Polyakov loops are set as

$$\begin{array}{c}\hfill \Phi =\left({\mathrm{Tr}}_{c}L\right)/N_c,{\Phi }^{*}=\left({\mathrm{Tr}}_c{L}^{†}\right)/N_c,\end{array}$$

(5)

and the Polyakov loop variables are defined as the thermal expectation values $\langle \Phi \rangle$ and $\langle {\Phi }^{*}\rangle$. For simplicity, we still use $\Phi$ and ${\Phi }^{*}$ for their arguments, and within a diagonal representation, $\Phi$ and ${\Phi }^{*}$ are complex scalar fields [11,12]. In a pure gauge theory, the Polyakov loop serves as an order parameter for deconfinement phase transition. It vanishes in the confining hadron phase where the quark has infinite free energy and becomes finite in the deconfined quark matter phase. However, in QCD with dynamical quarks, it should be taken as an approximate order parameter only. The situation is roughly analogous with the chiral condensate which is an exact order parameter in chiral limit for massless quark but only an approximate order parameter for massive quarks.

As the Polyakov loop has been included in the QM model, the dynamical quarks is now coupled with the uniform temporal background gauge field through the following replacement: $D_{\mu }={\partial }_{\mu }-i{A}_{\mu }$ and $A_{\mu }={\delta }_{\mu 0}{A}_0$. In the mean-field approximation, both quantum and thermal fluctuations of the meson and gauge fields are neglected, while the quarks and antiquarks are treated as quantum fields. After integration over the fermions, the quark/antiquark contribution to the grand thermodynamical potential density in the presence of the Polyakov loop can then be written as

$$\begin{array}{c}\hfill {\Omega }_{\overline{\psi }\psi }={\Omega }_{\overline{\psi }\psi }^{\mathrm{vac}}+{\Omega }_{\overline{\psi }\psi }^{\mathrm{th}}=-2N_f{N}_c\int \frac{d^3\overrightarrow{p}}{{\left(2\pi \right)}^3}{E}_q-2N_{f}T\int \frac{d^3\overrightarrow{p}}{{\left(2\pi \right)}^3}\left[\ln g_q^{+}+\ln g_q^{-}\right].\end{array}$$

(6)

Here, $N_f=2$, $N_c=3$, $E_q=\sqrt{{\overrightarrow{p}}^2+m_q^2}$, and the constituent quark masses for u and d quarks are given by $m_q=g{\sigma }_v$ together with ${\sigma }_v\equiv \langle \sigma \rangle$. In vacuum ${\sigma }_v=f_{\pi }=93.0\mathrm{MeV}$. The first term of Eq.(6) denotes the fermion vacuum one-loop contribution, which is divergent and can customarily be renormalized by using the dimensional regularization scheme. The second term of Eq.(6) $g_q^{+}$ and $g_q^{-}$ have been defined upon taking trace over color space

$$\begin{array}{ccc}\hfill g_q^{+}& =& \left[1+3(\Phi +{\Phi }^{*}{e}^{-(E_q-\mu )/T})\times e^{-(E_q-\mu )/T}+e^{-3(E_q-\mu )/T}\right],\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill g_q^{-}& =& \left[1+3({\Phi }^{*}+\Phi e^{-(E_q+\mu )/T})\times e^{-(E_q+\mu )/T}+e^{-3(E_q+\mu )/T}\right].\hfill \end{array}$$

(8)

The first term of Eq.(6) is just the fermion vacuum one-loop contribution to the effective potential at $T=\mu =0$ and the integral is ultraviolet divergent. To regularize this fermion vacuum term, we need to perform the dimensional regularization near three dimensions, $d=3-2ϵ$, in order to isolate the divergences. This yields a resulting potential up to zeroth order in $ϵ$ as given by [67,68]

$$\begin{array}{c}\hfill {\Omega }_{\overline{\psi }\psi }^{\mathrm{vac}}=\frac{N_f{N}_c}{16{\pi }^2}{m}_q^4\left[\frac{1}{ϵ}+ln\frac{{\Lambda }^2}{m_q^2}+ln\left(4\pi \right)-{\gamma }_E+\frac{3}{2}\right],\end{array}$$

(9)

where $\Lambda$ is an arbitrary renormalization scale parameter. Then the effective potential can be renormalized by adding a proper counter term to the Lagrangian as done for the two-flavor case in Refs. [69,70]. Up to some irrelevant constants, the renormalized fermion vacuum loop contribution reads

$$\begin{array}{c}\hfill {\Omega }_{\overline{\psi }\psi }^{\mathrm{vac}}={\Omega }_{\overline{\psi }\psi }^{\mathrm{reg}}=-\frac{N_c{N}_f}{8{\pi }^2}{m}_q^4\ln \left(\frac{m_q}{\Lambda }\right),\end{array}$$

(10)

It is worth pointing out that the parameters $\lambda$ and $\vartheta$ are technically dependent on the arbitrary renormalization scale $\Lambda$, however, as indicated in Refs. [69,70,71,72], the $\Lambda$ dependence can be neatly cancelled out by redefining the parameters in the model, and the physical quantities in the present model are therefore free from ambiguities related to the choice of the renormalization scale.

The phenomenological potential of the Polyakov loop $\mathcal{U}(\Phi ,{\Phi }^{*},T)$ as a function of T, $\Phi$ and ${\Phi }^{*}$ is constructed to reproduce some thermodynamical quantities of the pure-glue gauge theory calculated on the lattice QCD. The first one is that the terms of the potential $\mathcal{U}$ must be invariant under the $Z\left(3\right)$ center symmetry since the pure gauge QCD Lagrangian has center symmetric, and this underlying $Z\left(3\right)$ center symmetry is broken spontaneously in the pure gluonic case. Secondly, because there is no reason for any asymmetry between $\Phi$ and ${\Phi }^{*}$ in a pure gauge theory, the potential $\mathcal{U}$ should have additional requirement with the symmetry under the exchange between $\Phi$ and ${\Phi }^{*}$. Finally, as mentioned as above, $\Phi$ and ${\Phi }^{*}$ are treated as order parameters for the deconfinement phase transition, we thus require that the thermal expectation values $\Phi$ and ${\Phi }^{*}$ evaluated at the minimum of the potential $\mathcal{U}$ are $\Phi ={\Phi }^{*}=0$ when temperature is low, on the contrary, when the temperature is high, these values should become finite and asymptotically approach to a constant, e.g. the unity 1 for the logarithmic form of the Polyakov-loop potential. For the specific functional form of the potential $\mathcal{U}$, there still exist various possibilities in the literature with similar properties [14,15,73,74,75]. For the physical interest in the present work, we prefer to use the simplest polynomial potential based on a Ginzburg-Landau ansatz as proposed in Ref. [73]

$$\begin{array}{c}\hfill \frac{\mathcal{U}(\Phi ,{\Phi }^{*},T)}{T^4}=-\frac{b_2\left(T\right)}{4}{(|\Phi |}^2+|{\Phi }^{*}{|}^2)-\frac{b_3}{6}({\Phi }^3+{\Phi }^{*3})+\frac{b_4}{16}{\left(\right|\Phi |}^2+|{\Phi }^{*}{{|}^2)}^2,\end{array}$$

(11)

with a temperature dependent coefficient

$$\begin{array}{c}\hfill b_2\left(T\right)=a_0+a_1\left(\frac{T_0}{T}\right)+a_2{\left(\frac{T_0}{T}\right)}^2+a_3{\left(\frac{T_0}{T}\right)}^3.\end{array}$$

(12)

The potential parameters are adjusted to confront with the lattice data for the pure gauge theory thermodynamics so that the equation of state and the Polyakov-loop expectation values are reproduced on the behavior of the Polyakov loop as a function of temperature. This yields the following values

$$\begin{array}{c}\hfill a_0=6.75,a_1=-1.95,a_2=2.625,\\ \hfill a_3=-7.44,b_3=0.75,b_4=7.5.\end{array}$$

(13)

Moreover, for pure $SU\left(3\right)$ gauge theory, the deconfinement phase transition is of first-order at the critical temperature at $T_0=270$ MeV. However, we note that the potential presented above does not take into account the backreaction of the dynamical quarks on the gauge sector and therefore the influence of the quarks on the deconfinement phase transition. It is then natural to let $T_0=T_0\left(N_f\right)$ in order to guarantee that the behavior in the glue sector also depends on the number of quark flavors or even the baryon chemical potential [14,76]. In the present study we shall use $T_0=208$ MeV for two quark flavors [14,70].

Now the thermodynamic grand potential for the PQM model in the mean-field approximation with the fermion vacuum one-loop is given by

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{MF}}(T,\mu )=U(\sigma ,\overrightarrow{\pi })+\mathcal{U}(\Phi ,{\Phi }^{*},T)+{\Omega }_{\overline{\psi }\psi }.\end{array}$$

(14)

The equations of motion for the mean fields $\sigma$, $\Phi$ and ${\Phi }^{*}$ follow from the stationarity conditions

$$\begin{array}{c}\hfill \frac{\partial {\Omega }_{\mathrm{MF}}}{\partial {\sigma }_v}=0,\frac{\partial {\Omega }_{\mathrm{MF}}}{\partial \Phi }=0,\frac{\partial {\Omega }_{\mathrm{MF}}}{\partial {\Phi }^{*}}=0,\end{array}$$

(15)

which yield the temperature and chemical potential dependence of the chiral condensate $\sigma (T,\mu )$ as well as of the Polyakov loop expectation values $\Phi (T,\mu )$ and ${\Phi }^{*}(T,\mu )$. The parameters in the model are specified by the vacuum properties and we will take the values $m_{\pi }=138$ MeV, $m_{\sigma }=500$ MeV. The coupling constant g is usually fixed by the constituent quark mass in vacuum around $m_q=333$ MeV and $f_{\pi }=93$ MeV.

3. Phase Structure

The solutions of the above coupled gap equations (15) at zero and finite chemical potential yield the chiral condensate $\sigma$ and the Polyakov loop expectation values $\Phi$ and ${\Phi }^{*}$ as a function of the temperature. The temperature behaviors of these order parameters could be used to reveal the nature of chiral and deconfinement phase transition. In Figure 1 we study the temperature variation of the chiral condensate and Polyakov loop expectation values at $\mu =0$ MeV and $\mu =308$ MeV. In order to illustrate the subtle differences between these two types of QCD phase transition, we would like to split our studies into two parts, one is the chiral phase transition related to the chiral symmetry, another is the deconfinement phase transition that is directly relevant to the $Z\left(3\right)$ center symmetry.

Let us first discuss the restoration of chiral symmetry in high energy. The temperature behavior of the chiral condensate shows that the system experiences a quite smooth crossover transition at zero chemical potential. On the contrary, for a relatively larger chemical potential, i.g. $\mu =308$ MeV, there is a first-order phase transition due to the fact that the chiral condensate makes a sharp jump across the gap of the condensate when the temperature is nearby its critical value. Therefore, by compared two curves of the chiral condensate $\sigma$ in Figure 1, we find that the chemical potential would make the variation of chiral order parameter with temperatures sharp and sharper so that the crossover is to eventually change to a second-order phase transition at a critical endpoint (CEP), where the crossover turn into the first-order phase transition. The existence of the QCD critical endpoint guarantees the occurrence of the first-order phase transition, and vice versa. From Figure 1, for both a crossover and a first-order chiral phase transition, the temperature derivative of the chiral condensate $\sigma$ has only one typical peak at a specific temperature, which is usually defined as the critical temperature. Hence, for zero quark chemical potential, the chiral restoration occurs at $T_{\chi }^c\simeq 187$ MeV, whereas for a relatively larger quark chemical potential $\mu =308$ MeV, the critical temperature moves to the lower temperature around $T_{\chi }^c\simeq 41$ MeV.

This is unlike the chiral phase transition, in which the critical temperature is well defined at least for the case of a first-order phase transition. On the contrary, we are still ambiguous to conclusively identify the deconfinement phase transition through the Polyakov loop expectation values $\Phi$, ${\Phi }^{*}$ or their temperature derivatives [16,75], even though the chemical potential is very large. In order to explicate this problem, the temperature variations of the Polyakov loop expectation values $\Phi$ and ${\Phi }^{*}$ at $\mu =0$ MeV and $\mu =308$ MeV are demonstrated in Figure 1. From this figure, by taking the $\Phi$ as example, we can find that the evolution of the Polyakov loop expectation value $\Phi$ as a function of the temperature exhibits a rather smooth behaviour for zero chemical potential. For a larger chemical potential at $\mu =308$ MeV, $\Phi$ becomes more smoother and flatter except that there is a small jump coincident with the jump of the chiral order parameter $\sigma$. Because the temperature derivative ${\Phi }^{\prime }=d\Phi /dT$ generally has one more peak in calculations, especially when $\mu$ is high [16,70], it is obscure to determine the critical temperature for the deconfinement phase transition through the peaks of the Polyakov loop variation. The problem will arise again and become worse for chiral and deconfinement phase transition when the strange quark s is included in POM model [15,16].

In order to reveal more information about the natural properties of the chiral and deconfinement phase transitions, the effective potential itself is an important method and a powerful theoretical tool. Besides solving the equations of motion (15), based on the shape or geometry of an effective potential, we can directly search a global minimum of the effective potential ${\Omega }_{\mathrm{MF}}$ in Eq.(14) to get the expectation values of the chiral condensate and the Polyakov loop fields, or more interestingly, the local minima of the effective potential, which are corresponding to the false vacua during phase transitions. For convenience, we refer a study based on the geometry of effective potential as a geometric approach. For a single order parameter, such as the QM model, a geometric approach can precisely tell two types of chiral phase transition and simultaneously give out two minima of effective potential [55,77]. However, for the PQM model there are three order parameter variables $\sigma$, $\Phi$ and ${\Phi }^{*}$ in the grand canonical potential ${\Omega }_{\mathrm{MF}}$ in Eq.(14), so that it is extremely challenging to demonstrate and reveal the thermal effective potential via evolving these variables simultaneously in such a big parameter space. Inspired by previous studies in the PQM model [71,72], in order to simplify the problem and provide a more intuitive insight into the physics, the number of the active variable in the effective potential at finite temperature and chemical potential is to be set as one at a time. This means that the effective potential has been divided into two directions: (1) the gluonic direction, and (2) the mesonic direction. In the gluonic direction, the Polyakov loop $\Phi$ (or ${\Phi }^{*}$) is considered as the active variable while the $\sigma$ and ${\Phi }^{*}$ (or $\Phi$) fields remain on their expectation values all the time. By contrast, in the mesonic direction, the $\sigma$ field is taken as a active variable in the grand canonical potential ${\Omega }_{\mathrm{MF}}$ while fixing the Polyakov loop fields on their expectation values $\Phi$ and ${\Phi }^{*}$.

For simplicity, we firstly take the Polyakov loop $\Phi$ as a typical example in the following discussion, the qualitative result presented here, of course, can be applied to that of the ${\Phi }^{*}$ field too. The gluonic case is shown in Figure 2, where the left panel is the scaled grand canonical potentials ${\Omega }_{\mathrm{MF}}$ over $T^4$ as a function of the Polyakov loop $\Phi$ for zero chemical potential, and the right panel is the scaled grand canonical potentials ${\Omega }_{\mathrm{MF}}$ over $T^4$ as a function of the Polyakov loop $\Phi$ for $\mu =310$ MeV. From Figure 2 it is obvious that the grand canonical potentials for both $\mu =0$ MeV and $\mu =310$ MeV share similar behaviors: there is only one minimum for each of the grand canonical potentials. These minima, which are to be identified as the expectation values $\Phi$ of the Polyakov loop field, move to their higher values smoothly and slowly with the increase of the temperature. Even though the chemical potential is so large as $\mu =310$ MeV, when the chiral phase transition is of first-order, the grand canonical potential ${\Omega }_{\mathrm{MF}}$ still exhibits a trivial feature with a single minimum. This indicates that the disconnection of the order parameter $\Phi$ in Figure 1 for $\mu =308$ MeV and $T\simeq 41$ MeV should be driven by that of the chiral order parameter $\sigma$, such a sudden jump does not supported by the grand canonical potential of the Polyakov loop field itself. So that the peak in the Polyakov loop field derivative around the critical temperature at $T=T_{\chi }^c\simeq 41$ MeV is a fake signal for the definition of the deconfinement phase transition. This is also the reason why a critical temperature of the deconfinement conversion should be given by the another pseudo peak corresponding to $\Phi \left(T\right)/\Phi (T\to \infty )>1/2$ as discussed in Refs. [70,71].

Theoretically, the trivial behavior of the grand canonical potential as a function of the Polyakov loop can be interpreted qualitatively if we go back to the original pure-gluonic potential $\mathcal{U}$ in Eq. (11). It is easy to show that the potential $\mathcal{U}$ displays only one minimum at $\Phi =0$ as long as the temperature is less than the critical temperature at $T_0=208$ MeV. Only when T is larger than the critical value $T_0$, the potential $\mathcal{U}$ will develop a typical first-order potential with two local minima separated by a global maximum (or a barrier). It seems that the inclusion of the quark and meson fields in the grand canonical potential ${\Omega }_{\mathrm{MF}}$ does not change the nature of the deconfinement phase dramatically. As demonstrated in Figure 2, with the increase of the temperature, the global minimum of the potential moves continuously from its low value to a high one, but the characteristic shape of the potential remains intact as long as $T<T_0$.

Let us now investigate how the grand canonical potentials ${\Omega }_{\mathrm{MF}}$ evolving with the chiral order parameter $\sigma$ for different chemical potentials by fixing the Polyakov loops on their expectation values. The scaled grand canonical potentials are shown in Figure 3 as a function of $\sigma$ for $\mu =306$ MeV and $\mu =310$ MeV, respectively. From Figure 3, one clearly observes the characteristic pattern of a first-order phase transition potential when $\mu$ is very large: two minima corresponding to phases of restored and broken chiral symmetry are separated by a potential barrier and they will become degenerate at $T=T_{\chi }^c$. Chiral symmetry is approximately restored as $T>T_{\chi }^c$, where the previous local minimum at relatively low $\sigma$ value becomes the global/absolute minimum as shown in Figure 3. On the other hand, when the temperature is below the critical one $T_{\chi }^c$, the global minimum locates in a relatively larger sigma and the constituent quark becomes massive. In this time, the chiral symmetry is spontaneously broken. As a result, we can conclude that there is a firs-order chiral phase transition when chemical potentials are both at $\mu =306$ MeV and $\mu =310$ MeV. Furthermore, when the temperature is at $T=T_{\chi }^c$, the height of the potential barrier is about $0.66$$\mathrm{MeV}/{\mathrm{fm}}^3$ as $\mu =310$ MeV, but it will drop sharply to $0.30$$\mathrm{MeV}/{\mathrm{fm}}^3$ when the chemical potential reduces to $\mu =306$ MeV. So that with decreasing of the chemical potential, it is believed that the potential barrier combined two minima becomes flat and flatter. Consequently, the first order chiral phase transition turns into a second order phase transition at a specific point called as the CEP. In that moment the barrier disappears and the chiral phase transition exhibits a second-order quality.

For a large chemical potential, the grand canonical potentials ${\Omega }_{\mathrm{MF}}$ as a function of the chiral order parameter $\sigma$ in Figure 3 really guarantee the existence of the first-order chiral phase transition. But the nature of the phase transition revealed by the shapes of effective potentials is fairly distinct from each other. For $T>T_{\chi }^c$, both cases show that there have the critical spinodal temperatures $T_{\mathrm{sp}}$ since the local minima of effective potential are going to disappear at certain temperature. This suggests that the hadron-quark phase transition is a weak first-order phase transition if the system is beating up from low temperature to a high one. However, while the temperature T is below the critical value $T_{\chi }^c$, the evolutions of the potential with the order parameter $\sigma$ for these two chemical potentials exhibit quite different features. From the left panel in Figure 3, for temperatures very near the critical temperature $T_{\chi }^c\simeq 48.9$ MeV, the grand canonical potentials ${\Omega }_{\mathrm{MF}}$ for $\mu =306$ MeV indeed displays a local minimum which is separated by a barrier from another local minimum. But as the temperature is lowered, the local minimum at low $\sigma$ value approaches the intervening maximum. These two minima of the potential meet and form an inflection point at the spinodal temperature $T_{\mathrm{sp}}\simeq 37$ MeV. Accordingly, there has only one minimum in the grand canonical potentials ${\Omega }_{\mathrm{MF}}$ when $T\le T_{\mathrm{sp}}$. Because of the disappearance of the barrier between two local minima in the effective potential, this kind of the phase transition sometime can be denoted as a weak first-order conversion [63,66]. On the contrary, the increase of the chemical potential normally leads to a larger barrier between the two degenerate minima of the grand canonical potentials ${\Omega }_{\mathrm{MF}}$ when T is close to $T_{\chi }^c$. In other words, the first order phase transition would be strengthened with the growth of the chemical potential, and a larger barrier implies a larger latent heat. For the case of $\mu =310$ MeV, the grand canonical potentials ${\Omega }_{\mathrm{MF}}$ as a function of the chiral order parameter $\sigma$ for different temperatures are shown in the right panel in Figure 3. One significant difference from the case of $\mu =306$ MeV is that there always exists a potential barrier between two minima of the potential, no matter what the temperature is chosen. This indicates that the chiral phase transition for a quark-hadron phase transition is now a strong first-order phase transition when $\mu =310$ MeV [63,78].

The chiral phase transition of u and d quarks in the $T-\mu$ plane is shown in Figure 4. It is found that for two light flavors, there is a crossover in the low-density region and a first-order phase transition in high-density region, and the critical end point in the middle position is around $(T_{\mathrm{E}},{\mu }_{\mathrm{E}})\simeq (301.4\mathrm{MeV},62.1\mathrm{MeV})$. In the $T-\mu$ plane, a vast part of the QCD phase diagram is a crossover, only a very small window is left for the first order phase transition. As compared with previous studies in the QM model, the incorporation of effects of deconfinement in the language of the Polyakov-loop fields merely increase the critical temperature for chiral phase transition, while leaves the nature of the chiral phase transition unchanged. Since we are interested in the dynamics of the first order phase transition in the present study, the lines for the deconfinement crossover transitions are omitted in Figure 4.

In order to present a more detailed description about the chiral first-order phase transition, two additional spinodal lines are also included in Figure 4. These two spinodal lines are usually introduced to characterize the region of spinodal instability for a weak first-order phase transition. From the Figure 4, both the lower and upper spinodal lines rise up with the reduction of the chemical potential, but the distance between these two lines become small and small. In the end, the spinodal lines together with the coexistent critical line will bend forward to the CEP and terminate at that point. Furthermore, according to the low spinodal line as $T<T_{\chi }^c$, the structure of phase diagram below the coexistent critical line can be further classified into two categories: a week first-order phase transition and a strong one. For a weak first-order phase transition, it is usually characterized by the low barrier and the existence of the critical spinodal temperature $T_{\mathrm{sp}}$, whereas for a strong first-order phase transition, there has an effective potential with a zero-temperature potential barrier. Therefore, as shown in Figure 4, along horizontal axis, the corresponding chemical potential, when a week first-order phase transition turns into a strong one, is established as ${\mu }_{\mathrm{c}}=309$ MeV. For $\mu >{\mu }_{\mathrm{c}}$, it is a strong first-order phase transition and a potential barrier separated two minima never disappears with the decrease of the temperature, conversely, the local minimum aside from the global minimum will gradually disappear and there is no potential barrier for $T<T_{\mathrm{sp}}$. On the other side of the coexistent critical line, since there exists the critical spinodal temperature, the hadron-quark phase conversion should be always of weakly first-order.

4. Nucleation in the PQM Model

4.1. Homogeneous Thermal Nucleation

For a first-order phase transition, the hallmark of a potential is an energy barrier separating the two minima. For some potentials a barrier already exists at the classical tree-level, for others it could be induced by loop expansions or the integration of other dynamical fields. When the temperature is very close to its critical value $T_c$, these two minima are degenerate. However, as the temperature is deviating from $T_c$, one of the minimum becomes an absolute/global minimum (a true vacuum), the other is now changed to a local minimum (a false vacuum). The decay of a false vacuum consists of a first-order phase transition where the order parameter changes from a false vacuum to a true vacuum. The conversion may occur through nucleation of bubbles caused by quantum or thermal fluctuations, depending on how the ambient temperature compares to the energy barrier. In the present work, we were mostly concerned with the limit that thermal fluctuations dominate quantum fluctuations.

Based on the theory of homogeneous thermal nucleation, at finite temperature, we start with a four-dimensional Euclidean action for a scalar order parameter field, e.g. a scalar field $\varphi$,

$$S_E\left[\varphi \right]={\int }_0^{\beta }d\tau \int d{r}^3\left[\frac{1}{2}{\left(\frac{\partial \varphi }{\partial \tau }\right)}^2+\frac{1}{2}{\left(\nabla \varphi \right)}^2+V\left(\varphi \right)\right],$$

(16)

where the potential has a local minimum at $\varphi ={\varphi }_f$ and a global minimum at $\varphi ={\varphi }_t$. Note that a generalization of this problem to more than one field is rather straightforward by replacing the $\varphi$ field with $\Phi =({\varphi }_1,{\varphi }_2,\dots ,{\varphi }_N)$. The path integral expression for the partition function is given by

$$Z=\int \mathcal{D}\varphi e^{-S_E\left[\varphi \right]}.$$

(17)

In a saddle point approximation, the order parameter $\varphi$ is to expand around its classical equation of motion (EOM) as $\varphi =\overline{\varphi }+\phi$. Since $\overline{\varphi }$ is a solution of EOM, it will minimise the Euclidean action and dominate the path integral. Thus, in the vicinity of the saddle point, the Euclidean can be written in term of a Taylor series up to the second order as

$$S_E\left[\varphi \right[\simeq S_E\left[\overline{\varphi }\right]+\frac{1}{2}{S}_E^{\prime \prime }\left[\overline{\varphi }\right]{\phi }^2+\dots ,$$

(18)

where we note that the first-order variation is absent as the equations of motion. The second term in the above equation is a Gaussian integral, which can be evaluate exactly. After performing the functional Gaussian integration in the partition function in Eq.(17), one gets the following expression

$$Z\sim e^{-S_E\left[\overline{\varphi }\right]}{\left(\det S_E^{\prime \prime }\left[\overline{\varphi }\right]\right)}^{-1/2}.$$

(19)

Then our next ultimate task is to find the solution of the EOM and substitute back into the partition function to obtain the free energy of the system through the standard formula $\mathcal{F}=-TlnZ$. The nucleation rate $\Gamma$ could in turn be determined exactly from this Euclidean observable, the free energy, as $\Gamma =-2\mathrm{Im}\left(\mathcal{F}\right)$.

To compute the classical Euclidean field equation of motion, we extremize the Euclidean action

$$\frac{\delta S_E\left[\varphi \right]}{\delta \varphi }{|}_{\varphi =\overline{\varphi }}=0.$$

(20)

Then the classical equation of motion becomes

$$\left(\frac{{\partial }^2}{\partial {\tau }^2}+{\nabla }^2\right)\overline{\varphi }=V^{\prime }\left(\overline{\varphi }\right),$$

(21)

and solutions of this EOM with minimum energy are expected to be spherically symmetric [79], then the critical field configuration of the field $\varphi$ has $O\left(4\right)$ symmetry. However, for sufficiently high temperature, as argued by Linde [54], the problem becomes approximately three dimensional by performing the integration over the Euclidean time coordinate. Accordingly, the Euclidean action is now rewritten as

$$S_E\left[\overline{\varphi }\right]\equiv \frac{S_3\left[\overline{\varphi }\right]}{T},$$

(22)

where $S_3\left[\overline{\varphi }\right]$ is the rescaled action of a three-dimensional field theory. In this case, the saddle point will be given by the $O\left(3\right)$ symmetric solution of

$$\frac{d^2\overline{\varphi }\left(r\right)}{d{r}^2}+\frac{2}{r}\frac{d\overline{\varphi }\left(r\right)}{dr}=V^{\prime }\left(\overline{\varphi }\right).$$

(23)

This is a second-order partial differential equation and has many solutions that extremize the Euclidean action. In the semi-classical approximation for false vacuum decay, not all saddle points are taken into account, what we are really interested in is a nontrivial saddle-point solution which interpolated between the false vacuum ${\varphi }_f$ and some field value on the other side of the barrier. The field configuration ${\varphi }_b$, often called bounce or instanton, is associated with the nucleation bubble embedded in the homogeneous false vacuum. The interior of the bubble consists of the true vacuum, whereas, outside the bubble, the field $\overline{\varphi }$ should arrive at its false vacuum. Therefore, for the bounce, it is reasonable to employ the boundary condition of the false vacuum as $\underset{r\to \infty }{lim}\overline{\varphi }\left(r\right)={\varphi }_f$. Since the EOM in Eq.(23) is a second-order ordinary differential equation, another boundary condition $\frac{d\overline{\varphi }\left(r\right)}{dr}{|}_{r=0}=0$ should be imposed on this EOM by the requirement of finite energy at the origin.

Once we obtain the bounce ${\varphi }_b$, the explicit expression for the three-dimensional Euclidean action $S_3$ is given by

$$S_3=4\pi \int dr{r}^2\left[\frac{1}{2}{\left(\frac{d{\varphi }_b}{dr}\right)}^2+V\left({\varphi }_b\right)\right],$$

(24)

and the surface tension, an one-dimensional energy of the bounce, is evaluated accordingly as

$$\Sigma =\int dr\left[\frac{1}{2}{\left(\frac{d{\varphi }_b}{dr}\right)}^2+V\left({\varphi }_b\right)\right].$$

(25)

In the practical calculations, It is worth noting that if the false vacuum has a non-zero potential energy, an additional term $-V\left({\varphi }_f\right)$ should be included in the $S_3$ action and the surface tension $\Sigma$. Finally, based on a saddle-point approximation around the bounce solution, the false vacuum will decay to the true vacuum with the rate per unit time per unit volume given at finite temperature by [49,50,54]

$$\Gamma =\frac{{\omega }_{-}}{2\pi }{\left(\frac{S_3}{2\pi T}\right)}^{3/2}{\left[\frac{{\mathrm{Det}}^{\prime }(-{\nabla }^2+V^{\prime \prime }\left({\varphi }_b\right))}{\mathrm{Det}(-{\nabla }^2+V^{\prime \prime }\left({\varphi }_f\right))}\right]}^{-1/2}{e}^{\left(-\frac{S_3}{T}\right)}.$$

(26)

Here ${\det }^{\prime }$ signifies that the zero eigenvalues associated with the translation symmetry of the bubble are omitted, ${\omega }_{-}$ is the eigenvalue of the negative mode and ${\varphi }_b$ denotes the bounce solution of the EOM. Traditionally, the calculation and evaluation of the prefactor in Eq.(26) is a nontrivial matter, but for the most typical cases, the Euclidean action appearing in the exponent is large, and thereby the exponential is very small. Hence the bounce solution and its action $S_3$ play the dominant role in determining the nucleation rate, the prefactor can be roughly estimated by dimensional analysis and then approximately expressed as $T^4$ or $T_c^4$ for simplicity [59,80].

4.2. The Minima of the Effective Potential

To give an intuitive description of thermal phase transition in particle and nuclear physics, the classical potential $V\left(\varphi \right)$ should be generalized to the effective potential $V_{\mathrm{eff}}(\varphi ;T,\mu )$ in order to incorporated thermal and quantum corrections. In equilibrium, given the temperature and the chemical potential, the effective potential $V_{\mathrm{eff}}$ is then a function of the order parameter $\varphi$. The expected value or the condensate for the order parameter is determined by minimizing the effective potential with the order parameter field, $\partial V_{\mathrm{eff}}/\partial \varphi =0$. If there have many extrema, we need to calculate the second derivative to make sure it is a minimum point not a maximum point in the full parameter space. This is the standard procedure to search the minima of the effective potential for a single order parameter field.

However, if there are many order parameters ($N>1$) in the effective potential, the above method is to be extended to a more general but complicated procedure as follows. Firstly, the gap equations will be obtained by minimizing $V_{\mathrm{eff}}({\varphi }_1,{\varphi }_2,\dots ,{\varphi }_N)$ with respect to the order parameters ${\varphi }_i$ ($i=1,2,\dots ,N$),

$$\frac{\partial V_{\mathrm{eff}}}{\partial {\varphi }_1}=\frac{\partial V_{\mathrm{eff}}}{\partial {\varphi }_2}=\dots =0.$$

(27)

The solutions of this equation give the points where the effective potential has its extrema in the full order parameter space. Secondly, the $N\times N$ Hessian matrix,

$${\mathcal{H}}_{i,j}=\frac{{\partial }^2{V}_{\mathrm{eff}}}{\partial {\varphi }_i\partial {\varphi }_j},$$

(28)

has to be evaluated at the points where the gap equation (27) is satisfied. At last, as long as we get the Hessian matrix, the points corresponding to the minima of the effective potential can be selected by checking the signs of the eigenvalues of the Hessian matrix. If they are all positive, this point is a minimum, otherwise, it is either saddle point or a maximum. In a practical calculation, this method works very well while the effective potential is a "good” potential with a large gradient around the minimum point, but the method becomes less efficient or even fails to give out the minima when the effective potential becomes very flat, such as the effective potential nearby the CEP. Besides the above standard method to explore the minima of the effective potential [64], for the physical interest in the present study, we propose an alternative method to find the minima of the effective potential based on a geometric approach, especially the local minima existed in the potential.

For the first step, for a potential with more than one order parameters $V(\Phi )$, supposed we already know one of the minimum of the potential, namely $\overline{\Phi }=({\overline{\varphi }}_1,{\overline{\varphi }}_2,\dots ,{\overline{\varphi }}_N)$. Here, we suggest to chose the point at the global minimum of the potential because it is easy to locate by scanning the potential directly in the order parameter space. In what follows, we denote these values as the expectation values. The standard procedure is to describe as follows: (1) the N-dimensional potential can be extremely simplified by treating the ${\varphi }_1$ field as a variable while fixing other fields on their expectation values ${\overline{\varphi }}_i$ ($i=2,\dots ,N$). Now, the N-dimensional problem has been changed to an one-dimensional problem, and it is easy to dig up other local minima in the ${\varphi }_1$ field direction, i.e. ${\varphi }_1={\varphi }_1^{\alpha },{\varphi }_1^{\beta },\dots$. (2) the ${\varphi }_2$ field is considered as a variable whereas other fields maintain their expectation values all the time, then the local minima of the effective potential $V\left({\varphi }_2\right)$ can also be found as the case for a single order parameter field. These minima of the potential $V\left({\varphi }_2\right)$ are expressed as ${\varphi }_2={\varphi }_2^{\alpha },{\varphi }_2^{\beta },\dots$. (3) by using the same algorithm, we can find all minima in the direction of the ${\varphi }_3$ field as ${\varphi }_3={\varphi }_3^{\alpha },{\varphi }_3^{\beta },\dots$, and so on. (4) after obtaining the minima for the potential in every direction, the landscape of the potential can be roughly reconstructed by combining the potentials in different direction all together. In particular, if there is only one minimum in some directions, those parameters can be discard and the numbers of the dimension is to be reduced accordingly.

For the second step, the primary points obtained through the above algorithm in the first step are not necessary the true local minima of the potential. To trace a true local minimum of the potential, we have to increase the dimension of the parameter space or the number of the order parameter one by one. (1) starting from the initial point at $V=V({\varphi }_1^{\alpha \left(0\right)},{\varphi }_2,{\overline{\varphi }}_3,,\dots )$, by taking the ${\varphi }_2$ field as a variable, we can find the minima of the potential when other fields are fixed at ${\varphi }_1^{\alpha }={\varphi }_1^{\alpha \left(0\right)}$ and ${\varphi }_i={\overline{\varphi }}_i$ for $i=3,\dots ,N$. Let us assume one of the minima is ${\varphi }_2^{\alpha \left(0\right)}$. (2) we treat the ${\varphi }_1$ as a variable when fixing other fields at ${\varphi }_2={\varphi }_2^{\alpha \left(0\right)}$ and ${\varphi }_i={\overline{\varphi }}_i$ for $i=3,\dots ,N$. If one of the minima of the potential at ${\varphi }_1={\varphi }_1^{\alpha \left(1\right)}$ is equal to the original point at ${\varphi }_1={\varphi }_1^{\alpha \left(0\right)}$, we can conclude that the point (${\varphi }_1^{\alpha \left(0\right)}$, ${\varphi }_2^{\alpha \left(0\right)}$) at the (${\varphi }_1,{\varphi }_2$) plane is the true local minimum of the potential in the 2-dimensional parameter space. Otherwise, we need to repeat the first stage again by setting ${\varphi }_1={\varphi }_1^{\alpha \left(1\right)}$ in the potential in order to gain ${\varphi }_2^{\alpha \left(1\right)}$. After several iterations, we can luckily arrive at a right point when the value of the ${\varphi }_1$ satisfies ${\varphi }_1^{\alpha \left(k\right)}={\varphi }_1^{\alpha (k+1)}$, then the point (${\varphi }_1^{\alpha \left(k\right)}$, ${\varphi }_2^{\alpha \left(k\right)}$) is the exact local minimum what we are looking for. Otherwise, after several iterations, if the value of the ${\varphi }_1^{\alpha (k+1)}$ is far away from its previous value at ${\varphi }_1={\varphi }_1^{\alpha \left(k\right)}$ with the increase of the number of iteration k, the starting point (${\varphi }_1^{\alpha \left(0\right)}$, ${\varphi }_2^{\alpha \left(0\right)}$) is a saddle point, it should be discarded accordingly. (3) in the last stage, once we get the minimum point in the ${\varphi }_1$ and ${\varphi }_2$ plane, akin to the first and second stages, we can extract the local minimum of the potential in a 3-dimensional order parameter space by employing the field ${\varphi }_3$ as a variable. Ultimately, one by one, step by step, we can gain all minima of the potential in the full order parameter space.

4.3. The Tunneling between Two Minima

Generally, for a given potential, it is a tough job to find all minima in the N-dimensional order parameter space if $N\ge 2$. As discussed in above subsection, we can use two different algorithms to do calculations. For the first algorithm, it is straightforward and easy to apply. The advantage of this algorithm is that we can use it to find all minima of the potential very efficiently and quickly. But for some problems, since we don’t know the geometry of the potential, even though we have exact positions of two minima, it does not help to find the instanton connected with these two minima. For example, there is a two-dimensional potential in Fig.2 in Ref. [81] that does not have any instanton solution between these two local minima. On the other side, for the second algorithm developed by us, it is complicated or even awkward to be applied, due to the fact that we have to find a local minimum of the potential step by step. Nevertheless, the reward for such a sophisticated procedure is that we can change a multi-dimensional problem to a one-dimensional problem, which is usually easy to solve. More specifically, since the tunneling occurs in two nearby minima, in the practical calculation, we actually do not need to know all minima of the potential. These two minima which are close to each other and can be connected through the instanton are more important and crucial. Therefore, to study the problem of the false vacuum decay, we prefer to adopt the second algorithm to find two minima lived nearby.

To demonstrate the algorithm developed in the above section, we firstly take an interesting two-dimensional potential which has many minima as a typical example. The potential which is formerly proposed in Ref. [81] has the form

$$V({\varphi }_1,{\varphi }_2)=sin({\varphi }_1-{\varphi }_2)+\frac{1}{2}cos({\varphi }_1+{\varphi }_2)+cos3({\varphi }_1+{\varphi }_2)+2cos3({\varphi }_1-\frac{{\varphi }_2}{2}).$$

(29)

There are many minima in this potential, but we need not to know all of them in the parameter space. What we are interested are those they are close to each other. Assume we already know one of the minima at the point $p_1=(2.39338,2.82768)$, our next task is to find out the other minimum close to this one. Following the second algorithm developed in above, by taking ${\varphi }_1$ as a variable, the potential $V\left({\varphi }_1\right)$ has three local minima at $({\varphi }_1^{\alpha }=0.376583,{\varphi }_1^{\beta }=2.39338,{\varphi }_1^{\gamma }=4.56047)$ when fixed the second order parameter at ${\varphi }_2=2.82768$. This implies that there have three local minima in the direction of the ${\varphi }_1$ field. In the meantime, if we fix ${\varphi }_1=2.39338$, the potential $V\left({\varphi }_2\right)$ displays others three minima in the direction of the ${\varphi }_2$ field. So that there probably exist four local minima around the point $p_1$. It is worth noting that if we start from the initial point $p_1$, we can only find the four minimum points around this point, however, if we run over the minimum points in the potential one by one, the missing counting problem can be resolved and all missing minimal points can be recovered accordingly. Next, as mentioned in above, the initial point $({\varphi }_1^{\alpha \left(0\right)}=0.376583,{\varphi }_2^{\alpha \left(0\right)}=2.39338)$ could not be a true minimum in the potential. We need to take the ${\varphi }_2$ as a variable to search the local minima of the potential when fixing the ${\varphi }_1$ field at ${\varphi }_1^{\alpha \left(0\right)}=0.376583$. After a simple calculation, we can find one of the minimum of the potential is at ${\varphi }_2^{\alpha \left(1\right)}=2.74049$. Then by fixing the ${\varphi }_2$ field at ${\varphi }_2^{\alpha \left(1\right)}=2.74049$, we can obtain a minimum point at ${\varphi }_1^{\alpha \left(1\right)}=0.376074$. Through the iteration method, we can get ${\varphi }_2^{\alpha \left(2\right)}=2.74048$ and ${\varphi }_1^{\alpha \left(2\right)}=0.376074$. In this time, since we have ${\varphi }_1^{\alpha \left(2\right)}={\varphi }_1^{\alpha \left(1\right)}$, we can conclude the point at $p_2=({\varphi }_1^{\alpha \left(2\right)},{\varphi }_2^{\alpha \left(2\right)})=(0.376074,2.74048)$ is a true local minimum of the potential in the ${\varphi }_1-{\varphi }_2$ plane. Accordingly, the tunneling rate between these two minima $p_1$ and $p_2$ is described by the system of coupled ordinary differential equations

$$\begin{array}{ccc}\hfill \frac{d^2{\varphi }_1}{d{r}^2}+\frac{2}{r}\frac{d{\varphi }_1}{dr}& =& \frac{\partial V({\varphi }_1,{\varphi }_2)}{\partial {\varphi }_1},\hfill \end{array}$$

(30)

$$\begin{array}{ccc}\hfill \frac{d^2{\varphi }_2}{d{r}^2}+\frac{2}{r}\frac{d{\varphi }_2}{dr}& =& \frac{\partial V({\varphi }_1,{\varphi }_2)}{\partial {\varphi }_2},\hfill \end{array}$$

(31)

with the boundary conditions at $\underset{r\to \infty }{lim}\varphi \left(r\right)=p_1$ and $\frac{d\varphi \left(r\right)}{dr}{|}_{r=0}=0$ because the point $p_2$ has a lower potential than $p_1$ and is regarded as a true vacuum.

For a second example, we turn back to our present study, the effective potential ${\Omega }_{\mathrm{MF}}(T,\mu )$ in Eq.(14) now has three order parameters. Traditionally, it is hard to give out an intuitive landscape picture of the potential. Hence, it is necessary to simplify a multi-dimensional problem to an one-dimensional problem. According to the second algorithm, for a given temperature and chemical potential, we firstly find the global minimum of the effective potential by either solving the gap equations or scanning the effective potential directly. The exact values of the order parameters in the global minimum of the potential are called as the expectation values for the $\sigma$, $\Phi$ and ${\Phi }^{*}$ fields. But the big challenge is to find the possible local minima nearby this global one. Following the second algorithm, we can treat the $\sigma$ field as a variable, but fixed other parameters $\Phi$ and ${\Phi }^{*}$ at their expectation values, and the results have been depicted in Figure 3. On the contrary, the order parameter $\Phi$ can also be taken as a variable, but setting the $\sigma$ and ${\Phi }^{*}$ fields on their expectation values instead as shown in Figure 2. From these figures, it can be found that the local minima of the effective potential only exists in the direction of the $\sigma$ field and there is no minimum in both $\Phi$ and ${\Phi }^{*}$ directions. Therefore, we can dissociate one of the equation of motion for the $\sigma$ field from a system of three coupled ordinary differential equations very neatly. Since the effective potentials for the $\Phi$ and ${\Phi }^{*}$ fields always exhibit an “U” type, the equations of motion have two trivial solutions which are the expectation values.

Consequently, the bounce is then a solution to the equation of motion for the $\sigma$ field when taking the $\Phi$ and ${\Phi }^{*}$ fields on their expectation values

$$\frac{d^2\sigma }{d{r}^2}+\frac{2}{r}\frac{d\sigma }{dr}=\frac{\partial {\Omega }_{\mathrm{MF}}}{\partial \sigma },$$

(32)

with the boundary conditions $\underset{r\to \infty }{lim}\sigma \left(r\right)={\sigma }_F$ and $\frac{d\sigma \left(r\right)}{dr}{|}_{r=0}=0$. Here, ${\sigma }_F$ is the false vacuum or the local minimum of the effective potential. It is worth noting that our analyses are different from the previous works in Refs. [64,65], where the first-order hadron quark phase transition is also considered in the three-flavor PQM model. In their works, the jumps of the order parameters for the strange quark and the Polyakov-loops fields were treated as active variables when the temperature is close to a first-order coexistence line. Based on the basic picture of the thin-wall approximation, they have parameterized the bounce solution through the interpolation of two values of the order parameter corresponded to two minima of the effective potential for four order parameters, i.g ${\sigma }_x$ for ud quarks, ${\sigma }_y$ for s quark, $\Phi$ and $\Phi *$. However, according to our above analyses, the jumps of the order parameters for the Polyakov-loops fields are fake signal to be taken as a first-order phase transition. These sudden leaps are induced by the disconnection of the chiral order parameter $\sigma$ and do not support by the effective potential. Moreover, our analyses are also in agreement with the results in Refs. [16,70,82,83], where the deconfinement phase transition and the phase transition for the s quark are crossover. Therefore, we would like to treat the Polyakov-loops fields as a background fields rather than the active variables during the chiral phase transition.

5. Results and Discussion

In this section, a bounce is to be gained by numerically solving the equation of motion in Eq.(32) with the proper boundary conditions, $\sigma \to {\sigma }_F$ as $r\to \infty$ and ${\sigma }^{\prime }\left(0\right)=0$. In this case, a three-dimensional problem has been simplified to an one-dimensional problem when we constrain the order parameter space in the $\sigma$ field direction, thus the numerical package [81] and the discussions presented in previous studies in the quark meson model [63] can be applied straightforwardly.

For a typical first-order phase transition, the effective potential ${\Omega }_{\mathrm{MF}}$ in the $\sigma$ field direction displays three distinct extrema. One is the potential barrier, two of them are the local minima representing the quark and hadron phases, respectively. For convenience, the quark phase is denoted as ${\sigma }_l$, whereas the hadron phase is defined as ${\sigma }_h$. As $T=T_{\chi }^c$, these two minima are degenerate, however, with the increase of the temperature further, the quark phase at $\sigma ={\sigma }_l$ becomes the absolute minimum of the effective potential and chiral symmetry is approximately restored for $T>T_{\chi }^c$. So that the quarks lose most of their constituent masses and become almost massless in this phase. Simultaneously, the other minimum at $\sigma ={\sigma }_h$ is the local minimum, and it is taken as a false vacuum. On the other side, as the temperature is to drop down from its critical temperature $T_{\chi }^c$, the hadron phase at $\sigma ={\sigma }_h$ is the absolute minimum of the effective potential and should be considered as the true vacuum, whereas the quark phase becomes a metastable state and is treated as a false vacuum. Since the false vacuum and true vacuum will get flipped when the temperature goes across to its critical value $T_{\chi }^c$, ${\sigma }_F$ in the boundary conditions should be set as ${\sigma }_F={\sigma }_l$ for $T>T_{\chi }^c$ and ${\sigma }_F={\sigma }_h$ for $T<T_{\chi }^c$. Under these boundary conditions, in what follows, we will divide our studies into two categories: a weak first-order phase transition and a strong first-order phase transition in order to present an intuitional description of the false vacuum decay at finite temperature and density.

5.1. A Weak First-Order Phase Transition

In the region of the weak first-order phase transition when the chemical potential is less than ${\mu }_c\simeq 309$ MeV, as the temperature is lowered from its critical value at $T_{\chi }^c$, the local minimum at $\sigma ={\sigma }_l$ approaches the intervening maximum, then two extrema meet and form an inflection point at the spinodal temperature $T_{\mathrm{sp}}\equiv T_{\mathrm{c}1}$. Thus, there only exists one minimum in the effective potential for $T<T_{\mathrm{c}1}$. Similar results can also be applied to the case in another direction of the temperature. When the temperature is larger than the critical one, $T_{\chi }^c$, with the increase of the temperature, the local minimum at $\sigma ={\sigma }_h$ will move up to the intervening maximum and eventually these two extrema will merge together at another spinodal temperature $T_{\mathrm{sp}}\equiv T_{\mathrm{c}2}$, so that only the global minimum remains for $T>T_{c2}$. Correspondingly, when the temperature is among two spinodal temperatures $T_{\mathrm{c}1}$ and $T_{\mathrm{c}2}$, there has a false vacuum, therefore there could exist a bounce solution to connecting the false vacuum and the true vacuum, otherwise, we have only a trivial solution for the equation of motion, i.e. the expectation value of the order parameter.

For $\mu =306$ MeV, the exact numerical solutions of the equation of motion (32) are plotted in the left panel of Figure 5 as $T=40$, 41, 44, 46, 47 and 48 MeV when the temperature is below the coexistence line, $T<T_{\chi }^c\simeq 48.9$ MeV. Here, the specific boundary conditions are selected as $\sigma \to {\sigma }_l$ as $r\to \infty$ and $d\sigma \left(0\right)/dr=0$. Form this figure, as the temperature is very close to the critical value at $T=T_{\chi }^c$, the bounce solutions show an obvious “core" structure with $\sigma \simeq {\sigma }_h$ inside the bubble which is separated by a relatively thin wall from the outside false vacuum at $\sigma \simeq {\sigma }_l$. While the temperature approaches to another limit for the existence of the bounce solution, $T_{c1}\simeq 28$ MeV, the bubble profiles usually become a “coreless" structure due to the fact that the radii of the bubbles are comparable to the thickness of the walls, thereby the thin-wall approximation can not be applied. On the other side of the coexistence line as $T>T_c$. the bubble profiles are obtained by solving the equation of motion (32) with the boundary conditions $\sigma \to {\sigma }_h$ as $r\to \infty$ and $d\sigma \left(0\right)/dr=0$. The results are plotted in the right panel of Figure 5 when taking the temperatures at $T=49.5$, 50, 51 and 52 MeV for $\mu =306$ MeV. For the temperature between critical value at $T_{\chi }^c\simeq 48.9$ MeV and the up spinodal temperature at $T_{c2}\simeq 52.1$ MeV, the bubble profiles exhibit an uniform structure with the true vacuum inside separated by a relatively thin wall from the outside false vacuum, only if the temperature is very close to its critical one at $T=T_{\chi }^c$, the center of the bubble is to deviate from its true vacuum value at $\sigma ={\sigma }_l$ largely since the thickness of the bubble wall is of the same order as the radius.

To give a more specific description of the size of the bubble, the typical radius of the bounce can be roughly estimate by the maximal value of the first derivative of the $\sigma \left(r\right)$ field, which gives $R_c\equiv \left|{\sigma }^{\prime }\left(r\right)\right|$. It is worth pointing out that this kind of definition of $R_c$ is rather arbitrary and we should not confuse with the critical radius ${\mathcal{R}}_c$ defined in the thin-wall approximation. In the framework of the thin-wall approximation, when the radius is much larger than the thickness of the bubble wall, characterized by $\xi \simeq 1/m_{\sigma }$, the frictional force $2{\sigma }^{\prime }/r$ in the equation of motion (32) can be neglected and the three-dimensional Euclidean action $S_3$ is well approximated by the expression

$$S_3=4\pi r^2\Sigma -\frac{4}{3}\pi r^3\epsilon ,$$

(33)

where $\epsilon =\Omega \left({\sigma }_F\right)-\Omega \left({\sigma }_T\right)$ is the difference between the values of the effective potential at the false vacuum and at the true vacuum. From this equation (33), there is an energy competition between the energy gain from the phase conversion of the bulk and the energy cost from the creation of an interface of the bubble. If the bubble is small, the energy cost is higher than the energy gain, the bubble will shrink and disappear. On the contrary, a very large bubble usually represents a larger bulk energy gain, which is higher than the surface energy cost. Therefore, the energy competition implies the existence of a critical size of the bubble, which is determined by minimization of the action $S_3$ with respect to r and is given by ${\mathcal{R}}_c=2\Sigma /\epsilon$.

In Figure 6, typical radii of the bounces as a function of temperature T are obtained for $\mu =306$ MeV in the two metastable regions. From this figure, the typical radius goes to zero at the spinodal line for both cases when the temperature is above or below the coexistence line. The reason is because the false vacuum is to become unstable and the phase conversion occurs via the spinodal decomposition process rather than the bubble nucleation, hence there is no bounce solution anymore as $T\to T_{\mathrm{sp}}$. This is an obvious feature of the weak first-order phase transition and it seems this character does not appear in the thin-wall approximation, where the typical radius does not vanish but only becomes very small [59,84]. However, when the temperature closes in the critical line at $T=T_{\chi }^c$, the exact numerical calculation and the thin-wall approximation show a similar behavior: the radius of the bounce rise up sharply and become divergent as $T\to T_{\chi }^c$. This is due to the fact that the radius of the bounce incline to expand very quickly as the temperature T approaches to $T_{\chi }^c$, when the radius is large enough, the thin-wall approximation is applied, we can end up with the consequence, $R_c\sim {\mathcal{R}}_c$. From the definition of the critical radius $\mathcal{R}$, noting that $\epsilon \to 0$ as $T\to T_{\chi }^c$, the radius of the bounce should also become divergent when $T=T_{\chi }^c$.

The surface tension plays a central role in determining the process of the bubble nucleation because it measures the amount of energy per unit of area in the construction of a interface between the false vacuum and the true vacuum. Once the bounce solutions have been obtained, the surface tension of the nucleation bubble can be calculated by using the definition in the equation (25). In Figure 7, we show the surface tension $\sigma$ as a function of the temperature T for $\mu =306$ MeV as the temperature is among the up and down spinodal lines. For both cases, as the temperature leads to its spinodal critical value $T_{\mathrm{sp}}$, the surface tension reduces to zero since there is no potential barrier and we can only have a trivial solution to the equation of motion at that moment. Oppositely, the surface tension $\Sigma$ shows a quite different behavior in the hadron phase by comparison with that of the quark phase. When the temperature is below the coexistence line, with the increase of the temperature, the surface tension $\Sigma$ starts to grow from zero dramatically and then reaches a maximal value at a certain temperature, after that it will inflect and bend downwards. Such a nontrivial behavior of the surface tension was also reported for a first-order hadron-quark phase transition in the quark meson model [63,84] and the Friedberg-Lee model [78], where the exact numerical method has been applied. But the non-monotonic feature of the surface tension did not appear in the studies based on the thin-wall approximation [64,66,77], so that the inflection point of the surface tension is usually taken as a hint for the breakdown of the thin-wall approximation. On the contrary, when $T>T_{\chi }^c$, the surface tension $\Sigma$ demonstrates a trivial character with a simply monotonic function. With the decrease of the temperature, it will grow up almost linearly from zero to a largest value as $T\to T_{\chi }^c$ given that there is a biggest potential barrier and a smallest free energy difference between two vacua at this point.

To estimate the nucleation rate for homogeneous nucleation theory, the largest contribution to the result is the action $S_3$ evaluated at a saddle-point because it has the exponential form in the equation (26). While the prefactor in the nucleation rate that comes from fluctuations around the saddle-point solution is of negligible numerical contribution in comparison with the term $e^{-S_3/T}$ and it can then be estimated based on dimensional analysis. Especially, when the saddle point action $S_3/T$ is much larger than the unit one. In order to display the saddle-point action due to the appearance of the bounce, the $S_3/T$ exponent as a function of temperature T for $\mu =306$ MeV is plotted in Figure 8 when the temperature is between the up and down spinodal critical temperatures. From this figure, the action will begin at the zero value as $T=T_{\mathrm{sp}}$, then it will rise up very quickly and become divergent as $T\to T_{\chi }^c$. Therefore, in the most area of the two metastable regions between the spinodal lines, the action $S_3/T$ is larger than the unit one and the decay of the false vacuum is still exponentially suppressed by the saddle-point action. Unless the temperature is very close to the spinodal critical value $T_{\mathrm{sp}}$, the action $S_3/T$ is to go acrose below to the unit one. This indicates that the hadron phase can survive in the region above the critical temperature line until the temperature is close to the up spinodal line or the quark phase can also live in the region below the critical temperature down to the another spinodal line. Thus, we believe it is better to use the spinodal line to separate the phase boundary of the quark and hadron phase rather than the critical coexistence line for the first-order phase transition.

5.2. A Strong First-Order Phase Transition

As discussed in the above section, when the chemical potential $\mu$ is larger than ${\mu }_c\simeq 309$ MeV, there always exists a potential barrier even when the temperature goes to zero for $T<T_{\chi }^c\simeq 28.8$ MeV. Normally, we denote a phase conversion induced by an effective potential with a zero-temperature potential barrier as a strong first-order phase transition. Taking $\mu =310$ MeV as a typical example of a strong first-order hadron quark phase transition, we have numerically solved the equation of motion (32) with the proper boundary conditions and the results have been depicted in Figure 9. From the picture, as the temperature is close to the critical temperature $T_{\chi }^c$, the bounce solution usually exhibits a “core" structure with the true vacuum inside and the false vacuum outside and the two vacua has been separated by a thin wall. This implicates that the thin-wall approximation is applicable for the temperature nearby the critical temperature. On the contrary, when the system is far away from the critical coexistence line, the bounce solution will lose its core structure due to the fact that the “frictional’ term in the equation of motion (32) becomes more and more important and can not be discarded as the decrease of the radius r. Moreover, since there is no spinodal line for the strong first-order phase transition, we can inevitably have the bounce solution which is the continues saddle-point field profile connected two vacua when $T<T_{\chi }^c$. Whereas on the other side of the critical coexistence line as $T>T_{\chi }^c$, with the increase of temperature, the potential barrier begin to decline and finally disappear at $T_{\mathrm{sp}}\simeq 41.4$ MeV. Hence there has only one minimum in the effective potential, thereby we have no bounce solution anymore as $T\ge T_{\mathrm{sp}}$.

In order to reveal the size of the bounce evolving with the temperature for a strong first-order phase transition, we have drawn the radius $R_c$ of the bounce as a function of the temperature T both in the hadron phase and in the quark phase as depicted in Figure 10. In a traditional quark phase as $T>T_{\chi }^c$, while the temperature reduces from the spinodal line to a critical value, the radius of the bounce will increase sharply and become divergent at $T=T_{\chi }^c$ for the reason that two vacua get degenerate and the energy difference $\epsilon$ is zero at that moment. On the other hand, when the system is located at the hadron phase, with the ascent of the temperature, the radius of the bounce climbs very slowly and steadily until the temperature is near a specific temperature somewhere at $T\simeq 20$ MeV. After that temperature, it will raise up dramatically and diverge at $T=T_{\chi }^c$. This is an apparent difference between a strong first-order phase transition and a weak first-order phase transition in the hadron quark phase conversion.

Similar to the case of a weak first-order phase transition, the surface tension in the hadron phase also presents a non-monotonic behavior as shown in Figure 11. With the increase of the temperature, it will grow up slightly and arrive at a maximum at a certain temperature. After that it starts decreasing rapidly to a minimum $\Sigma \left(T\right)\sim 3$$\mathrm{MeV}/{\mathrm{fm}}^2$ at $T=T_{\chi }^c$. As mentioned in the previous discussion, the turning point of the surface in general suggests a limit for the application of the thin-wall approximation. However, from the Figure 11, there are two visible differences between a strong and a weak first-order phase transition. One is that the surface tension $\Sigma \left(T\right)$ as a function of the temperature does not goes to zero with the reduction of the temperature by reason of the existence of the bounce for a strong first-order phase transition as $T<T_{\chi }^c$. Meanwhile, according to Figure 10, since the radius of the bounce is almost a constant when the temperature is less than 18 MeV, the surface tension carries an analogous property, it does not change too much with the reduction of the temperature when the temperature is smaller than 18 MeV. Furthermore, the other obvious difference is the value of the surface tension at $\mu =310$ MeV is much larger than that of a weak first-order phase transition at $\mu =306$ MeV. As $T=T_{\chi }^c$, the former is about 3 $\mathrm{MeV}/{\mathrm{fm}}^2$, whereas the latter is only about half of that around $1.6$$\mathrm{MeV}/{\mathrm{fm}}^2$. Therefore, we can conclude that the surface tension will increase with the enhancement of the chemical potential due to the fact that a large chemical potential usually hints a large potential barrier.

In a quark phase as $T>T_{\chi }^c$, the surface tension as a function of temperature T has been demonstrated in the right panel of Figure 11. It is shown that the surface tension $\Sigma \left(T\right)$ monotonically goes down from a maximum $\sim 3$$\mathrm{MeV}/{\mathrm{fm}}^2$ to zero when the temperature is ascent up to the spinodal line at $T=T_{\mathrm{sp}}$, where there is no potential barrier so that there is no bounce solution any longer. This trivial behavior of the surface tension is one of the obvious character of a weak first-order phase transition.

When fixing the chemical potential at $\mu =310$ MeV, the resulting plot of the $S_3/T$ action as a function of the temperature T due to the appearance of the bounce are shown in Figure 12 for both cases of the hadron phase and quark phase. From the left panel of Figure 12, for a strong first-order phase transition, the $S_3/T$ action reduces firstly to a minimal point at $T\sim 18$ MeV, then it will rise up quickly and intend to diverge again as $T\to T_{\chi }^c$. The non-monotonic behavior of the $S_3/T$ action is one of the special character of a strong first-order phase transition, which has been also reported in recent studies on a strong cosmological first-order phase transition [85] and a strong hadron quark phase transition [63,78]. Moreover, since the $S_3/T$ action evaluated on the bounce persistently satisfies the requirement of $S_3/T\gg 1$ in the hadron phase, the nucleation rate of the hadron state inside the homogeneous quark phase has been exponentially suppressed in the whole region of the conversional hadron phase, even though the temperature is below the critical coexistence line at $T=T_{\chi }^c$. This is a meaningful result because if the system is cooling down from a very high energy, where the initial state is well prepared as a free quark, the system is more like to remain in a quark state rather than a hadron state in spite of the fact that the temperature is below the critical coexistence line at $T=T_{\chi }^c$ for $\mu >{\mu }_c$. In other words, for an area of a strong first-order phase transition in phase diagram, it is better to treat the system as a quark matter instead of a hadron matter when $T\le T_{\chi }^c$.

When the temperature is above the critical coexistence line at $T=T_{\chi }^c$, the hadron quark phase conversion has been changed to a weak first-order phase transition by reason of the disappearance of the potential barrier as $T\to T_{\mathrm{sp}}$. In the right panel of Figure 12, the $S_3/T$ action as a function of the temperature T due to the appearance of the bounce is illustrated when the chemical potential is fixed at $\mu =310$ MeV. Now it is a weak first-order phase conversion, the evolution of the $S_3/T$ action with the temperature exhibits a similar result as the prior cases at $\mu =306$ MeV. The $S_3/T$ action tends to become divergent when the temperature approaches to the critical one at $T=T_{\chi }^c$, however, in the opposite direction, as the temperature inclines to the spinodal line, it will descent to the zero dramatically. Unless the temperature is very close to the spinodal line, the action fulfils the a critical condition such that $S_3/T>1$. As a result the hadron quark phase transition dose not take place exactly at the critical temperature but up to a certain temperature nearby the spinodal line.

6. Summary

In this work, we have investigated the dynamics of a first-order phase transition via homogeneous thermal nucleation within the Polyakov quark-meson model at finite temperature and density. In the framework of the mean field approximation, the temperature and density-dependent effective potential in the presence of a fermionic vacuum term has been calculated and the phase diagram together with the two metastable regions has been obtained. It is found that at low density the chiral phase transition is a crossover, however it will turn to a first-order phase transition in the high-density region accompanied by a critical end point (CEP). By using a geometric method of the effective potential, we can precisely locate the position of the CEP at $(T_{\mathrm{E}},{\mu }_{\mathrm{E}})\simeq (301.4\mathrm{MeV},62.1\mathrm{MeV})$. Furthermore, in the region of a first-order phase transition, besides the coexistence line at $T=T_{\chi }^c$, two additional spinodal lines which are usually ignored in the standard phase diagram of QCD have been also presented explicitly. From phase diagram in Figure 4, we can find the distance of these two spinodal lines is to decline with the temperature descent. In the end, the spinodal lines will joint together with the coexistence line and terminate by the same point at the CEP, where the order of phase transition is of second order.

In order to study the instanton tunneling between two vacua, we need to search all minima of the effective potential. Aside from a standard way through the calculation of the Hessian matrix, we have provided an alternative method by constrained the effective potential to the mesonic field direction and the Polyakov-loop field direction in the parameter space for the sake of simplification. The advantage of the our geometric method is that it supplies an intuitive and vivid way to find two minima nearby which can actually have the tunneling solution. Moreover, we can change the complicated N-dimensional problem into a typical one-dimensional problem and this will make us much easy to solve the equation of motion to get the bounce solution both in an exactly numerical method and an analytically thin-wall approximation. Of course, the method in the present study have two main shortages. For the first one, the procedure of finding the minima is usually trouble and sophisticated, especially when the potential has many order parameter variables and many minima, since we have to explore all minima of the effective potential one by one. For the second one, there is an opportunity to miss some minima around a specific point of the effective potential, but the absent minima can be recovered again if we run over all the established minima. Fortunately, although there are three order parameters in the PQM model, the effective potential of the model has merely two minima and the instanton tunneling is more likely to happen in the $\sigma$ field direction rather than the Polyakov-loop field directions. Accordingly, the system of the three coupled equation of motions for the bounce has been dissociated, and what should be solved is an ordinary differential equation of the $\sigma$ field in Eq.(32). Thus the problem has been simplified completely and the method developed in a single field can be applied directly.

With the appropriate boundary conditions, the exact bubble profiles were numerically calculated for a first-order phase transition. For convenience, our discussions have been separated into a weak first-order hadron quark conversion and a strong one. For both cases, when the temperature is close to the critical coexistence line at $T=T_{\chi }^c$, the bubble profile exhibits a core structure. The kernel of the bubble stays in a true vacuum at $\sigma ={\sigma }_T$, while it is separated from a homogeneous false vacuum at $\sigma ={\sigma }_F$ by a thin wall. On the other side, when the temperature leaves far away from the critical one, the bubble profile tends to display a coreless structure since the thickness of the bubble wall has the same order of its radius. To give a more detailed information about the size of the bounce, the typical radius has been roughly estimated by the maximal value of the first derivative of the bubble profile. For a weak first-order phase transition, as $T\to T_{\chi }^c$, the radius $R_c$ goes to the infinite, whereas it will drop down to zero very quickly as the temperature approaches to the spinodal line. For a strong first-order phase transition as $\mu >{\mu }_c$ and $T<T_{\chi }^c$, the radius becomes divergent as $T\to T_{\chi }^c$, but it remains almost a constant with the descent of the temperature. This is one of special features of a strong first-order phase transition.

The surface tension plays an important role during the procedure of bubble nucleation because it is an amount of the energy cost per unit area in production of the interface between two vacua. In the scenario of a weak first-order phase transition, the surface tension reaches to a relative larger value at $T=T_{\chi }^c$, then it will continuously reduce to zero as $T\to T_{\mathrm{sp}}$. on the other hand, in the context of a strong first-order phase transition, the surface tension is likely to keep its constant value as $T\to 0$ since there is no spinodal temperature any longer. Similarly, the $S_3/T$ action evaluated on the bounce solution shows a common characteristic for a weak first-order phase transition. It is divergent as the critical temperature closes in, while in the opposite direction, it will fall down dramatically to zero as $T\to T_{\mathrm{sp}}$. Considered the exponential dependence of the nucleation rate on the $S_3/T$ action, the decay of the false vacuum is still exponentially suppressed when $S_3/T>1$. Hence the false vacuum could survive and exist as a metastable state for a relatively long time so long as the system lies between the up and low spinodal lines. Only when the temperature is very close to the spinodal critical line, the $S_3/T$ will go through the unity 1. As a result, the phase boundary for a weak first-order phase transition should be resized accordingly with the results in the present work. Moreover, for a strong first-order phase transition, the situation will become more worse, since there is no down spinodal line, the nucleation rate of the true vacuum is always exponentially suppressed as $S_3/T\gg 1$. In other words, a “conventional” hadron matter below the critical coexistence line for $\mu >{\mu }_c$ should be potentially treated as a quark matter either, such a conclusion together with the low values of the surface tension will favor a more complicated structure of strong interaction matter in high density.