{'fraction_non_alphanumeric': 0.05265748031496063, 'fraction_numerical': 0.0204361790302528, 'mean_word_length': 4.370566142718041, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Extensions and further applications of the nonlocal Polyakov--Nambu--Jona-Lasinio model

The nonlocal Polyakov-loop-extended Nambu--Jona-Lasinio (PNJL) model is further improved by including momentum-dependent wave-function renormalization in the quark quasiparticle propagator. Both two- and three-flavor versions of this improved PNJL model are discussed, the latter with inclusion of the (nonlocal) 't Hooft-Kobayashi-Maskawa determinant interaction in order to account for the axial U(1) anomaly. Thermodynamics and phases are investigated and compared with recent lattice-QCD results.


Introduction
The Polyakov-loop-extended Nambu and Jona-Lasinio (PNJL) model has become a widely used "minimal" approach to deal with important aspects of the QCD phase diagram, in close contact and comparison with lattice-QCD results whenever available. Its basis is the NJL model [1][2][3][4][5], generalized by incorporating Polyakov-loop dynamics that is controlled by a suitably chosen effective potential [6][7][8][9][10]. Once the dynamical input is fixed in terms of a few parameters by reproducing basic facts from hadron physics in vacuum and QCD thermodynamics, the PNJL model permits to explore, at least schematically, regions of the phase diagram that are not accessible to lattice QCD but nonetheless of principal interest. A particularly attractive feature of this model is its ability, through the coupling of the Polyakov loop to the quarks as dynamical quasiparticles, to study the correlated pattern of the deconfinement and chiral crossover transitions.
In recent work [11][12][13][14][15][16] the PNJL model has been generalized and extended by introducing nonlocal effective interactions of the quark fields. This was an important step in order to avoid the traditional sharp momentum-space cutoff characteristic of former (local) NJL models. The nonlocal version now permits to establish contacts with QCD at the level of quark quasiparticle properties. In Refs. [11,12] the first part of such a program has been accomplished by generating momentum-dependent dynamical quark masses from a generalized gap equation, in close correspondence with (Landau gauge) Dyson-Schwinger equations or (extrapolated) lattice-QCD results. At this stage, however, quark wave-function renormalization effects were not yet incorporated. The residue of the quark propagator at its quasiparticle pole was strictly Z ≡ 1. The commonly used chiral quark interactions are not explicitly momentum dependent and, consequently, do not provide a mechanism for reducing the wave-function renormalization factor Z(p) from unity (at least not at mean-field level).
In the present work we generalize the PNJL model further in order to improve on this issue. A vector-type derivative coupling between quarks will be introduced such as to reproduce results for Z(p) from lattice QCD. The two-flavor version will be dealt with in Sect. 2, followed by the three-flavor case with inclusion of the axial U(1) anomaly in Sect. 3. The thermodynamics based on this extended PNJL model will be derived in Sect. 4. The discussion of chiral and Polyakov-loop crossover transitions at zero quark chemical potential will be updated, also in view of recent lattice results. Finally, extrapolations to nonzero quark chemical potentials and the resulting phase diagram will be explored. We present a summary and an outlook in Sect. 5.

.1 Quark propagator and wave-function renormalization
Consider as a starting point the full, renormalized Euclidean quark propagator S(p), the inverse of which is represented by the Dyson-Schwinger (DS) equation It involves the renormalized (dressed) gluon propagator D µν ab and the full quark-gluon vertex function Γ ν , together with renormalization constants Z 1 , Z 2 . The solution of the DS equation (2.1) has the general form . (2. 2) The functions A(p 2 ) and B(p 2 ) depend on the renormalization scale and encode the Lorentzvector and -scalar dressings of the propagating quark by its interaction with the full gluon field. The wave-function renormalization factor, Z(p 2 ) = A −1 (p 2 ), accounts for the deviation of the quark's quasiparticle strength from unity, induced by momentum-dependent vector-type interactions. The dynamical quark mass function, M (p 2 ) = B(p 2 )Z(p 2 ), is generated primarily by nonperturbative mechanisms involving the scalar quark condensate ψ ψ that acts as an order parameter of spontaneously broken chiral symmetry. In the chiral limit (m 0 = 0), the mass M (p 2 ) remains strictly zero at all orders of QCD perturbation theory.

Nonlocal Nambu-Jona-Lasinio model with two quark flavors
In Refs. [11,12] we have already shown how a nonlocal generalization of the Nambu and Jona-Lasinio (NJL) model [1] can be constructed starting from the SD equation (2.1) in the case of Z(p 2 ) ≡ 1. We are now going to write down the extended nonlocal action, S E , for the two-flavor NJL model including a momentum-dependent renormalization function Z(p 2 ) (following Refs. [14,15,[17][18][19]). Assuming separability as in Refs. [11,12], i. e, g 2 D µν ab (p − q) → G 2 C(p) C(q) δ µν δ ab , where G is a coupling strength of mass dimension −2, we can write in Euclidean space-time: with the quark fields ψ(x). We work in the isospin limit with equal current quark masses for up and down quarks, m u = m d ≡ m 0 . The nonlocal currents j a (x), J(x) are given by Here the usual chiral scalar-isoscalar and pseudoscalar-isovector NJL operators, Γ a∈{0,...,3} = (1, iγ 5 τ ) are now supplemented by vector currents with derivative couplings,ψ(x ) The nonlocality distribution C(z) is the same as in Refs. [11,12]. In momentum space: normalized as C(p = 0) = 1. The high-momentum behavior, C(p 2 ) ∝ α s (p 2 )/p 2 , with the QCD running coupling α s (p 2 ), is known from the operator product expansion for the dynamically generated quark mass [20]. It is matched at a scale λ to the nonperturbative infrared sector of C(p 2 ) that is conveniently parametrized as a Gaussian. Its width reflects a typical instanton size or the characteristic gluonic field strength correlation length, d C ∼ (0.3-0.4) fm. The matching scale λ is fixed by requiring continuity and differentiability of C(p 2 ). The detailed low-momentum form of (2.6) is actually not essential because all relevant integrals involve C(p 2 ) multiplied by powers of p 3 . All four quark currents j a (x) with a = 0, 1, 2, 3 have the same distribution C(z) to ensure chiral invariance. The new element in this generalized nonlocal approach is the vector-type derivative current J(x) of Eq. (2.5) that is separately invariant under global chiral transformations. We follow Refs. [18,19] at this point. The derivative coupling proportional to γ µ ∂ µ in J(x) induces the momentum-dependent wave-function renormalization in Eq. (2.2). Apart from the overall coupling G that has been factored out in Eq. (2.3), the coupling strength now features an additional independent parameter κ with mass dimension 1 in the denominator at each vertex. The nonlocality distribution F(z) associated with J(x) is again normalized in momentum space as F(p = 0) = 1. It is given a Gaussian parametrization, with d F to be adjusted such that the momentum-dependent renormalization factor Z(p) is well reproduced. Note that this ansatz is actually consistent with the one for C(p) in Eq. (2.6): according to Ref. [21] there is no modification of Z(p 2 ) = 1 for p 2 > λ 2 coming from a Landaugauge operator product expansion up to first order in the current quark mass. Given the action (2.3), we apply the standard bosonization procedure. Consider the partition function Z = Dψ Dψ e −S E (2.8) and introduce five bosonic fields σ, v and π k plus five auxiliary fields S, V, P k (k ∈ {1, 2, 3}). Next, insert a "one" in terms of (Euclidean) delta functions on the right-hand side of Eq. (2.8), and write the partition function as (2.9) The inserted delta functions imply a replacement of j 0 , J and j k by S, V and P k , respectively. The first exponential, after re-insertion of Eqs. (2.4) and (2.5) for the currents, has terms of the formψÂ ψ. The path integration over the fermionic fieldsψ and ψ is carried out in the standard way: (2.10) The integration over the auxiliary fields S, V, P k is readily performed by completing the square, leading to where detÂ is the fermion determinant. In momentum space one finds after Fourier transfor- (2.12) Finally, comparison of Eq. (2.11) with the definition (2.8) of the partition function gives the bosonized Euclidean action, where we have used integration in momentum space 2 (for later convenience).
We omit the "hats" on the Fourier transforms and set henceforth C(p) ≡Ĉ(p), φ(p) ≡φ(p) , etc. 2 The arguments of the quadratic expressions in σ, v, π follow immediately from the Plancherel theorem under the assumption that those fields are real-valued in coordinate space.
with q ± = q ± p/2. The fields σ and v obviously mix. In the present context this mixing is of no relevance and we will not discuss it further. More important in the present work is the inverse pion propagator, G π (p 2 ), that is used for the determination of the pion mass: (2.20) Note that no summation is implied in this second derivative: with tr{τ i · τ k } = 2δ ik only the diagonal elements in isospin space survive.

Mean-field approximation and chiral condensates
As mentioned earlier (see Eq. (2.15)), the mean-field valuesσ andv are determined by imposing a minimal action S MF . From Eq. (2.16) we obtain the following two gap equations:  where S F (x, y) is the full quark propagator while S F (x, y) denotes its perturbative part. For nonvanishing current quark mass, m 0 = 0, the perturbative contribution needs to be subtracted, leading to (2.23)

Pion mass and pion decay constant
From the inverse pion propagator G π (p 2 ) the pion mass is determined by the condition Furthermore, the square of the quark-pion coupling constant, g 2 πqq , is defined as the residue of the pseudoscalar-isovector quark-antiquark amplitude, hence g −2 πqq = dG π (p 2 ) dp 2 (2.26) Finally, comparing the definition of the (inverse) pion propagator, Eq. (2.20), with Eq. (2.25) the renormalized one-pion state,π k , is defined as such that δ 2 S bos /(δπ(p) δπ(−p)) = p 2 + m 2 π . The derivation of the pion decay constant in the presence of the wave-function renormalization factor, Z(p), is considerably more involved than the simpler case discussed in Refs. [11,12]. The decay constant f π is defined as the matrix element of the axial current between the vacuum and the renormalized one-pion state at the pion pole, In order to calculate this matrix element one has to gauge the nonlocal action (2.13) by introducing a set of axial gauge fields A i µ (x), i ∈ {1, 2, 3}. Since we are dealing with a nonlocal theory this gauging requires a generalization of the parallel transport. The "connection" that links operators at two different space-time points x, y is given by a Wilson line, W(x, y), defined as where s runs over an arbitrary path connecting x with y. This procedure implies the replace- In the case of local operators it gives the well-known minimal gauge-coupling rule Starting from Eq. (2.9) the replacements (2.30) enter first directly in the kinetic term and, secondly, in the currents j 0 , j k and J of Eqs. (2.4) and (2.5). The fermion determinant of the gauged theory reads, in coordinate space, 7 (2.31) Choosing a straight line connecting the points x and y in the Wilson line W(x, y), we have and its derivative with respect to the gauge field A i The matrix element (2.28) becomes 8 : , Furthermore, the mass function M (p) is given by Eq. (2.18) and includes the wave-function renormalization factor Z(p). With Eq. (2.34) the pion decay constant is now given as evaluated at the pion pole p 2 = −m 2 π . The pion-field component with index i is singled out in this expression and no summation is implied.
Finally Taylor-series expansions of the inverse pion propagator, Eq. (2.20), and of the pion decay constant, Eq. (2.36), recover the Goldberger-Treiman relation, and the Gell-Mann-Oakes-Renner relation, at leading order in the quark mass m 0 . This means that the nonlocal NJL model including wavefunction renormalization preserves chiral low-energy theorems and current-algebra relations. m 0 G κ 3.0 MeV 2.37 fm 2 5.21 GeV Table 1: Current quark mass m 0 and coupling-strength parameters G, κ. Table 2: Output resulting from the parameters given in Table 1. All numbers are given in GeV.

Parameter fixing
The two-flavor case discussed in this section involves primarily the coupling-strength parameters G, κ and the (current) quark mass m 0 . These parameters are chosen such that the pion mass, m π 0.14 GeV, and the pion decay constant, f π 0.09 GeV, are reproduced. At the same time, the normalized nonlocality distributions C(p) and F(p), associated with the chiral scalarpseudoscalar and vector-derivative interactions, are constrained by the available information from lattice QCD or Dyson-Schwinger-equation results for the dynamical quark mass M (p) and the renormalization factor Z(p) (see Figs. 1, 2). The lattice results taken here for orientation are Landau-gauge simulations with two dynamical quark flavors for M (p) extrapolated to the chiral limit [22,23], and with 2 + 1 dynamical quarks for Z(p) (with renormalization scale set at 3 GeV) [24].
With these constraints the nonlocality length scales d C and d F of the Gaussian distributions in Eqs. (2.6) and (2.7) become well within the expected range of gluonic correlation lengths. The matching scale at which C(p) turns over to its high-momentum expression proportional to α s (p 2 )/p 2 is then λ = 0.83 GeV.
The input parameters are summarized in Table 1. The resulting output is collected in Table 2. Note that the small value of the current quark mass, m 0 3 MeV, together with the large magnitude of the chiral condensate qq , is consistent with the renormalization scale µ ∼ 3 GeV characteristic of the lattice-QCD results for M (p) and Z(p) that have been used to constrain the quark quasiparticle propagator. The scale-independent product m 0 qq enters in the Gell-Mann-Oakes-Renner relation (2.38).
The parameter κ, associated with the vector-type derivative coupling, comes in the combination √ Gp/κ wherep is a typical momentum (p ∼ 0.5 GeV) appearing in the relevant amplitudes and integrals. This coupling strength is effectively about one order of magnitude smaller than √ G itself, indicating a convergent hierarchy in a derivative expansion of the quark quasiparticle interactions.
The dynamical quark mass M (p) in Fig. 2, extrapolated towards the low-momentum region, is subject to uncertainties concerning its limit M (p = 0), commonly interpreted as a constituent quark mass. Our standard choice is the solid line in   .17)). Lattice data of Z(p) are taken from Ref. [24] for different current quark masses as indicated.   .18)). Diamonds show unquenched lattice results for the two-flavor case, extrapolated to the chiral limit [22]. The open triangles and circles correspond to extrapolated two-flavor lattice results [23] for quenched and unquenched simulations, as indicated. The solid curve is calculated using the input of Table 1  3 Nonlocal NJL model with three quark flavors 3

.1 Three-flavor action
The action for the improved nonlocal three-flavor model is determined by the action given in Ref. [12], plus the extra term introducing the vector derivative interaction, Eq. (2.5): with U(1) flavor matrices λ α (α = 0, 1, . . . , 8), and The functions C(p) and F(p) are given as in the two-flavor case; again we haveψ(x ) The term in the second line of Eq. (3.1) derives from the 't Hooft-Kobayashi-Maskawa determinant that generates the anomalous breaking of the axial U(1) A symmetry. This term produces an η meson mass about twice as large as the mass of the η meson, once the strength H of this determinant interaction is properly chosen.
The bosonization procedure of the action has already been outlined in the context of the two-flavor model (Sect. 2.2), so we restrict ourselves to the result for the partition function: (3.5) Here we have introduced 18 scalar, pseudoscalar bosonic fields σ α , π α together with a vectorlike field v (α ∈ {0, . . . , 8}) and, additionally, 18 auxiliary fields S α , P α necessary to deal with the six-fermion interactions induced by the 't Hooft term plus one auxiliary field V . In momentum space, the fermion determinant detÂ reads (3.6) The major difference compared with the two-flavor partition function (2.10) are the cubic expressions in the auxiliary fields S α , P α . As a consequence an analytic evaluation of the functional integral over the auxiliary fields is not possible anymore. A suitable method to proceed at this point is the stationary phase approximation (SPA): choose the fields S α , P α so as to minimize the integrand in the bosonized partition function (3.5). A necessary condition imposed on the fields is therefore: where S α , P α are now to be considered as (implicit) functions of σ α , π α . Performing, in addition, the integration over the vectorlike fields V the bosonized action can thus be written: (3.8) From here on we can apply the methods developed for the two-flavor case (Sects. 2.2-2.5) in order to reproduce the pseudoscalar meson spectrum within the nonlocal three-flavor NJL model.

Mean-field approximation, gap equations and chiral condensates
Starting from the action S bos E , Eq. (3.8), a power series expansion is performed around the expectation values of the fields σ α , π α (α = 0, 1, . . . , 8) and Given the constraint imposed by charge conservation, one needs to consider only the fields σ 0 , σ 3 , σ 8 . (In the isospin limit studied later, one has the additional constraint that σ 3 also vanishes.) It is useful to introduce and analogously, S = diag(S u , S d , S s ) = S 0 λ 0 + S 3 λ 3 + S 8 λ 8 . Since π α = P α = 0 to leading order, the action in mean-field approximation reads  (3.7)). One finds: . (3.13c) Finally, the chiral condensate qq is calculated in a way analogous to the two-flavor case using the definition (2.22). This leads to (3.14) now with q = u, d, s. Note that M q (p) → m q for large p. As in Eq. (2.23), the subtraction of the free quark propagator makes sure that no perturbative artifacts are left in qq for m q = 0.

Second-order terms and meson masses
We focus now on the masses of the (pseudoscalar) meson octet, plus the η meson with its large mass related to the anomalous breaking of the U(1) A symmetry. Consider once again second-order corrections to the mean-field action, extracted from a functional Taylor expansion, where the second derivatives, , are evaluated at the mean-field values σ α (x) =σ α , etc. Mixed second-order derivatives including one differentiation with respect to π α (x) vanish because of tr{γ 5 } = tr{γ µ γ 5 } = tr{γ µ γ ν γ 5 } = 0.
A basis change with π ij = 1 √ 2 (λ α π α ) ij gives the standard representation of the pseudoscalar octet: Defining analogously a corresponding matrixσ for the scalar boson octet, the fermion determinant, Eq. (3.6), can be written as The resulting second-order contributions to the action are given by where

18)
q ± = q ± p 2 , and (r ± ) −1 is given as the solution of the system For completeness we have also determined the (inverse) propagators associated with the vectortype field v: .

Renormalization and decay constants
As in the two-flavor case, renormalized fields 9 ,φ(p) = Z −1/2 ϕ ϕ(p), are introduced so that the quadratic part of the Lagrangian can be written in the standard form Given the pole structure of the propagators G −1 P (p 2 ), one arrives at an explicit expression for the renormalization constants: The pseudoscalar meson decay constants in the three-flavor case are defined as where J µ A,α (x) =ψ(x)γ µ γ 5 λα 2 ψ(x) denotes the axial-vector current. As described in Sect. 2.5 we have to gauge the nonlocal action in Eq. (3.1) connecting the fields by a Wilson line and introducing a set of axial gauge fields A α µ (α ∈ {0 . . . , 8}). The only modifications compared to the two-flavor calculation are the replacements τ · ¡ A → λ α ¡ A α , σ → λ α σ α and τ · π → λ α π α in the gauged fermion determinant (2.31). Following the lines of the two-flavor calculation, we find for the matrix elements (3.31): 10 , (3.32) 9 Here ϕ(p) stands generically for any of the fields πα(p), . . . 10 As in the two-flavor case, the integration over α can be performed analytically [25].  with (3.33) The decay constants are derived from the expressions (3.32) and their definitions, Eq. (3.31), by contraction with p µ : evaluated at the corresponding mass p 2 = −m 2 β . Owing to the properties of the Gell-Mann matrices and assuming isospin symmetry (i. e., m u = m d ) one has 11 f αβ = δ αβ f π for α ∈ {1, 2, 3} and f αβ = δ αβ f K for α ∈ {4, 5, 6, 7}.

Parameters and results
The enlarged set of parameters in the three-flavor case now includes, in addition, the strange current quark mass m s and the coupling strength H of the U(1) A -breaking (nonlocal) 't Hooft-Kobayashi-Maskawa interaction. The profiles of the nonlocality distributions, C(p 2 ) and F(p 2 ), are left unchanged, assuming that these distributions are governed by (flavor-independent) gluon dynamics, so flavor-symmetry breaking effects can be ignored in these quantities. With the parameters collected in Table 3, the results for chiral condensates, dynamical quark masses at p 2 = 0 and properties of the pseudoscalar meson nonet including octet-singlet mixing are summarized in Table 4. The input current quark masses, m u,d 3 MeV and m s 70 MeV, are once again compatible with MS masses at a renormalization scale µ ∼ 3 GeV and correspondingly large magnitudes of quark condensates at that scale. The resulting pseudoscalar meson-octet masses are generally within less than 2 % of their measured values, while the pion decay constant comes out just slightly smaller than its value f (0) π 86 MeV in the chiral limit, about 10 % below its empirical value. The result for the η-η mixing angle, θ η = −29.9 • , agrees well with the empirical value from Ref. [26], θ = −29.0 • . 12 Fig. 3 shows the dynamically generated quasiparticle masses of the u-and s-quarks, M u (p) = Z(p)(m u +σ u C(p)) and M s (p) = Z(p)(m s +σ s C(p)), where the flavor-independent distributions C(p) and F(p) enter. Note the interplay of the chiral condensates Eq. (3.14), the scalar fields in Eqs. (3.13) and the wave-function renormalization factor Z(p) in determining those dynamical quark masses. 11 This follows from the fact, that the summands in Eq. (3.32) can be written as Assuming mu = m d one gets the stated result. Note, in particular, that for the two-flavor case, Aij is independent of i, j, hence leading the anti-commutators of the Pauli matrices, as given in Eq. (2.34). 12 Note the different definitions of the η-η mixing angle in this work and in Ref. [26]. The cited number θ = −29.0 • has, however, already been translated to the definition, Eq.   Table 4: Output for condensates, dynamical quark masses, pseudoscalar meson masses, decay constants (all given in GeV) and η-η mixing angles resulting from the parameters listed in Table 3.  4 Thermodynamics: Polyakov loop and nonlocal PNJL model 4

.1 Preparations
The nonlocal NJL approach provides a basic understanding for the mechanisms governing spontaneous chiral symmetry breaking and the dynamical generation of massive quasiparticles from almost massless current quarks, in close contact with QCD itself. But the thermodynamics based on the NJL model goes wrong: colored quasiparticles would exist in regions of temperature T and baryon chemical potential µ B where they should be suppressed by confinement. This is the point where the Polyakov loop enters. It acts as an order parameter for the confinementdeconfinement transition in pure-gauge QCD. 13 The coupling of the dynamical quarks to the Polyakov loop promotes the NJL scheme to the PNJL model [6][7][8][9][10][11][12][13]16]. Once this coupling is active, it prevents colored objects (quarks, diquarks and color-octet qq pairs) from propagating in the hadronic low-temperature phase, while color-singlet (mesonic) quark-antiquark pairs and three-quark compounds are permitted (although the formation of localized baryons is not yet realized at this stage).
In the PNJL context the Polyakov loop Φ is expressed in terms of homogeneous (constant) background fields, A 3 4 and A 8 4 , representing the Euclidean temporal gauge-field components associated with the color-diagonal generators of the SU(3) gauge group: Here N c = 3 is the number of colors, tr c indicates the trace over color space only, β := 1/T is the inverse temperature, and λ 3 , λ 8 are the Cartan-algebra elements of the color SU(3) Lie group. One can derive an effective potential, the Polyakov-Fukushima potential U(Φ, Φ * , T ), that describes the first-order confinement-deconfinement transition in the pure-gauge case (here we adopt the version proposed in Refs. [9,10]): This effective potential satisfies the condition of being symmetric under the center Z(3) of the gauge group: Φ → exp [2πin/3] Φ, with n = 1, 2, 3. The temperature-dependent coefficients b 2,4 (T ) are given in Ref. [9], parametrized such that the first-order deconfinement transition in pure-gauge lattice QCD [28] is accurately reproduced together with its transition temperature T 0 = 270 MeV. Although U(Φ, Φ * , T ) is parametrized in terms of the longitudinal components A 3 4 , A 8 4 of the Cartan algebra, it implicitly includes effects of the non-Cartan elements of color SU(3), as well as effects of the transversal components of the gluon fields A i µ , as pointed out in Refs. [13,29,30].
The PNJL model introduces the coupling of the Polyakov loop to the quarks in the standard gauge-covariant derivative and adds the effective potential U to the mean-field action (2.16). The thermodynamics is described using the well-known Matsubara formalism. The Euclidean time component p 4 of the quark four-momentum is replaced by with Matsubara frequencies ω n = (2n + 1)πT, n ∈ Z and quark chemical potential µ. For antiquarks, the corresponding replacement is Four-momentum integrations are replaced as

Two-flavor nonlocal PNJL model
The PNJL thermodynamic potential with N f = 2 quark flavors becomes 14 HereS −1 is the inverse quark propagator expressed in Nambu-Gor'kov space, where the momentum-dependent mass matrixM is diagonal in color space, andẐ is the wave-function renormalization matrix, also diagonal in color space, In Eqs. (4.5a) und (4.5b), M (p 4 , p ) is the dynamically generated mass, Eq. (2.18), and Z(p 4 , p ) is the wave-function renormalization factor, Eq. (2.17), already encountered in the zero-temperature treatment. 15 The generalized Matsubara frequencies appearing in M (p 4 , p ) and Z(p 4 , p ) are: with ω n = (2n + 1)πT, n ∈ Z the standard fermionic Matsubara frequencies.
In the limit of vanishing quark chemical potential, µ = 0, the thermodynamical potential (4.3) for N f = 2 can be written in the simpler form where the generalized Matsubara frequencies are again given in Eq. (4.5c). The finite-temperature generalization of the gap equations (2.21) are obtained by the requirement of a minimal grandcanonical potential Ω, leading to the necessary conditions According to Refs. [9,10] one has Φ = Φ * in mean-field approximation and, consequently, A 8 4 = 0. In the following we first discuss mean-field solutions at µ = 0. Eq. (4.7) determines the expectation values of the scalar and vector fields,σ(T ) andv(T ), together with the Polyakov loop Φ(T ). The mean fieldsσ andv enter in the calculation of the thermal condensate, qq (T ) = −∂Ω/∂m 0 , via the dynamical quark mass M (ω i n , p ) and the renormalization factor Z(ω i n , p ), when taking the derivative of the pressure P = −Ω with respect to the quark mass m 0 .
A key result is shown in Fig. 4. As in our previous work [11], the entanglement of the chiral and deconfinement transitions is a characteristic result of the PNJL approach, a feature that persists in the nonlocal version including wave-function renormalization effects. The crossover transition temperature in this example is T c = 205 MeV, to be compared with the critical temperature, T 0 = 270 MeV, of the first-order deconfinement transition in pure-gauge lattice QCD that has been used to fix the temperature scale in the Polyakov-loop effective potential (4.2). In Fig. 4 no change has been made in this potential when introducing the coupling to dynamical quarks. In Ref. [31] it has been argued that the transition temperature, T 0 , for the deconfinement transition depends on the number of dynamical quark flavors: T 0 (N f = 2) = (208 ± 30) MeV instead of T 0 (N f = 0) = 270 MeV mentioned earlier. Using this shift in T 0 as input in the Polyakov-loop effective potential U, the chiral transition temperature is lowered to T c 170 MeV shown in Fig. 5.
Compared to the nonlocal PNJL model without wave-function renormalization, there is now the condensatev that becomes temperature dependent. The result is shown in Fig. 6. It is evident, thatv stays constant up to a temperature that is about 1.5T c . Above this temperature the absolute value of the condensatev decreases. The reason for this can easily be understood: the wave-function renormalization factor approaches one at high momenta (see the lattice data in the right picture of Fig. 1); a similar effect must appear when the temperature is increased.

Thermodynamical potential
The thermodynamics of three quark flavors introduces the strange quark, with its heavier mass m s well separated from the masses m u,d of the light quarks. It also introduces the axial U(1) anomaly as a mechanism that generates additional nonlinear dynamics in the resulting gap equations [11,12]. By analogy with the two-flavor case, the stationary phase approximation gives the grandcanonical mean-field thermodynamical potential 16 in terms of the fields (3.13): ). The solid line shows our PNJL model calculation with T 0 = 270 MeV as input in the effective potential U. The dashed curve corresponds to the pure gauge case, compared to lattice data from Ref. [34].
At this point our findings appear to be at variance with considerations in Ref. [31] suggesting that the feedback induced by the dynamical quarks should change the Polyakov-loop effective potential such that T 0 (N f = 0) = 270 MeV is shifted downward to T 0 (N f = 2 + 1) = (187 ± 30) MeV. Such a shift in U would run the crossover profile for the Polyakov loop (solid curve in Fig. 8) down to much lower temperatures in the present approach, out of scale for comparison with the N f = 2 + 1 lattice-QCD data. On the other hand, a more sophisticated treatment of the effective potential U using renormalization-group methods [35,36] might help clarifying the situation. Fig. 9 shows the temperature dependence of the chiral condensates ūu = d d and ss resulting from the present calculation. The softening of the strange quark condensate clearly reflects the stronger explicit chiral symmetry breaking by the s-quark mass m s . One notes that in these condensates, the influence of the quark wave-function renormalization factor Z(p) through the gap equations (3.13) and (3.14) is only marginal in comparison with previous calculations [11,12] that used Z(p) ≡ 1.
Recent (2+1)-flavor lattice-QCD computations of the chiral condensate have been performed using physical light quark masses by the Budapest-Marseille-Wuppertal (BMW) group [33,37] and by the "hotQCD" collaboration [32,38] (the latter employing a highly improved staggered quark (HISQ) action). These lattice data suggest a relatively smooth chiral crossover with transition temperature in the range T c ∼ (160-170) MeV. Our mean-field nonlocal PNJL calculation, using a Polyakov-loop effective potential U with T 0 (N f = 0) = 270 MeV as before, is not able to reproduce these updated lattice results (see dashed curve in Fig. 10). Running the T 0 scale in U down to T 0 (N f = 2 + 1) = 187 MeV and adding pionic contributions from chiral perturbation theory [39] give the qualitatively correct pattern for qq (T ) (see solid curve in Fig. 10). However, as already mentioned, the Polyakov loop would now be systematically shifted to lower temperatures, out of scale with the corresponding lattice data. In addition, such small values of T 0 enhance the inference of the Polyakov loop such that the chiral transition tends to become first order, in contradiction with the lattice results. Again, improved Polyakov-loop effective potentials with systematic implementation of the backreaction from light quarks [35,36] are under discussion in this context.

Finite densities and phase diagram for 2 + 1 flavors
Consider now the phase diagram emerging from the nonlocal PNJL model at nonzero quark chemical potential. The starting point is the (2 + 1)-flavor thermodynamical potential (4.11) treated in mean-field approximation. For calculational simplicity we choose the symmetric combination µ ≡ µ u = µ d = µ s . The case µ s = 0 together with µ u = µ d has been studied earlier in Ref. [12]. The flavor-symmetric choice of the quark chemical potential shifts the µ-scale just slightly upward as compared to the µ s = 0 case. The resulting (T, µ) phase diagram, Fig. 11, is the generic one for PNJL type models at mean-field level: a chiral crossover transition ending in a critical point followed by a first-order transition line down to µ ∼ 0.3 GeV at T = 0. The behavior of the Polyakov loop, detached from the chiral first-order-transition boundary has led to an interpretation in terms of an intermediate "quarkyonic" phase [40]. These features are not changed by the additional effects of quark wavefunction renormalization in the nonlocal PNJL model.
However, recent two-flavor investigations beyond mean field [35,36], including both the quantum fluctuations in the quark quasiparticle sector and the backreaction from the matter fluctuations to the Polyakov-loop effective potential, suggest that the coincidence of chiral and deconfinement crossover transitions continues all the way down to a critical end point that now appears at low temperature, T CEP ∼ 25 MeV, at µ ∼ 0.3 GeV. Similar tendencies have been found in Ref. [41]. That is then the terrain of nuclear physics, however, where the PNJL approach is not supposed to work and different methods (e. g., chiral effective field theory with baryons) must be applied [42]. Ref. [39]). Dashed curve: mean-field result with T 0 = 270 MeV (same as solid curve in Fig. 9). Lattice-QCD data from Refs. [32,33,37,38].
Along the chiral crossover transition at small chemical potential µ, it is convenient to introduce a transition band, T c (µ), specified by the region in which the chiral condensate has dropped to about half of its vacuum value: 0.4 < ψ ψ / ψ ψ 0 < 0.6. This is also the domain in which the chiral susceptibility develops a pronounced maximum. In this region an expansion where the error reflects the uncertainty in locating the transition line in the crossover region. This should be compared to recent values found in lattice QCD: K lat 0.06 from Ref. [43] (analysis of the scaling properties of the chiral condensate and its susceptibilities), and K lat 0.08 from Ref. [44] (Taylor-series expansion), both with an uncertainty of about 10 %.

Summary and outlook
In this work the two-and three-flavor PNJL models have been further extended with inclusion of quark wave-function renormalization effects. This brings the PNJL approach in close contact with Dyson-Schwinger calculations and lattice-QCD results for the Landau-gauge quark quasiparticle propagator, with the wave-function renormalization factor Z(p) of the quark quasiparticle propagator introduced in addition to the momentum-dependent dynamical fermion mass M (p). In the following we summarize the main results of the present studies and point out improvements compared to previous work: (i) Including wave-function-renormalization effects requires a careful re-assessment of chiral low-energy theorems. Pseudoscalar meson masses and corresponding decay constants at zero temperature have been re-derived. The results clearly show that the formalism incorporates fundamental chiral relations such as the Gell-Mann-Oakes-Renner and Goldberger-Treiman relations. In the three-flavor case, the inclusion of the 't Hooft-Kobayashi-Maskawa interaction leads to the correct mass splitting between the η and the η meson.
(ii) The PNJL thermodynamics has now been developed with systematic inclusion of the quark quasiparticle renormalization factor Z(p). The temperature dependence of the chiral condensate and of the Polyakov loop has been calculated, indicating chiral and deconfinement crossover transitions. We have compared our results with recent lattice-QCD computations. Finally, a quark chemical potential has been introduced that enables extensions to the finite-density region of the QCD phase diagram.
(iii) The impact of the wave-function renormalization factor Z(p) compared to previous calculations [11,12] setting Z(p) ≡ 1 is generally quite small over the whole relevant momentum range. This can be understood considering the gap equations at zero temperature: since Z(p) deviates significantly from unity only in the momentum range p 1 GeV, its effect does not contribute much to the relevant integrals because of its suppression by the integration measure.
(iv) With inclusion of Z(p), the chiral and deconfinement crossover transitions tend to become smoother compared to our previous investigations.
(v) The flavor dependence, T 0 (N f ), of the deconfinement temperature scale is an important issue in lowering chiral-transition temperatures in accordance with the tendency recently observed in lattice-QCD computations. We have investigated the impact of different values for T 0 in the Polyakov-loop effective potential and found that such modifications do not relax the strong entanglement between chiral condensate and Polyakov loop, at least at mean-field level.
(vi) Concerning the behavior of the phase diagram at finite quark (or baryon) chemical potentials µ, the backreaction of the matter fluctuations on the Polyakov-loop sector is an important effect, as pointed out in Refs. [35,36], that calls for further studies. Finally, the challenge remains to construct a phase diagram at finite baryon densities in accordance with realistic constraints from nuclear physics.