QCD phase structure from functional methods

We discuss the QCD phase structure at finite temperature and chemical potential for $2$-flavour and $2+1$-flavour QCD. The results are achieved by computing QCD correlation functions within a generalised functional approach that combines Dyson-Schwinger equations (DSE) and the functional renormalisation group (fRG). In this setup fRG precision data from arXiv:1706.06326 for the vacuum quark-gluon vertex and gluon propagator of $2$-flavour QCD are used as input, and the respective DSEs are expanded about this input. While the vacuum results for other correlation functions serve as a self-consistency check for functional approaches, the results at finite temperature and density are computed, for the first time, without the need of phenomenological infrared parameters.


I. INTRODUCTION
A detailed understanding of the QCD phase structure at finite temperature and density are essential for our understanding of the formation of matter and the evolution of the universe. It is also necessary for interpreting data collected at running and planned heavy-ion experiments as well as further predictions, for reviews see [2][3][4][5][6][7][8][9][10][11].
Functional investigations of the gluon sector have revealed a very interesting structure. While a quantitative access to the gluon propagator in Yang-Mills theory without any phenomenological input requires very elaborate approximations, see [20], its modification in the presence of dynamical quarks converge already in rough approximations to the matter fluctuations: the dominating effect of the matter fluctuations by far is the change in the momentum scale running in the ultraviolet which is already captured very well by perturbation theory. This property is also reflected in the fact that the gluon propagator does not change significantly if changing the pion mass from the chiral limit to masses of about 400 MeV, see [1]. Moreover, the light flavour quark-gluon vertex does not change significantly in the presence of a heavy strange quark: without flavour changing parts the presence of the strange quark in first order only changes the momen-tum running of the gluon propagator. Note that this approximate decoupling or rather separation of fluctuations changes in the chiral limit with a vanishing strange quark mass.
In combination this provides us with an optimised expansion scheme of DSEs for QCD, based on the quantitative fRG results in [1]: All correlation functions except the quark propagator are expanded about the 2-flavour case. This allows us to access correlation functions and the phase structure of QCD at baryon-chemical potential up to µ B /T < ∼ 3 with quantitative reliability, but the results extend beyond this regime. In the vacuum, at finite temperature, and at not too high baryon-chemical potential quantitative functional and lattice results serve as benchmark tests for our method.

II. FRG-ASSISTED DYSON-SCHWINGER EQUATIONS
Functional relations for QCD such as functional renormalisation group equations or Dyson-Schwinger equation are formulated as master equations for the effective action Γ[Φ] with the superfield Φ = (A µ , c,c, q,q). The derivatives of Γ with respect to the fields, Γ (n) Φi 1 ···Φi n (p 1 , ..., p n ) = are the one-particle-irreducible (1PI) correlation functions and are derived by taking field derivatives of the fRG or DSE for the effective action. This leads to oneloop (fRG) or two-loop (DSE) exact relations between full 1PI correlation functions of QCD. It is worth emphasising that the notion one-or twoloop has nothing to do with perturbation theory. The latter actually can be rederived by a recursive solution of the relations about the classical correlation functions. While being relations for the same set of correlation functions, in non-perturbative approximations they differ due to the different non-perturbative loop structure. Accordingly, an important respective self-consistency check is provided by agreeing results obtained in given approximations to different functional hierarchies. Moreover, their structure of closed one-or two loop relations between correlation functions allows to use correlation functions from external input such as other functional computations or lattice simulations, see e.g. [10, 28-31, 35-44, 58, 59]. By now, the lattice approach to Landau gauge QCD has produced quantitatively reliable vacuum ghost, gluon and quark propagators with small statistical errors as well as gluon propagators at finite temperature. Results for full dynamical QCD at finite temperature and physical quark masses are still lacking. Moreover, while impressive progress has been made in the computation of vertices, so far these computations lack the statistical precision required for a quantitative input.
In turn, functional computations in vacuum QCD have by now reached a quantitative level of precision not only for the propagators but also for vertices such as the quark-gluon vertex. In the present work we use the quark-gluon vertex and gluon propagator fRG data of vacuum 2-flavour QCD at physical quark masses from [1] as an external input in the DSEs for propagators and vertices in 2-and 2 + 1-flavour QCD at finite temperature and density.
We pursue slightly different strategies for the correlation functions involved. The quark propagator is computed from its DSE, the gap equation, see Figure 1 for all temperatures and densities. For the gluon propagator and the quark-gluon vertex we expand the respective DSEs about the fRG results of vacuum 2-flavour QCD from [1]: the DSE for the inverse gluon propagator Γ (2) AA (p) and the quark-gluon vertex Γ (3) qqA (p 1 , p 2 ) are reformulated as one for the difference between the vacuum two-flavour gluon propagator and quark-gluon vertex and the full gluon propagator and quark-gluon vertex. For the gluon propagator this DSE is solved, while we employ an ansatz for the strange, thermal and density fluctuations for the quark-gluon vertex. The latter is built on regularity relations between the different tensor structures in the vertex as well as the Slavnov-Taylor identities (STI).

A. Quark gap equation and chiral phase transition
The quark gap equation relates the inverse quark propagator Γ (2) qq to its classical counter part S (2) qq , to the quark and gluon propagators, and the classical and full quark-gluon vertex, see Figure 1. In the vacuum we write using the notation of [1]. In (2) the Dirac tensor structure is proportional to the quark wave function renormalisation Z q (p). Its RG-scale or momentum derivative is the quark anomalous dimension. In turn, the scalar dressing M q (p) is the quark mass function and carries the information of chiral symmetry breaking. Finite temperature and density singles out a rest frame and we parameterise, with the wave function renormalisation for the mode parallel to the rest frame, Z q . The overall factor Z q = Z ⊥ q is the wave function renormalisation for the transverse modes perpendicular to the rest frame. We have pulled out the transverse wave function renormalisation since the transverse modes carry more weight in the DSE loop integrals: within our construction of the thermal and density fluctuations of the quark-gluon vertex used in the DSE, see Section II C, we shall use a uniform wave function renormalisation Z q . This implies the approximation Z q /Z q ≈ 1 in the vertices. The quark momentump includes the chemical potential viap 0 = p 0 − i/3µ B . Eq. (3) straightforwardly extends to general correlation functions (1) withp ij with j = 1, ..., n that carry the chemical potential of the respective field Φ ij , for example, The factor s Φi is the baryon number of the component field Φ i of the superfield Φ = (A µ , c,c, q,q). We have s q = 1/3, sq = −1/3, and s Φi = 0 for Φ i = q,q. Note also that the correlation functions at finite density and/or finite temperature also carry a genuine dependence on the baryon-chemical potential. Eq. (3) signifies the mass function M q (p) as the 'pole mass' function of the quark as the full mass function is Z q M q . Alternatively, we could have defined the full mass function Z q M q . Then M q would be the 'screening mass' In (5) the quark and gluon propagators are given by G qq = (1/Γ (2) ) qq and G AA = (1/Γ (2) ) AA respectively, and all momenta are counted as incoming. At finite temperature the frequency integrals turn into Matsubara sums with A µ , c,c : ω n = 2πT n , q,q : ω n = 2πT n + 1 2 Gluon DSE for the difference ∆Γ (2) AA between the full gluon propagator and the vacuum 2-flavour gluon propagator. Lines and vertices with black blobs are the full propagators and vertices at finite temperature and density. Lines and vertices with grey blobs are full vacuum propagators and vertices for N f = 2. The square bracket contains the temperature and density fluctuations and does not require renormalisation. We have dropped the part of the classical gluon propagator that carries the renormalisation of the strange loop.
Eq. (5) can be solved with the knowledge of the gluon propagator G A and the quark-gluon vertex Γ (3) qqA . In Section II B and Section II C we detail how these correlation functions at finite temperature and density are computed on the basis of the vacuum two flavour fRG-data from [1].
B. DSE for the strange quark, temperature, and density corrections to the gluon propagator The DSE for the inverse gluon propagator for general flavours N f at finite temperature and density is expanded about that in the vacuum for two-flavour QCD. This reads schematically and leaves us with a DSE for the difference ∆Γ (2) AA between the full inverse gluon propagator and that of vacuum 2-flavour QCD, depicted in Figure 2. At finite temperature and density the gluon two-point function has color-electric and color-magnetic components. Then, the two-point function Γ where we have separated the gauge-fixing part p µ p ν /ξ. In (8) we introduced the color-magnetic wave function renormalisation Z M A and the color-electric one, Z E A . The projection operators onto the color-electric and colormagnetic directions are given by In the present work we only consider the Landau gauge, ξ = 0. The loop integrals in the DSE are dominated by the color-magnetic components and in the following we will use the O(4)-symmetric approximation Z E A (p) = Z M A (p). In the quark gap equation this approximation is reflected in the identification of the color-magnetic and electric running couplings α E s (p) = α M s (p). Now we discuss the advantages of the split (7) as well as generic strategies for its numerical solution. Since the propagator shows a subleading dependence on both the additional quark flavour (for N f = 2 + 1) as well as temperature and density, the DSE for ∆Γ (2) AA can e.g. be iterated starting with ∆Γ (2) AA = 0. The subdominant behaviour of changes in the quark look has been evaluated in [1] with a study of the missing current-quark mass dependence in two-flavour QCD: the strange quark loop is quantitatively captured with its induced perturbative momentum and RG-scale running of the gluon propagator. The latter behaviour is already captured in simple approximations. This argument also carries over to the mild temperature and in particular very mild chemical potential dependence, see e.g. [10,23,35]. These studies also provide further numerical justification for the approximation Z E A (p) = Z M A (p). Note also that the expansion in (7) reduces the renormalisation subtleties in the vacuum. Moreover, if iterating about the respective vacuum result for N f = 2 and N f = 2 + 1, the temperature-and density independence of the renormalisation procedure is apparent in the finiteness of ∆Γ (2) AA . A further numerical stabilisation is achieved by separating the thermal sums and integrals in ∆Γ (2) AA as follows with q 0 = ω n at T = 0 and q 0 = ω at T = 0, and loop(q, p) stands for the loops in the second line of Figure  2. Then, the second line in (10) encodes the direct effect of thermal and density fluctuations, namely substituting the q 0 -integral by the Matsubara sum as well as introducing the chemical potential in the quark loop while keeping the vacuum correlation functions. Evidently this is finite and numerically stable. Moreover, it can be computed from the knowledge of the vacuum correlation functions only. The third line vanishes for ∆Γ (2) AA = 0, and the difference of loops in the square bracket decays faster than loop(q, p) itself. This results in a numerically stable second line which is only triggered by the direct effects of thermal and density fluctuations in the third line. Moreover, it is subleading for sufficiently small T, µ B .
In the present work we utilise the quantitative smallness of the correction ∆Γ (2) AA , and resolve it in a simplified approximation: First of all the temperature dependence of ghost-gluon correlation function has been shown to be subdominant for the temperatures considered, both with functional methods and on the lattice , e.g. [60,61]. The dependence on the chemical potential is even more suppressed as it enters the DSEs for ghost-gluon correlations only via the chemical-potential dependence of the quark loops in gluonic correlations. Consequently, we approximate ghost propagator and the ghost-gluon vertex at finite temperature and density on the right hand side of (10) with their vacuum counterparts in [1].
Moreover, we have approximated the vacuum threegluon vertex with the classical one. In this approximation ghost-gluon contributions to the third line in (10) vanish. This approximation is trustworthy as long as ∆Γ (2) AA remains a small perturbation. We have monitored this property in our explicit computation. Due to the finiteness of the momentum integrals this approximation is easily lifted which will be considered elsewhere.

C. Quark-gluon vertex at finite T and µB
In this Section we derive an Ansatz for the strange quark, thermal and density corrections ∆Γ (3) qqA of the quark-gluon vertex for 2-and 2 + 1-flavour QCD based on the N f = 2 flavour quark-gluon vertex in the vacuum computed in [1], see also [42,44,62,63]. The strange-quark vertex in the vacuum is identified with the two-flavour vertex: for perturbative momentum scales with p 2 /m 2 s 1 the mass-dependence is sub-leading and for infrared momenta the symmetry-breaking part of the mass function dominates. Hence, similarly to the parameterisation of the gluon two-point function (5) we expand the vertex about the two-flavour vertex in the vacuum, for both the light quarks, q = l, and the strange quark, q = s. Note that in the present work we assume isospin symmetry with m l = m u = m d , see (30) for twoflavours and (31) for 2+1 flavours. At finite temperature and density the vertex has color-electric and magnetic components. We identify them within our O(4)symmetric approximation, which is consistent with the O(4)-symmetric approximation of the gluon two-point function.
Eq. (11) leaves us with the task of solving the Dyson-Schwinger equation for the difference of flavours and temperature/chemical potential, depicted in Figure 3. Evidently, for small corrections ∆Γ In the present work we refrain from computing the corrections ∆Γ (3) qqA , their computation will be presented elsewhere. Here we follow a strategy that reduces the computational costs significantly: We construct vertices by utilising firstly non-trivial relations between the dressings of the different tensor structures of the quark-gluon vertex, secondly its regularity at momenta p 2 > ∼ 1 GeV, and thirdly the STIs for the dressings of the different tensor structures. Apart from the reduced numerical costs this construction guarantees the quantitative compatibility of the running of the different 'avatars' of the strong couplings derived from the purely gluonic, ghost-gluon and quark-gluon vertices, see [1]: small deviations in the perturbative and semi-perturbative momentum runnings of the respective couplings have a huge effect on the chiral symmetry breaking scale.
Within this construction the vertex dressings of all tensor structures of the quark-gluon vertex can be expressed in terms of quark, gluon and ghost dressings and the dressing of the ghost-gluon vertex. This allows us to make use of the hierarchy for the latter dressings in terms of their T, µ B and s-quark dependence: While ghostgluon correlations are effectively independent of T, µ B and s-quark fluctuations, pure gluon correlations show a mild T dependence. The dependence of pure gluon correlations on the strange current-quark mass m 0 s is negligible. This argument extends to the µ B -dependence of pure gluon correlation, which is also subleading. Consequently, only the quark dressings show a significant T, µ B -dependence, that cannot be neglected in the computation of the quark-gluon vertex. The above hierarchy also entails, that to leading order differences between the light-quark-gluon vertex, Γ llA , and the strange-quarkgluon vertex Γ (3) ssA are solely induced by differences in the quark dressings. This supports the self-consistent and quantitative nature of the present approximation scheme.
The idea behind this construction can be nicely illustrated by discussing the quark-gluon coupling α qqA (p) and the ghost-gluon coupling α ccA (p) given by Here, λ qqA , λ (1) ccA are the dressings of the classical tensor structures of the respective vertices. For the classical action we have λ ccA . Furthermore, we identify the wave function renormalisations of all quark flavours, Z q = Z l = Z s due to their very mild dependence on the quark mass. This implies that in the present approximation differences in the vertex dressing for Γ (3) llA and Γ (3) ssA are only triggered by the quark mass functions consistent with our discussion above. In (12) the vertices are evaluated at the symmetric point where p is the incoming quark (ghost) momentum, q is the incoming anti-quark (anti-ghost) momentum and k + is the incoming gluon momentum. Moreover, Z A (p) and Z c (p) are gluon and ghost dressing function, for more details on the notation see [1]. The two couplings in (12) are related with STIs and two-loop universality. They agree to a quantitative degree due to rather trivial quark-gluon scattering kernels for perturbative and semi-perturbative momentap 2 > ∼ 3 GeV. There we have see Figure 14 in Appendix A. Eq. (14) is a rather interesting relation: the ghost-gluon dressing λc cA and the ghost wave function renormalisation Z c (p) are effectively independent of the strange quark, temperature and in particular of density fluctuations. For T and µ this reads with the strange current-quark mass m 0 s . The very mild temperature and current-quark mass dependence has been checked explicitly with functional methods and lattice simulations. In turn, the chemical potential dependence of ghost correlation functions is only triggered very indirectly: The explicit chemical potential dependence of the quark propagator triggers a very mild one in gluon correlations which then feeds into the ghost correlations.
In summary, for the temperatures and chemical potentials relevant for the phase structure we can rely on vacuum ghost-gluon correlations. As an important consequence the whole strange-quark, temperature and chemical potential dependence of the quark-gluon dressing λ (1) qqA of the classical tensor structure in the perturbative and semi-perturbative regime stems from the factor Z q (p) in (14). Using (11) for the dressing λ (1) qqA leads us to the general parametrisation where the superscript λ (in) indicates the input data, here the 2-flavour QCD dressing in the vacuum. With (14) the correction ∆λ (1) qqA is given by withλ Hence, the strange quark, thermal and density dependence of λ (1) qqA is entirely carried by Z q . Note that λ (1) qqA carries a non-trivial momentum dependence also for p > ∼ 1 GeV as we refrained from dividing out the ratio λ (1) ccA /Z c . As discussed above, this ratio has a negligible s, T, µ B -dependence, see (15), and in the present approximation we simply use the vacuum values. Then it is a global prefactor in the difference and can be absorbed iñ λ (1) qqA as done in (17). The above derivation sets the stage for the discussion of the full quark-gluon vertex. We discuss regularity and STI-constraints for the different dressings, which finally provide us with a vertex ansatz. To that end we expand the full quark-gluon vertex in a full tensor basis which contains twelve elements. In [1] the full transversely projected quark-gluon vertex has been written as follows, with k + = p + q, see (13), and the transverse and longitudinal projection operators The tensor basis T i qqA in (19) is given by λ (1) qqA is the dressing of the classical tensor structure. The basis (21) can be derived from quark-anti-quark and gluon derivatives of the gauge invariant operators qqA , see [1,15]. The tensors T In [1,15] is has been also observed that for perturbative and semi-perturbative momenta the different dressings λ (i) qqA of those tensor structures obtained from the same operator are indeed related at the symmetric point. Moreover, in [1,15] it has been found that the dressings λ (i) qqA only show a mild angular dependence. For similar considerations in a basis with defined C-parity see [62][63][64][65]. This allows us to employ angular averages of the input dressings, qqA , while λ (3,8) qqA ≈ 0. Moreover, for momentap > ∼ 1 GeV we have which leaves us with the task of writing ∆λ (1,3,4,7,8) qqA in terms of quark, gluon and ghost dressings and the dressing of the ghost-gluon vertex with the help of the STIs.
The validity of the latter procedure is checked in the vacuum, see Appendix A. We are led to ∆λ (1,5,6,7) qqA (p, q) =λ (1,5,6,7) qqA In (24) and (25) we have used the abbreviations for dressings factors that originate from the quark legs of the vertex. The dressingsλ ccA , are depicted in Figure 15 and Figure 16 in Appendix A at the symmetric point. Within our approximation with identical dressings for all quarks, Z q = Z l = Z s , the quark-gluon vertex dressings λ (1,5,6,7) are flavour-independent while λ (2,3,4,8) carry a -mildflavour-dependence due to ∆ ZqMq in ∆λ (2,3,4,8) , see (24). This concludes our construction of an STI-compatible quark-gluon vertex at finite temperature and density on the basis of the 2-flavour quark-gluon vertex in the vacuum. We emphasise that the construction is general, and is easily adapted to generic input data not only from the fRG but also from other functional approaches as well as the lattice.
with the pion decay constant as a first prediction. Also, the quark mass function M q (p) can be compared with the respective fRG and lattice results. We emphasise that no phenomenological infrared parameter such as the infrared vertex strength is adjusted, the current setup has the full predictive power of a first principle approach.
with the normalisation N π of the Bethe-Salpeter wave function of the pion, where Z q (p 2 ) = ∂ p 2 Z q (p 2 ) and Z q (p 2 ) = ∂ 2 The relations in (28) are reductions of the exact relations with additional form factor that are set to unity in (28). In particular, the approximation (28) leads to an underestimation of the pion decay constant. In summary, with (28) we can compute the pion mass and decay constant from the quark , for more details see e.g. [1,23,67]. The double lines stand for scalarpseudoscalar degrees of freedom introduced with dynamical hadronisation, [1,15,16,23,[67][68][69][70].
propagator. This allows us to determine the respective current-quark masses m 0 u/d = m 0 l . The pion mass m π and decay constant f π can be computed from the quark propagator, see e.g. [71,72], for reviews see [65,73,74]. We use the Pagels-Stoker approximation for the pion decay constant, and the Gell-Mann-Oakes-Renner relation for the pion mass, This concludes our setup for fRG-assisted DSE computations described in Section II and the present section Section II A. Now we determine the 2-flavour currentquark masses with (27), using (28). This leads us to the quark propagator depicted in Figure 5 and With N π = 88 MeV the difference between the pion decay constant and the normalisation of the pion wave function is rather small: (f π − N π )/f π = 0.043. For large differences the approximations used in (28) lack reliability, while a small difference serves as a self-consistency check.
An important consistency check is given by the momentum dependence of the quark mass function M q (p) and, to a lesser extent, by the quark wave function renormalisation Z q (p). Note that the latter has a relatively strong dependence on the renormalisation scheme. The result is shown in Figure 5. Our result is consistent with the lattice data and the quark propagator from the fRG computation. In particular we find that the constituent quark mass at vanishing momentum, M (p 2 = 0) = 387 MeV in Figure 5, is slightly smaller than that in the fRG computation, and is rather close to the result from lattice QCD.
In this context we emphasise that the resummation schemes for the DSE and the fRG are working rather differently. The former equation is depicted in Figure 1, the latter one is depicted in Figure 6, for more details see [1]. In comparison, the fRG equation consists of more diagrams as well as a different vertex structure: First of all, all vertices are dressed and in particular the flow diagram with quark-gluon vertices has two dressed quarkgluon vertices. Moreover, the flow equation for the quark  [1]. We also display 2+1-flavour vacuum results from the fRG, [23], and the lattice, [75], as well as fRG results at T = 1.3 Tc from [23].
two-point function contains tadpole terms with the full two-quark-two-gluon scattering vertex, the two-quarktwo-ghost scattering vertex and the quark four-point vertex that are absent in the respective DSE. In particular the latter tadpole is important for the strength of chiral symmetry breaking. It carries in particular resonant meson exchanges of the σ-mode and pions. In summary, the small fRG-DSE deviations are within the expected systematic error of the approximations used both in the fRG and the DSE, which is a respective non-trivial consistency check for both functional approaches.
We proceed with the 2 + 1-flavour case. The results for the vacuum gluon propagator and quark propagator are depicted in Figure 7 and Figure 8 together with finite temperature results that are discussed in the next Section IV. As in the 2-flavour case we determine the current-quark masses with (27), using (28), With N π = 86 MeV the difference between the pion decay constant and the normalisation of the pion wave function is rather small: (f π − N π )/f π = 0.034. As in the 2flavour case this serves as a self-consistency check of the approximations in (28). The small reduction of the pion decay constant in comparison to the full one with 93 MeV is a well-known artefact of the approximations in (28). This completes the setup. We emphasise that in the current approach there is no phenomenological infrared parameter, the only input are the fundamental mass parameters of QCD: the light current-quark masses fixed with the pion pole mass, and, in the case of 2+1 flavours, the ratio of the light and strange current-quark masses. Accordingly, the quark propagator and in particular the pion decay constant f π , the mass function M q (p) and the constituent quark mass M q (0) are already predictions within the current setup. Note also that we could have fixed the physical point with the pion decay constant, then the pion mass would have been a prediction. In summary this amounts to a fully first principle setup, we only have to determine the fundamental parameters and all observables are predictions in the current approach.

IV. RESULTS AND DISCUSSIONS
In this Section we present our results for 2-and 2 + 1flavour QCD at vanishing and finite temperature and density. We depict our results for the physical case of 2 + 1-flavour QCD, the 2-flavour results show a similar behaviour. In Figure 7 and Figure 8 we show the 2+1flavour gluon and quark propagators in the vacuum and at finite temperature. Again the results are in quantitative agreement with the respective lattice results in the vacuum (gluon: [75,[77][78][79], quark: [80,81]) and the functional results at vanishing and finite temperature (gluon and quark: fRG [23], DSE [10,30]). At finite temperature the gapping of the gluon propagator increases, consistent with the functional results in [10,23,30], reflecting the increase of the thermal screening masses.
These results readily extend to finite baryon-chemical potential. We are interested in particular in the chiral phase structure at finite temperature and density. In the present work we employ the definition of the chiral transition temperature T c (µ B ) given by the thermal susceptibility of the renormalised light chiral condensate ∆ l,R . Up to renormalisation terms the quark condensate ∆ qi is The lattice data at vanishing chemical potential are from [82]. The crossover steepens, finally leading to 2nd order CEP and a first order regime. A comparison with the fRG results in [23] is done in Appendix B.
defined as with q i = u, d, s. The renormalised light chiral condensate comprises the thermal and density part of the chiral condensate. In particular, the renormalised light chiral condensate is given by Where N R is a convenient normalization which leads ∆ l,R dimensionless, and here we choose N R = m 4 π . The subtraction of the vacuum condensate eliminates the necessity of explicitly discussing the renormalisation of (32). In Figure 9 we show the renormalised light chiral condensate obtained from the current computation in comparison to lattice data at vanishing chemical potential, [82]. Figure 9 also contains results for different baryonchemical potentials. The peak of the thermal susceptibility, ∂ T ∆ l,R , provides the chiral transition temperature. The susceptibilities are shown in Figure 10 for different baryon-chemical potentials. The crossover temperatures for 2-and 2 + 1-flavour QCD at vanishing chemical potential are computed as The chiral transition temperatures in (34) are in quantitative agreement with the lattice, [50,82,83], and fRG results, [23]. In Appendix B we also compare the temperature dependence of the renormalised light chiral condensate computed in the present work and that from [23]  for different chemical potential. For small chemical potential they are in quantitative agreement. In turn, they show significant differences for µ B /T > ∼ 3, triggered by the different locations of the CEP in both computations, see Figure 13. Now we proceed to the full phase structure at finite density. Before we present the respective results, we briefly discuss the self-consistency of our present approximation at finite baryon-chemical potential. At vanishing temperature and below the baryonic onset µ B < µ * B , no correlation function has an explicit dependence on µ B , Γ (n) Φ1···Φn (p 1 , ..., p n ; µ B ) = Γ  µ B -dependence of a correlation function is carried by the frequency arguments. This is the Silver Blaze property for correlation functions as discussed in [35,84]. In the present approximation we only consider the explicit µ Bdependence in γ 0p0 and ignore the above subtleties in the dressing functions M q (p 2 ) ≈ M q (p 2 ) and Z q (p 2 ) ≈ Z q (p 2 ). In Figure 11 we show the mass function M q (p 2 = 0), that is atp 2 = −µ 2 B /9. For T = 0 the µ B -dependence is rather small and simply reflects the mild momentum dependence of the quark propagator for 0 ≤ p < ∼ 0.5 GeV for momenta p 2 < 0. Accordingly, the violation of the Silver Blaze property is very small.
This leads us to one of our main results, the phase structure of 2-and 2+1-flavour QCD summarised in Figure 12. The transition temperature T c (µ B ) in Figure 12 decreases with increasing baryon-chemical potential. The curvature of the transition temperature line T c (µ B ) at vanishing chemical potential is determined within the expansion about µ B = 0, with T c = T c (µ B = 0). From the results for the thermal susceptibility we are led to in the two flavour case, and Both results are compatible with the respective fRG results in [23]. While the 2-flavour lattice results show reference curvature κ N f = 2 N f = 2 + 1 fRG-DSE: this work 0.0179(8) 0.0150 (7) fRG: [23] 0.0176(1) 0.0142 (2) lattice: [57] 0.0153 (18) lattice: [52] 0.0144 (26) lattice: [54] 0.015 (4) DSE: [32] 0.038 DSE: [10,30,85] 0.0456 0.0238 fRG: [25] 0.0051 lattice: [86] 0.0078 (39) lattice: [87] 0.0056(6) a wide spread, the 2 + 1-flavour lattice results for the curvature are in quantitative agreement with the present results. A comparison of all results is provided in Table I. We now move to the discussion of the phase structure at large baryon-chemical potential. This requires a discussion of the reliability bounds of the current approximation: the quark gap equation depends on the gluon propagator and the quark-gluon vertex. As expected, the former correlation function shows little dependence on the baryon-chemical potential in agreement with respective results from the fRG, [23] and other DSE studies [10,35]. Accordingly, the respective systematic error is supposedly small. In turn, the quark-gluon vertex inherits a strong dependence on the baryon-chemical potential related to that of the quark propagator via the STI. Moreover, strongly resonant baryon and meson channels are potentially not well-captured by the present STIconstruction of the vertex. Based on this analysis and the results in [23] on non-trivial momentum dependencies for µ B /T > ∼ 3 we conclude that the present approximation has to be upgraded towards solving the DSE for ∆Γ For 2+1-flavour QCD we find a critical end point with Both CEPs are at large ratios µ B /T , As already mentioned above, in this regime the current approximation has lost its quantitative reliability. Its systematic improvement is under way: we solve the DSEs for the vertex dressings ∆λ (i) with i = 1 − 8 in ∆Γ (3) qqA . This analysis also includes scalar-pseudoscalar meson and baryon exchange channel, see [31]. Note that the current fRG-DSE setup allows to utilise fRG results on the scalar, pseudoscalar and density four-quark channels as well as higher order scatterings of resonant channels. We hope to report on respective results soon.
We close this section with a comparison of our main result, the phase structure of N F = 2 + 1 flavour QCD with other theoretical results as well as freeze-out data from different groups, see Figure 13. We would like to emphasise again that the current computation relies only on the determination of the fundamental parameters of QCD. The light current-quark masses are fixed with the pion mass, and the strange current-quark mass is fixed with the ratio of m 0 s /m 0 l ≈ 27. All correlation functions, and the renormalised light chiral condensate as well as the critical temperatures T c (µ B ) are predictions in first principle functional QCD and compare well to the fRG results and the lattice results at small densities. In summary, the current combined functional setup and respective DSE and fRG studies provide a promising starting point for getting to quantitative predictions in the large density regime including the potential CEP in the next few years.

V. SUMMARY
In this work we have evaluated the QCD phase structure at finite temperature and chemical potential with a combination of the functional renormalisation group (fRG) and Dyson-Schwinger (DSE) equations set-up in Section II: fRG results from [1] for the 2-flavour gluon propagator and quark-gluon vertex in the vacuum have been used as an input for the DSEs of 2-and 2+1-flavour QCD at finite temperature and density. We have solved the coupled set of DSEs for the quark propagator and for the thermal and density corrections of the gluon propagator for 2-and 2+1-flavour QCD. In the latter case we also have computed the strange quark contributions for all temperatures and densities. The respective corrections of the quark-gluon vertex have been deduced from regularity properties of the vertex for momenta p > ∼ 1 GeV, relations between the dressings of the different tensor structures, and the Slavnov-Taylor identities for the dressing  functions, see Section II C. No phenomenological IR-scale has to be tuned and the only parameters are the fundamental ones in QCD, the current-quark masses, see Section III. The latter are determined with the pion pole mass m π = 140 MeV, and, in the 2+1 flavour case, with the ratio of light and strange current-quark masses. The approach has been benchmarked in the vacuum and at finite temperature: the results for gluon propagator, quark propagator and the renormalised light chiral condensate are consistent with lattice results, e.g. [75,76,94], fRG results, [23], and DSE results, e.g. [10] and references therein. The respective comparisons are depicted in Figures 5,7,8,9. After these successful benchmarks the setup has been used for the computation of the phase structure for 2-and 2 + 1-flavour QCD at finite temperature and density, see Figure 12. For 2 + 1-flavour QCD the results have been compared with other theoretical results as well as freeze out data, see Figure 13 and the discussion there. At small chemical potential the present result agrees quantitatively with lattice and fRG results. For µ B /T < ∼ 3, the present approach can be readily used for the computation of observables such as net baryon number fluctuations relevant for the access to the freeze-out curve as well as for the physics of the QCD epoch in the early universe, see e.g. [95].
For larger chemical potential, µ B /T > ∼ 3, lattice simulations are obstructed by the sign problem. In this regime the present fRG-assisted DSE results agree quantitatively with the fRG results in [23]. This non-trivial agreement enhances the reliability of the respective results, given the very different resummation schemes of DSE and fRG used in particular in the matter sector. Still, the approximation used in FRG computation [23] and the present approximation lack quantitative reliability in this regime: the current construction of the thermal and density corrections of the quark-gluon vertex can only be trusted quantitatively as long as it still a correction. This calls for a self-consistent computation of the thermal and density fluctuations that contribute to the quark-gluon vertex in particular in the presence of resonant interactions triggered by the finite density as well as scaling phenomenon in the vicinity of a potential critical end point. To a smaller degree this reliability analysis also applies to the present linear feedback of temperature and density corrections in the gluon DSE. With the important caveats in mind, the current computation shows a critical endpoint at (T CEP , µ B CEP ) = (93, 672) MeV, see (40).
We emphasise that the approximations used here do not pose conceptual or computational challenges and can be resolved in an extension of the current work. In order to consolidate the large density results, the approximation is currently systematically improved, as well that used in the respective fRG computation from [23]. These combined works should finally allow for a quantitative access to the large density regime including the position of the potential CEP.
The full quark-gluon vertex can be expanded in a tensor basis {T i qqA } with i = 1, ..., 12. A tensor basis for the transverse part of the vertex, (19), has been obtained from the (21) with i = 1, ..., 8 via transverse projections. The four remaining tensors are obtained from the longitudinal projections of these tensors. These projections only contain four linearly independent tensors, and we choose T (9,10,11,12) qqA µ (p, q) = Π µν (k + ) T (1,2,6,8) qqA ν (p, q) . (A1) The dressing functions of these tensor structures can be determined via the STI for the quark-gluon vertex, which relates the dressings to the dressings of the quark, gluon and ghost propagators as well as scattering kernels. A very detailed analysis of the respective STI-relation for the longitudinal dressings λ (9−12) qqA is given in [1], Appendix D.
Regularity for perturbative and semi-perturbative momenta relates in particular λ (1) qqA and λ (9) qqA with For deriving (A2) we first remark that λ (1) qqA is the dressing of the transverse projection of T 1 qqA , while λ (9) qqA is the dressing of its longitudinal projection. The full relation also involves λ (5,7) qqA and can be found in [1], Appendix D. Then (A2) follows with a further relation, The cancellation between λ (5) qqA and λ (7) qqA is suggested by the relation of the respective tensor structures: both can be derived from the gauge invariant quark-gluon interaction term,q related toqD / 3 q, see [1,15]. There, (A3) has also been checked numerically. For the last of these tensor structures related toqD / 3 q it has been found numerically, Note that such relations also hold for the dressings of the tensor structures λ (2,4) qqA derived fromqD / 2 q, λ (2) qqA (p, q) 2 ≈ λ This closes the discussion of relations between dressing functions from gauge invariant quark-gluon interaction terms and regularity.

STI-construction
Now we utilise the STIs for the dressings λ  qqA also for the transverse dressing λ (1) qqA . This leads us tõ forp > ∼ 1 GeV. Here, λ ccA is the dressing of the transverse part of the ghost-gluon vertex, and the differences ∆ X and sums Σ X of the quark dressings defined in (26) carry the momentum and RG-behaviour of the quark two-point functions. Moreover, in the present work the superscript (in) indicate the vacuum N f = 2-flavour data. The dressingλ (1) qqA (p, q) is nothing but the full scattering kernel in the STI (up to the convenient rescaling with Z c (k + )).
The ratios in (A7) are approximately constant for p > ∼ 1 GeV in line with regularity. In turn, forp < ∼ 1 GeV regularity fails, see [1], which is reflected in a non-trivial momentum dependence of the ratios, see Figure 14. Importantly, the STI links the dressing λ (1) qqA of the classical tensor structure to ghost-gluon correlation functions and the quark wave function renormalisation. It is wellknown that Z c and λ (1) ccA are nearly insensitive to finite temperature and small chemical potentials, as well as the strange quark.
Note also that the non-regularity of vertices is in oneto-one correspondence to the occurrence of a gluon mass gap which signals confinement, for a detailed discussion see [1,20]. Accordingly, the irregularity in the quarkgluon vertex is triggered by the necessary one in the glue system. In summary this suggests the following parametrisation of the vertex: the non-trivial deviation from the STI-construction based on the regularity assumption is triggered by the gluon sector and hence is well-captured by the N f = 2 flavour vertex. In turn, the strange quark and finite temperature and density fluctuations are regular, potential corrections to the irregularities only come indirectly via the gluon correlation functions in the diagrams for strange, T and µ B fluctuations. These subleading effects are discarded here, but are evaluated in a forthcoming publication.
Importantly, under the assumption of regularity of the vertex in the vanishing gluon momentum limit k + → 0 one can relate combinations of the transverse dressings λ to the longitudinal ones, hence fixing the former by the STIs. It has also been shown in [1,15], that the regularity of the quark-gluon vertex holds for momenta k 2 + > ∼ 1 GeV 2 , but fails below. Moreover, also the scattering kernels only deviate from unity in the infrared. In Figure 15 this is shown for the most relevant dressings λ (1,4,7) qqA at the symmetric point (13). Given the small angular dependence, the momentum behaviour at the symmetric point gives a very good estimate of the validity of the STI-regularity relations in the loops. This will provide us with the dressing functions ∆λ where the momentum dependence has been suppressed, and the superscript (in) indicates the input, here 2-flavour vacuum QCD data. We start with the generalised version of (17) for λ (1) qqA , and are led to ∆λ (1) qqA (p, q) =λ with Σ X and ∆ X as defined in (26). At the symmetric point (A9) reduces to (17). As an estimate for the systematic error we shall use both approximations forλ qqA is subleading, the other two are relevant, and are depicted in Figure 4.
We first discuss the dressing λ (4) qqA of the chirally breaking tensor structure T λ ( 1 ) q q A Z c / (Σ Z q λ ( 1 ) c c A ) λ ( 4 ) q q A p / ( Z 1/2 Α ∆ Z q M q ) λ ( 7 ) q q A p 3 Z c / Σ Z q M q ≡ 0 and a consistent extension to T, µ B = 0 should reflect this property. In the momentum regime with regularity it is related to the longitudinal dressing λ (10) qqA but none of the other longitudinal ones. This suggests a construction with ∆ ZqMq . Moreover, RG-consistency either enforces a dressing with Z 1/2 A or with 1/Z c , the difference being the momentum running of the (ghost-gluon propagator-) coupling. This leads us tõ , (A11) depicted in Figure 15.
We also emphasise that while the other STI relations used in the present work related to λ (1) 12 , the scalar tensor structure in the ghost-quark scattering kernel, that for λ ( 1 ) q q A Z c / (Σ Z q λ ( 1 ) c c A ) λ ( 2 ) q q A p / ( Z 1/2 Α ∆ Z q M q ) λ ( 3 ) q q A p 2 / (∆ Z q M q λ ( 1 ) c c A ) λ ( 4 ) q q A p / ( Z 1/2 Α ∆ Z q M q ) λ ( 5 ) q q A p 3 Z c / Σ Z q λ ( 6 ) q q A p 3 / Σ Z q λ ( 7 ) q q A p 3 Z c / Σ Z q λ ( 8 ) q λ (12) qqA also depends on λ 12 proportional to σ µν . This may explain the lack of Z c in the relation above. In any case it simply hints at a more complicated relation for λ (8) qqA .
Here we refrain from a more detailed discussion as λ qqA (p, q) =λ (7) qqA (p, q) Σ Zq − Σ (in) Zq (p, q) , ∆λ (8) qqA (p, q) =λ The respective ratios related to theλ (i) qqA with i = 1, ..., 8 are depicted in Figure 16, and the complete list of T, µ Bdependent dressing is summarised in (24). µ B /T < ∼ 3 or µ B < ∼ 400 MeV provide an estimate of the systematic error of the respective works. They can be traced back to the different approximation schemes in DSE and fRG. For larger baryon-chemical potential, µ B /T > ∼ 3, we approach the CEP in the fRG computation at µ B = 635 MeV, and the CEP in the present work is approached for larger chemical potential, µ B = 672 MeV. Hence, while the different locations of the CEPs also provide systematic errors for the prediction of the CEP, the larger deviations of the condensates does not indicate a qualitative failure of the approximations. Still, for µ B /T > ∼ 3 the respective approximations have to be systematically improved for narrowing down the location (and even presence) of the CEP.