Equation of state of a quark-meson mixture in the improved-PNJL model at finite chemical potential

TL;DR

This paper investigates the equation of state of QCD using an improved three-flavor PNJL model that includes unquenched quark effects and mesonic contributions, comparing results with lattice QCD and exploring finite chemical potential.

Contribution

It introduces an improved PNJL model incorporating unquenched quarks and mesonic effects, providing detailed thermodynamic calculations and phase diagram analysis.

Findings

Pressure and entropy density agree with lattice QCD at zero chemical potential

Phase diagram and meson masses are explored at finite chemical potential

Thermodynamical pressure shows consistency with theoretical expectations

Abstract

We study the equation of state of QCD using an improved version of the three-flavor Polyakov-Nambu-Jona-Lasinio model beyond the mean-field approximation. It incorporates the effects of unquenched quarks into the Polyakov-loop effective potential, as well as mesonic contributions to the grand-canonical potential. We study in full detail the calculation of the thermodynamical potential in this approach and compare the resulting pressure and entropy density with the most-recent lattice-QCD calculations at zero baryochemical potential. Finally, we present some exploratory results at finite chemical potential which include the phase diagram of the model, the quark and meson masses, and finally, the thermodynamical pressure.

Figures

Tables

Table 1. Table 1: Parameters for 𝒰 YM ​ ( T ; Φ , Φ ¯ ) subscript 𝒰 YM 𝑇 Φ ¯ Φ {\cal U}_{\rm{YM}}(T;\Phi,\bar{\Phi}) in Eq. ( 5 ), taken from Ref.  Ratti:2005jh .

Parameter	$a_0$	$a_1$	$a_2$	$a_3$	$b_3$	$b_4$	$T_0$
Value	6.75	-1.95	2.625	-7.44	0.75	7.5	270 MeV ($N_f=0$)
							190 MeV ($N_f=2+1$)

Equations

$${\mathcal{L}}_{PNJL}$$

$${\text{tr}}_f({\tau }^a{\tau }^b)=2{\delta }^{ab},$$

$$L(\mathbf{x})=\mathcal{P}\exp \left(i{\int }_0^{1/T}d\tau A_4(\tau ,\mathbf{x})\right),\Phi (\mathbf{x})=\frac{1}{N_c}{\mathrm{tr}}_{c}L(\mathbf{x}),$$

$${\frac{\partial \mathcal{U}(T;\Phi ,\bar{\Phi })}{\partial \Phi }|}_{\Phi =⟨\Phi ⟩(T),\bar{\Phi }=⟨\bar{\Phi }⟩(T)}=0,{\frac{\partial \mathcal{U}(T;\Phi ,\bar{\Phi })}{\partial \bar{\Phi }}|}_{\Phi =⟨\Phi ⟩(T),\bar{\Phi }=⟨\bar{\Phi }⟩(T)}=0.$$

$$\frac{{\mathcal{U}}_{\mathrm{YM}}(T;\Phi ,\bar{\Phi })}{T^4}=-\frac{b_2(T)}{2}\bar{\Phi }\Phi -\frac{b_3}{6}\left({\Phi }^3+{\bar{\Phi }}^3\right)+\frac{b_4}{4}{\left(\bar{\Phi }\Phi \right)}^2,$$

$$b_2(T)=a_0+a_1\frac{T_0}{T}+a_2{\left(\frac{T_0}{T}\right)}^2+a_3{\left(\frac{T_0}{T}\right)}^3,$$

$${\Omega }_{\mathrm{PNJL}}(T,{\mu }_i;⟨{\bar{\psi }}_i{\psi }_i⟩,\Phi ,\bar{\Phi })={\Omega }_q(T,{\mu }_i;⟨{\bar{\psi }}_i{\psi }_i⟩,\Phi ,\bar{\Phi })+{\mathcal{U}}_{YM}(T;\Phi ,\bar{\Phi }),$$

$${\Omega }_q={\Omega }_q^{(-1)}+{\Omega }_q^{(0)}+\cdots ,$$

$${\Omega }_q^{(-1)}(T,{\mu }_i;⟨{\bar{\psi }}_i{\psi }_i⟩,\Phi ,\bar{\Phi })$$

$$m_i=m_{i0}-4G⟨{\bar{\psi }}_i{\psi }_i⟩+2H⟨{\bar{\psi }}_j{\psi }_j⟩⟨{\bar{\psi }}_k{\psi }_k⟩,j,k\ne i;j\ne k.$$

$$L=\text{diag }(e^{i\varphi },e^{i{\varphi }^{\prime }},e^{-i(\varphi +{\varphi }^{\prime })}).$$

$${\text{ tr}}_c\log \left(1+L{e}^{-(E_i-{\mu }_i)/T}\right)=\sum _{j=1}^3\log \left(1+e^{-(E_i-{\mu }_i)/T}{L}_{jj}\right)=\log \prod _{j=1}^3\left(1+e^{-(E_i-{\mu }_i)/T}{L}_{jj}\right).$$

$$\Phi =\frac{1}{3}{\text{ tr}}_{c}L=\frac{1}{3}\left(e^{i\varphi }+e^{i{\varphi }^{\prime }}+e^{-i(\varphi +{\varphi }^{\prime })}\right),$$

$${\text{ tr}}_c\log \left(1+L{e}^{-(E_i-{\mu }_i)/T}\right)$$

$${\text{ tr}}_c\log \left(1+L^{†}{e}^{-(E_i+{\mu }_i)/T}\right)$$

$$\frac{\partial {\Omega }_{PNJL}}{\partial ⟨{\bar{\psi }}_i{\psi }_i⟩}=0,\frac{\partial {\Omega }_{PNJL}}{\partial \Phi }=0,\frac{\partial {\Omega }_{PNJL}}{\partial \bar{\Phi }}=0,$$

$$⟨⟨{\bar{\psi }}_i{\psi }_i⟩⟩=-6\int \frac{d^{3}k}{(2\pi {)}^3}\frac{m_i}{E_i}\left[1-f_{\Phi }^{+}(E_i-{\mu }_i)-f_{\Phi }^{-}(E_i+{\mu }_i)\right],$$

$$f_{\Phi }^{+}(E_i-{\mu }_i)$$

$$f_{\Phi }^{-}(E_i+{\mu }_i)$$

$$\frac{{\mathcal{U}}_{glue}}{T^4}(T;\Phi ,\bar{\Phi })=\frac{{\mathcal{U}}_{YM}}{T_{YM}^4}(T_{YM};\Phi ,\bar{\Phi }),$$

$$\frac{T_{YM}-T_{YM}^{cr}}{T_{YM}^{cr}}=0.57\frac{T-T_{glue}^{cr}}{T_{glue}^{cr}},$$

$${\Delta }_{l,s}(T)=\frac{⟨⟨{\bar{\psi }}_q{\psi }_q⟩⟩(T)-\frac{m_{q0}}{m_{s0}}⟨⟨{\bar{\psi }}_s{\psi }_s⟩⟩(T)}{⟨⟨{\bar{\psi }}_q{\psi }_q⟩⟩(0)-\frac{m_{q0}}{m_{s0}}⟨⟨{\bar{\psi }}_s{\psi }_s⟩⟩(0)}.$$

$$P(T)=-[{\Omega }_{\mathrm{PNJL}}(T,{\mu }_q=0)-{\Omega }_{\mathrm{PNJL}}(T=0,{\mu }_q=0)],$$

$${\mathcal{K}}_M=G-\frac{H}{2}⟨⟨{\bar{\psi }}_s{\psi }_s⟩⟩.$$

$${\Omega }_q^{(0)}(T,{\mu }_i)=\frac{1}{2}{\int }_0^1\frac{d\lambda }{\lambda }\int \frac{d^{3}p}{(2\pi {)}^3}T\sum _n\text{ Tr }[S^{\lambda }(i{\nu }_n,\mathbf{p}){\Sigma }^{\lambda }(i{\nu }_n,\mathbf{p})],$$

$$S^{\lambda }(i{\nu }_n,\mathbf{p})=S_H(i{\nu }_n,\mathbf{p}).$$

$${\Sigma }^{\lambda }(i{\nu }_n,\mathbf{p})=\sum _{M}T\sum _m\int \frac{d^{3}q}{(2\pi {)}^3}\Omega S_H(i{\omega }_m,\mathbf{q})\bar{\Omega }{t}_M^{\lambda }(i{\nu }_n-i{\omega }_m,\mathbf{p}-\mathbf{q}),$$

$$t_M^{\lambda }(i{\nu }_n-i{\omega }_m,\mathbf{p}-\mathbf{q})=\frac{2\lambda {\mathcal{K}}_M}{1-2\lambda {\mathcal{K}}_M\Pi (i{\nu }_n-i{\omega }_m,\mathbf{p}-\mathbf{q})},$$

$${\Omega }_q^{(0)}(T,{\mu }_i)=\sum _{M\in J^{\pi }=\{0^{+},0^{-}\}}{\Omega }_M^{(0)}(T,{\mu }_M({\mu }_i)),$$

$${\Omega }_M^{(0)}(T,{\mu }_M({\mu }_i))$$

$${\Omega }_M^{(0)}(T,{\mu }_M)$$

Full text

1. I Introduction

2. II Thermodynamics at $\mu_{q}=0$: the mean field approximation

3. II.1 PNJL Lagrangian and grand-canonical potential at ${\cal O}(N_{c})$

4. II.2 Effect of quarks into the Polyakov-loop effective potential

5. III Beyond mean field: Mesonic fluctuations at $\mu_{q}=0$

6. III.1 Grand-canonical potential at ${\cal O}((1/N_{c})^{0})$

7. III.2 Pressure of a quark—Polyakov loop—meson mixture at $\mu_{q}=0$

8. IV Equation of state at finite density

9. IV.1 Phase diagram and boundaries: Maxwell construction

10. IV.2 Properties at finite $\mu_{q}$

11. IV.2.1 Approximate determination of meson masses

12. IV.2.2 Pressure of the PNJL model at ${\cal O}((1/N_{c})^{0})$ and $\mu_{q}\neq 0$

13. V Conclusions and extensions

14. VI Acknowledgements

15. A In-medium quark-antiquark propagator

16. B Scattering phase shift

I Introduction

One of the most widely used low-energy realizations of quantum chromodynamics (QCD) for the quark degrees of freedom is the Nambu–Jona-Lasinio (NJL) model Nambu:1961tp ; Vogl:1991qt ; Klevansky:1992qe ; Hatsuda:1994pi ; Alkofer:1995mv ; Buballa:2003qv . In this effective model one only considers dynamical quarks which exchange small momenta in their interaction, whereas gluons are integrated out. This model shares the global symmetries of QCD, and it is able to describe the phenomenon of spontaneous symmetry breaking and its high-temperature/density restoration. In particular, the local version of the model (where quarks interact locally in space through a four-point interaction) has the property of simplicity and transparency, as opposed to other more complicated approximations to QCD. In addition, thanks to its straightforward extension to finite temperatures and densities, it allows for the study of the QCD phase diagram.

Static properties of gluons can be implemented in the model by a coupling to a Polyakov-loop effective potential Fukushima:2003fw ; Megias:2004hj ; Ratti:2005jh ; Hansen:2006ee ; Fukushima:2008wg . Such static properties, like the equation of state, are nicely described when the effective potential is conveniently parametrized by using lattice-QCD calculations at finite temperature. However, as opposed to QCD on a lattice, the Polyakov-Nambu-Jona-Lasinio (PNJL) model allows for a direct analysis of the phase diagram in the whole temperature ($T$)-baryochemical potential ($\mu_{B}$) plane. At low baryonic densities one can compare to lattice-QCD calculations that use a Taylor expansion around $\mu_{B}=0$ Ratti:2005jh .

Old tensions in the equation of state at $\mu_{B}=0$, computed by two different lattice-QCD collaborations, have been solved by newer results and a nice agreement has been achieved. It has also been confirmed that for physical $N_{f}=2+1$ quark masses the equation of state at vanishing chemical potential presents a crossover phase transition Aoki:2006we , with a “critical temperature” (the temperature at the inflection point of the chiral quark condensate) around $T_{c}=155$ MeV Borsanyi:2010bp ; Bazavov:2011nk . In addition, the temperature-dependence of the QCD pressure (as well as other thermodynamical quantities) turns out to be very similar for both groups within a wide range of temperatures Borsanyi:2013bia ; Bazavov:2014pvz .

Given this robust result for the equation of state on the lattice, one would like to check whether the thermodynamics of the PNJL model is consistent with it. Unfortunately, when one compares the lattice-QCD results with the classical mean-field calculations of the pressure in the PNJL model Ratti:2005jh one realizes several deficiencies. At low temperature the pressure is well below the lattice-QCD results, and above $T_{c}$ it increases too fast towards the ideal gas expectations, as opposed to what is seen in lattice results (see Ref. Marty:2013ita for an illustration).

In this paper we revisit the thermodynamics of the PNJL model and show that these two discrepancies with respect to the lattice-QCD results have two natural origins, and can be cured using already-established approaches. First, at low temperature, where hadronic degrees of freedom dominate, the pressure is mainly produced by the lightest mesonic states. The calculation of the thermodynamical potential beyond mean field includes meson-like fluctuations Quack:1993ie ; Hufner:1994ma and brings the hadronic pressure much closer to the lattice data. Second, at moderate temperatures the interaction of quarks in the gluonic environment modifies the Polyakov loop effective potential. Introducing an effective potential which includes the quark back reaction to the gluons Haas:2013qwp we show that quarks are able to control the pressure of the glue sector (in a dual way in which the Polyakov loop regulates the quark pressure below the critical temperature) and this interaction lowers the pressure around $T/T_{c}\simeq 1.5-2$ to lattice-QCD values.

Therefore, we implement for the first time the quark back reaction into the Polyakov-loop effective potential in a NJL-like model (our results at $\mu_{B}=0$ have been presented in a short note in Torres-Rincon:2016ahl ), determining the grand canonical potential—and the phase boundary—beyond mean field. To the best of our knowledge, this is the first time that these techniques have been jointly implemented into a NJL model at finite density.

In Sec. II we introduce the PNJL model with scalar and pseudoscalar interactions, and compute the grand-canonical potential at mean-field level. We introduce the “glue” potential that accounts for the back reaction of quarks onto the Polyakov loop effective potential. We will show that at high temperatures the pressure is consistent with lattice-QCD data. In Sec. III we add mesonic fluctuations to the thermodynamic potential which enter as collective modes of the propagating quarks. This improves the description of the low-temperature pressure, bringing it closer to the lattice-QCD results. After achieving a reasonable agreement with lattice-QCD calculations at $\mu_{B}=0$, we will compute results at finite chemical potential in Sec. IV. We will determine the phase diagram of the model as a function of $T$ and $\mu_{B}$, show some examples of the behavior of quark and meson masses as functions of the baryochemical potential, and finally, compute the pressure of the model at finite $\mu_{B}$. In Sec. V we present our conclusions and describe the possible extensions of this work to reach a finer description of the QCD thermodynamics.

II Thermodynamics at $\mu_{q}=0$: the mean field approximation

In this section we will review the effective Lagrangian and the basic properties of the PNJL model at finite temperature and finite chemical potential Klevansky:1992qe ; Buballa:2003qv ; Fukushima:2003fw ; Fukushima:2008wg ; Costa:2010zw . We will detail the grand-canonical potential $\Omega(T,\mu_{q})$ at mean-field level and explain the role of the effective potential for the Polyakov loop. Then, we will argue for the need to account for modifications of the Polyakov loop effective potential due to the presence of quarks, and adopt the findings of Ref. Haas:2013qwp to improve the description of the gluonic sector of the theory. Finally, we will compute the equation of state of the improved PNJL model at mean-field level and $\mu_{q}=0$.

II.1 PNJL Lagrangian and grand-canonical potential at ${\cal O}(N_{c})$

We consider the Lagrangian of the PNJL model Fukushima:2003fw ; Megias:2004hj ; Ratti:2005jh ; Hansen:2006ee ; Fukushima:2008wg ; Torres-Rincon:2015rma with (color neutral) pseudoscalar and scalar interactions (neglecting the vector and axial-vector vertices for simplicity),

[TABLE]

where the flavor indices $i,j,k,l=1,2,3$ and $\tau^{a}$ ($a=1,...,8$) being the $N_{f}=3$ flavor generators with normalization

[TABLE]

with $\textrm{tr}_{f}$ denoting the trace in flavor space.

In the Lagrangian (1) the bare quark masses are represented by $m_{0i}$ and their chemical potential by $\mu_{i}$. The covariant derivative in the Polyakov gauge reads $D^{\mu}=\partial^{\mu}-i\delta^{\mu 0}A^{0}$, with $A^{0}=-iA_{4}$ being the temporal component of the gluon field in Euclidean space (we denote $A^{\mu}=g_{s}A_{a}^{\mu}T_{a}$). The coupling constant for the scalar and pseudoscalar interaction $G$ is taken as a free parameter (fixed e.g. by the pion mass in vacuum).

The third term of Eq. (1) is the so-called ’t Hooft Lagrangian. It mimics the effect of the axial $U(1)$ anomaly, accounting for the physical splitting between the $\eta$ and the $\eta^{\prime}$ mesons. $H$ is a coupling constant (fixed by the value of $m_{\eta^{\prime}}-m_{\eta}$) and $\mathbb{I}$ is the identity matrix in Dirac space.

Finally, ${\cal U}(T;\Phi,\bar{\Phi})$ is the so-called Polyakov-loop effective potential used to account for static gluonic contributions to the pressure. The Polyakov line and the Polyakov loop are respectively defined as

[TABLE]

where ${\cal P}$ is the path-integral ordering operator, and the trace ${\rm tr}_{c}$ is taken in the color space.

In the limit of infinitely heavy quarks —i.e. for a Yang-Mills (YM) theory— the expectation value (minimum of the effective potential) of the Polyakov loop (EVPL) $\langle\Phi\rangle(T)$ is a true order parameter of the deconfinement transition. That means that the EVPL is strictly zero (nonzero) in the confined (deconfined) phase (cf. right panel of Fig. 1). In the standard approach, the EVPL is taken to be a real function, but independent of its complex conjugated Dumitru:2005ng .

Following Ratti:2005jh we take an homogeneous Polyakov loop field $\Phi({\bf x})=\Phi=$const., and calculate the expectation values $\langle\Phi\rangle(T),\langle\bar{\Phi}\rangle(T)$ that minimize the effective potential ${\cal U}(T;\Phi,\bar{\Phi})$ at a given temperature

[TABLE]

At this point it is convenient to remind that for a given temperature, although the effective potential is a function of the Polyakov loop (and its conjugate), the only physical information is carried by the EVPL, where the effective potential takes its minimum. The reasoning is in the spirit of the Landau-Ginzburg functional for an homogeneous order parameter. In addition, the value of the effective potential at the EVPL gives the physical pressure of the YM theory, which is a function of the temperature only, $P(T)=-{\cal U}(T;\langle\Phi\rangle(T),\langle\bar{\Phi}\rangle(T))$. The EVPL and the pressure of the $SU(3)$ YM system can be computed by lattice QCD as functions of the temperature Borsanyi:2012ve . Therefore, lattice-QCD results can be used to constrain the minima of ${\cal U}$ as a function of the temperature; and eventually, to parametrize the effective potential as done in Ref. Ratti:2005jh .

However it is important to notice that the knowledge of $\langle\Phi\rangle(T)$ and $P(T)$ can only help to fix the effective potential at its minimum. The complete dependence of the effective potential on $\Phi$ (and ${\bar{\Phi}}$) is not determined by these functions. For this reason, several ansätze for ${\cal U}$ have been introduced in the literature (compare for example the different parametrizations considered in Ref. Dutra:2013lya ), with the requirement of being compatible with the breaking pattern of the ${\cal Z}_{3}$ center symmetry of QCD.

Following our previous work Torres-Rincon:2015rma we will take the so-called “polynomial parametrization” for the Yang-Mills potential, introduced in Ref. Ratti:2005jh ,

[TABLE]

where

[TABLE]

and the parameters $a_{0},a_{1},a_{2},a_{3},b_{3},b_{4}$ and $T_{0}$ are given in Ref. Ratti:2005jh (see Table 1). The effective potential is plotted in the left panel of Fig. 1 for several temperatures around $T_{0}$.

The EVPL is obtained by performing the minimization (4) of the (Yang-Mills) effective potential in (5). The phase transition is of first order with a transition temperature of $T_{0}=270$ MeV (see right panel of Fig. 1).

In the presence of dynamical quarks, the effective potential ${\cal U}_{\rm{YM}}(T;\Phi,\bar{\Phi})$ complements the Lagrangian of the quark sector in Eq. (1). The associated grand-canonical potential $\Omega_{\rm{PNJL}}$ accounts for the thermodynamical properties of a gas of quarks and “static” gluons,

[TABLE]

where $\langle\bar{\psi}_{i}\psi_{i}\rangle$ is the scalar quark bilinear for each flavor $i$. It should not be confused with its expectation value (to be defined later), for which we will use a double bracket. The quark bilinear is merely a variable of the grand-canonical potential, not yet fixed to its physical value. We use a single bracket notation to distinguish it from the $\bar{\psi}_{i}(x)\psi_{i}(x)$ operator.

The quark potential will be eventually expanded in inverse powers of $N_{c}$,

[TABLE]

where the exponent $(a)$ denotes the term at order ${\cal O}(1/N_{c}^{a})$.

For the time being we consider the mean-field result Hufner:1994ma , which is equivalent to the leading-order term ${\cal O}(N_{c})$ and reads Hansen:2006ee ; Torres-Rincon:2015rma

[TABLE]

where $E_{i}=\sqrt{k^{2}+m_{i}^{2}}$. Notice that $G\sim{\cal O}(1/N_{c}),H\sim{\cal O}(1/N_{c}^{2})$ and $\langle\bar{\psi}_{i}\psi_{i}\rangle\sim{\cal O}(N_{c})$, making all terms in Eq. (9) to be ${\cal O}(N_{c})$.

In Eq. (9) it is important to notice that due to interactions, the constituent quark mass $m_{i}$ depends on the quark bilinears $\langle{\bar{\psi}}_{j}\psi_{j}\rangle$ via a “gap equation”. In a diagrammatic way, this equation is sketched in Fig. 2, where we have added to the 4-point and 6-point interactions a finite spatial range (in dashed line) to detail the topology of the diagrams.

In the mean-field approximation the gap equation reads Buballa:2003qv ; Hansen:2006ee

[TABLE]

In addition, the color traces ($\textrm{tr}_{c}$) of Eq. (9) can be performed for $N_{c}=3$ using a mean-field approximation for the Polyakov line $L=\exp\left(iA_{4}/T\right)\in SU(3)$ Fukushima:2003fw . In the Polyakov gauge, it can be represented by Fukushima:2003fw

[TABLE]

Using this diagonal representation one finds

[TABLE]

Expanding the last product, and writing the result in terms of the Polyakov loop (and its conjugate):

[TABLE]

it is straightforward to get

[TABLE]

Similarly to Eq. (4), the grand-canonical potential (7) should be minimized with respect to $\langle\bar{\psi}_{i}\psi_{i}\rangle,\Phi$ and $\bar{\Phi}$,

[TABLE]

to get the expectation value of the quark bilinear—the quark condensate—, the EVPL and its conjugate, denoted respectively as $\langle\langle\bar{\psi}_{i}\psi_{i}\rangle\rangle$, $\langle\Phi\rangle$ and $\langle\bar{\Phi}\rangle$.

The first equation of Eq. (16) provides the expression for the quark condensate (recall that the derivative with respect to the quark bilinear also affects the quark constituent mass)

[TABLE]

which can be written as $\langle\langle\bar{\psi}_{i}\psi_{i}\rangle\rangle=3m_{i}A_{i}/4\pi^{2}$ in terms of the function $A_{q}$ defined in Eq. (54) of App. A.

The Fermi-Dirac distribution functions are modified by the Polyakov loop,

[TABLE]

For $N_{f}=3$ the gap equation (10), together with Eq. (17), and the last two equations in (16) form a system of five coupled equations. We solve this system numerically to obtain the value of $\langle\langle\bar{\psi}_{i}\psi_{i}\rangle\rangle$ (or equivalently $m_{i}$), $\langle\Phi\rangle$ and $\langle\bar{\Phi}\rangle$. In the isospin limit, two gap equations are degenerate giving $m_{u}=m_{d}$ (and $\langle\langle{\bar{\psi}}_{u}\psi_{u}\rangle\rangle=\langle\langle\bar{\psi}_{d}\psi_{d}\rangle\rangle$). In addition, at vanishing chemical potential one has $\Phi=\bar{\Phi}$, which is evident from our equations.

When minimizing (7) using the last two equations in (16), it is clear that the presence of quarks will affect the determination of the EVPL (simply because $\Omega_{q}$ also depends on $\Phi,\bar{\Phi}$). The quark sector is affected as well by the EVPL, cf. Eq. (17). The most drastic effect is the transformation of the deconfinement transition from a first-order transition in the pure gauge case to a crossover in the unquenched theory, with an approximate transition temperature at $T_{0}=270$ MeV.

As argued in Ref. Schaefer:2007pw on the basis of renormalization group arguments, the presence of dynamical quarks should also be accompanied by a reduction of the $T_{0}$ parameter of the gluon potential (6). This induces a decrease of the transition temperature. Following our last work Torres-Rincon:2015rma we assume a value of $T_{0}=190$ MeV for $N_{f}=2+1$ flavors. The EVPL as a function of the temperature can be seen in the left panel of Fig. 3. In particular note that $\Phi$ becomes larger than one at high temperatures, as already noticed in Ref. Schaefer:2007pw . When the Polyakov-loop effective potential is coupled to the quark sector, the restriction $\Phi,{\bar{\Phi}}\leq 1$ in the Yang-Mills case does not necessarily hold anymore.

It is important to note that the change of $T_{0}$ does not modify the functional form of the Polyakov-loop effective potential. Therefore, at a given temperature, the EVPL does not coincide anymore with the minimum of the effective potential, simply because $\langle\Phi\rangle(T)$ is obtained by the minimization of the grand canonical potential, and not of ${\cal U}_{YM}$. We illustrate this in the right panel of Fig. 3 where we plot ${\cal U}_{\rm{YM}}/T^{4}$ (with $T_{0}=190$ MeV) for three different temperatures $T=170,190,210$ MeV, together with the location (marked with arrows) of the EVPL for each temperature. One can see that the EVPL might be considerably away from the minimum of the effective potential.

For this reason, to obtain the “gluonic” contribution to the pressure $P_{\rm{gluon}}=-{\cal U}_{YM}(T;\langle\Phi\rangle(T),\langle\bar{\Phi}\rangle(T))$, one should evaluate the effective potential in a region not modified by the presence of dynamical quarks, i.e. away from its minimum (points marked with arrows in Fig. 3). In view of this situation, it is legitimate to ask for a possible effect of quarks on the Polyakov-loop effective potential out of its minima. Any modification of ${\cal U}_{\textrm{YM}}$ due to quarks would result critical for the correct computation of the gluonic pressure.

II.2 Effect of quarks into the Polyakov-loop effective potential

To pursue such a study, we introduce the quark back reaction onto ${\cal U}$, not limited to their effect at the minimum of the potential. This has been considered in Ref. Haas:2013qwp in the context of the functional renormalization group (FRG) equations applied to QCD.

In Ref. Haas:2013qwp the pure gauge effective potential ${\cal U}_{YM}$ is compared to the “glue” potential ${\cal U}_{glue}$, where quark-antiquark excitations are allowed in the gluon propagator. The technical details of the computation within the framework of the FRG are explained in detail in that reference. The authors found remarkable differences between the two effective potentials. However, the study shows that the glue potential turns out to be very approximately related to the original ${\cal U}_{YM}$ but with a rescaling of the temperature. The glue potential is related to the YM one as Haas:2013qwp :

[TABLE]

where the relation between $T_{YM}$ and the temperature $T$ is given by

[TABLE]

with $T_{YM}^{cr}$ playing the role of the deconfinement temperature in the YM case (and fixed to $T_{YM}^{cr}=270$ MeV), and $T_{glue}^{cr}$ the transition temperature in the unquenched case. The numerical coefficient 0.57 is the outcome from the comparison of the two effective potentials.

The value $T^{cr}_{glue}=210$ MeV is used in Ref. Haas:2013qwp in the context of the quark-meson model, fixed to reproduce lattice-QCD data. In our case, we simply fix this value to $T_{glue}^{cr}=190$ MeV, reminiscent of the $T_{0}$ parameter that we used in Ref. Torres-Rincon:2015rma . As we will check later, this value will eventually provide a reasonable result when comparing to lattice-QCD data. In addition, this choice lies within the uncertainties quoted in Haas:2013qwp .

Before presenting our results we will comment the choice of parameters and the conventions used. For the NJL sector we will use the same parameters of our Ref. Torres-Rincon:2015rma , chosen to reproduce experimental results of meson and baryons masses, as well as the pion decay constant and the value of the quark condensate in vacuum. As we do not treat baryons here, nor vector mesons, we only have 5 parameters, viz. $m_{q0}$, $m_{s0}$, $\Lambda$, $G$ and $H$. They are fixed in vacuum with the help of $m_{\pi}$, $m_{K}$, $\langle\langle\bar{\psi}_{q}\psi_{q}\rangle\rangle$, $f_{\pi}$ and $m_{\eta}-m_{\eta^{\prime}}$. At finite temperature we parametrize the YM Polyakov loop effective potential as described in Sec. II.1 with the parameters given in Table 1.

As a low-energy realization of the QCD dynamics, the PNJL Lagrangian is not renormalizable. We use an ultraviolet (UV) regulator $\Lambda$ that cuts off the momentum integrals. This is consistent with the idea that the model is only valid at low momenta (below $\Lambda$). Nevertheless, at finite temperature the high-momentum modes are automatically suppressed by the distribution function inside the integrals. For this reason, we regularize our momentum integrals in such a way that $\Lambda$ is set to infinity in all UV-finite integrals, but we keep it finite in the divergent integrals (like vacuum terms). With this prescription one can reproduce the Stefan-Boltzmann limit of the thermodynamic functions in the asymptotic $T\rightarrow\infty$ limit Zhuang:1994dw ; Costa:2009ae ; Costa:2010zw .

Our first observation is that the use of ${\cal U}_{\rm{glue}}$ does not appreciably change the mass of quarks as compared to the YM potential. We compare the two sets of results in the left panel of Fig. 4. As being directly related to $m_{i}$, the quark condensates do not present sensible differences. However, in the right panel of Fig. 4 we see some changes in the EVPL, $\langle\Phi\rangle$, calculated by Eq. (16). The glue effective potential provides a reduction of the EVPL at high temperatures to values smaller than 1 for the temperatures considered here. A similar reduction in the Polyakov loop is seen in the context of the quark-meson model in Ref. Haas:2013qwp .

The chiral transition temperature $T_{c}$ can be defined as the inflection point of the quark condensate. Using the results in Borsanyi:2010bp (with a $T_{c}=157\pm 6$ MeV) we compute the renormalized chiral condensate $\Delta_{l,s}$ defined as

[TABLE]

In Fig. 5 we show the close correspondence of $\Delta_{l,s}$ with lattice results when plotted as a function of $T/T_{c}$. In the PNJL model the transition temperature is typically shifted to higher temperatures. In this case we get $T_{c,l}=209$ MeV for the light sector, and $T_{c,s}=208$ MeV for the strange one. For this reason, to compare with lattice-QCD results we opt to plot this condensate as a function of $T/T_{c}$.

Finally, we present the result for the mean-field pressure at $\mu_{q}=0$. It is computed as,

[TABLE]

where we subtract the vacuum pressure to have $P(T=0)=0$. We remind that at mean-field level the grand-canonical potential of the model is $\Omega_{\rm{PNJL}}=\Omega_{q}^{(-1)}+{\cal U}$.

In Fig. 6 we plot the pressure of the PNJL model at mean-field level using ${\cal U}={\cal U}_{YM}$ with $T_{0}=190$ MeV (dashed line), and using ${\cal U}={\cal U}_{glue}$ (solid line). The YM potential provides a pressure which increases fast to the gluon SB limit, whereas the glue potential limits this increase and gives a lower pressure at high temperatures. We also include the lattice-QCD results from Refs. Borsanyi:2013bia (Set 1) and Bazavov:2014pvz (Set 2). We make note that the fast overshooting of the lattice-QCD results by the YM potential cannot be cured by modifying the parameters of the PNJL model alone, as we have checked with several parametrizations and regularization procedures in the literature. It is natural that the PNJL model cannot explain the QCD pressure at some intermediate range of temperatures because it lacks from attractive gluon interactions, which should contribute to the pressure above $T\simeq 2-2.5T_{c}$ Ratti:2005jh (lattice-QCD data at higher temperatures can be described by the hard-thermal-loop perturbation theory Haque:2014rua ). At asymptotically large temperatures the Polyakov-loop effective potential approaches the SB limit by construction.

We end this section by pointing out the interesting dual role of the coupling between the NJL Lagrangian and the glue effective potential. On the one hand, it is known that the Polyakov loop potential helps to suppress non-physical quark pressure at low temperature, contributing to the effective confinement of the system. On the other hand, we have seen that quarks also play a similar role at high temperature, helping to reduce the gluonic contribution to the pressure. We find this observation particularly interesting. It should be mentioned that the glue potential has been recently used in the context of the $N_{f}=2+1$ quark-meson model in Refs. Kovacs:2016juc ; Stiele:2016cfs providing an analogous gluon pressure reduction, and a similar comparison to lattice-QCD data.

In the next sections we will fix ${\cal U}={\cal U}_{glue}$ given by Eq. (20), and analyze the low-temperature regime by computing a systematic correction to the pressure beyond the mean-field level.

III Beyond mean field: Mesonic fluctuations at $\mu_{q}=0$

The modification of the Polyakov-effective potential in the PNJL model provides a reasonable approximation to the pressure computed in lattice QCD at moderate temperatures around $T\sim(1.5-2)T_{c}$, as can be seen in Fig. 6. However, at low temperatures the pressure of the system is underestimated by the mean-field calculation.

Hadronic states should also contribute to the pressure of the system, especially at $T<T_{c}$, where they represent the main degrees of freedom (quarks and gluon remain confined at low temperatures and their pressure should be suppressed). They should be responsible for the pressure deficit below $T_{c}$ as compared to lattice data in Fig. 6. On the other hand, the addition of the hadronic states cannot spoil the current agreement at high temperature. One expects that hadrons deconfine around $T_{c}$ and become part of the partonic medium.

To include the effect of the lightest hadrons on the thermodynamic potential, we need to go beyond the mean-field approximation, by accounting for next-order diagrams in the large-$N_{c}$ expansion.

We will illustrate the large-$N_{c}$ expansion of the thermodynamic potential by the diagrammatic picture adapted from Refs. Quack:1993ie ; Hufner:1994ma . This expansion is also sketched in Ref. Zhuang:1994dw . In practice—due to technical simplifications abrikosov ; fetter —the $1/N_{c}$-corrections are computed, not directly on the grand-canonical potential, but at the level of the quark self-energy, as we will detail later. We believe, however, that a more pedagogical understanding of the large-$N_{c}$ expansion is provided in terms of the grand-canonical potential.

III.1 Grand-canonical potential at ${\cal O}((1/N_{c})^{0})$

At claimed in Sec. II.1, the mean-field potential corresponds to the leading order in a $1/N_{c}$ expansion, $\Omega_{q}^{(-1)}(T,\mu_{i})$ Quack:1993ie ; Hufner:1994ma . This term contains the $O(N_{c})$ noninteracting, Hartree and ’t Hooft diagrams Zhuang:1994dw ; Hufner:1994ma , which we plot in Fig. 7. We have given a spatial range to the contact interaction to distinguish the different topology between the ${O}(N_{c})$ Hartree (second diagram in Fig. 7) and the ${O}((1/N_{c})^{0})$ Fock (first diagram in Fig. 8) interactions.

This set of diagrams, computed at mean-field level in the imaginary-time formalism, and including the Polyakov loop corrections, provide the terms contained in Eq. (9): The expression of the noninteracting diagram can be found in any textbook Kapusta:2006pm with the standard Fermi-Dirac factors modified by the coupling to the Polyakov loop, Eqs. (18,19). At mean-field level, the Hartree diagram contains two powers of the quark condensate and a factor $G$ coming from the interaction. The ’t Hooft diagram contains one quark condensate of each flavor and a factor $H$, cf. Eq. (1).

After incorporating the Hartree diagram into the thermodynamic potential, one is tempted to include the Fock diagram as well (first diagram in Fig. 8). In the NJL limit—when the dashed lines are contracted into a point—the Hartree and the Fock diagrams are indistinguishable. Nevertheless, in the large-$N_{c}$ limit the Fock diagram is explicitly $O((1/N_{c})^{0})$, and it is only one of the infinite number of terms that contribute at this order. These are depicted in Fig. 8 and denoted as “ring diagrams”. A similar resummation of ring diagrams appears in the random-phase approximation fetter of a dense electron gas where electron-hole loops are summed. When the interaction can be taken as local like in the NJL model, it becomes equivalent to the present case.

One has to consider all these diagrams together to be consistent with $N_{c}$ counting.

We now explain the meaning of the filled squares at the vertices of the diagrams in Fig. 8. Using mean field we can simplify the original ’t Hooft interaction to an effective 4-point vertex which is, in addition, combined with the 4-point interaction $G$. The resulting quartic coupling ${\cal K}_{M}$ depends on the particular flavor channel (denoted by $M$), i.e. depends on the flavor of the incoming and outgoing quarks. For example, if all the interacting quarks are light, then

[TABLE]

The combinations for the other flavor channels are given in Ref. Klevansky:1992qe ; Rehberg:1995kh ; Torres-Rincon:2015rma .

Notice that this combination do not spoil $N_{c}$-counting. The effective vertices ${\cal K}_{M}(T,\mu_{i})$ of the diagrams depicted in Fig. 8 are denoted with a filled square.

A direct summation of the ring diagram is difficult due to the combinatorial factors in the perturbative expansion of the potential abrikosov ; fetter . In practice, it is more convenient to use the so-called “coupling-constant integration” abrikosov ; fetter ; Klevansky:1992qe , which is equivalent to the diagrammatic summation. One computes the ${\cal O}((1/N_{c})^{0})$ potential as,

[TABLE]

where the fermion self-energy $\Sigma^{\lambda}(i\nu_{n},{\bf p})$ and the quark propagator $S^{\lambda}(i\nu_{n},{\bf p})$ admit an expansion in powers of $1/N_{c}$ and the coupling constant is rescaled by ${\cal K}_{M}\rightarrow\lambda{\cal K}_{M}$. Finally the trace Tr is taken in color, flavor and spin spaces.

The full quark propagator $S^{\lambda}$ contains self energy insertions to all orders in $N_{c}$. In addition, the quark self-energy $\Sigma^{\lambda}$ is also computed order by order using the full quark propagator. Instead of solving such self-consistent problem, we attach ourselves to a strict $N_{c}$ expansion keeping the thermodynamic potential at ${\cal O}((1/N_{c})^{0})$. We follow the same procedure proposed in Refs. Klevansky:1992qe ; Hufner:1994ma .

The quark propagator is kept at leading order ${\cal O}(1)$, i.e. the Hartree propagator shown in Fig. 2,

[TABLE]

Then, the fermion self-energy $\Sigma^{\lambda}(i\nu_{n},{\bf p})$ is written in terms of the scattering amplitude, $t_{M}$.

[TABLE]

with $\Omega=\mathbb{I}_{c}\otimes\tau^{a}\otimes\{\mathbb{I},i\gamma_{5}\}$ being factors in color, flavor and spin spaces Torres-Rincon:2015rma depending on the meson channel $M$. In this work the meson channel can be scalar or pseudoscalar, $J^{\pi}=\{0^{+},0^{-}\}$.

The scattering amplitude is strictly ${\cal O}(1/N_{c})$,

[TABLE]

where $\Pi$ is the quark-antiquark correlation function given in Eq. (50) which is ${\cal O}(N_{c})$ when the Hartree propagators are used.

Then, the fermion self-energy in Eq. (27) is also ${\cal O}(1/N_{c})$, and the thermodynamical potential is, as desired, ${\cal O}((1/N_{c})^{0})$ (the color trace introduce one more power of $N_{c}$). Had we kept higher powers of $1/N_{c}$ in the quark propagator of Eq. (26), then we would have generated ${\cal O}(1/N_{c})$ terms in the thermodynamic potential, which are higher-order to the terms considered in this work.

How mesons appear in this calculation? It is convenient to mention that the same resummation as performed in Eq. (28) for the scattering amplitude has been done for the solution of the Bethe-Salpeter equation for a $q{\bar{q}}$ scattering Klevansky:1992qe ; Rehberg:1995kh ; Buballa:2003qv ; Torres-Rincon:2015rma . These works show that mesons are generated as poles of the $q\bar{q}$-scattering amplitude, with masses and decay widths given by the location of the poles in the complex-energy plane. Therefore, the presence of mesons is reflected in the increase of the scattering amplitude (28), subsequently affecting the thermodynamic potential in Eq. (25).

Combining Eqs.(25), (26), (27) and (28) we express the thermodynamic potential as

[TABLE]

where $M$ denotes the contribution of each mesonic spin-flavor channel. The individual contribution of a mesonic channel is

[TABLE]

where $g_{M}$ is the spin-isospin degeneracy of the meson state. Notice that the meson chemical potential is a function of the quark chemical potential, $\mu_{M}=\mu_{i}-\mu_{i^{\prime}}$, with $i$ and $i^{\prime}$ denoting the valence quark and antiquark of the meson.

We perform a change of variables $i\omega_{m}\rightarrow i\nu_{n}-i\omega_{m},{\bf q}={\bf p}-{\bf q}$ and get

[TABLE]

where in the second line we recognize the quark-antiquark correlation function $\Pi(i\omega_{m},{\bf q})$ (cf. Eq. (50)). Therefore,

[TABLE]

The integration over $\lambda$ can be performed exactly Zhuang:1994dw ; Hufner:1994ma ; Blaschke:2013zaa ; Wergieluk:2012gd

[TABLE]

Now we perform the remaining Matsubara summation in Eq. (33) by the complex integration around Matsubara frequencies Kapusta:2006pm

[TABLE]

We deform the contour of the integral Hufner:1994ma to a contour around the real axis (where the two-quark propagator has analytical cuts, cf. App. A). We can distinguish two integrations

[TABLE]

We split the integrals in Eq. (35) ($\int_{-\infty}^{+\infty}=\int_{-\infty}^{0}+\int_{0}^{+\infty}$) and perform the variable shift $\omega\rightarrow\omega-\mu_{M}$,

[TABLE]

Now we perform the change of variable $\omega\rightarrow-\omega$ in the second integral, and use the property of the polarization function

[TABLE]

where we have used the fact that the polarization function $\Pi(\omega,{\bf p})$ is a function of the energy always through the combination $\omega+\mu_{M}$ (cf. Eqs. (53, A)). We obtain

[TABLE]

The integration limits can be again restored to $(0,+\infty)$ by adding and subtracting the same quantity

[TABLE]

Then, we use the relation

[TABLE]

to get

[TABLE]

Finally, one has

[TABLE]

It is interesting to relate this equation to the so-called Beth-Uhlenbeck approach, exploiting the use of the scattering phase shift. We provide a discussion of the scattering phase shift in App. B, and refer the reader to Refs. Yamazaki:2012ux ; Yamazaki:2013yua ; Blaschke:2013zaa for more information.

Using Eq. (61) we can express $\Omega^{(0)}_{M}$ in Eq. (42) as an integration over the scattering phase shift,

[TABLE]

where we have neglected the first term within the brackets of Eq. (42) (this is usually called “no-sea” approximation).

An equivalent formula is sometimes used after integrating by parts

[TABLE]

which resembles the Beth-Uhlenbeck approach to the thermodynamics of an interacting gas Blaschke:2013zaa ; Beth:1936 ; Beth:1937zz .

A simplifying assumption is that the argument of the phase shifts is approximately Lorentz invariant (we check the validity of this claim in Appendix B). Then, we can introduce the Mandelstam variable $s=\omega^{2}-p^{2}$, and approximate

[TABLE]

Introducing $s$ in Eq. (43), we obtain the meson part of the grand-canonical potential,

[TABLE]

where the phase shift is computed using Eq. (62).

The limits of integration must be chosen in accordance with the UV cutoff of our theory. The momentum integral is limited to $(0,\Lambda)$. However the UV cutoff is eventually removed in all convergent integrals, for consistency with the mean-field calculation. The integral in $s$ has a lower limit at [math] (the quark-antiquark propagator is evaluated at real values of $\sqrt{s}$). The upper limit is taken at $(\Lambda_{1}+\Lambda_{2})^{2}=(\sqrt{\Lambda^{2}+m_{1}^{2}}+\sqrt{\Lambda^{2}+m_{2}^{2}})^{2}$, where the unitary cut ends (cf. Appendix B for details of terminology). For low temperatures the phase shift present several discontinuities as a function of the energy (see Fig. 23). To avoid systematic numerical uncertainties we use the known limits of the Landau and unitary cuts to perform the integration over $s$ in a support where the function is nonzero (cf. Fig 21).

III.2 Pressure of a quark—Polyakov loop—meson mixture at $\mu_{q}=0$

In this section we finally present our results for the pressure of the PNJL model computed beyond the mean-field approximation at ${\cal O}((1/N_{c})^{0})$. For the time being we fix the quark chemical potential to zero and plot the pressure as a function of the temperature,

[TABLE]

where the grand-canonical potential taken at ${\cal O}((1/N_{c})^{0})$,

[TABLE]

with $\Omega^{(-1)},\Omega^{(0)}_{M}$ and ${\cal U}_{glue}$ given by Eqs. (9), (46), and (20), respectively. Notice that all parameters have been fixed at mean-field level, so that the mesonic contribution is parameter free for every temperature.

We show in Fig. 9 the total pressure of the system when we add more and more mesonic states on top of the mean-field (MF) result. For the sake of clarity, in this plot we have only included pseudoscalar and scalar mesons with a vacuum mass below 1 GeV. As expected, the dominant contribution at low temperatures (up to $T\simeq 100$ MeV) comes from the pions. Then, more massive mesons (in particular the four kaons) start to affect the thermodynamics. Notice that heavier states make only tiny contribution to the pressure as they are Boltzmann suppressed due to their higher masses. At high temperature, the mean-field pressure dominates, and the mesons contribute as resonant states as they melt above their Mott temperature. This temperature is defined as the one at which the meson pole acquires an imaginary part, thus generating a nonzero decay width. Above this temperature the meson has nonzero probability to decay into a quark-antiquark pair.

In Fig. 10 we compare the total pressure and entropy density of the system with the most recent lattice-QCD calculations. Set 1 corresponds to Ref. Borsanyi:2013bia (Wuppertal-Budapest group) and Set 2 to Ref. Bazavov:2014pvz (HotQCD collaboration). In the left panel we include the same mesonic states as accounted in Fig. 9, i.e., those with a vacuum mass below 1 GeV. We already see a large improvement with respect to the pressure at mean field (cf. Fig. 6). The mesonic contribution enhances the low-temperature pressure towards the lattice-QCD data. At intermediate temperatures there is still a small deviation which is naturally associated to the lack of more hadrons. This area can be systematically improved adding more and more states. For example, if we include all scalar and pseudoscalar states ($\pi,K,\bar{K},\eta,\eta^{\prime},a_{0},K_{0}^{*},\bar{K}_{0}^{*},f_{0},f_{0}^{\prime}$) we observe in the right panel of Fig. 10 that the pressure is closer to the lattice-QCD data. Of course, other missing mesons (and baryons) also contribute in this region.

Finally, we observe that the mesonic contribution overestimates the pressure of lattice-QCD, which already was correctly reproduced at high temperatures (cf. Fig 6).

Mesons contribute as scattering states to the pressure even at high temperatures. Although the contribution of each degree of freedom is tiny (of the order of 1 % of the total pressure), the addition of many states makes a non negligible contribution at large temperatures. This mesonic pressure is nonphysical at such temperatures, and one has to go beyond the present approach to cure this problem. One expects that a reduction of the mesonic pressure at high temperatures, the total pressure flattens, thus reducing the entropy density as well.

A way to remove the non-physical meson pressure at high temperatures might be to incorporate the effect of mesonic fluctuations into the quark propagator. This idea has been suggested for the NJL model in Blaschke:2015bxa , and applied recently to a simplified version of the PNJL model in Blaschke:2016hzu ; Blaschke:2016fdh . In this approach the meson phase-shift in Eq. (43) is corrected by an extra factor $-\frac{1}{2}\sin^{2}\delta(\omega,{\bf p};T,\mu_{M})$, which drastically reduces the mesonic pressure at high temperatures (see Fig. 3 in Ref. Blaschke:2016fdh ).

In Blaschke:2016hzu ; Blaschke:2016fdh a (negative) perturbative correction to the pressure from QCD Kapusta:2006pm is included to lower the pressure at high temperature. We believe that this correction is not appropriate within this scheme. It takes into account contributions arising from the diagrams in Fig. 11 at ${\cal O}(\alpha_{s})$ in QCD for the momentum modes between $\Lambda$ and $\infty$ (not included in the original NJL model). Even assuming that perturbative QCD can be applied to temperatures $T\lesssim 300$ MeV, these diagrams have been already accounted for in the calculation. The diagram involving a quark propagator becomes part of the Hartree and Fock diagrams. Once removing the UV cutoff in our integrations (to reproduce the Stefan-Boltzmann asymptotic limit) the momentum modes beyond $\Lambda$ have been already integrated. The diagrams involving gluons and ghosts belong to the Yang-Mills sector, whose pressure is already reproduced by construction thanks to the parametrization of the Yang-Mills effective potential ${\cal U}_{YM}$. Therefore, these diagrams already participate in the effective description, where gluons with momentum between 0 and $\infty$ are integrated out. As a result, the inclusion of these diagrams would involve a double counting of UV degrees of freedom.

From the results of this manuscript we claim that these corrections are not only not appropriate, but also not needed for the suppression of the mean-field pressure. The back reaction of quarks into the gluon effective potential already takes care of this reduction.

IV Equation of state at finite density

We present now some results at finite baryochemical potential by turning on the quark chemical potential. We set the strange chemical potential to zero, $\mu_{s}=0$, whereas the light quark chemical potential is taken in the isospin limit, $\mu_{q}=\mu_{u}=\mu_{d}$. In terms of baryon, charge and strangeness chemical potentials we use $\mu_{S}=\mu_{Q}=0$, $\mu_{B}=3\mu_{q}$.

The extension to a small chemical potential —where the transition is still a crossover— is straightforward. However, at high densities one faces the problem of multiple solutions in the gap equation (10) for a given temperature. This is the signature of a coexistence of different phases and the presence of a first-order phase transition. It is important to determine the precise location of the phase transition by means of the Maxwell construction. In the next section, we detail how to compute the phase coexistence region and the phase boundary of the model.

IV.1 Phase diagram and boundaries: Maxwell construction

We now describe the procedure to determine the phase boundary of the model. To simplify the discussion, we will work with the mean-field thermodynamic potential at ${\cal O}(N_{c})$ in this section. However, we will also incorporate the next-to-leading order terms in our final results.

For a temperature of $T=1$ MeV, we find a region of the phase transition where several solutions for the quark mass exist for a unique value of the chemical potential. In the left panel of Fig. 12 we show the solutions for the quark mass within a range of $\mu_{q}$.

Between $\mu_{q}\simeq 380$ MeV and $\mu_{q}\simeq 490$ MeV, the gap equation has three solutions for the quark mass. At lower chemical potentials we only find the chirally-broken solution (with a quark mass of $m_{q}=480$ MeV), whereas at higher chemical potentials we have the chirally-symmetric phase (with a quark mass close to the bare one, $m_{q0}$). In between, the three solutions are depicted with different lines styles.

The phase transition is determined by an analogue of the Maxwell construction. In the presence of different solutions, the energetically favorable phase is the one with the minimum value of $\Omega(T,\mu_{i})$. In the right panel of Fig. 12 we plot the value of (minus) the thermodynamic potential without vacuum subtraction. At a given $\mu_{q}$ the solution with highest $-\Omega(T,\mu_{q})$ represents the stable phase. In this example the chirally-broken phase is the stable one up to $\mu_{q}=424$ MeV. At this point the quark mass, and hence the order parameter of the chiral transition, present a discontinuity. For larger $\mu_{q}$, the stable phase is the chirally-restored one (with small quark masses).

Repeating this procedure for other temperatures we map the first-order critical boundary on the phase diagram, up to the critical end-point in which the three solutions degenerate into a single one. This signals a second-order critical point, and the beginning of the crossover regime. The phase boundary is shown in the left panel of Fig. 13.

When expressed in terms of the baryochemical potential ($\mu_{B}=3\mu_{q}$) the critical point is given by the coordinates $(\mu_{B,c},T_{c})=(733,165)$ MeV. For the sake of comparison, the critical point using the ${\cal U}_{YM}$ with $T_{0}=190$ MeV provides $(\mu_{B,c},T_{c})=(735,169)$ MeV. We also show the area where the presence of a critical point is unlikely according to a recent lattice-QCD calculation Bazavov:2017dus . Notice that vector interactions among quarks—not considered in this work—can modify these values at mean field, and lower the critical temperature Fukushima:2008wg .

Quark masses are computed as functions of $T$ and $\mu_{q}$ using Eq. (10). We present them in the right panel of Fig. 13 for several quark chemical potentials. Notice the discontinuity of the masses for $\mu_{q}>244$ MeV, as a consequence of the first-order phase transition.

Finally in Fig. 14 we show the mean-field pressure of the system (23) as a function of $T$ and $\mu_{q}$. We recall that the glue potential ${\cal U}_{\rm{glue}}$ has no explicit dependence on $\mu_{q}$ in our scheme. It only depends on $\mu_{q}$ through the coupling of the quark field to the EVPL.

Notice that the pressure is always a monotonically increasing function of the temperature, and continuous at the phase transition by construction, but with a discontinuous first derivative for the case of a first-order phase transition (see lines with $\mu_{q}>244$ MeV).

IV.2 Properties at finite $\mu_{q}$

Following the scheme outlined for $\mu_{q}=0$ to include mesonic fluctuations in the thermodynamical potential, we will compute the PNJL model pressure at ${\cal O}((1/N_{c})^{0})$ at finite chemical potential. From the results at $\mu_{q}=0$ we can assume that the main contribution to the pressure comes from pions —lightest mesons, hence easy to create thermally— and from kaons/antikaons —more massive, and therefore Boltzmann suppressed, but with an isospin degeneracy of 4. For simplicity, we do not consider other mesons due to their expected small contribution.

IV.2.1 Approximate determination of meson masses

In order to understand the influence of mesons, we show the dependence of the masses of $\pi$, $K$ and $\bar{K}$ on the temperature and chemical potential. The dynamically generated mesons appear as poles of the $q\bar{q}$-scattering amplitude (28),

[TABLE]

in the appropriate flavor channel $M$ Rehberg:1995kh ; Torres-Rincon:2015rma . A standard determination of the meson mass (and possible decay width) comes from the pole position $1-2{\cal K}_{M}\Pi(z=M-i\Gamma/2,{\bf p}={\bf 0})=0$ in the complex-energy plane, when the scattering amplitude is analytically continued to the appropriate Riemann sheet. The polarization function $\Pi(z,{\bf p})$ is shown in Eq. (53).

If the pole is not very far from the real axis, one can simplify the search by neglecting the imaginary part of the solution in the argument of the $B_{0}$ function in Eq. (53). In this case, one does not need to extend the $B_{0}$ function to the second Riemann sheet Rehberg:1995kh . To illustrate the meson masses we follow this path, even if at large temperatures/chemical potentials it might not be suitable (the decay width increases fast after the Mott temperature, and the pole moves away from the real axis). In our calculation of the pressure we make no such approximation; it is only used in this section to extract the meson masses.

In the left panel of Fig. 15 we show the pion masses, and in the right panel the $K$ and $\bar{K}$ masses, for a given set of chemical potentials.

The kaon masses show a peculiar behavior due to the fact that isospin degeneracy is broken by the chemical potential. As we take $\mu_{s}=0$, the $K$ ($K^{+}$ and $K^{0}$ isodoublet) will carry a meson chemical potential $\mu_{K}=\mu_{q}$, whereas the ${\bar{K}}$ ($K^{-}$, $\bar{K}^{0}$ doublet) will have $\mu_{{\bar{K}}}=-\mu_{q}$. This difference produces non equal masses for $K$ and $\bar{K}$ at high temperatures and high baryochemical potentials. In the right panel of Fig. 15, we kept a small temperature range for clarity between the curves.

In Fig. 16 we summarize the Mott temperatures $T^{\textrm{Mott}}$ for the lightest mesons versus the quark chemical potential. It might be illustrative to known the temperature at which the different mesons acquire a decay width as a function of the chemical potential, and how this temperature decreases as long as $\mu_{q}$ increases.

It is interesting to see that the ${\bar{K}}$ states have a slightly larger Mott temperature that the $K$ mesons. This indicates that ${\bar{K}}$ is more stable in the plasma than the $K$, especially at high chemical potential.

IV.2.2 Pressure of the PNJL model at ${\cal O}((1/N_{c})^{0})$ and $\mu_{q}\neq 0$

Using Eqs. (47,48) we now present the pressure of the PNJL model for several values of the chemical potential including mesonic fluctuations. The pressure created by pions and kaons (those states with strangeness = 1) is plotted in Fig. 17, as a function of temperature for several chemical potentials. The Mott temperature is also indicated in the plot for each chemical potential.

We obtain a jump of the pressure when crossing the first-order phase transition ($\mu_{q}=300$ MeV). It is more pronounced in the kaon case. Although across a first-order transition the total pressure should be continuous, the partial pressure of a given meson is in general discontinuous.

In Fig. 18 we show the pressure of antikaons (those with strangeness = 1) in the left panel, and the total pressure of the system, i.e. a mixture of quarks, the dominant mesons ($\pi,K,\bar{K}$), and the effects of the Polyakov loop.

It is important to remind that in this section we have determined the location of the first-order phase transition by minimizing the grand-canonical potential at ${\cal O}((1/N_{c})^{0})$. Had we used the potential at ${\cal O}(N_{c})$, we would have obtained a continuous mean-field pressure (as in Fig. 14), but a non-physical discontinuous total pressure. In other words, the phase boundary of the model computed at mean field is corrected by the mesonic fluctuations. For example, at $\mu_{q}=300$ MeV the temperature at which the first-order transition takes place is modified by $\sim 5$ MeV.

Finally, we present a calculation of the pressure compared to recent lattice-QCD results at finite chemical potential Bazavov:2017dus . To be consistent with lattice-QCD data we have used $\mu_{s}=\mu_{q}=\mu_{B}/3$, so we have recomputed all our results using a nonzero strange quark chemical potential. The outcome is shown in Fig. 19.

The results follow the trend of our Figs. 6 and 10. The mean-field results (dashed lines) provide a good description of lattice-QCD data at high temperature thanks to the use of the glue effective potential, but with a poor description of pressure at low temperatures. Addition of meson fluctuations (pions, kaons and antikaons) considerably improves the low-temperature part, still not close to the data due to the absence of additional meson states.

V Conclusions and extensions

In this work we have revisited the thermodynamic properties of the $N_{f}=2+1$ PNJL model including ${\cal O}((1/N_{c})^{0})$ corrections to the grand-canonical potential to account for physical mesonic excitations in the system. In addition, we have taken into account the modification of the Yang-Mills Polyakov-loop effective potential by the presence of quarks.

At zero chemical potential we have obtained a pressure which is close to the most recent common lattice-QCD results from several groups in a fairly wide range of temperatures. Then, we have computed the phase diagram of the model, and explored the equation of state at finite chemical potential. The results provide thermodynamical information of regions in the phase diagram away from the $\mu_{q}=0$ axis. This is the relevant equation of state for the quark-gluon-plasma created in heavy-ion collisions with center-of-mass energies of a couple of GeV, like the low-energy collisions at the Beam Energy Scan at RHIC collider in Brookhaven (U.S.A.), the current NA61/SHINE experiment from the LHC accelerator at CERN at Geneva, the future CBM experiment at FAIR/GSI facility in Darmstadt (Germany) and NICA experiment at Dubna (Russia). For a final implementation of the PNJL model into numerical simulations Marty:2012vs ; Marty:2014zka ; Marty:2015iwa at these energies, one needs to extract masses and cross sections of the different species of particles (including quark, mesons, diquarks and baryons) at finite temperatures and densities. This is currently under progress.

Including the pseudoscalar and scalar mesons, this approach describes reasonably well the low-temperature side of the phase diagram, which is dominated by mesons. We can also describe the phase diagram at moderate temperatures due to the combined contribution of quarks and the Polyakov-loop effective potential. At higher temperatures it is evident that dynamical gluons cannot be neglected and the description if not suitable (although the Stefan-Boltzmann limit is asymptotically recovered). However, the PNJL model is valid for the typical temperatures (and densities) that can be reached presently in heavy-ion reactions (the temperature extracted from direct thermal photons in ALICE Pb+Pb collisions at $\sqrt{s}_{NN}=2.76$ TeV is $T\sim 300$ MeV Adam:2015lda ).

This work can be improved in several ways. A straightforward extension is to consider more mesonic states like vector and axial vector states, by enlarging the effective Lagrangian to include these interaction channels. Diquarks can be effectively introduced as described in Ref. Blaschke:2014zsa . Their contribution to the pressure is suppressed at low temperatures by the Polyakov loop in a similar way than for quarks. Baryons Blanquier:2011zz ; Torres-Rincon:2015rma should also be introduced, especially at finite chemical potential, where the contribution of nucleons can largely affect the thermodynamics. However, it is not evident how to perform this extension at the level of the thermodynamic potential. An alternative effective description has been used in Ref. Blaschke:2015zwa .

As discussed in the text we believe that $\alpha_{s}$ corrections are not justified to improve the mean-field sector Blaschke:2016fdh ; Blaschke:2016hzu . They are not even required for a good matching of lattice-QCD data, because the back reaction of quarks into the gluon effective potential Haas:2013qwp already corrects the mean-field pressure at moderate temperatures without the additional hypothesis of the validity of perturbative expansion.

An undesirable effect beyond mean field is the non negligible contribution of mesons to the pressure around $T\sim 300$ MeV. A promising path to suppress this contribution is sketched in Ref. Blaschke:2015bxa , where the mesonic correction of the quark propagator produces a negative correction to the mesonic pressure. However, this approach cannot be accommodated into a strict $N_{c}$ counting, like the one presented here. An alternative solution might be to perform a full self consistent approach, for example, based in a two-particle-irreducible potential like the one in Ref. Muller:2010am . Such a study involves a big technical effort and it is beyond the scope of this work.

One could also adjust the parameters of the Polyakov-loop effective potential after having included the mesonic contributions. This would allow us to suppress the total pressure around $T=300$ MeV, and reach lattice-QCD data. However, this would hide the problem of having an unphysical meson pressure at high temperatures. Moreover, if more and more mesons are included in the calculation, then a new parametrization must be fixed each time to reduce the pressure at high temperatures. In our approach the parameter set is fixed at mean field once for all.

VI Acknowledgements

We thank D. Blaschke for useful discussions and K. Yamazaki for email exchange. JMTR also acknowledges D. Blaschke and the Institute of Theoretical Physics at University of Wrocław for hospitality during the realization of this work. This work has been funded by the programme TOGETHER from Région Pays de la Loire, a Helmholtz Young Investigator Group VH-NG-822 from the Helmholtz Association and GSI, and EU Integrated Infrastructure Initiative HadronPhysics3 Project under Grant Agreement n. 283286. JMTR also thank support from the Spanish Ministerio de Ciencia e Innovación under contract FPA2013-43425-P.

Appendix A In-medium quark-antiquark propagator

For a particular meson channel $M$, the quark-antiquark 2-propagator $\Pi^{M}_{qq^{\prime}}$ at finite temperature and density is of key importance for the derivation of the grand-canonical potential at ${\cal O}((1/N_{c})^{0})$ in the large-$N_{c}$ expansion. In the imaginary-time formalism it is expressed as Rehberg:1995nr ; Torres-Rincon:2015rma

[TABLE]

where $S_{q}(i\nu_{n},{\bf k})$ is the propagator of quark of flavor $q$ in the Hartree approximation. The trace is to be taken in color, flavor and spin spaces ($\textrm{Tr}=\textrm{tr}_{c}\textrm{ tr}_{f}\textrm{ tr}_{\gamma}$). The factor $\Omega_{q}$ is

[TABLE]

where $\mathbb{I}_{c}$ is the unit in color space, and $\Gamma_{M}=\{\mathbb{I},i\gamma_{5}\}$ for scalar and pseudoscalar channels, respectively. The flavor trace provides a factor 2 for all mesons (already factorized out in the expressions of Ref. Rehberg:1995nr ), and the color trace gives a trivial factor $N_{c}$. Then, we have

[TABLE]

After performing the Dirac trace and the Matsubara summation one obtains Rehberg:1995kh ; Rehberg:1995nr

[TABLE]

where we have analytically continued the bosonic Matsubara frequency $i\omega_{n}$ to an arbitrary complex variable ($i\omega_{n}\rightarrow z\in\mathbb{C}$). The sign $\mp$ corresponds to pseudoscalar and scalar channels, respectively.

The $A_{q}$ function reads

[TABLE]

with the generalized Fermi-Dirac distribution functions of Eqs. (18,19) accounting for the Polyakov loop effects. The $B_{0}$ function is defined as Rehberg:1995kh ; Rehberg:1995nr

[TABLE]

where $E_{q}=\sqrt{k^{2}+m_{q}^{2}}$ and $E_{q^{\prime}}=\sqrt{({\bf p}-{\bf k})^{2}+m_{q^{\prime}}^{2}}$. Notice that the $B_{0}$ function is a complex function, whose first argument, $z$, is complex in general.

To compute meson masses and widths one should look for possible poles of the $q\bar{q}$ scattering amplitude in the complex energy plane, taking ${\bf p}=0$ (meson at rest). The scattering amplitude in a certain meson channel $M$ reads Torres-Rincon:2015rma

[TABLE]

where ${\cal K}_{M}$ is the effective four-fermion coupling introduced in Sec. III.1.

Following the standard results of scattering theory, any pole of this amplitude along the real axis below any two-body threshold (physical Riemann sheet) is identified with a bound mesonic state. Alternatively, a pole in the lower half complex plane above the two-body threshold (in the second or superior Riemann sheets) is associated to a resonant mesonic state with certain decay width (related to the imaginary part of the pole position). The first situation occurs for temperatures below the Mott temperature, and the later above the Mott temperature. Therefore, the Mott effect represents the effective realization of deconfinement in the context of the NJL/PNJL models.

Appendix B Scattering phase shift

The Jost representation of the S-matrix element $s_{M}$ is taylor1972scattering

[TABLE]

where the denominator plays the role of the Jost function taylor1972scattering . Due to the Schwarz reflection property

[TABLE]

for arbitrary $z\in{\mathbb{C}}$, one can deduce that $s_{M}$ is a pure phase. This phase is what we call the scattering phase shift, defined as in taylor1972scattering

[TABLE]

Equating Eq. (57) and Eq. (59) one obtains

[TABLE]

The structure in the right-hand side appears in the last line of Eq. (42). The logarithm can be simply expressed as

[TABLE]

to obtain Eq. (43).

On the other hand, the scattering phase shift is computed with the help of Eq. (60) and Eq. (58),

[TABLE]

To exploit Eq. (62) one needs to evaluate the quark-antiquark propagator (App. A) for real energies at $z=\omega-\mu_{M}+i\epsilon$ with nonzero momentum ${\bf p}\neq 0$. For the particular case ${\bf p}=0$ one can use the formulas presented in Ref. Rehberg:1995nr to compute meson and diquark masses, as we already did in our previous work Torres-Rincon:2015rma . However, at ${\bf p}\neq 0$ we have found that the method in Ref. Rehberg:1995nr does not provide consistent results in the so-called Landau cuts, perhaps due to an inconsistent change of integration variables in that reference. Therefore, to study the case ${\bf p}\neq 0$ we have evaluated Eq. (62) by numerical complex integration and taking thereafter the limit $\epsilon\rightarrow 0$.

We have numerically checked whether the scattering phase shift at finite momentum behaves approximately in a Lorentz invariant way. This assumption is used in some of the previous works Zhuang:1994dw ; Hufner:1994ma ; Blaschke:2013zaa , arguing that the phase shift can be taken as a function of the Mandelstam variable $s=\omega^{2}-p^{2}$ only. We have computed the scattering phase shifts at finite momentum ${\bf p}\neq 0$ for the pion case. In Fig. 20 we present the phase shift as a function of $p$ and $s=\omega^{2}-p^{2}$ for two different temperatures: $T=1$ MeV (left panel) and $T=300$ MeV (right panel), and $\mu_{q}=0$. If the phase shift is Lorentz invariant, then it should be independent on the momentum $p$ when it is written as a function of $s$.

In the left panel we observe that this invariance is slightly broken at high $s$. However, this region is expected to be thermally suppressed by the distribution function when the phase shift is integrated in Eq. (46). In the right panel, we observe an additional violation due to the presence of the Landau cut ($s<0$) at finite momentum, not present at $T=0$. However, the numerical value of the phase shift at this cut is small, so one can expect that the contribution of that region to the thermodynamical potential is negligible. In conclusion, we can safely assume that the assumption of Lorentz-invariance is approximately valid.

When the quark masses are different, even at ${\bf p}=0$, a Landau cut at finite temperature appears due to an additional cut along the real axis of the $B_{0}$ function. This is important, for instance, for the kaon/antikaon case, where one propagates a light quark and a strange antiquark (or vice versa). Therefore, it is convenient to know the general analytical structure of the quark-antiquark propagator (or the $B_{0}$ function). For the case ${\bf p}=0$, with a finite UV cutoff $\Lambda$, and two different quark masses $m_{1}\neq m_{2}$, the branch points for the different cuts are shown in Fig. 21. Without loss of generality we take $m_{1}>m_{2}$ to avoid a crossing of the Landau cuts, which would make the pictorial representation difficult.

In Fig. 21 we have denoted $\Lambda_{1,2}=\sqrt{m_{1,2}^{2}+\Lambda^{2}}$. We keep the UV cutoff $\Lambda$ in all the momentum integrals in this Appendix. Hence, we have been able to numerically check the limits described in Fig. 21. However, as argued in the text, in our actual calculations we remove the cutoff in all convergent integrals at finite temperature.

To grasp the meaning of this branch points let us consider the case of a pion at rest. At any temperature below the Mott temperature, the pion appears as a bound state in the scattering amplitude of a light quark and a light antiquark. It is unable to decay into a quark-antiquark pair because its mass is smaller than the mass of the quark-antiquark pair Torres-Rincon:2015rma , unless it has enough kinetic energy to do so. In this case, the pion pole lies on the real energy axis at $\omega_{pole}=m_{\pi}$, below the two quark threshold $\omega_{th}=2m_{q}$. As a consequence, the decay width (imaginary part of the pole position) is zero. The polarization function $\Pi_{\pi}(\omega+i\epsilon,0)$ is real for $\omega<\omega_{th}$.

We find that the function $1-2{\cal K}_{\pi}\Pi_{\pi}(\omega+i\epsilon,0)$ is positive for $\omega<\omega_{pole}=m_{\pi}$. Therefore from Eq. (62) we expect $\delta_{\pi}=0$. At $\omega_{pole}$ this function becomes zero by definition, and changes sign for $\omega>\omega_{pole}$. The phase shift presents a sudden jump from $\delta_{\pi}=0$ to $\delta_{\pi}=\pi$. This jump signals the existence of the bound state at $m_{\pi}$.

The phase shift keeps its constant value $\delta_{\pi}=\pi$ until $\omega_{th}=2m_{q}>m_{\pi}$, where the propagating quark-antiquark pair has enough energy to dissolve. This is the beginning of the so-called unitary cut, producing an imaginary part of the polarization function, which becomes discontinuous along the real axis. The nonzero imaginary part makes the phase shift to decrease according to Eq. (62), evolving non-trivially with $\omega$. At the upper limit of the unitary cut $\omega_{\Lambda}=\Lambda_{1}+\Lambda_{2}=2\sqrt{\Lambda^{2}+m_{q}^{2}}$, the polarization function takes again real values, and the phase shift equals again $\delta_{\pi}=0$, in accordance with the Levinson theorem taylor1972scattering ,

[TABLE]

We plot the phase shift for $T=1$ MeV in Fig. 22, where one can clearly observe the jump from 0 to $\pi$ at $\omega_{pole}=m_{\pi}$. Then, the phase shift starts to decrease from $\omega_{th}=2m_{q}\simeq 970$ MeV until $\omega_{\Lambda}=\Lambda_{1}+\Lambda_{2}=2\sqrt{\Lambda^{2}+m_{q}^{2}}\simeq 1495$ MeV.

At temperatures above the Mott transition the pion is a resonant state, whose pole position contains a negative imaginary part. The polarization function is a complex function for all energies and the phase shift is nonzero everywhere. If the decay width is small, then the position of the resonance can still be captured by a rapid increase of $\Delta\delta_{\pi}\sim\pi$ in the scattering phase shift.

As anticipated, an additional complication occurs when the two quark masses are different, like for kaons and antikaons. In this case the additional Landau cut provides an extra contribution at $T\neq 0$ to the kaon phase shift Yamazaki:2013yua .

In Fig. 23 we summarize all these considerations by plotting the scattering phase shifts as functions of the energy and temperature for ${\bf p}=0$. We show the pion case in the left panel and the kaon one in the right panel. The Landau cut is visible in the kaonic case at low energies. We have checked that its range coincides with the theoretical limits $(\sqrt{\Lambda^{2}+m_{s}^{2}}-\sqrt{\Lambda^{2}+m_{q}^{2}},m_{s}-m_{u})$. This is the kinematic range where a kaon can decay by absorbing a quark from the thermal bath, leaving another quark in the final state.