The phase structure of the Polyakov--quark-meson model beyond mean field

The Polyakov-extended quark-meson model (PQM) is investigated beyond mean-field. This represents an important step towards a fully dynamical QCD computation. Both the quantum fluctuations to the matter sector and the back-reaction of the matter fluctuations to the QCD Yang-Mills sector are included. Results on the chiral and confinement-deconfinement crossover/phase transition lines and the location of a possible critical endpoint are presented. Moreover, thermodynamic quantities such as the pressure and the quark density are discussed.


I. INTRODUCTION
Driven by the heavy-ion programs at GSI, CERN SPS, RHIC and LHC there is strong interest in the properties of strongly interacting matter at extreme temperatures and baryon densities. The access to the phase diagram of QCD is hampered by the fact that lattice simulations at finite density suffer from the sign problem [1]. However, in recent years impressive progress has been made in the understanding of the phase diagram of QCD within QCD effective models such as the Polyakov-loop extended Nambu-Jona-Lasinio (PNJL), see e.g. [3,4], and the Polyakov-loop extended quark-meson (PQM) model, see [6]. In these models the information about the confining glue sector of the theory is incorporated in form of a Polyakov-loop potential that is extracted from pure Yang-Mills lattice simulations. The matter sector of these models has been studied in detail also beyond the mean-field level by taking into account the quark-meson quantum fluctuations mostly within a functional renormalization group (FRG) approach, for recent studies see [7,8].
However, the most difficult problem within these models is the question of how to embed the quantum back-reaction of the matter sector to the gluonic sector. This problem has effectively been resolved in [6] where the change of Λ QCD , and hence the confinementdeconfinement transition temperature T 0 , in the presence of dynamical quarks has been computed within hard thermal and/or hard dense loop (HTL/HDL) perturbation theory. This leads to a flavor and chemical potential dependence of the transition temperature T 0 , and is a qualitatively viable procedure: the change of the non-perturbative scale Λ QCD can be very well estimated within perturbation theory as has been shown and confirmed in many computations at zero-temperature and finite-temperature QCD. Moreover, this perturbative estimate in [6] has been confirmed by first principle QCD computations with the functional renormalization group at real and imaginary chemical potential in [9,10]. On the other hand, very recently the µ-dependence of T 0 has been estimated by constraining PNJL results with those in the statistical model [11]. This links the model parameter T 0 (N f , µ) under certain weak assumptions to experimental data about the chemical freeze-out curve and nicely confirms the theoretical prediction in [6].
This shows that even though being QCD effective models, they provide valuable information about the phase diagram, and in particular help to exclude certain scenarios. Moreover, the models can be understood as specific controlled approximations to full dynamical QCD. Most directly this is realized within the above-mentioned FRG approach put forward in [9,10,[12][13][14]. This link enables us to systematically extend the PNJL/PQM models towards full dynamical QCD.
In the present work we take an important step towards full dynamical QCD and present the first computation with fully dynamical matter sector in the PQM model at finite density. An interesting first FRG computation at vanishing density has been put forward in [15]. There, however, the back-coupling to the glue sector has been neglected. Here, we include the back-reaction of the matter sector as described above. The uncertainty can be estimated by comparing it to the QCD computation in [9,10]. Within the present PQM approach beyond mean field we provide results on the phase boundaries for the confinement-deconfinement and the chiral transitions as well as on thermodynamic quantities such as the pressure, quark number density and quark number susceptibility.
The outline of the paper is as follows: in the next section we discuss the back-reaction of the matter fluctuations to the QCD Yang-Mills sector in more detail. In Sec. III, the Polyakov-loop and its potential is introduced and coupled to the quark-meson model which then defines the Polyakov-quark-meson model. At first, the grand potential of this model is evaluated in mean-field approximation. Afterwards, the flow equation for the grand potential of this model is derived in Sec. IV where also the choice of model parameters is discussed. Sec. V is devoted to the phase structure and some thermodynamical applications, in particular we evaluate the pressure, quark number density and quark number susceptibility. The phase structure, i.e., the chiral and confinementdeconfinement phase transition is explored in detail. Subsequently, the influence of the Polyakov loop on the thermodynamics is investigated and in Sec. VI concluding remarks are drawn.

II. THE PQM MODEL BEYOND MEAN FIELD
The Polyakov-quark-meson (PQM) model was studied in a mean-field approximation for the matter sector in [6]. In this model, the coupling between the pure glue sector, realized by an effective Polyakov-loop potential U, and the quark-meson matter sector is provided by the fermionic determinant. However, in [6] an important step beyond the mean-field approximation was already introduced. More specifically, the pure glue potential U was adjusted on the basis of phenomenological arguments: the phase transition temperature T 0 in the potential U relates to the dynamical scale in QCD, Λ QCD . If the dynamics of the quarks is switched on, this scale is lowered. In [6] we have provided a flavor and quark chemical potential dependent transition temperature T 0 (N f , µ) which was estimated from hard thermal and hard dense loop considerations. Here, we shall elaborate this argument in more detail, also in order to show that no double counting is involved in the computation and to tighten the link to first principle QCD. We also stress that the above phenomenological argument works well at large chemical potential µ as well as at small chemical potential. For vanishing chemical potential we also compare our results with a recent dynamical QCD calculation [9]. This allows us to qualitatively adjust the phase transition lines at small chemical potential and small temperature, and gives us access to information about a possible critical endpoint as well as the size of a possible quarkyonic matter region at large chemical potential in the QCD phase diagram as suggested in [16].
What is then left is to include the quantum and thermal fluctuations of the matter sector in the presence of the non-trivial Polyakov loop. This is achieved with functional renormalization group (FRG) methods and is done in Sec. IV, for reviews see e.g. [17,18]. In combination, this allows for a Polyakov-loop effective model computation which already takes into account the full dynamics of QCD phenomenologically.
We proceed with the explanation of the N f and µ dependence of the transition temperature T 0 . In [9] a full two-flavor QCD computation was put forward for vanishing and imaginary chemical potential. The related functional RG equation for the effective action is provided in a simple diagrammatic form in Fig. 1. The first two loops in Fig. 1 stand for quantum fluctuations (gluons and ghosts) in the pure glue sector that, e.g., generate the Polyakov-loop potential. The third loop represents the quark fluctuations, and the last term encodes mesonic fluctuations generated by dynamical hadronization [18- 20]. The crosses denote the cut-off insertion which restricts the loop momenta to that about the cut-off scale, and φ stands for all fields. For more details see [9,10,12].
We would like to emphasize that the above equation for the effective action is fully coupled. An important example is the gluonic (and indirectly the ghost) propagation that is modified in the presence of dynamical quarks. This aspect is visualized in Fig. 2 where the quark contri- amounts to dropping the matter loops in Fig. 1. The pure Yang-Mills computation has been done in [21], and the related T 0 agrees quantitatively with lattice predictions which are used to adjust the T 0 parameter in the Polyakov-loop potential U. These findings also elucidate that T 0 is nothing else but the dynamical YM-scale Λ YM inherent in the ghost and gluon propagators. This supports the quantitative accuracy of such an approach.
In turn, the quantum dynamics of the quark-meson model is obtained by dropping the Yang-Mills diagrams in Fig. 1, that is the gluon and ghost loop. This has been studied extensively in the literature, for recent works see e.g. [7,8]. In summary, it is this simple additive structure for the different contributions that allows to systematically improve quark-meson models towards full QCD within an FRG approach. Moreover, it also allows us to directly use full QCD information within the PQM/PNJL models. This reduces the parameterdependence of these models, and qualitatively enhances the predictive power.
An important example for this structure is the Polyakov-loop potential in these models: the Polyakovloop potential U in full dynamical QCD still comes from the Yang-Mills contributions in Fig. 1. The related Polyakov-loop potential in QCD can be computed from the full dynamical QCD computation in [9], and receives a flavor and chemical potential dependence via Fig. 2. It is clear that the quark contributions to the gluon propagation (and ghost propagation) shown in Fig. 2 lower the dynamical scale in the gluon and ghost propagation, and hence T 0 . In turn, with constant T 0 without a flavor and µ-dependence, this contribution is estimated by the pure Yang-Mills results. We conclude that the hard thermal and hard dense loop adjustment of T 0 is necessary for the full inclusion of quark effects to the Polyakov-loop potential. In particular, the above contributions are not included in the fermionic determinant and thus no double counting is involved.
In summary, the inclusion of the N f and µ-dependence in T 0 as well as the dynamics of the quark-meson sector of the model in the presence of a non-trivial Polyakovloop expectation value gives us a good qualitative control about the full QCD dynamics, and the link is tightened by comparing T 0 (N f , µ) with QCD results in [9,10]. The direct use of the full QCD Polyakov-loop potential will be discussed elsewhere.

III. THE POLYAKOV-QUARK-MESON MODEL
The Polyakov-quark-meson model as put forward in [6] already provides a good approximation of full QCD at low energies and temperatures and not too high densities. Its classical Euclidean action with N f = 2 light quarks q = (u, d) and N c = 3 color degrees of freedom reads where D / (Φ) = γ µ ∂ µ − i gγ 0 A 0 (Φ) with gauge coupling g and h is the Yukawa coupling between mesons and quarks. The temporal component of the gauge field is linked to the order parameter of quark confinement (in pure Yang-Mills), the Polyakov-loop variable Φ, where the color trace tr is in the fundamental representation, P denotes path ordering, and β = 1/T , the inverse temperature. The purely mesonic potential is defined as The isoscalar-scalar σ field and the three isovectorpseudoscalar pion fields π together form a chiral vector field φ. Without the explicit symmetry breaking term c in the mesonic potential the Lagrangian is invariant under global chiral SU (2) L × SU (2) R rotations.
The parameters of the model defined in Eq. (1) are fixed at the low energy physics at vanishing temperature.

A. Polyakov-loop potential
In recent years, several Polyakov-loop potentials have been proposed [22,23], for comparison see [24]. Their functional form is motivated by the underlying QCD symmetries in the pure gauge limit and they all reproduce a first-order transition at T c ∼ 270 MeV for N c = 3 colors in this limit.
The above potential reproduces well the equation of state and the Polyakov-loop expectation value, in particular around the transition [24]. The parameter T 0 = 270 MeV corresponds to the transition temperature in the pure YM theory. As motivated in the previous section, it is left to fix T 0 in the presence of dynamical quarks. Since it is directly linked to the dynamical mass-scale Λ QCD this parameter necessarily has a flavor and chemical potential dependence in full dynamical QCD and T 0 → T 0 (N f , µ).
In the present work we use perturbative relations for fixing the relative scales [6,25]. The results compare well to the full N f = 2 QCD computation in the chiral limit in [9,10]. The latter thus allows for an error estimate of the present procedure. The one-loop β-function of QCD with massless quarks is given by Here, we have assumed a RG scheme that minimizes (part of) the higher-order effects. At leading order the corresponding gauge coupling is given by with α 0 = α(Λ) at some UV-scale Λ. The scale Λ QCD = Λ exp(−1/(α 0 b)) corresponds to the Landau pole of Eq. (10). The temperature dependence of the coupling is also governed by Eq. (10) with the identification p ∼ T . This yields the relation [6] T whereT and α 0 are free parameters. Eq. (11) allows us to determine the N f -dependence of the critical temperature T 0 (N f ). Analogously to [6] we chooseT to be the τ -scale,T = T τ = 1.77 GeV. This constitutes a reasonable UV scale for the mean-field model. Then the pure Yang-Mills input, T 0 (N f = 0) = 270 MeV, leads to α 0 = 0.304. In the present work we shall stick to these values. In addition to the arguments given in [6], the ratio T 0 /T χ in the chiral limit compares well with that computed in the full two-flavor QCD calculation in [9]. Table I summarizes the N f -dependent critical temperature T 0 in the Polyakov-loop potential for massless flavors:  Massive flavors lead to suppression factors of the order T 2 0 /(T 2 0 + m 2 ) in the β-function. For 2 + 1 flavors and a current strange quark mass m s ≈ 150 MeV we obtain T 0 (2 + 1) = 187 MeV. We estimate the systematic error for T 0 (N f ) being of the order +15 −20 MeV related to the scale matching of the present PQM computation with the QCD computation in the chiral limit in [9]. Note, however, that the link to QCD qualitatively improves the error estimate in comparison to the estimate done in [6].
As argued in the last section, in addition to the flavor dependence of T 0 we introduce a chemical potential dependence via a µ-dependent running coupling b, which should push the confinement-deconfinement transition temperature down close to the chiral transition line. This can be achieved by defining with The factorγ is a parameter governing the curvature of T 0 (µ) and b µ ≃ 16 π N f as in [6]. As for the N f -dependence the µ-dependence in Eq. (12) compares well to that found in QCD [9,10]. Based on the results there we estimate the systematic error with 0.7 γ 1, and we shall investigate theγ-dependence of our results in Sec. V.

B. Grand Potential in Mean-Field Approximation
All thermodynamic properties of the PQM model follow from the grand potential. It is a function of the temperature and one quark chemical potential since we consider the SU (2) f -symmetric case in this work and set µ ≡ µ u = µ d .
In the mean-field approximation certain quantum and thermal fluctuations in the path integral representation of the grand potential are neglected. The mesonic quantum fields are replaced by their corresponding classical expectation values and only the integration over the quark loop is performed which is modified by constant gluon background fields in the PQM model [26]. The final potential in mean-field approximation reads and consists of the quark contribution including the Polyakov-loop variables with the quark/antiquark single-quasiparticle energies E p = p 2 + m 2 q and the constituent quark mass m q = hσ. The purely mesonic potential U is given by Eq. (3) and the effective Polyakov-loop potential U, e.g., by Eq. (4). Details of the potential derivation can be found in [6]. The quark contribution involves a divergent vacuum term which can be regularized. As shown in [27,28] this term is important and modifies the underlying thermodynamics. Since this term upgrades the standard mean-field approximation it is neglected here whereas it is included in the full RG approach.
The solution of the corresponding equations of motion are obtained by minimizing the thermodynamic potential with respect to the three mean fields σ, Φ andΦ, i.e., The solutions to Eq. (16) provide the chiral σ and Polyakov-loop expectation values Φ and Φ as functions of the temperature and quark chemical potential.

IV. FLOW EQUATION
The non-perturbative FRG method has a wide range of applicability. In the context of equilibrium statistical physics it represents a very efficient way to describe critical phenomena and in particular phase transitions. One particular formulation of RG flows is based on the concept of the effective average action Γ k where k denotes a RG momentum scale, for reviews see e.g. [17,18]. For a system with bosonic (ϕ) and fermionic fields (ψ), the variation of Γ k with the RG scale (t = ln k) is governed by the flow equation where Γ (i,j) k denotes the ith (jth) derivative of Γ k with respect to the ϕ (ψ) fields. The trace involves a ddimensional momentum integration and a summation over all inner spaces (flavor, color and/or Dirac). The regulator functions R k,i and their derivatives implement Wilson's idea of integrating successively over narrow momentum shells and ensure that the flow equation is both infrared and ultraviolet finite. We employ the optimized regulator functions [29] which depend only on the spatial components of the momentum and for bosonic fields is given by whereas the fermionic regulator for the PQM model reads Further details and the finite temperature and density generalization of the flow equation can be found, e.g., in Refs. [8,[30][31][32][33][34].
In the PQM model with the classical action, Eq. (1), the flow Eq. (17) for the quark-meson sector also depends on the Polyakov loop via the Dirac term. Here we will take into account the full dependence of the quark and meson fluctuations on a general constant Φ,Φ background. In summary, this corresponds to the truncation Γ k = d 4 x ψ (D / + µγ 0 + ih(σ + iγ 5 τ π)) ψ (20) where Ω k [σ, π, Φ,Φ] is the full effective or grand potential in a general background Φ andΦ. It is the quantum analogue of the mean-field potential Ω MF , given in Eq. (14) and, in the present truncation, carries the full k-dependence of the effective action Γ k . The approximation Eq. (20) has also been used for the matter sector in [9] and in the first PQM study beyond mean field at vanishing density [15]. Finally, we are interested in Ω k=0 evaluated at the solutions of the quantum equations of motion. This corresponds to the thermodynamic potential where from all temperature, chemical potential and field derivatives follow. The spatial gauge field is set to zero while we keep the temporal component as a constant mean-field background. With this truncation the flow equation for the effective potential in leading order derivative expansion reads with the Polyakov-loop enhanced quark/antiquark occupation numbers N q (T, µ; Φ,Φ) = 1 + 2Φe (Eq−µ)/T + Φe 2(Eq−µ)/T 1 + 3Φe (Eq−µ)/T + 3Φe 2(Eq−µ)/T + e 3(Eq−µ)/T and Nq(T, µ; Φ,Φ) ≡ N q (T, −µ;Φ, Φ) , (22) see also the QCD study [9] at vanishing and imaginary chemical potential, and the PQM study [15] at vanishing chemical potential. In the present work we explore the full phase diagram at real chemical potential. The number of internal quark degrees of freedom is denoted by ν q = 2N c N f = 12. This equation describes the flow of the full quark and mesonic subsystem modified by a constant background field. The grand potential also depends on the Polyakov-loop variables and the expectation value of the square of the chiral 4-component field φ 2 = σ 2 + π 2 which coincides with σ 2 since π 2 = 0.
The quasi-particle energies for i = q, σ, π are given by with the corresponding squared constituent quark and meson masses respectively Primes denote the φ 2 -derivative of the grand potential, i.e., Ω ′ k := ∂Ω k /∂φ 2 . In the limit of vanishing background fields, i.e., when Φ,Φ → 1, the extended occupation numbers simplify to the usual Fermi-Dirac distribution functions for quarks and antiquarks N q (T, µ; 1, 1) = Nq(T, µ; 1, 1) = 1 1 + exp((E q + µ)/T ) , and the flow of the quark-meson model is recovered [8,35]. The flow equation (21) constitutes a set of coupled, highly non-linear partial differential equations that cannot be solved with analytical methods since the right hand side of (21) depends on derivatives of the unknown potential Ω k . For the sake of full quantitative precision and in order to facilitate the access to the potential firstorder phase transition region in the QCD phase diagram we do not resort to Taylor expansions but rather compute the full numerical solution of the effective potential Ω k [σ, π, Φ,Φ]. This extends the analyses of the full effective potential in the quark-meson model, see [8], to the present case. The solution of Eq. (21) yields the thermodynamic potential Ω eff [σ, π, Φ,Φ] ≡ Ω k=0 [σ, π, Φ,Φ] for the Polyakov-extended quark-meson model in the infrared. The expectation values of the fields are determined by the quantum equations of motion, in analogy to the mean-field analysis in Eq. (16).
A full study including algorithmic differentiation (AD) techniques [36] will be presented elsewhere [37]. In the present work we shall simply evaluate the effective potential Ω k on the solutions Φ(σ),Φ(σ) of the mean field EoMs. We shall argue that this is already a quantitatively reliable approximation to the full solution: note that the present truncation introduces a cut-off, and hence implicitly a momentum-dependent fermion propagator. For large temperatures T in the deconfined regime the fermion propagator is cut-off by the Matsubara mass πT . At small temperatures the theory is chirally broken and the fermion propagator exhibits a mass of the order of the chiral scale ∼ 360 MeV. Moreover, at sufficiently large densities the fermionic fluctuations are already well-described by the one loop determinant. In summary, the fluctuation dependence of the Polyakov loop can be treated as a perturbation in the whole phase diagram. Accordingly, the back-reaction of the Polyakov loop beyond mean field to the matter fluctuations is suppressed for all T and µ. Indeed, in comparison to the full solution the minimum for σ varies within ±3 MeV whereas Φ,Φ are naturally more sensitive to this approximation and vary within ±20 MeV. Note, that the above structure provides further non-trivial reliability for the Polyakov-loop extended models. Seen as an expansion towards full QCD they partially allow for perturbative arguments about the mean-field analysis. The full computation and a more detailed analysis of this structure will be provided in [37]. The above structure also emphasizes the crucial input: the back-reaction of the matter sector to the Polyakov loop via T 0 (N f , µ).
It remains to determine the initial effective action Γ Λ , or, more precisely, the effective potential Ω Λ at the initial scale Λ = 950 MeV in the UV. First of all, consistency with the flow Eq. (21) leads to a (relevant) term originated in the fermionic loop, Ω ∞ Λ [σ, π, Φ,Φ], see Eq. (28). This term is relevant for the correct thermodynamics and also includes fermionic vacuum fluctuations. The fielddependent part of Ω Λ consists of a sum of the quarkmeson potential Eq. (3) and the external glue input, the Yang-Mills Polyakov-loop potential U. This yields It is left to determine the model parameters in the quark-meson sector. As in the mean-field analysis, they are fixed to reproduce the physical quantities in the vacuum such as the pion decay constant, f π = 93 MeV, pion and sigma meson masses, m π = 138 MeV, m σ = 500 MeV, and the constituent quark mass m q = 298 MeV in the infrared. The explicit symmetry breaking term c = m 2 π f π is an external field and is scale independent.

V. PHASE STRUCTURE
We first discuss the order parameters at vanishing chemical potential and compare our findings with lattice results [38,39] as well as continuum QCD results [9]. We evaluate the validity of our approach, in particular, the need for the inclusion of the back-reaction of the matter sector via a N f and µ-dependence of the dynamical transition parameter T 0 in the Polyakov-loop potential. Fig. 4 summarizes our findings for µ = 0 and the dynamical transition parameter T 0 = 208 MeV. It shows in the left panel the chiral condensate, normalized with the vacuum value and the (degenerated) Polyakov-loop variable as a function of the temperature. In the right panel, the corresponding temperature derivatives of the order parameters are displayed. As one can see from Fig. 4, both widths in the temperature derivatives are comparable. For the chiral transition a slightly larger critical temperature, T c ∼ 190 MeV as for the deconfinement transition, T c ∼ 175 MeV, is found. This relates to the standard scenario that chiral symmetry restoration occurs at higher temperatures than deconfinement. We remark that the arguments for the standard scenario are strictly valid only for phase transitions rather than for smooth crossovers as in the present case. Moreover, the chiral and confinement-deconfinement crossover temperatures agree within the respective widths. The above findings compare well with a two flavor QCD computation in the chiral limit [9], and also with recent lattice results [38,39]. For the comparison with the latter one has to bear in mind that the definitions of the chiral order parameter are ambiguous which influences the crossover temperature.
If we switch off the back-reaction of the matter sector and stick to the pure Yang-Mills value of the transition parameter, T 0 = 270 MeV the situation changes quantitatively as can be seen from Fig. 5: the chiral and Polyakov derivatives peak at similar critical temperatures about T c ∼ 220 MeV but the width of the Polyakov derivative is very broad suggesting a peak substructure. This broad structure does not compare well with the lattice predictions nor does the high crossover temperatures. This demonstrates the importance of the dynamics of T 0 already at vanishing chemical potential. Now we extend our analysis to finite chemical potential. Our findings are summarized in Fig. 6. The left panel shows, for comparison, the phase diagram without back-reaction of the matter sector to the glue sector. The right panel shows the phase diagram with the full dynamics: the quark-meson dynamics in the presence of the Polyakov-loop background is included with the flow Eq. (21) whereas the back-reaction of the matter sector to the glue sector is encoded in T 0 (N f , µ) as given in Eq. (12). In the right panel of Fig. 6 we have usedγ = 0.85 and T 0 (2, 0) = 208 MeV which compares well to the QCD results in [10] for small densities. The shaded band corresponds to the width of dΦ/dT at 80% of its peak height. We observe, that the width of this band shrinks with increasing µ. Thus, the deconfinement transition gets sharper at higher µ. Furthermore, we can unambiguously distinguish the peak positions of the chiral and deconfinement transition in the temperature derivatives. Even when we vary the parameterγ almost no differences emerge. Only for value ofγ close to one an ambiguity in the peak positions of the temperature derivatives arises, similar to the constant T 0 -case discussed below.
We see that the deconfinement transition line stays close to the chiral phase boundary. It has been speculated that the quarkyonic phase [16] is signaled by a region with confinement and chiral symmetry. Mean field computations in the PNJL and PQM model with constant T 0 have supported this scenario. As has been shown in [6], this does not hold true if the dynamics of the transition parameter T 0 is taken into account. Here we have confirmed that in the fully dynamical PQM model the prediction in [6] holds.
The chiral first-order transition line arising at small temperatures and large chemical potentials terminates in a CEP which is located at µ CEP ∼ 292 MeV, T CEP ∼ 23 MeV. As discussed in [40] the precise location of the CEP depends on the chosen parameters in the vacuum, in particular on the sigma meson mass. The back-bending of the first-order transition line in the phase diagram towards smaller chemical potential is typical for a FRG calculation [8]. But very close to the µ-axis the slope of the first-order line tends back to infinity similar to a mean-field treatment (not shown in the figure). This behavior is also in agreement with the Clausius-Clapeyron relation according to which the transition line should hit the µ-axis perpendicular.
If we switch off the back-reaction of the matter fluctuations to the Yang-Mills sector and choose a constant T 0 = 208 MeV, the above picture changes drastically. For finite chemical potential the Polyakov loops Φ andΦ start to deviate and the widths of their temperature derivatives increase over the whole phase diagram. The resulting phase diagram for a constant T 0 is summarized in the left panel of Fig. 6. In the vicinity of the intersection point of the chiral transition and the deconfinement transition which are both smooth crossovers around µ ∼ 180 MeV a double peak structure in the corresponding temperature derivatives of the Polyakov-loop variables occur. This hampers a unique identification of the transition point. In order to clarify this behavior we plot the maximum of the peak location of the T -derivatives in the phase diagram together with a band around 80% of the maximum value. For larger chemical potential the peak locations deviate strongly and a coincidence of the chiral and deconfinement transitions can be excluded. This brings back the chirally symmetric and confined region which has been connected to the quarkyonic phase. In summary, we have shown that this signature is very sensitive to the correct implementation of the back-reaction of the matter sector to the glue sector, and is most likely not present.
At high chemical potential and small temperatures, a first-order chiral phase transition takes place which ends in a critical endpoint of second order located at µ CEP = 293 MeV, T CEP ∼ 32 MeV.

A. Thermodynamics
In general, the initial action given at the ultraviolet cutoff is independent of the temperature and chemical potential which limits the reliable calculation of thermodynamic quantities to a certain temperature and density region. Due to the ultraviolet cutoff high thermal and quantum fluctuations are suppressed. As a consequence, cutoff-independent predictions can be obtained for temperatures T Λ/8 which in our case yields an upper temperature bound of T ∼ 120 MeV. In order to cure this cutoff remnant at high temperature one has to combine the FRG with the expected perturbative results. The inclusion of the missing high-momentum modes can be achieved in an effective way by adding to the original flow, Eq. (21), a flow equation for an interacting Polyakov-loop quark system for scales k > Λ. Note, that an explicit gluon contribution to the flow equation is neglected here because the effective Polyakov-loop potential is fitted to reproduce the Stefan Boltzmann (SB) limit at high temperatures. In this way we employ the equation with E q = k 2 + m 2 q as previously defined. This equation is integrated from k = ∞ to k = Λ and yields the UV-contribution Ω ∞ Λ (T, µ), which is then added to the grand potential resulting from the solution of the RG flow equation (21). The divergent zero mode contribution in Eq. (28) is neglected here. However, for vanishing quark masses this term represents an unobservable shift in the grand potential.
In Fig. 7 the thermodynamic pressure normalized to the Stefan-Boltzmann pressure for three different quark chemical potentials is shown as a function of temperature. In the left panel of this figure, a fixed T 0 = 208 MeV has been used while in the right panel the µ-corrections are taken into account. The inset displays the pressure for µ = 290 MeV which is close to the critical endpoint in the phase diagram. For small temperatures the pressure decreases and has a kink at the critical temperature due to the first-order transition. Without the back-reaction of the matter fluctuations to the Yang-Mills sector a similar behavior in all three curves is observed. In the vicinity of the chiral transition the pressure increases due to the melting of the quark masses and saturates at about 80% of the corresponding ideal gas limit which reads for N f massless quarks and (N 2 c − 1) massless gluons (29) Including the back-reaction of the matter sector via the inclusion of T 0 (µ) changes the thermodynamics at larger chemical potential. The pressure increases much faster and saturates at µ = 290 MeV at 95% of the Stefan-Boltzmann limit. A similar trend is seen in the entropy, Fig. 8, and quark number density, Fig. 9, if the µ-corrections are taken into account. The entropy density decreases for small temperatures at µ = 290 MeV since the number of active degrees of freedom decreases when approaching the first-order transition from below. At the transition the entropy jumps. The bump around T ∼ 90 MeV (left panel) is a remnant from the smooth chiral crossover transition. This effect is completely washed out when the µ-corrections are included (right panel). Similar to the findings for the pressure these corrections become more significant at larger chemical potential. This also appears in the quark number density n q = −∂Ω/∂µ which is plotted in Fig. 9. For comparison the corresponding SB-limits (dashed lines) are also shown in this figure. The quark density approaches the SBlimit always from below. Without the µ-corrections the Polyakov loop suppresses the quark densities for chemical potential larger than the intersection point of the chiral and deconfinement transition in the phase diagram Fig. 6. With the T 0 (µ) corrections both transitions coin-cide over the whole phase diagram and as a consequence the quark number density approaches much faster the SB-limit (right panel of Fig. 9).
In Fig. 10 the scaled quark number density (left panel) and the corresponding scaled quark number susceptibility (right panel) for three different temperature slices around the critical endpoint (T CEP , T CEP ±5 MeV) as a function of the quark chemical potential are collected. In this figure the µ-corrections in T 0 are omitted while in Fig. 11 they are taken into account. Due to the chiral critical endpoint which is a second-order transition the susceptibility diverges with a certain power law [41]. There are no strong modifications in the structure of the susceptibility divergence if the back-reaction of the matter sector is taken into account or not. As a consequence it seems that the size of the critical region around the CEP is not strongly modified by these fluctuations. The only difference is that including the µ-corrections the peak height of the susceptibility is more pronounced towards the CEP.

VI. SUMMARY AND CONCLUSIONS
In the present work we have studied the Polyakovextended quark-meson model beyond mean-field approximation. The quark-meson fluctuations to the matter sector are included within a functional renormalization group approach. In turn, the quark-meson fluctuations to the gluonic sector are estimated by a perturbative computation as suggested in [6]. Interestingly, the latter procedure has been recently confirmed by a full dynamical QCD computation [9,10].
The validity range of the present model includes the temperature regime about the crossover temperature T c ≈ 200 MeV up to medium quark chemical potential of µ q ≈ 100−200 MeV. At low temperature and large chemical potential the model suffers from missing inclusion of baryonic degrees of freedom and further resonances.
However, the quantum fluctuations lead to strong modifications of the phase structure. Such modifications have already been observed in FRG studies of the quark-meson model. One of the prominent effects is a shrinking of the size of the critical region around the critical endpoint in the phase diagram [42]. Furthermore, quantum fluctua-tions push the potential critical endpoint towards lower temperatures and larger chemical potential, see, e.g., [27].
Our computation entails a small likelihood for a critical point with µ B /T ≈ 1 − 2 as predicted by some recent lattice studies [43]. Indeed, this is in accordance with the arguments of [44] as well as with recent developments concerning the convergence of the Taylor expansion about vanishing chemical potential [45]. The latter also complicates the extraction of the location of the critical point from the analysis of higher moments such as investigated in [46]. Thus, an extension of the present computations beyond mean field towards the low temperature high density regime will provide valuable information.
Further results concern thermodynamical quantities such as the pressure and the density. At vanishing density these quantities agree well with related lattice predictions and, in particular, show the correct behavior in the transition regime from mesonic degrees of freedom to the quark-gluon plasma regime.
In a forthcoming publication [37] we will present a full solution including novel algorithmic differentiation techniques [36] for the flow Eq. (21) which also allows us to discuss the systematics of the present approach towards QCD, as well as studying further thermodynamic quantities. We also plan to extend the present study to the 2 + 1 flavor case, as well as tightening the relation to full dynamical QCD.