Propagators and phase structure of Nf=2 and Nf=2+1 QCD

We investigate the phase structure of QCD at finite temperature and chemical potential by solving a coupled set of truncated Dyson-Schwinger equations for the quark and gluon propagator. In contrast to previous calculations we take into account the full back-reaction of the quarks onto the Yang-Mills sector and we include the effects of strange quarks. We discuss the resulting thermal mass of the unquenched gluon propagator and extract order parameters for the chiral and deconfinement transition from the quarks. Our result for the temperature dependence of the quark condensate at zero chemical potential agrees well with corresponding lattice calculations. We determine the phase diagram at finite chemical potential and find a potential critical endpoint at (mu_q,T) = (190,100) MeV.


I. INTRODUCTION
The existence of a high temperature and/or density phase with quarks and gluons as thermodynamically active degrees of freedom has been a major prediction since the early days of QCD. Over the years this topic has received a lot of attention from both, the theoretical and experimental sides. Today, high quality results are available from lattice simulations at zero chemical potential, see e.g. [1,2] and references therein. As it turned out, the low and high temperatures phases of QCD are not separated by a phase transition but rather continuously connected. By now, lattice QCD has firmly established the notion of a crossover at physical quark masses and zero chemical potential.
Unfortunately the situation is much less clear at (real) finite chemical potential, where lattice calculations are hampered by the notorious fermion sign problem. Although various extrapolation methods agree with each other at rather small chemical potential [3][4][5][6][7], regions in the (T, µ)-plane with µ q /T > 1 are hardly accessible. Therefore, despite these efforts the basic structure of the phase diagram of QCD is not at all settled yet, see e.g. [8,9] and references therein.
An alternative approach to the QCD phase diagram is functional methods, i.e. the functional renormalization group (FRG) and Dyson-Schwinger equations (DSEs). Naturally, these approaches have to rely on approximations in most cases 1 , which are less controlled than those in lattice Monte Carlo simulations. These approximations, however, are controlled by constraints such as symmetries and conservations laws and also by 1 There are, however, limits where exact results are possible, see e.g. Ref. [10].
comparison of the results with corresponding ones from lattice calculations at zero and imaginary chemical potential. If this procedure is followed carefully, one may hope to obtain meaningful results also for the finite chemical potential region with µ q /T > 1.
One of the potential advantages of the FRG and DSE-approach to QCD as compared with models like the Nambu-Jona-Lasinio model (NJL) [11], its Polyakov loop extended versions [12][13][14] and the Polyakov loop extended quark-meson model (PQM) [15][16][17] is the direct accessibility of the Yang-Mills sector. In the models, gluons are not active degrees of freedom and their reaction on the medium can thus neither be studied nor directly taken into account. This is possible within the functional approaches to QCD. With functional renormalization group equations the gluons of quenched QCD where studied in Ref. [18], and QCD at zero and imaginary chemical potential in Refs. [19,20]. With Dyson-Schwinger equations, lattice data for the quenched gluon where used as an input for the quark DSE to study chiral and deconfinement transitions as well as quark spectral functions in quenched QCD [21,22]. These calculations were generalized to the two-flavor case and to finite real chemical potential in Ref. [23].
In the latter study the quarks have been backcoupled to the Yang-Mills sector using a hard thermal loop (HTL) like approximation for the quark loop in the gluon DSE. In this work we will lift this restriction and consider the back-reaction with dressed quarks, as discussed in section II and the appendix. We will therefore be able for the first time to calculate the unquenched two plus one flavor quark and gluon propagators at finite temperature and density. We also determine the chiral and deconfinement transition regions from the quark condensate and the dressed Polyakov loop [21,[24][25][26]. Our results for the QCD phase diagram with N f = 2 and N f = 2 + 1 quark flavors are presented and discussed in section III. In the concluding section IV we summarize our results and discuss remaining uncertainties due to our truncation scheme.

II. IN-MEDIUM PROPAGATORS AND ORDER PARAMETERS
The objects of interest in this work are the propagators of the quark and the gluon, from which we will extract their thermal properties and order parameters for chiral symmetry breaking and confinement. In the Landau gauge and in the medium they are given by The vector dressing functions A and C and the scalar dressing function B of the quark propagator depend on the momentum p = (ω n , p) and, implicitly, on temperature and chemical potential. The same is true for the gluon dressing functions Z L and Z T with longitudinal and transversal orientation with respect to the heat bath. The corresponding projectors are given by The Matsubara modes ω n are subject to antiperiodic boundary conditions for fermions, ω n = πT (2n + 1), and periodic ones for the gauge boson, ω n = πT 2n. The quark dressing-functions will be determined self-consistently from the quark DSE, while the gluon dressing-functions will be given by a combination of lattice results and the gluon DSE as discussed in section II A below. From the quark propagator, Eq. (1), one can calculate the condensate where Z 2 is the quark wave function renormalization constant. For finite bare quark masses the resulting quantity is quadratically divergent and can be regularized according to where the divergent part mΛ 2 from the strange quark condensate ψ ψ s cancels the corresponding part of the light quark condensate ψ ψ l when multiplied with the ratio m l /m s of bare light to strange quark masses. The quark condensate is a strict order parameter for chiral symmetry breaking in the chiral limit and serves as an indicator for the chiral crossover at physical quark masses.
In [23] we presented a first calculation of the dressed Polyakov loop at finite real chemical potential in two-flavor QCD. This object is defined by [24][25][26] where ψ ψ ϕ is the quark condensate evaluated at generalized, U (1)-valued boundary conditions to act as order parameters for deconfinement it is mandatory to implement the generalized, U (1)valued boundary conditions only on the level of observables, but not in the partition function itself. All closed quark loops therefore maintain the physical value ϕ = π, whereas ϕ ∈ [0, 2π[ in the quark-DSE. This procedure breaks the Roberge-Weiss symmetry, a necessary condition for the dual condensate to act as an order parameter for center symmetry breaking [20,27].
In the presence of a crossover for chiral symmetry breaking and confinement, we determine the maxima of the derivative with respect to the quark mass of the order parameters to obtain the pseudocritical temperatures.
The main difficulty in the Dyson-Schwinger framework is to find a truncation scheme that correctly describes the relevant physics. In the case of QCD thermodynamics this means to cut off the infinite tower of DSEs in such a way that the temperature and density dependence of the n-point functions that are not determined self-consistently is carefully approximated. Fig. 1 shows the DSE for the quark propagator, which depends on the fully dressed gluon propagator and quark-gluon vertex. For the quark-gluon vertex we will use a carefully designed expression constructed along its Slavnov-Taylor identity and matched to provide the correct renormalization group running in the quark and gluon DSEs. For the gluon we will use temperature dependent lattice data for the quenched propagator, and unquench it by invoking the quark part of the gluon DSE. In particular this means that the gluon becomes sensitive to the chiral dynamics of the quark sector.  [28,29], and are mostly restricted to large and infinite temperatures. Another possibility to access the finite-temperature gluon is to use the functional renormalization group, where interesting first results have been reported recently [18]. Although the systematic errors are still considerable, probably the most reliable resource, however, is lattice QCD, where a number of quenched studies are available [22,28,[30][31][32]. These errors are most pronounced around the deconfinement transition temperature of the pure gauge theory, but fortunately have not much impact on our unquenched calculation [33].
Here, we will continue to use the lattice data from Ref. [22] as also done in our previous work [23].
In terms of diagrams, quenched QCD means that all quark loops are neglected, i.e. the last diagram of Fig. 2. One then recovers the first (second) order deconfinement transition of SU(3) (SU(2)) Yang-Mills theory at the correct transition temperatures. Using the dressed Polyakov loop this transition has been studied in Ref. [22], a corresponding study for the Polyakov loop potential using FRGs is reported in Ref. [34]. In order to investigate the phase diagram at physical quark masses, however, one needs to take the quark loop explicitly into account. In fact, as we will see a careful treatment of the back-reaction of the quarks onto the gluon sector is mandatory for a meaningful study of the phase diagram. Consequently, this also introduces an implicit chemical-potential dependence of the gluon dressing functions which is most pronounced in their thermal masses, see section III B.
Technically, we implement this dependence by substituting the Yang-Mills part of the gluon DSE by the quenched propagator from the lattice, and merely add the quark loop. The resulting equation is shown diagrammatically in Fig. 3, where we already distinguished between the (isospinsymmetric) light quark and the strange quark contribution. This procedure of merely adding the quark loops does neglect all quark loop effects inside the gluon or ghost loops and is therefore only approximate. However, we explicitly checked that this approximation works well in the vacuum [35], where differences to the fully self-consistent solution are below the five percent level. This may still be the case at finite temperature and chemical potential.
In [23] we furthermore approximated the quark loop by taking only bare quarks into account, which amounts to using a hard-thermal loop (HTL) like approximation with only the quarkgluon vertex being dressed (for a review on HTL results see e.g. [36]). This approximation, introduced for technical reasons, is much more severe and can only be justified for temperatures well above the critical one. In this work we therefore improve upon this situation and calculate the quark loop explicitly with the quark dressing functions obtained from the quark-DSE. By doing so, the quark and gluon DSE become coupled and have to be solved simultaneously. This means also, that the gluon becomes sensitive to the chiral dynamics in the quark sector and in particular to the chiral transition. The coupled system of DSEs in Figs. 1 and Fig. 3 can be written as where q = p − l, S f is the quark propagator for one flavor f = u, d, s, C F = 4 3 is the Casimir operator, Γ ν the dressed quark-gluon vertex, Z 1F and Z 2 are the vertex and wave function renormalization constants, and the Matsubara sum as well as the integration over the loop three-momentum l is represented by l = T n d 3 l (2π) 3 . For the momenta we often abbreviate generically p = ( p, ω p ) and p 2 = p 2 + ω 2 p . With the transversal and longitudinal projectors from Eqs. (3,4), the quark loop can be decomposed as The medium effects in the gluon-DSE manifest themselves predominantly in contributions to the thermal mass of the gluon. As a result, the quark loop can be split into a finite part that is similar to the vacuum quark loop and an infrared divergent thermal part that is proportional to 1/p 2 : with the infra-red regular part Π reg. and thermal masses A technical problem associated with the use of a sharp momentum cut-off Λ in the quark loop is the appearance of quadratic divergencies. This problem is particularly troublesome in the medium, since the Λ 2 term appears like a squared thermal mass and has to be carefully removed without spoiling the T 2 and µ 2 terms which carry the physics. We do this using the Brown-Pennington projection method [37], as detailed in Appendix A.
After removing the quadratic divergences and projecting on longitudinal and transversal parts the quark loop in the medium is given by with q = p − l and D q (p) = p 2 A 2 (p) +ω 2 p C 2 (p) + B 2 (p) for one quark flavor and with the abbreviationω n = ω n + iµ. The quark-gluon vertex has been split in three scalar functions, Γ, Γ s and Γ 4 which will be explained below. For the limit p → 0 we find the thermal masses. Of course, the expressions (10) together with (11), (12) also have the correct high-temperature limit, where the HTL results m T th.
Some technical details of the evaluation of this object are given in App. A. It is noteworthy that the quark loop does not contribute to the transversal thermal masses, in contrast to the Yang-Mills part of the gluon self-energy.

B. The quark-gluon vertex
An important quantity that appears in both, the quark and gluon DSEs, is the dressed quarkgluon vertex. At zero temperature the nonperturbative structure of this vertex has been explored to some degree in the past [38][39][40], however not much is known about the behavior at finite temperature and density. Therefore in order to close the system of DSEs we have to resort to an educated guess for the form of the quark-gluon interaction. Fortunately, there are constraints which serve as important guides in the construction of such a guess. On the one hand, at large momenta temperature and density effects are exponentially suppressed and one can entirely rely on the vacuum structure of the vertex. It is then well known, that the vertex has to combine with the gluon dressing functions such that the resulting combination of dressing functions runs like the running coupling in the ultraviolet. In addition, the vertex has to satisfy its Slavnov-Taylor identity. This identity has not been solved yet in QCD, however important guidance can be obtained from the corresponding Abelian Ward identity, which serves to express part of the vertex dressing in terms of the quark dressing functions. All these constraints have already been taken into account in previous works [21][22][23]. We therefore employ the same construction which uses the first term of the Ball-Chiu vertex, satisfying the Abelian Ward-Takahashi identity, multiplied with an infrared enhanced function Γ(p 2 , k 2 , q 2 ) that accounts for the non-Abelian dressing effects and the correct ultraviolet running of the vertex. The resulting expression reads where p and k are the fermion momenta, q is the gluon momentum and d 1 as well as the scales d 2 and Λ are parameters. Whereas d 2 and Λ control the renormalization group running of the vertex function from the large into the low momentum region, d 1 controls the strength of the quark-gluon interaction at small momenta and therefore the amount of quark mass generation in the hadronic phase. In the ultra-violet δ = −9 N c 44Nc−8N f is the anomalous dimension of the vertex and β 0 = 11Nc−2N f 3 . In Eqs. (11,12) we have used the abbreviations Γ s (p, q) = A(p)+A(q) 2 and Γ 4 (p, q) = C(p)+C(q) 2 . The squared momentum variable x is chosen to be q 2 in the quark DSE, but p 2 + k 2 in the quark loop. This change in the momentum dependence is well justified by the need to maintain multiplicative renormalizability in the gluon-DSE [35]. Note that, since the Ball-Chiu part is determined from the quark dressing functions our vertex does include effects from temperature and chemical potential.

C. Strange quarks and unquenching effects
The evaluation of the quark loop with dressed quarks allows us to add strange quarks in a nontrivial way. In general the light and strange quarks couple via unquenching effects in QCD, as shown in Fig. 3. In the vacuum, results for the propagators of N f = 2 + 1 QCD have been discussed in Ref. [35]. Here we present the first finite temperature and chemical potential study of the propagators in the unquenched theory. However, we also wish to emphasize that not all the unquenching effects are included yet. In addition to the ones in the gluon propagator there are also corresponding effects in the quark-gluon vertex. Diagrammatically, part of these can be written as hadronic contributions like pion and kaon exchange [39]. In Ref. [41] it has been argued that these effects are necessary to obtain critical scaling beyond mean field at the chiral phase transition temperature in the chiral limit. In our case of physical quark masses, however, we expect their influence on the critical temperatures to be small. Furthermore, since mesons have zero quark number the chemical potential dependence of the meson exchange is expected to be small as well. In this work we therefore do not include these contributions explicitly, but rather absorb their effects in our vertex ansatz (13). While this strategy may work rather well for meson contributions, it presumably fails to account for diquark and baryon contributions at large densities and low temperatures. This is the region of the phase diagram where nuclear effects are expected to play a major role [8] and results in our present truncation scheme should be interpreted with great care. We will come back to this point below.
In our calculation we only vary the chemical potential of the up/down quarks and set the strange quark chemical potential to µ s = 0, since the net number of strange quarks in heavy ion collisions vanishes. The light quark mass can be fixed by calculating the pion mass via the Gell-Mann-Oakes-Renner relation and the Pagel-Stoker approximation for the pion decay constant. This gives a mass of m u/d ≈ 2 MeV at a renormalization point of 80 GeV. We furthermore fix the strange quark mass by its ratio to the light quark mass by choosing ms m u/d = 27. For the vertex parameters we take d 1 = 7.5 GeV 2 , d 2 = 0.5 GeV 2 and Λ = 1.4 GeV. For the coupling constant we take α(µ) = 0.3. With this set of parameters and with two massless and one strange quark, the pion decay constant in the vacuum is f π ≈ 88 MeV and at the physical point of quark masses at vanishing density, T c ≈ 156 MeV. Those numbers agree very well with expectations, which is in contrast to simpler truncation schemes where usually the critical temperature comes out rather small [42].
Note again, however, that unquenching effects in the vertex are taken into account indirectly via the quark dressing functions A and C that appear explicitly in the vertex construction (13). We emphasize that it is nontrivial that this vertex construction leads to correct quenched [22] as well as unquenched transition temperatures as will be detailed in the next section.

A. Condensates, thermal gluon mass and
phase diagram for N f = 2 We first discuss the effect of coupling the twoflavor quark and gluon DSEs together as explained above. Here we will encounter some qualitative effects, which are still valid in the N f = 2 + 1 case discussed below.
The quark condensate and dressed Polyakov loop which are shown in Fig. 4 for two light quark flavors at vanishing chemical potential exhibit a crossover behavior around T c ≈ 200 MeV, with the transition region for chiral symmetry restoration and deconfinement being almost equal. We clearly identify the two most important effects from including the quark loop: a reduction of T c and a change from a first order phase transition to a crossover. Fig. 4 also shows the contribution from the quark loop to the electric part of the thermal gluon-mass as defined in Eqs. (10), (12). The thermal mass is small in the hadronic phase, shows a sharp rise in a temperature range around the transition temperature and converges to the T 2behavior expected from HTL-calculations in the high temperature phase. This behavior can be expected, since the quark loop has an inverse dependence on the quark mass-function. It is therefore suppressed by dynamical chiral symmetry breaking in the hadronic phase and becomes large in the high temperature phase, where the quark condensate becomes small. As compared to a truncation scheme where the quarks in the loop are not dressed, i.e. the HTL-like approximation we used in [23], the quark condensate shows a steeper transition around T c . This behavior originates from the growing thermal gluon mass due to the quark loop, which in turn decreases the interaction strength in the quark DSE. Therefore, the back coupling of the quark and gluon DSEs accelerates the chiral phase transition.
In Fig. 5 we present our results for the chiral and deconfinement transition lines determined in this work and compared to our previous results of Ref. [23], where we used the HTL approximation for the quark loop. Note that the deconfinement transition is determined solely from Σ −1 ; since the corresponding calculation of Σ +1 is plagued by considerable numerical uncertainties we refrain from showing the results. We first observe that the pseudo-critical temperatures for restoration of chiral symmetry and deconfinement are now closer together, in comparison to the HTL results where a difference of up to 10 MeV was observed. This closer matching of the critical temperatures can be explained by the gluon becoming sensitive to the chiral phase transition due to the explicit appearance of the quark mass functions in the full quark loop calculation. Moreover, the critical end-point has moved to a smaller quark chemical potential and a larger temperature. We attribute this to the accelerated chiral transition due to the improved back-coupling of the quarks onto the gluon, as explained above. The effect is substantial and underlines once more the general importance of details in the back-coupling effects from the quarks onto the Yang-Mills sector of QCD. As a result we obtain a ratio of critical quark chemical potential to critical temperature at the end point of µ q EP /T EP ≈ 1.1. Thus for N f = 2 QCD this ratio stays above one (although not much) in agreement with the claim made in Ref. [23]. We will see below that the addition of the strange quark further increases this ratio.

B. Condensates and phase diagram of
In addition to the two light quark flavors we now also add the heavier strange quark to the system. Recall, however, that following the conditions in a heavy ion collision we do not include a chemical potential for the strange quark in our calculations. We discuss the effect of this choice below. Let us first have a look at the resulting unquenched gluon propagator. The main effect of the temperature and chemical potential dependent quark loop onto the gluon is the change of the thermal electric mass. We plotted this quantity, normalized by its asymptotic value, as a function of temperature and up/down-quark chemical potential in Fig. 6. One clearly sees the continuous change for small chemical potential similar to the two-flavor case in Fig. 4. For larger chemical potential the transition becomes steeper until it becomes discontinuous in the vicinity of the critical endpoint of the chiral transition, discussed below. This behavior underlines again the fact that the chiral transition in the quark sector of the unquenched theory  also drives the changes in the Yang-Mills sector of QCD. In the asymptotic region at large temperature and/or chemical potential we finally recover the HTL-result m 2 th ∼ T 2 + 3µ 2 /π 2 for the thermal gluon mass.
Let us now focus on the quark sector of QCD and come back to the behavior of the quark condensate at zero chemical potential. Since we have two plus one quark flavors at our disposal we can directly compare with corresponding lattice data at vanishing chemical potential. In Fig. 7, we compare the renormalized quark condensate ∆ l,s defined in Eq. (6) with lattice data from Ref. [1]. The agreement is good except for the region above the phase transition, where the DSE results are larger than the lattice data. We attribute this discrepancy to the truncation of our quark-gluon vertex, which employs medium effects only in the Ball-Chiu part, while the ansatz function Γ is temperature independent. This is in contrast to [43], where it has been suggested that the parameter d 1 of the vertex infra-red part is reduced above T c . Therefore, neglecting the T -dependence of the vertex infra-red part leads to an over-estimation of the interaction strength in the chirally symmetric phase, and thus to too large condensates. Also the scalar parts in the vertex are neglected, which are known to be important when it comes to the details of dynamical chiral symmetry breaking and restoration [38]. Additionally, the (neglected) temperature dependence of explicit hadronic parts of the vertex may play a role [41]. Nevertheless, we find our present truncation sufficiently accurate in particular with respect to the steepness of the crossover to justify its extension into the finite chemical potential domain.
In Fig. 8 we show our results for the light and strange quark condensates in the temperature and chemical-potential plane. The condensates are not renormalized and shown on different scales for the sake of clarity. In the background of the figure we see again the crossover behavior of the light quark condensate at small values of the chemical potential. This turns into a first order transition at large chemical potential, visible in the foreground of the figure. The strange quark condensate shows a behavior similar to the light quark condensate around the phase transition, but continues to melt at larger temperatures. Certainly, this behavior is caused by the larger explicit breaking of chiral symmetry for the strange quark. At densities where the light quark shows a first order phase transition, the strange quark also shows a discontinuity at the same temperature. This is inherited from the discontinuity in the gluon propagator, see Fig. 6, and demonstrates how the gluon couples light and strange quarks.
In Fig. 9 we show the resulting phase diagram for 2 + 1 flavor QCD in our approximation scheme. As compared to the two-flavor case, Fig. 5, we find a chiral critical endpoint at slightly larger values of chemical potential. Thus, although the additional strange quark lowers the transition temperatures by roughly 45 MeV, the melting of the condensate is slightly less steep as for the N f = 2 case. As a consequence, one hits the first order transition region at larger values of the chemical potential. At µ = 0 we furthermore observe a slightly wider spread between the chiral and the deconfinement transition for the N f = 2 + 1 case as compared to N f = 2. This effect may be explained by the different strength of the back-reaction of the quarks onto the gluon sector. For N f = 2 + 1 the backreaction is stronger, resulting in a wider crossover transition region which in turn allows for a wider spread between the transition temperatures. On the other hand, the quenched calculation with no back-reaction is entirely driven by deconfinement and both transitions happen at the same temperatures [22].
So far, we set the strange quark chemical potential to zero. Setting it equal to the light quark chemical potential has very little effect on the results.
From the chiral transition line one can extract the curvature κ, which is defined by as the coefficient of a power expansion around µ = 0. This quantity can be measured on the lattice despite the fermion sign problem, and therefore serves as a comparison of the chemical potential dependence of the transition line. Two different lattice groups determined the curvature recently and found κ ≈ 0.059(2)(4) [5] and κ ≈ 0.059(20) [6] for 2+1 flavors. In comparison with the N f = 2 result κ ≈ 0.051(3) of Ref. [44], there is a general trend of increasing curvature when more flavors are taken into account in agreement with a large N canalysis [45]. In our case, we extract κ by fitting Eq. (14)  of the HTL-like quark loop. We also checked that the order (µ/T ) 4 is negligible in these intervals.
Our results for the curvature κ from the HTLtruncation of Ref. [23] and the two results of this work are summarized in Tab. I together with our findings for the location of the critical end-point (CEP). We find that the back-coupling increases the curvature compared to the HTL-like truncation. This is easily understood from the shape of the transition lines of Fig. 5, where the CEP from the fully back-coupled results is almost on top of the transition line of the HTL-approximated result, while at the same time the transition temperature at µ = 0 is significantly larger. In general, our values for κ agree in size with QM and PQM calculations [46,47] but are much larger than the lattice results. Part of this discrepancy may be attributed to volume effects in the lattice calculations [46], another part of it may be due to limitations of our present approximation scheme. This is particularly apparent in the N f -dependence of the curvature. As discussed above, the lattice results indicate a significant increase of curvature with N f , whereas our results indicate the opposite. The reason for this discrepancy may be found in the incomplete treatment of unquenching effects in our present approach. While we took into account quark loop contributions in the gluon polarization, the inclusion of corresponding effects in the quarkgluon vertex is numerically demanding and will be investigated in future work.

IV. CONCLUSION AND OUTLOOK
In this work we have determined the phase boundaries of N f = 2 + 1 QCD in the temperature and (up/down) quark chemical potential plane. We have solved a coupled set of Dyson-Schwinger equations for the Landau gauge quark and gluon propagators. For the gluon we worked with lattice results for the temperature dependence of the electric and magnetic part of the quenched propagator and added quark loop effects determined in the DSE-framework. For the quark-gluon interaction we worked with a construction which takes into account the correct renormalization group running of the vertex combined with leading order temperature and chemical potential effects as encoded in its Slavnov-Taylor identity.
From the temperature and chemical potential dependence of the quark propagator we extracted information on the chiral as well as the deconfinement transition encoded in the quark condensate and the dressed Polyakov loop. As a result we find a crossover behavior at small chemical potential in quantitative agreement with corresponding lattice results. We also find a critical endpoint at (T EP , µ EP q ) ≈ (130, 145) MeV followed by a first order transition. We therefore confirm our previous claim of the nonexistence of a CEP for µ EP q /T EP < 1 [23] in agreement with extrapolated results of lattice QCD [4,6].
The calculations presented in this work are, to our knowledge, the first effort based on continuum QCD to determine the phase diagram with realistic up/down and strange quark masses combined with a dynamical gluon sector. In this respect we wish to emphasize again the importance of the backcoupling effects of the quarks onto the Yang-Mills sector of the theory. These effects are responsible for the change of the order of the transition from first order in the quenched theory to the crossover behavior of full QCD and the corresponding reduction of the (pseudo-)critical temperatures. Both phenomena are nicely reproduced in our present truncation scheme.
Nevertheless, this is not the end of the story and from the present approximation scheme no firm conclusions on the precise location of the CEP nor on its very existence can be drawn. The reason is to be found in our present expression for the quarkgluon interaction. Although it does take into account the above-mentioned constraints from QCD itself, it lacks in detail when it comes to the precise dependence on temperature and chemical potential encoded in hadron back-reaction effects onto the quark propagator, as discussed in section II B. In the DSE-framework, these have only just begun to be explored [41]. In particular, baryon effects may be important when it comes to large chemical potential, as has been emphasized in [8] and explicitly checked in a two-color framework in Ref. [48]. Furthermore, at large chemical potential inhomogeneous phases may be energetically favored over homogeneous ones, see e.g. [49,50] and references therein. Whether these effects show up in addition to the existence of a CEP or instead of the CEP is an important open question that needs to be investigated further in the future.

V. ACKNOWLEDGEMENTS
We thank Michael Buballa, Jens Braun, Leo Fister, Daniel Müller, Jan Pawlowski, Bernd-Jochen Schaefer, Lorenz von Smekal and Jochen Wambach for fruitful discussions. This work has been supported by the Helmholtz Young Investigator Grant VH-NG-332 and the Helmholtz International Center for FAIR within the LOEWE program of the State of Hesse.