The QCD phase structure at finite temperature and density

We discuss the phase structure of QCD for $N_f=2$ and $N_f=2+1$ dynamical quark flavours at finite temperature and baryon chemical potential. It emerges dynamically from the underlying fundamental interactions between quarks and gluons in our work. To this end, starting from the perturbative high-energy regime, we systematically integrate-out quantum fluctuations towards low energies by using the functional renormalisation group. By dynamically hadronising the dominant interaction channels responsible for the formation of light mesons and quark condensates, we are able to extract the phase diagram for $\mu_B/T \lesssim 6$. We find a critical endpoint at $(T_\text{CEP},{\mu_B}_{\text{CEP}})=(107, 635)\,\text{MeV}$. The curvature of the phase boundary at small chemical potential is $\kappa=0.0142(2)$, computed from the renormalised light chiral condensate $\Delta_{l,R}$. Furthermore, we find indications for an inhomogeneous regime in the vicinity and above the chiral transition for $\mu_B\gtrsim 417$ MeV. Where applicable, our results are in very good agreement with the most recent lattice results. We also compare to results from other functional methods and phenomenological freeze-out data. This indicates that a consistent picture of the phase structure at finite baryon chemical potential is beginning to emerge. The systematic uncertainty of our results grows large in the density regime around the critical endpoint and we discuss necessary improvements of our current approximation towards a quantitatively precise determination of QCD phase diagram.


I. INTRODUCTION
The detailed understanding of the QCD phase structure at finite temperature and density is not only essential for our understanding of the formation of matter, but also for the interpretation and prediction of the wealth of data collected at running and planned heavy-ion experiments. For overviews on both, experimental measurements and theoretical studies, see e.g. the reviews [1][2][3][4][5][6][7][8] and references therein. While our experimental and theoretical understanding of the phase structure at small densities has advanced rapidly in the past decade, at large densities theoretical and experimental investigations have been hampered by several intricacies. They range, e.g., from the sign problem of lattice gauge theory [9] to the influence of finite detector efficiency on signatures of the phase transition [10]. This leaves us with many highly relevant open questions regarding the phase diagram, in particular the existence and location of a critical end point (CEP) and the phase structure at small temperatures and large densities. The relevance of a CEP derives from the fact that the phase transition is of second order at this point. The resulting critical long-range correlations can potentially be observed, e.g., in particle number correlations measured in heavy-ion experiments, see e.g. [11]. Within the Beam Energy Scan (BES) Program at RHIC, significant measurements have been performed in this direction [2,[12][13][14][15]. This will be extended in BES phase II. Experimental studies of the QCD phase structure are also planned or run at other facilities with different collision energies and luminosities, such as CBM at FAIR [16] and HADES at GSI [17] in Germany, NA61/SHINE at CERN [18], the NICA/MPD in Russia [19], J-PARC-HI in Japan [20], and HIAF in China [21], see also [22][23][24].
In the present work we evaluate the phase structure of N f = 2 and N f = 2 + 1 flavour QCD within the fRG approach as introduced in [52], see also [53,54]. For QCDrelated reviews see e.g. [55][56][57][58][59][60][61]. We built upon the fRGresults for Yang-Mills theory in the vacuum, [31], and at finite temperature, [33], as well as vacuum QCD results in the quenched approximation, [27], and in full unquenched QCD, [28,29,32]. This work is extended to unquenched QCD at finite temperature and density, which gives us access to the chiral and confinement-deconfinement phase structure of QCD in terms of QCD correlation functions.
Our study cumulates in a prediction for the QCD phase diagram for µ B /T 6. For N f = 2 + 1 it is presented in Figure 1 on the next page together with a survey of other theoretical predictions as well as a compilation of freezeout points from different incarnations of the hadron resonance gas. We define the chiral phase boundary through the renormalised light chiral condensate ∆ l,R . At small and intermediate baryon chemical potential the transition is a crossover. At µ B = 0 the pseudocritical temper- Phase diagram for N f = 2 + 1 flavour QCD in comparison to other theoretical approaches and phenomenological freeze-out data. Our result for the chiral crossover is depicted by the black dashed line. The crossover temperature has been determined through the peak position of the thermal susceptibility of the renormalised light chiral condensate, ∂T ∆ l,R , at fixed baryon chemical potential µB. For more details see Section V A, and in particular Figure 10. We show dotted black lines for µB/T = 2, 3 to indicate the reliability bounds for the lattice and functional methods. The phase boundary globally agrees well with recent lattice results. In particular the curvature of the phase boundary for small chemical potential, κ = 0.0142 (2), is consistent with recent lattice results, κ = 0.0149 (21) in [44], κ = 0.0144 (26) in [47], and κ = 0.015 (4) in [49], for an overview see [62]. We find a critical end point at (TCEP, µB CEP ) = (107, 635) MeV. Indications for an inhomogeneous regime close to the chiral phase transition for µB 420 MeV are depicted by the hatched red area. For quantitative statements in this area the current approximation has to be upgraded systematically. Accordingly the hatched red area also serves as a reliability bound for the current approximation. For more details see Section V B and Figure 21.
Other theoretical results: lattice QCD based on an analytic continuation from the imaginary chemical potential [44] (WB), lattice QCD based on a Taylor expansion in chemical potential [49] (HotQCD), DSE approach with backcoupled quarks and a dressed vertex [37] (Fischer et al.), and DSE calculations with a gluon model [63] (Gao et al.). Freeze-out data: [2] (STAR), [ ature is T c = 156 MeV. The curvature of the chiral phase boundary at small chemical potential is κ = 0.0142 (2). With increasing µ B , the crossover becomes sharper and we find a critical endpoint at (T CEP , µ B CEP ) = (107, 635) MeV .
Our results for the chiral phase boundary are depicted by the black dashed line in Figure 1.
In addition to a CEP, we also find indications for an inhomogeneous regime for µ B 420 MeV in the vicinity and above the chiral phase boundary. It is given by the region in the phase diagram where mesonic dispersion relation develop a minimum at nonvanishing spatial momentum, for more details we refer to Section V B. This indicates a potential instability towards the formation of an inhomogeneous quark condensate. The region where this regime has significant overlap with the homogeneous chiral condensate is shown by the red hatched area in Figure 1. Within this area, a competition between homogeneous and inhomogeneous quark condensation has to be taken into account. Hence, this already suggests that the systematic error of the present approximation grows large for µ B /T 3.
In Figure 1 also we compare our results to recent predictions of lattice gauge theory for the phase structure at small µ B /T from the Wuppertal-Budapest Collaboration [44] (WB) and the HotQCD Collaboration [49] (HotQCD). Our result for the pseudocritical temperature and the curvature of the phase boundary agree very well with the lattice. We also show predictions of the DSE approach from different groups, [37] (Fischer et al.) and [63] (Gao et al.). Finally we included the freeze-out data from [2] (STAR), [64] (Alba et al.), [3] (Andronic et al.), [65] (Becattini et al.), and [66] (Vovchenko et al.). The freeze-out points are surprisingly close to our result for the chiral phase boundary, even at larger µ B . All in all, we see that a consistent picture of the QCD phase boundary at finite density starts to emerge form a culmination of results from different sources.
In order to discuss the implication for CEP searches, it is instructive to convert µ B to the center-of-mass beam energy per nucleon, √ s. Assuming the connection between these quantities is captured by the statistical hadronisation scenario, one finds to a very good approximation for central collisions the relation √ s = (a/µ B − 1)/b with a = 1307.5 MeV and b = 0.288 GeV −1 [3]. This yields for our prediction of the location of the CEP the beam energy √ s CEP ≈ 3.7 GeV .
This is clearly below the smallest beam energy of current BES measurements of √ s = 7.7 GeV, but well within reach of future experiments such as FAIR's SIS100 [16], NICA MPD [19], J-PARC HI [20], and STAR's Fixed-Target (FXT) program [68], see also [22][23][24]. Our results therefore provide a strong motivation for CEP searches at these future experiments. Furthermore, the inhomogeneous regime appears to be also within reach of heavyion collisions at small beam-energies. Hence, looking for experimental signatures of this regime might be a worthwhile endeavour. This work is organised as follows. In Section II we introduce the functional renormalisation group approach to QCD. In Section III, IV we discuss in detail the underlying systematic truncation scheme, and specify the flows for correlation functions including the propagators, vertices, and the effective potential. In Section V we present the numerical results and discuss them in detail. In Section VI we analyse the systematic error in the current approximation, also in relation to that of other functional approaches. In Section VII we close with a short summary. Many technical details of our calculations are deferred to the appendices.

A. QCD with dynamical hadronisation
The fRG approach is based on an infrared regularisation for momentum modes p 2 k 2 of the theory at hand, where k is the infrared (IR) cutoff scale. This is achieved by adding a momentum-dependent mass term to the action which vanishes for momenta p 2 k 2 . Accordingly, the ultraviolet quantum fluctuations of the theory with momenta p 2 k 2 are untouched. In turn, in the presence of the infrared cutoff, quantum fluctuations with momenta p 2 k 2 carry a mass proportional to k and are therefore suppressed. The presence of this regulator leads to a scale dependent effective action Γ k , which includes quantum fluctuations of momentum modes p 2 k 2 .
Thus, in the ultraviolet (UV) at large cutoff scales, k 2 /Λ 2 QCD 1, the effective action Γ k tends towards the bare QCD action, which is well under control with perturbation theory. This can be used as the starting point of a renormalisation group (RG) evolution of Γ k , where the IR cutoff scale k plays the rôle of an RG scale. By evolving Γ k from the UV to the IR, quantum fluctuations at a given momentum scale p 2 ∝ k 2 are included successively, and the full quantum theory of QCD is resolved as k → 0. Γ k=0 is the full quantum effective action of QCD. Its flow k∂ k Γ k , is given by the Wetterich equation [52], see also [53,54]. It is evident from the discussion above that the fRG is a practical realisation of the Wilson RG.
At low energies or momenta the dynamical degrees of freedom of QCD are hadrons rather than quarks and gluons. The offshell dynamics relevant for the flow equation is then dominated by the lightest hadrons. At small and intermediate densities, these are the pseudo-Goldstone bosons of spontaneous chiral symmetry breaking, in particular the pseudoscalar pions π, and the scalar resonance f 0 (500), or σ. Formulated in terms of the fundamental degrees of freedom, this requires taking care of the scalar-pseudoscalar four-quark interaction channels where these mesons emerge as resonances, as well as their scatterings. This is done conveniently by introducing a composite field where q is a Dirac spinor with N f flavours and N c colours. T a f , a = 1, . . . , N 2 f − 1 are the generators of the SU (N f ) flavour group and T 0 f = 1/ 2N f . For example, in the two-flavour case with up and down quarks, q = (u, d) we have the tensor structure (iγ 5 σ/2), where σ = (σ 1 , σ 2 , σ 3 ) are the Pauli matrices. This part of φ then corresponds to the three pions. The full field φ reflects the underlying SU (N f ) flavour symmetry for N f = 2, which translates into a O(4) symmetry for φ in case of isospin symmetric matter. In the physically more relevant case of N f = 2 + 1 with the light l = (u, d) quarks, the heavy s quark and q = (l, s) we simply embed (3) accordingly. This is discussed later, see (16).
The systematic introduction of the composite field φ in the fRG approach is done via dynamical hadronisation, see Refs. [56,[124][125][126]. The present formulation follows [27,28,32,56]: we introduce a scale-dependent composite fieldφ k (q,q) and φ = φ k through a source x J φφk and a respective regulator term. This field carries scalar and pseudoscalar quantum numbers. The effective action is then given by a modified Legendre transformation of the Schwinger functional W k [J] w.r.t. the current J = (J A , J c , Jc, J q , Jq, J φ ) , (4) including that of the composite fieldφ k . This leads us to with the currents J[Φ] defined by the equations of motion (EoM), and the cutoff term The cutoff term implements the IR regularisation though a momentum dependent mass-like term, as discussed above. The components of the superfield Φ in (5), (7) are gluons, ghosts, quarks and the scalar-pseudoscalar mesonic field φ, Φ = (Φ f , φ) , Φ f = (A, c,c, q,q) , φ = (σ, π) , (8) and q = l with l = (u, d) for N f = 2, and q = (l, s) for N f = 2 + 1. Note also that the cutoff term includes an infrared regularisation with R φ for the composite field.
This leads us to the regulator matrix The regulator R k (p) specifies the momentum dependence of the mass-like cutoff. It is chosen such that it suppresses quantum fluctuations with momenta smaller than the cutoff scale, p 2 k 2 , while it leaves the UV physics with momenta p 2 k 2 unaffected. Furthermore, to recover the full quantum effective action at k = 0, we demand R k→0 → 0.
This setup is closely related to a two-particle irreducible (2PI) or rather two-particle point irreducible (2PPI) formulation. If also considering the density channelqγ 0 q it resembles density functional theory, for a more detailed discussion of these relations for generic composite fields see [56].
The scale evolution ofφ k , or rather its expectation value ∂ tφk , can be chosen freely. Its choice corresponds to a reparameterisation of the theory. We emphasise that the mean field φ = φ k in the effective action is k-independent. Note also that on the EoS of the auxiliary field the effective action reduces to the standard effective action of QCD in terms of the fundamental fields, Γ QCD [Φ f ]: The EoM of the composite field including the cutoff term entail a vanishing current, Since the composite field is introduced only trough the source J φ and its regulator R φ in the first place, the EoM (10) removesφ k from the path integral at vanishing cutoff, reducing it to the standard gauge fixed path integral of QCD. We are led to At finite cutoff scale the composite field can also be eliminated. There, vanishing of the current J φ is obtained for δ(Γ k + ∆S k )/δφ = 0, and the infrared regularised path integral at J φ = 0 still depends on the cutoff term of the composite field. This amounts to inserting a UVirrelevant four-quark interaction in the classical QCD action. This procedure does not spoil the renormalisability as a pointlike NJL-type interaction does, but solely provides an IR regularisation of the respective resonant four-quark interaction.
In summary this setup encodes the full information of the QCD correlation functions but also allows for a simple access to bound state information such as the Bethe-Salpeter wave functions, see [27,28,33,127]. Note also that in general the QCD correlation functions now involve derivatives of the composite field. As an important example we consider general four-quark vertices. They are given by functional derivatives of the QCD effective action Γ QCD on the EoM Φ f,EoM = 0 at T . At finite temperature Φ f,EoM contains a nonvanishing temporal gluon background field, A 0,EoM , which carries the information of confinement, see [69,[128][129][130].
In (12) it is understood that φ = φ EoM [Φ f ], and we have restricted ourselves to the T = 0 case with Φ f,EoM = 0.
In (12) we have also introduced our notation for n-point correlation functions or vertices, If we choose the composite field such that it completely absorbs a given momentum channel in the four-quark scattering, the first term on the right hand side in (12) vanishes and we are left with exchange terms of the composite field. For example, for the pseudoscalar channel the second and third line in (12) comprise terms with a pion propagator (1/Γ (2) ) ππ with two Bethe-Salpeter wave functions Γ (3) qqπ attached. Note also that nontrivial terms such as in the second and third line of (12) only occur in correlation functions of the fundamental fields Φ f for more than one quarkantiquark pair. For example, for the quark two-point function or inverse propagator we find schematically Γ (2) QCD,q1q1 = Γ (2) 0,q1q1 + δ 2 φ δq 1 δq 1 Γ 0,φ = Γ (2) 0,q1q1 , (14) with φ = φ EoM [Φ f ]. Then the term proportional to Γ (1) 0,φ vanishes as the latter is the EoM for φ. Finally we use Φ f = 0 as in (12). For the dynamical hadronisation in the present work we use the option of absorbing the dominant four-quark interaction channel completely. This is achieved by choosing the flow ofφ k such that the scalar-pseudoscalar channel in the effective four-quark interaction of the light quarks u, d is eliminated. Formally, this can be viewed as a successive bosonisation in terms of a Hubbard-Stratonovich transformation of this interaction channel at each RG scale k. For more details see Section III E. This puts forwards the general hadronisation relation with τ = (T 0 2 , iγ 5 T 2 ) and T a f defined in (3). The vector e σ points in the σ-direction, with the convenient normalisationê σ = 1/ 2N f matching that of τ . Since we only consider the dynamical hadronisation of the π-σ channel of the light quarks, we embed the SU (2) scalarpseudoscalar field (σ, π) into SU (3). This implies in (15), in a slight abuse of notation. The hadronisation functionṡ A k ,Ḃ k andĊ k can be chosen freely. The hadronisation functionȦ k controls the overlap of φ with the scalarpseudoscalar channel, whileḂ k simply changes the wave function renormalisation Z φ of the composite field. The latter function can e.g. be used to conveniently choose Z φ = 1. As it does not change the parameterisation of the theory we discard it in the following, usingḂ k = 0. Finally, the hadronisation functionĊ k introduces a shift in the scalar field σ. In general, this shift can be used to entirely absorb the quark mass of the respective quark flavour, leading to chiral quarks. However, this comes at the expense of diffusing the symmetry properties of the effective action. For example, if we start with a Z 2 -symmetric effective potential V (σ) = V (−σ) necessary (but not sufficient) for chiral symmetry, a shift in σ introduces odd powers in the σ field in the effective potential. Of course, chiral symmetry is not lost, but the respective symmetry transformation is not simply σ → −σ anymore. A prominent example is the σmass term 1/2m 2 σ σ 2 . A shift in σ with σ → σ − c σ /m 2 σ leads to a linear term in the effective action, while still keeping the original σ-mass term. In the presence of higher powers of σ, more σ-odd terms are generated as well as changing the coefficients of the σ-even terms. Consequently we expect thatĊ k cannot be chosen freely if we restrict ourselves to dynamical hadronisation setups that keep chiral symmetry apparent. Indeed this constraint leads us toĊ k ≡ 0, see Appendix B. The linear term (17) plays a special rôle in the effective action. To see this more clearly, we solve the Legendre transform (5) for the Schwinger functional and use the EoM for the current leading to J[Φ] defined in (6). With these definitions (17) implies that J[Φ] contains the term −c σ . Then, the Schwinger functional reads, where J = J[Φ]. Evidently the c σ -dependences in the first two terms on the right hand side of (18) cancel each other. We conclude already from here that the flow equation of the effective action with dynamical hadronisation should be explicitly c σ -independent. Moreover, necessarily the left hand side of (18) is also c σ -independent.
Consequently, the c σ -dependence of the current in the Schwinger functional is canceled by that of W k [0]. This leaves us with and J = J[Φ] satisfies the EoM (6). The current is shifted in the J σ -direction,Ĵ[Φ] = J[φ] + c σ δ Φσ , in components: The shifted currentĴ [Φ] in (20) does not depend on c σ , that is ∂ cσĴ ≡ 0. With these relation we finally arrive at a convenient form of the effective action, We would like to elucidate the above relations within a simple example which is relevant in the present context: Consider a composite σ-field that is just proportional to the scalar quark bilinear, with a cutoff dependent prefactor D k . Note that the hats indicate that this relation holds on the level of the fluctuation fields in the path integral. We also emphasise that our example (22) has an overlap with the first term in the dynamical hadronisation flow (15) . Importantly, the linear term − x c σ σ in the effective action (17) is in one to one correspondence to the quark (current) mass term in the classical action, S QCD ∝ m 0 q qq . Obviously, the latter term can be absorbed into a shift of the source term forσ with We conclude that the Schwinger functional with a source for the composite fieldσ k in (22) and the classical QCD action with a quark current mass term, S QCD ∝ m 0 q qq , is that in the chiral limit with m 0 q = 0 and a shifted current for the composite field. However, m 0 q = 0 simply is c σ = 0, see (24), and we arrive at Note that (25) entails that the full dynamics of QCD in the presence of explicit chiral symmetry breaking is that of the theory with a composite condensate fieldσ with full chiral symmetry: The explicit chiral symmetry breaking is completely absorbed in a shift of the external current J →Ĵ.
Flow of the effective action of QCD. The first three diagrams arise from the gluon, ghost, and the quark degrees of freedom respectively. The last diagram is that of the mesonic contribution. The double line with the up-down arrows indicates the nature of the mesons as quark-antiquark composites. The crossed circles indicate the regulator insertion in the flow equation.
In summary this leads us to the flow equation for QCD with dynamical hadronisation, [56,131], with the RG time t = ln(k/Λ), and The shift in the dynamical hadronisation term in the first line in (26) subtracts the explicit chiral symmetry breaking term in Γ (1) σ . The t-derivative ∂ t Γ k in (26) is taken at fixed c σ as theċ σ terms from ∂ t Γ k cancel with that from The (∂ t J sym )-terms are cancel by those from the Schwinger functional in (21). We emphasise that (26) is explicitly c σ -independent, and is a novel representation of dynamical hadronisation.
Note also that (26) with (15) entails that we do not have to specifyφ k but only the expectation value of its flow, ∂ tφk , as done in (15). The second term in the first line of (26) accounts for the cutoff dependence of the reparameterisation that originates in the cutoff dependence ofφ k . It can be understood as a generalised anomalous dimension of φ. Indeed, for this term carries nothing but the anomalous rescaling of the composite field φ. However, we emphasise again that ∂ t is a derivative at fixed mean superfield Φ, and in particular at fixed φ. The second term in the second line also accounts for the cutoff dependence of the reparameterisation via ∂ tφk . The first term in the second line of (26) is the standard functional flow. The trace sums over momenta, internal indices and species of fields including the composite field.
More explicitly it reads where all propagators are Φ-dependent. We show the flow equation (26) with ∂ tφ = 0 diagrammatically in Figure 2 for QCD. The first gluon and ghost diagrams constitute the glue contributions, and the last two arise from the matter sector, which are the quark and meson degrees of freedom, respectively. Note that the scalar and pseudoscalar mesonic loops are present here because we introduced them explicitly with dynamical hadronisation. This is nothing but a convenient parametrisation of resonant interaction channels. We emphasise that the absense of explicit loops for other hadrons, e.g. other mesons and baryons, does not entail that their dynamics is not included.

III. TRUNCATION SCHEME
In this section we discuss in detail the expansion scheme for the effective action and the flow equations used in the present work. This includes in particular a discussion of the approximations used, and their quantitative or qualitative validity.

A. Vertex expansion
In QCD, (26) can only be solved within given truncations to the effective action. The systematic expansion scheme behind such approximations is the vertex expansion. It is an expansion of the effective action in powers of the field, the expansion coefficients being the n-point functions Γ (n) (p 1 , ..., p n ). Typically, this expansion shows a rather good convergence if no divergent exchange processes or couplings are present. In Landau gauge QCD the strong couplings related to gluon exchange grow towards the infrared but tend to zero for momenta or cutoff scales below approximately one GeV, reflecting the QCD mass gap, see [27-29, 32, 131]. This behaviour supports the apparent convergence of the vertex expansion. Another source for divergent exchange processes are resonant interactions that potentially spoil or at least slow down the convergence of the vertex expansion. In particular for this reason dynamical hadronisation is crucial as it controls the resonant interaction channels and allows to take into account multi-scatterings of resonances in a technically accessible way via the respective effective potential of the composite fields. We can rely on quantitative results for vacuum Yang-Mills theory, [31] and N f = 2 flavour QCD, [27,32], and finite temperature Yang-Mills theory, [33], as input as well as benchmark tests for our current computations. State-of-the-art truncations which involve a large set of correlation functions have been used in these works.
Furthermore, it is well-known that mild momentumdependences of vertex functions and propagators are well captured by scale-dependent dressing functions; for investigations in QCD and low energy effective models see e.g. [28,70,73,95,131]. This suggests an approximation scheme in which only the dominant nontrivial momentum dependences are taken into account explicitly, while the rest of the dependences is approximated by scaledependent dressings.
This approach has been successfully applied to vacuum QCD in [28,29,131] and for QCD at finite temperature and imaginary chemical potential in [26]. In these works QCD flows have been expanded about the full gluon and ghost two-point functions of Yang-Mills theory. All other correlation functions considered there have been approximated with scale-dependent dressings. Schematically we decompose the two-point functions as Γ (2) Z S (2) kin + (additional tensor structures). For example, for the quark, S (2) kin carries the Dirac term, while the mass term is buried in the additional tensor structures. For the gluon this entails for the two-point function is Γ AA (p), but we have to consider different Lorentz tensors. The nontrivial momentumand scale-dependence of the kinetic term of a given field is fully captured by the anomalous dimensions, Note that (30) only schematically provides the anomalous dimension, we have neither specified the projection procedure nor the momentum evaluation. We further exemplify the setup with the relevant case of the gluon two-point function. There, one also utilises the observation that the anomalous dimension η A of the gluon can be parameterised in terms of the running coupling α YM k in Yang-Mills theory, with the Yang-Mills running coupling α YM k and the transversal gluon mass parameter, where Π ⊥ µν is the transversal projection operator, see (N2). The renormalised gluon mass parameterm 2 A differs from m 2 A in (32) by an appropriate renormalisation.
This renormalisation has to be applied to all coupling parameters and fields and is discussed in detail in the next Section III C. In (31) the Yang-Mills gluon mass parameter (m YM A,k ) 2 has been used. η glue A stands for all diagrams involving gluons and ghosts that contribute to the gluon anomalous dimension. This parameterisation of the glue part of the anomalous dimension allows for a simple representation of the anomalous dimension of full QCD, with the full QCD coupling α k and renormalised mass parameterm 2 A,k = m 2 A,k /Z A , respectively. Importantly, this setup also takes care of the backreaction of the quark fluctuations to the pure glue and ghost diagrams, for a detailed description see [28,131]. A benchmark test concerns the full gluon propagator in unquenched twoflavour QCD, which has been shown to be in quantitative agreement with the lattice results in [131], more details can be found in Appendix D.
A simpler approximation is given by the 'direct sum' With (34) the backreaction of the quark loop on the pure glue couplings and propagators in η glue A is neglected. While such an approximation also lacks full RG invariance, it works well for initial scales of Λ ≈ 10 − 20 GeV. Alternatively a readjustment of the scales can be done on top of the simple approximation (34). A respective benchmark test for the full gluon propagator in unquenched two-flavour QCD has been done in [131], for more details see Appendix D. The direct-sum approximation has been used very successfully in many functional applications to the phase structure of QCD, notably in fRG application to QCD, [26], as well as the DSE applications [7,36,37,[132][133][134][135][136][137][138]. We also emphasise that the approximation may get even better for large chemical potential, as the latter lifts the fermi sea and successively more quark fluctuations are buried in it. This effect reduces the impact of the quark loop on the purely gluonic correlation functions. However, this reasoning does not apply straightforwardly to the quark-gluon vertex, and has to be taken with a grain of salt beyond the onset of baryon density. Still, in summary the approximation is well-founded and well-tested.
Here we apply this promising approximation scheme to the unquenched QCD gluon propagator at vanishing temperature as well as the glue contribution at finite temperature. Let us describe this procedure with an expansion of the finite temperature and density theory about the vacuum. Schematically this is described by the following separation of the flow for a given correlation function at finite T and µ, where we suppressed the subscript k, and ∆Γ (n) is de-fined by, Its flow is given by the difference of the flow diagrams for the correlation function at finite T, µ and the vacuum, with Flow (n) While the first term on the right hand side of (37) depends on ∆Γ In summary this allows for the use of scale-dependent dressings for ∆Γ (n) with a mild momentum dependence at finite temperature and density, while still maintaining the quantitative nature of the approximation: Quantitative results for vacuum QCD and finite temperature Yang-Mills theory [27,[31][32][33] are used for the nontrivial momentum dependences of this part of the correlation functions. Note that this setup also allows for the use of quantitative results from other approaches such as lattice simulations and DSE computations.
In the present work we apply this approach to the gluon two-point function Γ (2) AA . Its density corrections and thermal quark fluctuations are computed with the input of the quantitatively reliable vacuum QCD results from [32] and the finite temperature Yang-Mills results in [33]. For the ghost two-point function Γ (2) cc we only use input from vacuum QCD from [32]. The details can be found in Section IV A 2.

C. Truncation for the effective action
Having captured the nontrivial momentumdependence of the gluon propagation, we adopt the following truncation for the Euclidean scale-dependent effective action Γ k for both, N f = 2 and N f = 2 + 1 flavours, (3). Note that all couplings depend on the RG scale k, which in most cases is omitted for the sake of clarity. The effective potential carries a dependence on the mesonic field ρ = (σ 2 + π 2 )/2, and the temporal gauge field A 0 , which is related to the Polyakov loop L[A 0 ]. The expectation values of these fields are approximate order parameters for the chiral phase transition (ρ) and the confinement-deconfinement phase transition (A 0 or L[A 0 ]), for more details see Section III F. In (40) we separated the contribution from the gluonic loop in (26), see Figure 2, that is V glue,k and that of the quark and meson loops, V mat,k . We also neglected the subleading ρ-dependences of the gluon loop, The gluon two-point function Γ AA has been already discussed in the previous Section III B. The subtraction with Z A p 2 in the second line removes the A 2 -contribution coming from the F 2 µν -term for avoiding double counting. In the present fRG approach the theory is fully determined by the strong running coupling α s,Λ at the initial cutoff scale k = Λ in the UV, the light current quark mass m 0 l for N f = 2 and additionally the strange current quark m 0 s for N f = 2 + 1 with the choice With dynamical hadronisation the current quark masses are encoded in c σ and c σs . Their values in the physical case are fixed by the pion mass and the kaon mass. Alternatively, the pion and kaon decay constants can be used to pin down these parameters. This is explained below. The values for an initial cutoff Λ = 20 GeV can be found in Table I at the end of this Section. The field strength tensor F a µν and the covariant derivatives in the ghost and quark kinetic terms also carry a cutoff dependence, and are defined and discussed in detail in Section III D, (60) and (57) respectively. The Z Φi 's are the cutoff dependent wave function renormalisations of the respective component fields Φ i of the superfield defined in (8). Note that the effective action is renormalisation group invariant (but k-dependent) w.r.t. the underlying renormalisation group scale µ of the theory at k = 0, for more details see [56]. This suggests the introduction of renormalised fields as well as renormalised coupling parameters, for examplē Since we perform our computations in Euclidean spacetime, the renormalised masses in (43) are curvature masses in contradistinction to the physical pole masses m Φ,pol . Within the fRG setup this difference has been discussed in detail in [95]. For example, the constituent quark masses are curvature masses. In turn, the pion pole mass m π,pol is used for fixing the value of explicit chiral symmetry breaking, see Table I. For more details on the definition of the mesonic pole masses see Section IV A 1. Similar relations hold for the other coupling parameters and are discussed later. Note that the epithet renormalised is related to the k-dependence. Note also that the renormalised fields carry classical dispersions if we use fully momentum-dependent Z Φ 's in (42). Accordingly in pole renormalisation schemes the respective masses are pole masses. Moreover, observables are provided in terms of the renormalised fields and couplings, and we shall discuss physics in terms of these objects.
The quark chemical potential matrixμ considered here is understood to couple to all flavours equally,μ = diag(µ q , µ q , µ q ). Thus, it is directly related to the baryon chemical potential µ B = 3µ q . The constituent quark massm s also has to be read as a diagonal flavour matrix diag(0, 0,m s ).
In the present work we do not consider dynamical hadronisation of the full N f = 3 flavour multiplet T a N f (σ a + i γ 5 π a ), with the U (N f ) generators T a N f , a = 0, ..., N 2 f − 1. In the singlet-octet basis they are given by T 0 Instead, we only consider part of the embedded N f = 2 multiplet: the scalar σ as well as the pseudoscalar pions. The four-quark interaction that gives rise to the corresponding resonances is given by λ q . The other mesons of the scalar and pseudoscalar nonets are too heavy for playing a significant rôle in the offshell dynamics. Of the strange part of the multiplet we only consider the scalar σ s . The two scalars σ, σ s are related to the ones in the singlet-octet basis via for more details see e.g. [34,73,89,93,120,139]. Both, light and strange quark constituent masses are given in terms of condensates, with the renormalised Yukawa coupling and the renormalised σ, σ s -expectation valueσ = Z 1/2 φ σ andσ s = Z 1/2 φ σ s respectively, see (42) and (43). Their values are determined via the explicit breaking terms c σ σ and 1/ √ 2 c σs σ s , as well as dynamical chiral symmetry breaking. Both phenomena lead to nonvanishing σ EoM and σ s,EoM and hence to nonvanishing renormalised expectation values (σ,σ s ) EoM . The explicit breaking coefficient c σ in the light sector is fixed with the physical pion mass. Alternatively one could use the ratio of the pion decay constant f π with that in the chiral limit, f π,χ for fixing c σ , that is f π /f π,χ ≈ 93/88, [140]. Here, we choose the pion mass instead of this ratio as it is more easily accessible. Furthermore, small errors in the determination of the pion mass only propagate to small errors in other observables. In turn, the explicit breaking coefficient c σs may be either determined by e.g. the kaon mass or the ratio f K /f π ≈ 111/93, [140]. We refer to Appendix E 3 for a more detailed discussion of the scale setting.
The relative size of the explicit breaking coefficients, c σs /c σ may also be determined by their relation to the current quark masses m 0 l and m 0 s . For large momentum scales of the order of the electro-weak scale ∼ 90 GeV, that is in the absence of any chiral dynamics, the constituent quark masses (45) reduce to the current quark masses. Note also that the renormalised mesonic quantities tend towards bare ones for large momentum or cutoff scales as Z φ → 1. Moreover, Z q → 1 due to the Landau gauge. Accordingly,σ,σ s → σ, σ s andh,m q → h,m q . This limit entails, where m 2 φ is the unique mesonic mass function for large momentum scales. Equation (46) entails that the ratio of the current quark masses agrees with that of the explicit breaking parameters, in N f = 2 + 1 flavour computations. The estimate comes from lattice computations, see [141]. Equation (46) and (47) also enable us to relate the chiral condensates to c σand c σs -derivatives, for more details see Appendix A. Finally, one may also adjust the constituent strange quark massm s or the difference to the constituent light quark massm l , on the basis of quantitative functional or lattice results in the Landau gauge.
In the present work we compute observables in the light quark sector. We use that offshell flavour-mixing terms are small, as they always involve the propagator of a heavy mesonic state. This is in stark contradistinction to onshell flavour-mixing terms, which are e.g. maximal in the pseudoscalar sector due to the axial anomaly. Accordingly, light quark and gluon correlation functions are sensitive to strange quark fluctuations only via the gluon propagator or rather the gluon dressing. The latter carries the momentum and RG running of the strangequark-induced vacuum polarisation, and it is well-known from respective quantitative N f = 2 flavour computations that the gluon propagator is almost insensitive to changes of the quark mass. This has been studied in [33], where the pion mass (and hence the light quark mass) has been changed from very light quarks to heavy ones in the range 60 MeV m π 300 MeV. This change had no effect on the gluon propagator within the systematic error bars of the result. This analysis carries over readily to the present N f = 2 + 1 flavour computation. Indeed, the influence of the strange quark mass is even smaller due to the smaller relative importance of the strange quark and its more effective decoupling in the infrared due to the larger explicit chiral symmetry breaking.
Consequently we use a simple approximation to the strange sector. The chiral offshell dynamics are dominated by the pions, and we approximatē With (49) the chiral dynamics in the strange quark sector are the same as in the light quark sector. All these determinations have to be taken with a grain of salt due to the rough nature of our approximation of the strange quark sector. While this approximation has to be improved for an access to observables with strangeness, the observables considered here depend only very mildly on the difference between the constituent light and strange quark masses, (48), For the determination of ∆m sl in the present work we use the ratio of the decay constants, which relates the current strange quark mass directly to observables. In the mean field approximation in low energy effective theories we typically have f π ≈σ EoM and f K ≈ 1/2σ EoM + 1/ √ 2σ s,EoM . We emphasise that these relations do not hold true in QCD as the determination of the decay constant requires the full momentum dependence of the quark mass functions; for a detailed discussion of the pion decay constant see Appendix E 3. However, the cutoff scale (and momentum) dependence of both the light and the strange quark masses are very similar (after being rescaled by their value at vanishing momentum). Hence we conclude that the mean field relations should hold even quantitatively for the ratios of the decay constants. Using (49) in the mean field approximation for the ratio of kaon decay constant and pion decay constant, this is well adjusted with ∆m sl = 120 MeV, Equation (51) is in good agreement with the actual value f K fπ ≈ 1.19, see [140]. More accurately we may derive the mass difference with a comparison to N f = 2 + 1 lattice QCD (with isospin symmetry). The ratio of the decay constant is determined with f K /f π = 1.194 (5), [141], and provides ∆m sl ≈ 135 MeV. For our choice ∆m sl = 120 MeV, the constituent strange quark mass is found to bem s = 467 MeV, in good agreement with quark model values. Note however, that the Landau gauge constituent strange quark mass is considerably higher,m s 500 MeV, see e.g. [7], which would amount to a current strange quark mass of ∆m sl 150 MeV in the present approach. This is further discussed in Appendix A. We emphasise again that a variation of ∆m sl within this range does not influence our results for light quark observables.
Next we discuss the treatment of the four-quark sector in (39): We only consider the four-quark interaction λ q of the N f = 2 scalar-pseudoscalar multiplet as well as the corresponding Yukawa interaction h and the mesonic composite fields. In turn, the Dirac term depends on all quark fields, so either u, d in the N f = 2 flavour case or, u, d, s in the N f = 2 + 1 flavour case. Accordingly, the mesonic field φ = (σ, π) in (39) is in the O(4)-representation with ρ = φ 2 /2 for both, N f = 2 and N f = 2+1. The linear term −c σ σ breaks the two flavour chiral symmetry explicitly, and leads to the current quark masses for u, d quarks. We assume light isospin symmetry here, so they have identical masses. They are absorbed in a shift of σ, leading to current quark mass terms from the Yukawa term. In the present setup the c σ -term is generated from the u, d-quark mass term via dynamical hadronisation with an appropriate choice oḟ C k in (15), for more details see Section III E. It is now apparent that the shifted σ-current coupled to the dynamical hadronisation flowφ k is chirally symmetric, as it does not depend on c σ , In turn, we do not consider offshell fluctuations from the four-quark interactions with strangeness: In (39) neither four-quark terms nor mesonic fields with strangeness are included. This approximation is based on the observation that at T, µ = 0 even the two-flavour scalar-pseudoscalar terms produce negligible contributions for cutoff scales k 500 MeV above the onset of chiral symmetry breaking, [27,33]. For cutoff scales in the vicinity of the onset of chiral symmetry breaking and below, k 500 MeV the other two-flavour channels and even more so the s-quark channels are not dynamical anymore due to their large mass scales. Indeed, sizable contributions at T, µ = 0 are only triggered by the pion channel. This is in line and supports chiral perturbation theory. the strong coupling αs,Λ, the pion pole mass m π,pol via cσ, and the ratio fK /fπ via the constituent quark mass difference ∆m sl =ms −m l for the present N f = 2 and N f = 2 + 1 flavour computations.
Middle part: IR enhancement of quark-gluon coupling αqAq → (1 + a) αqAq below k ≈ b. The value of a is adjusted with the constituent light quark massm l , for more details see Appendix E 2. This phenomenological IR-enhancement effectively accounts for the effect of nonclassical tensor structures in the quark-gluon vertex which are missing in the present approximation. If taking the full quark-gluon vertex into account, this is not necessary andm l is a prediction, see [27,32].
Lower part: predictions in vacuum QCD: the strange constituent quark massesms and the σ-massmσ. Also fπ is a prediction as we have fixed the pion pole mass m π,pol instead of fπ. Fixing the latter relative to the pion decay constant fπ,χ in the chiral limit (in the present work fπ,χ = 88.6 MeV for N f = 2 + 1) would have been the more physical but less accessible choice, for a detailed discussion see Appendix E 3.
Note however, that e.g. in the vicinity of the phase transition, kaons and the eta meson may play a rôle. Neglecting this is part of our current systematic error. We also emphasise that at large densities we expect relevant offshell contributions from diquark and/or density channels. The importance of the additional N f = 2 flavour channels has been investigated thoroughly in effective theories in [74,76,115], leading to a semiquantitative agreement of both approximations up to large densities after an appropriate rescaling including the critical region found in the present work for N f = 2 flavour QCD, see Section V D and Figure 20. This estimate is fully confirmed in a QCD study with the fRG in [142]. Note that the observed dominance of the scalar-pseudoscalar channel for N f = 2 flavours in [74,76,142] translates to µ B /T 7 in the present N f = 2 + 1 flavour study, and includes the respective critical region, see Figure 1.
Accordingly, the N f = 2 and N f = 2 + 1 theories differ by the Dirac term in the effective action, and hence by the  [144][145][146]. The momentum dependence of the fRG results is given by p 2 = k 2 , in accordance with the same identification of the N f = 2 flavour input data; for more details see Section IV A 2.
respective additional s-quark loops. This is important for purely gluonic correlation functions and amounts to a relative change in the physics scale Λ QCD as well as the respective β-functions. This is very similar to respective DSE computations, a difference being the backreaction onto the purely gluonic diagrams which is partly taken into account in the present work.
In summary the couplings α s,Λ , c σ and ∆m sl (or the ratio of the current quark masses m 0 s /m 0 l = c σs /c σ , see (47)) with Z Φ,Λ = 1 in the initial action Γ Λ at the initial UV scale Λ = 20 GeV are fixed by the fundamental parameters of QCD; for the strong coupling and the current quark masses, as well as phenomenological infrared enhancement parameters a, b for the quark-gluon coupling, see Table I. The latter phenomenological parameters, (a, b), effectively account for the infrared effects of the missing nonclassical tensor structures of the quark-gluon vertex, see Appendix E 2.
Strictly speaking, absolute momentum scales should be measured in the pion decay constant in the chiral limit. In the present approximation it is computed as see (E16). The details of the scale setting procedure are discussed in Appendix E 3. Apart from the quantitative precision of the predictions for observables in Table I in particular for N f = 2 + 1, the N f = 2 + 1 gluon dressing is a prediction within the current setup. It agrees quantitatively with the respective lattice results, see Figure 3.
Note that the gluon propagator shows the expected deviations in the deep infrared: the N f = 2 flavour propagator in [32] is of the scaling type, related to a different infrared closure of the Landau gauge as that implicitly defined by the lattice gauge fixing. This has been discussed at length in the literature, see [32] and references therein. Importantly, as part of the gauge fixing, it has no impact on observables.
D. Strong couplings αs from gluon, ghost-gluon and quark-gluon vertices The quark-gluon and purely gluonic correlation functions in the current approximation give rise to different 'avatars' of the strong couplings α Φ1···Φn with Φ i = A µ , c,c, q,q and n = 3, 4. For details see [27,[31][32][33], for similar considerations in quantum gravity see [147,148]. They are defined by the respective vertex dressings λ Φ1···Φn , divided by the wave function renormalisations Z

1/2
Φi of the attached legs evaluated at vanishing momentum, see e.g. [32]. Here we consider the couplings of the purely gluonic sector, and that of the matter sector, αq Aq = (αl Al , αs Ac ) with two entries for the light quark-gluon and strange quarkgluon couplings respectively, Note that the strong couplings (54), (55) are natural definitions of gluon exchange couplings. For example, the quark-gluon couplings involve two quark-gluon vertex functions and one gluon dressing 1/Z A , related to the one-gluon scattering of two quark currents, see Figure 4. For perturbative and semiperturbative scales k 1 − 3 GeV these couplings are related by modified Slavnov-Taylor identities (mSTIs). This property also holds true in the quantitative approximation [27,[31][32][33] for the purely gluonic couplings, for more details see in particular [32]. We emphasise that in the present setup gauge invariance is encoded in these mSTIs. They arise from BRST-variations of the gauge fixed action in the presence of the regulator terms, see e.g. [56,[149][150][151][152][153][154][155][156][157][158]. For smaller momenta and cutoff scales k 1 − 3 GeV they start to differ for two reasons: first of all the gluon mass scale starts to influence the running. More importantly the mSTIs only relate the longitudinal couplings, while the flow depends on the transversal couplings. It can be shown that some of the transversal α's cannot be identified with their longitudinal counterparts as otherwise the gluon mass scale would be absent, see [31,[159][160][161][162]. However, in the absence of the gluon mass scale in the propagator the theory lacks confinement, [69,128,129].
Here we compute the flows of the three-gluon coupling, α A 3 , and the quark-gluon couplings αq Aq . From the vacuum results for the couplings in [28,32] we infer that the four-gluon coupling runs closely to α A 3 , and the ghost-gluon coupling runs closely to αq Aq in the relevant momentum regime k 1 GeV. In the deep infrared all these couplings. However, except of the ghostgluon coupling, all the avatars of the strong coupling decay for these scales, see Figure 17 in Section IV B 3, and [28,32]. The ghost-gluon coupling only enters in diagrams that are suppressed except for the deep infrared with k 100 MeV, a more detailed discussion of this fact is provided in Section IV B 2. In summary the differences lead to negligible effects and we use For gluons and quarks this leads us to covariant derivatives in the fundamental and adjoint representation of SU (N c ) respectively, where The covariant derivative in the ghost terms is an adjoint one but carries gq Aq . In (57), f abc are the structure con- The terms in the first two lines on the r.h.s. of (39) are built from the operators in the classical QCD action. The gluonic field strength tensor reads and the Landau gauge ξ = 0 is chosen in this work. A similar truncation has been found to be successful in describing the transition from the quark-gluon regime at high energy to the hadronic one at low energy in the vacuum [28]. The truncation (39) for the effective action takes into account all order scatterings of the resonant scalar and pseudoscalar channels of the four-fermi interaction through the scale-dependent effective potential V k (ρ, L,L). It is complementary to taking into account a Fierz-complete four-fermi basis as done in [27,32,74,76,122]. A respective two flavour study with a Fierz complete basis and all order interactions of resonant channels is work in progress, [163].
E. Dynamical hadronisation in the σ − π channel As discussed in Section II, we take into account the resonant interactions in the σ − π channel with the help its where it is, are still open to be answered. For- of dynamical hadronisation. We use (15) in the flow equation (26). Our strategy is to choose the coefficients in (61) such that the running of the fourquark interaction λ q in (39) is exactly cancelled. This amounts to a scale dependent bosonisation of this channel. We emphasise that this only transfers the information carried by four-quark interaction from the quarkgluon sector to the meson sector. It is a practical way to enter the symmetry broken phase in the presence of resonant channels and still retain the full information of the underlying dynamics of the fundamental degrees of freedom. In case of the scalar-pseudoscalar channel considered here, a resonance at vanishing momentum indicates the formation of a quark condensate qq and therefore chiral symmetry breaking.
To this end, we project the flow equation on the flow of the scalar-pseudoscalar four-quark vertex, λ q in (39) at vanishing momentum. We define the respective dimensionless, renormalised four-quark coupling, Yukawa coupling and hadronisation function as where we suppressed the subscript k, and the renormalised Yukawa coupling has already been introduced in (43). With these definitions the flow of the renormalised four-quark function reads where Flow (n) are the dimensionless, renormalised flow diagrams. The projection on the scalar-pseudoscalar four-quark channel in (64) is done along the lines of [60] with an expansion of the flow in terms of quark-bilinears. This is indicated by the subscript (qτ q)(qτ q) . Flow (n) is defined in (38), and d λΦ i 1 ···Φ in is the canonical momentum dimension of the vertex dressing λ Φi 1 ···Φi n . Now we resort to fully hadronised flows in the σ-π channel by demandinḡ With this, all diagrams with four-quark verticesλ q are absent and the flow is depicted by This choice entails that the quantum fluctuations contributing to the quark scattering in the σ-π channel are transferred completely to effective hadronic degrees of freedom, the σ and π fields, their propagation, selfscattering and scattering with quark-antiquark pairs. The original four-quark interaction is described with instantaneous meson propagators and a Yukawa interaction of the σ, π with a quark-antiquark pair. Multiscatterings of the resonant channels are encoded in the effective potential V k . The flow of the respective couplings and propagators is discussed in Section IV.

F. Thermodynamic potential & order parameters
The equilibrium thermodynamic potential density Ω is given by where V is the spatial volume. Ω[Φ; T, µ] is nothing but V k=0 (ρ, L,L) in (40). It depends on the given background Φ, the temperature T and the quark chemical potential µ. In the present case we only consider homogeneous background fields and the infinite volume limit. Equation (67) can be obtained from the evolution of the effective action Γ k with (26). To simplify the calculations, we split Γ k into two parts, where the glue sector corresponds to the first two loops in Figure 2, and the matter sector to the latter two. The thermodynamic potentials related to the two parts are dealt with differently in this work; that is, we do not evolve the flow of Γ glue,k , but instead replace it with the QCD-enhanced glue potential [88,89], to wit, In (69) we have introduced the traced Polyakov loops L,L with and Here, P on the r.h.s. standing for the path ordering. In (71) the gauge field is the fluctuation field and A 0 = Â 0 is the temporal mean gauge field in Φ. Accordingly, the expectation values (70) are nontrivial functions of the mean field A 0 and This has been discussed at length in [69,[128][129][130]164], for related work see e.g. [165][166][167][168][169][170]. However, in [130] it has been shown that the difference in (72) can largely be attributed to a temperaturedependent rescaling.
In the present work this difference is ignored, it will be discussed elsewhere. We approximate L[A 0 ] ≈ 1 Nc Tr P[A 0 ] which allows us to utilise pure glue lattice results for the expectation value L and the correlations of the Polyakov loop for the construction of the Polyakov loop potential. The glue potential V glue employed in this work is that computed in [171], where the quadratic fluctuations of the Polyakov loop, P[A 0 ](x)P[A 0 ](y) , are taken into account. We specify the glue potential in Appendix G. Finally, one obtains the thermodynamic potential density as follows The potential (73) allows us to access both, the confinement-deconfinement phase transition or crossover via the Polyakov loop expectation value L, where Φ EoM now includes L EoM ,L EoM instead of A 0,EoM in a slight abuse of notation. In turn, explicit and spontaneous chiral symmetry breaking is encoded in the expectation value σ := σ EoM of the σ-field, defined with the respective EoM, The expectation value σ EoM is an order parameter for chiral symmetry breaking similar to the magnetisation in the Ising model. It is closely related to the chiral condensate of the light quarks, ūu +dd /2, which is a function of the constituent quark mass of the light flavours, The potential (73) also gives us access to the thermodynamics of the system which will be discussed elsewhere.

IV. CORRELATION FUNCTIONS
In this section we discuss the correlation functions presented in (39), and derive their flow equations. This includes the propagators of all fields and the respective anomalous dimensions (Section IV A), the strong couplings related to pure glue, ghost-gluon and quark-gluon vertices (Sections IV B 2, IV B 1), the four-quark scattering and the Yukawa coupling between pions, σ and the quarks due to dynamical hadronisation (Section IV C), and the flow of the effective potential (Section IV D).

A. Propagators and anomalous dimensions
The flow equation for the (inverse) propagators, ∂ t Γ (2) k (p) is obtained by taking the second derivative w.r.t. the fields of the Wetterich equation (26). They are depicted in Figure 5 for the quark, gluon, and meson fields, respectively. We are led to the anomalous dimensions η Φi with Φ i = (A, q,q, φ). While the η A comprises the full flow of the inverse gluon propagator, the mass terms of quarks and mesons are field dependent. They are captured by the Yukawa term and the second derivative of the effective potential respectively.

Mesons and quarks
The full quark and meson two-point functions (at vanishing pion fields π = 0 and vanishing quark chemical potential µ q = 0) are given by where both the wave function renormalisations and the mass terms also depend on the chosen mesonic field φ.
The notation in (77) is that in [27,32], where the full momentum dependence of the meson and quark two-point function of two-flavour QCD in the vacuum has been computed. The running quark mass parameter in the present work is given bym l,k = M l,k (0), and this relation holds on a quantitative level, see the Figures 26 and 28 in Appendix E 3. Evidently, the fully momentumdependent wave function renormalisations Z q (p) and mass functions M q (p) of the quarks are uniquely defined as the prefactor of the Dirac tensor structure and the scalar tensor structure respectively, where tr stands for the Dirac and colour trace, and we do not sum over the flavour index i = 1, ..., 3. The different tensor structures for the kinetic term and the mass term provide us also with unique projections of the flow equation on the flows ∂ t Z q (p) and ∂ t M q (p) respectively. The latter flow is related to the flow of the Yukawa coupling and is discussed in Chapter IV C. The flow of the wave function renormalisation is encoded in the anomalous dimension (30) of the quarks. In the present work we do not introduce a thermal splitting of the wave function renormalisation in parts longitudinal and transversal to the heat bath, but use the transversal part throughout. This leads us to where no sum over the flavours i is taken. It is wellknown that the anomalous dimension of the quark carries a mild momentum dependence, for results in the present fRG setting as well as respective DSE and lattice results see [7,27,32,172]. In particular, from [27,32] we infer a very small momentum dependence of η q in the regime p 2 ≤ k 2 . Hence, we use p = 0 which facilitates the explicit expressions, On the side of theories, lattice QCD simulation is a nonperturbative first-principle approach, which has provided us with lots of remarkable insights and understanding in recent years, such as the determination of freezeour parameters through the confrontation of lattice re- ]. In this work, we would like to perform the first-principle QCD calculations at finite temperature and baryon chemical potential within the FRG approach. QCD phase transitions including the chiral phase transition and the color deconfinement phase transition will be investigated. We will also study the QCD correlation functions and their dependence on the external parameters. Moreover, a T µB phase diagram will be presented based on our computation. This paper is organized as follows. In Sec. ?? we introduce the functional renormalisation group approach to QCD. In Sec. ?? correlation functions including propagators, the strong couplings, and the dynamical hadronization are discussed. In Sec. ?? numerical results and related discussions are presented, and a summary and outlook is given in Sec. ??. Details about our calculations, such as the flow equations for the e↵ective potential and couplings, anomalous dimensions, the glue potential, numerical setup, and the threshold functions, are presented FIG. 5. Flow equations for the inverse propagators of quarks, gluons, mesons respectively, in the present approximation; in particular the tadpole contribution from two-quark-twomeson as well as the remnant four-quark scatterings are missing, as is the ghost-gluon tadpole. The derivative∂t is defined in the caption of Figure 4. The left hand sides stand for, in a slight abuse of notation, ∂tΓ where p 0 is a small, but nonvanishing, frequency, for more details see Appendix J. Its nontrivial choice is also related to fact that the expression in the square bracket is complex for nonvanishing chemical potential. This is well-understood and originates in the Silver-Blaze property: At vanishing temperature correlation functions involving a quark field q(ω) below density onset are real functions of the complex frequency variable ω − i µ q , for a detailed discussion in the fRG and 2PI approaches see [72,96,173]. Therefore we project the flow on its real part, denoted by Re. The trace in (80) projects on the Dirac tensor structure in the quark two-point function, that is its kinetic part. The explicit expression of η q,k is also given in Appendix J.
It follows from the parameterisation (77) that the mesonic wave function renormalisations and masses are given by In contradistinction to the quark wave function renormalisations and masses, that of the mesonic degrees of freedom do not allow for unique definitions as they cannot be distinguished by their tensor structures. Accordingly, one can easily shift momentum dependence from the wave function renormalisation Z φ to the mass m 2 φ in (77). Furthermore, for given φ 0 = (σ 0 , π = 0), the σ and π parts of Γ (2) φφ are different, The σ and pion masses are given by the curvature of the mesonic potential V k in the respective field directions, and are hence called curvature masses, see e.g. [95]. Evidently, the pion mass vanishes in the chiral limit as it is simply given by the EoM for c σ = 0. The two masses differ by a higher ρ-derivative term which is only present for ρ 0 = 0, that is σ 0 = 0 as π ≡ 0. The wave function renormalisations follow similarly, In the present approximation we do not consider the dependence of the wave function renormalisations on the mesonic fields. These terms amount to momentumdependent mesonic self-scatterings that vanish at vanishing momentum. These terms are suppressed at low energies which is the only regime where the offshell mesonic fluctuations can play a rôle in the first place. This leaves us with a uniform wave function renormalisation which we define via the pion two-point function, A natural definition of Z π comes from the definition of the pion pole mass see e.g. [95]. In the present Euclidean setup we are restricted to p 2 0 ≥ 0. The approximation Z π (im π,pol , 0) ≈ Z π (0, 0) ≡ Z π (0) is therefore the optimal choice for the wave function renormalisation in (87), as it is closest to the pion pole. Since this pole, in turn, is close to the origin in the region where mesons are dynamically relevant, this approximation is even quantitative if Z π (p) is only mildly momentum dependent in the small region |p 2 | m 2 π . Indeed, this has been demonstrated to be true for cutoffs 1 GeV within low energy effective theories in [95]. Hence, we arrive at which also implies that the pole mass and the (renormalised) curvature mass of the pion agree. Therefore we consider explicitly the pion wave function renormalisation, and with (85) we arrive at For momenta close to the pion pole this agrees well with the pole renormalisation. In the present work we consider p 2 ≥ 0 and hence the optimal choice is p 2 = 0 and p 0 = 0, where the latter limit is taken first. This leads us to with the pion pole mass (88). As for the quarks we have neglected the thermal splitting of the wave function renormalisation in parts longitudinal and transversal to the heat bath. Instead we use the transversal part throughout. The validity of this approximation has been checked explicitly in e.g. [121].
In the present work we also approximate the full momentum dependence of the quark and meson propagators in the diagrams by cutoff-scale-dependent wave function renormalisations and masses. The cutoff scale dependence carries the -averaged-momentum dependence at p 2 ≈ k 2 . Accordingly, this approximation provides even semiquantitative results in the flow if the momentum dependence of the propagators is small for momenta p 2 k 2 . Moreover, the spatial momentum q loop integrals of flows of correlation functions at vanishing external spatial momenta p i = 0 are peaked at q 2 ≈ k 2 for generic cutoffs, and q 2 = k 2 for the present cutoff choice. Both properties originate in the spatial momentum measure with I k (q 2 > k 2 ) = 0 for the present cutoff choice (N1). We conclude that if approximating the full Z φ (p) with (89) at p = k, the flow diagrams with meson propagators are well-approximated. This leads us to the following approximation of the full momentum dependence of the meson two point function in the right hand side of the flow equations, where (0, k) stands for (p 0 = 0, p 2 = k 2 ). In summary, best use of these consideration is made ifZ φ with (92) is used in the flow diagrams, while Z φ (0) defined in (90) is used for determining the current quark masses via the pole mass of the pion in the vacuum.
It is left to derive the flow equations for the wave function renormalisationsZ φ and Z φ (0). Both can be read of from the anomalous dimension (30) computed from the k-derivative of (89). We use that lead to the flow ofZ φ with Note that (93) and (94) imply thatZ φ = Z φ (0, k). This can be seen from the t-derivative of Z φ (0, k), where the second term comes from the spatial momentum |p| = k in Z φ (0, k). We emphasise again that the present approximation is based on the mild momentum dependence of the mesonic wave function renormalisation for p 2 k 2 of both low energy effective theories as well as full QCD. In the vacuum this has been shown in [27,32,95]. Note, however, that at finite density a nontrivial momentum dependence of the meson wave function renormalisation can be induced by a modulated spatial structure, e.g. due to an inhomogeneous quark condensate. As discussed below, we find indications for such a regime here. For the determination of the pion mass in the vacuum we also require the wave function renormalisation at vanishing momentum, Z φ (0). Its flow is read-off from (90), The explicit expression for (93) and (96) are deferred to Appendix I. Note that Z φ (0) and Z φ (0, k) agree at k = 0, whileZ φ does not. This deviation can be used for a systematic error estimate, since the wave function renormalisations only enter via η φ . Accordingly, it is the difference η φ (0, k) + ∂ t ln Z φ (0, k), that is the last term on the right hand side of (95), which is missing in the flows. We have checked numerically that this difference does not play any rôle.
In the current approximation, the pion pole mass in the vacuum is given by Similarly we approximate However, we emphasise that (98) carries a qualitatively larger systematic error in comparison to (97): Apart from the fact that the physical meaning of the σ-resonance in the scalar four-quark channel is still debated, [174,175], it requires a small momentum behaviour for a larger momentum range |p 2 | m 2 σ,pol . Nonetheless, in an extension of the current setup to Minkowski frequencies the respective observable provides the position of the lowest resonance in the scalar channel.
At finite chemical potential the situation changes significantly: The mesonic dispersion develops a minimum for nonvanishing momentum. This entails that Z φ (0, k) and Z φ (0) may and do differ qualitatively. Moreover, such a behaviour indicates an inhomogeneous regime, for more details see Section V B and Figure 21. Still, we have checked that this does not affect the purely fermionic observables.

Gluons and ghosts
It is left to specify the gluon and ghost anomalous dimensions. The ghost propagator has been shown to be rather insensitive to quark contributions for N f = 2 and N f = 2 + 1 flavours as well as finite temperatures T 1 GeV, see e.g. [33] and references therein. For this reason we use the ghost anomalous dimension at vanishing temperature from [32], Note that (99) only enters explicitly the ghost triangle in the vacuum flow of the three-gluon vertex depicted in Figure 6, see Section IV B 2. Implicitly, η c is also present via the input from [32,33] and in the temperature and density fluctuations of purely gluonic two-and three-point functions.
For the gluon anomalous dimension we decompose the flow for the inverse gluon propagator into the pure glue diagrams and the quark loop, see Figure 5. As discussed in Section III B, we utilise quantitative results for vacuum two-flavour QCD, [32], and finite temperature Yang-Mills theory, [33]. Then, the gluon anomalous dimension η A,k at finite temperature and density is decomposed into a sum of three parts, The first term on the r.h.s. of the equation above accounts for the vacuum contribution to the gluon anomalous dimension. In turn, ∆η A A and ∆η q A account for the medium contributions to the gluon anomalous dimension, from gluon loops and quark loops respectively.
For N f = 2 we infer it directly from the corresponding gluon dressing function Z QCD A,k=0 (p) from [32] with Full RG-invariance of the procedure is then obtained by rewriting the anomalous dimension as a function of the running coupling α s,k as done in [28,131]. This is described further in Appendix D. In the present work we simply adjust the input coupling α s,Λ consistent with its implicit value given in η A : We choose α s,Λ such that the quark-gluon and purely gluonic couplings show the same ultraviolet running. Their running is proportional to 1/2η A and 3/2η A respectively, the rest of the β-functions is given by diagrams proportional to α s,k . Only for the consistent initial coupling both runnings can agree as they should. For more details see Section IV B 3. We also emphasise that instead of the fRG data from [32] we also could have taken lattice data or other data from other functional approaches such as the DSE. As mentioned before, it is an important feature of the current setup that results obtained within other approaches can be systematically included. This allows for systematic improvements and hence enhances the reliability of the current approach.
For N f = 2 + 1 we include the contribution of the s-quark via its flow. Then we use with The transversal projection operator Π ⊥ in (103) is defined in (N2). The superscript (s) denotes the strange quark contribution. As in the two-flavour case the initial coupling α s,Λ is adjusted by RG-consistency. For more details see Section IV B 3. Now we proceed to the gluon anomalous dimension at finite temperature and density. The difference to the vacuum anomalous dimension is comprised in the second and third term on the right hand side of (100). Here, ∆η q A,T in (100) denotes the contribution of the quark loop at finite temperature and quark chemical potential µ q . With (64) it reads where the superscript (q) indicates the contribution of the different quark flavours. The vacuum contribution is subtracted in (104), and we have used the transverse magnetic tensor in the projection, Moreover, the contraction of Lorentz and group indices is implicitly understood in (104). The explicit expression for ∆η q A can be found in (H1).
The thermal contribution of the glue sector is encoded in the second term on the r.h.s. of (100), i.e., ∆η glue A . It has been discussed in [33] that it is well-captured by a thermal screening mass, its inclusion is described in detail in Appendix H. Note that the present setup does not take into account the full backcoupling of the thermal and density fluctuations of the quark in the pure glue diagrams. A respective error estimate is discussed in Appendix D.
There it is shown that this approximation, the 'direct sum' of the contributions, already works semiquantitatively for the much bigger backcoupling effects related to adding the quarks to pure glue in the vacuum. In the present work we take the vacuum backcoupling fully into account, and only apply the direct sum approximation to the much smaller thermal and density contributions. This suggests a quantitative nature of the current setup.

B. 'Avatars' of the strong coupling αs
In the present truncation we consider different 'avatars' of the strong coupling, introduced in Section III D, as well as the Yukawa coupling between quarks and the scalarpseudoscalar mesons, introduced via dynamical hadronisation in Section III E.

Quark-gluon couplings
The strong couplings and their running with the RG scale, especially the quark-gluon couplings gq Aq = (αl Al , αs As ), play a significant rôle in dynamical chiral symmetry breaking, see e.g., [27,28,125,176] for more detailed discussions. As we discuss in III D, different strong couplings are identical in the perturbative regime due to the mSTIs, however, they deviate from each other in the nonperturbative or even semiperturbative regime [32]. In Figure 6 we show the flow equations for the quark-gluon and three-gluon couplings. The first and third diagrams on the r.h.s. of the equation in the first line are the usual QCD contributions, while the second one arises from the mesonic fluctuations. Thus the relevant flow equation for the quark-gluon coupling is given by the projection of the flow of Γ Aqq on the classical tensor structure T Aqq /g. This leads us to where the trace sums over Dirac indices and the fundamental representation of gauge group. The classical tensor structure is given by where the coupling has been divided out. The contracted flow in (106) is normalised with the trace of the classical tensor structure squared, In (106), {p} denotes the set of external momenta for the three-point vertex. The momentum for the gluon is chosen to be vanishing, and that for the quark as illustrated in Appendix J. We divide the term in the second line of (106) into two parts, i.e., ∂ t gq Aq A and ∂ t gq Aq φ , which correspond to the relevant contributions from the quarkgluon and quark-meson couplings, respectively. Their expressions are given in (M2) and (M3).
It is left to discuss the effect of the current approximation with only the classical tensor structure. In quantitatively reliable approximations to QCD the initial conditions at a perturbative initial scale Λ 10 GeV, as described in Appendix E 1, lead to QCD physics at vanishing cutoff scale k = 0. In particular, this holds for the physics of spontaneous chiral symmetry breaking, see [27,32] for quenched and unquenched N f = 2 flavour results respectively. From these computations we also know that not only the classical tensor structure of the quark-gluon vertex carries the fluctuations important for chiral symmetry breaking. In turn this implies, that an approximation with only the classical vector tensor structure may lack interaction strength. This typically calls for a phenomenological infrared enhancement of the quark-gluon coupling αq Aq that compensates for the missing tensor structures. This is common to all functional approaches within approximations that lack full quark gluon vertices, for related discussions in DSEs see Section VI and and e.g. [35-37, 39, 63, 132-135, 172], for a recent review see [7].
In the present work we follow the approach proposed in [28], where the infrared enhancement is directly introduced in the infrared part of the flow. More details can be found in Appendix E 2. Here, we simply summarise the results: In the present setup it turns out that the required infrared enhancement is rather small, only for cutoff scales smaller than 2 GeV the flow is enhanced with a factor 1.034, see (E9) and (E10) in Appendix E 2. While being small, its effect is significant as the constituent quark mass of the light quarks would be m q = 217 MeV without the enhancement. Thus, there is only a fine line between 'too little' (or no) and 'too much' chiral symmetry breaking in QCD. This has already been observed in [27,32], where a quantitative approximation to two-flavour QCD in the vacuum was studied. See also [177] for a related two-flavour DSE study with N c = 2. This peculiar behaviour of first principle QCD which may also has significant consequences for the interplay of confinement and On the side of theories, lattice QCD simulation is a nonperturbative first-principle approach, which has provided us with lots of remarkable insights and understanding in recent years, such as the determination of freezeour parameters through the confrontation of lattice re- ]. In this work, we w like to perform the first-principle QCD calculation finite temperature and baryon chemical potential wi the FRG approach. QCD phase transitions including chiral phase transition and the color deconfinement p transition will be investigated. We will also study FIG. 6. Flow equations for the three-point functions computed here: quark-gluon, three-gluon and Yukawa couplings. In the present approximation, the diagrams with higher order vertices and the remnant four quark vertices are dropped.∂t is defined in the caption of Figure 4. chiral symmetry breaking at finite temperature and density, can be phrased as: QCD is living on the edge. In [27,32], in particular all tensor structures and momentum dependences of the quark-gluon vertex have been included selfconsistently. There it has been observed that the nonclassical tensor structures, in particular the chiral symmetry breaking ones, contribute significantly to the size of dynamical chiral symmetry breaking.
Still, we interpret the smallness of the enhancement as a sign of robustness and quantitative reliability of the present setup: At finite µ B the coupling decreases and we are sensitive to the functional form of how the infrared enhancement is switched of at larger momentum scales. Note that the latter scale also includes temperature and chemical potential. Accordingly, the smaller the infrared enhancement is the less influence this can have. In turn, in the presense of a significant infrared enhancement this calls for a systematic study of its shape.

Gluonic couplings
The gluonic couplings includes the three-and fourgluon couplings g A 3 , g A 4 , and the ghost-gluon coupling gc Ac . In the vacuum we adopt the approximations: g A 4 = g A 3 and gc Ac = gl Al , see also (56). This is consistent with the results obtained in [28,32] for cutoff and momentum scales k 1 GeV. For smaller cutoff scales k 1 GeV the purely gluonic couplings, g A 3 , g A 4 , decay rapidly and the relative difference is subleading. In turn, the ghost-gluon coupling gc Ac approaches a sizable value (αc Ac (k = 0) ≈ 3). However, in the present setup, the ghost-gluon coupling only contributes explicitly to the vacuum part of the flow of the three-gluon vertex in the ghost triangle in Figure 6. This is a subleading diagram except in the deep infrared for k 100 MeV. In this regime it leads to the zero-crossing of the three-gluon vertex and a small negative IR value of the coupling. This infrared limit is not reproduced in the present approximation, but this deep-IR property has no effect on the present results.
This leaves us with the quark-gluon couplings αq Aq with the flows (106) at finite temperature and density, and the three-gluon coupling α A 3 . We strive for RGconsistency also at finite temperature and density and identify the thermal contributions of the quark-gluon and ghost-gluon couplings. For the computations we use that of the quark-gluon coupling, (106), which is the most relevant for the chiral dynamics of the theory. This entails for the medium contribution ∂ t ∆g The vacuum part has been computed in [28] and we use the respective flows The momentum configuration of Flow A 3 (p, −p) is chosen such that the three external gluons carry incoming momenta p, −p and 0 respectively. The contraction with δ µν p ρ f abc and the normalisation comes from the projection on the classical tensor structure T (1) As in (107), the strong coupling has been divided out. For the diagrammatic representation of the quark-, gluon-, ghost-triangle diagrams and the tadpole diagram from the four-gluon coupling, see Figure 6.

RG-consistency of the 'avatars' of the strong coupling
The flows computed in the present work utilise input from the quantitative N f = 2 flavour vacuum QCD [32] and finite temperature YM computations [33]: Parts of the gluon and ghost anomalous dimensions η A and η c are taken from these works, see Section IV A 2. The input is given in terms of the running scale k in either N f = 2 flavour vacuum QCD or finite temperature Yang-Mills theory. The current setup differs in several aspects. First of all, we also consider N f = 2 + 1 flavours. Moreover, the regulators chosen in the present work differ from that in [32,33]. These differences are fully accounted for if one resolves the anomalous dimension as a function of α s and further parameters, that is η A,k = η A (α s,k ). This has been put forward in [28,131], see also Appendix D.
This scale matching entails in particular that the strong coupling at the initial scale has to be chosen consistently with that implicit in η A (α s ). In the present work we guarantee the above RG-consistency in the following sense: In the perturbative regime the β-functions of the different 'avatars' of the strong coupling, i.e. the three-gluon coupling (54) and the quark-gluon coupling (55) have to agree. Their flows are given by (110) and (106) respectively. In the perturbative regime we have g A3 = gq Aq = g s , which entails schematically where 'diagrams(α s )' stands for the difference of the diagrams for the three-gluon and quark-gluon vertices depicted in Figure 6. The constraint (112) is the requirement of RG-consistency of the couplings. For a given initial cutoff scale k = Λ it defines the consistent initial condition for the strong coupling. In particular, it adjusts for the relative RG-scaling in N f = 2 and N f = 2 + 1 flavour QCD. The couplings αq Aq , αs As , α A 3 determined by this procedure agree well with each other for cutoff scales k 5 GeV, and hence are RG-consistent, see Figure 17. For smaller cutoff scales, k 5 GeV the different couplings start to separate which is also seen in [32]. This is expected and necessary, as for smaller scales we enter the nonperturbative regime which is influenced by the confining gluon gap. Moreover, the initial coupling α s,Λ derived from (112) agrees quantitatively with that in [32]. Both properties constitute nontrivial reliability tests of the current approach. Moreover, in Appendix D it is shown that the procedure adapted here only leads to small quantitative effects even if going from Yang-Mills theory to N f = 2 flavour QCD, which requires a far larger scale correction. In summary the combined checks suggest a small systematic error of the current procedure, the attraction of which is its simplicity.

C. Yukawa coupling
In the full effective action without approximation the Yukawa coupling in (39) is field-dependent, h = h(ρ), and the ρ-dependence takes into account higher scattering processes of quarks with the scalar-pseudoscalar mesonic channel [70]. This already entails that strictly speaking we have to deal with two Yukawa couplings at the expansion point φ 0 = (σ 0 , π = 0), Both vertices are derived from h(ρ) and are present in the flow equations for correlation functions. In the present work we consider the ρ-dependence of the Yukawa coupling as subleading and identify h σ = h π = h. The corresponding error is minimised by identifying h = h π as we have three pions and only one radial mode σ. Moreover, the pions are lighter except in the scaling regime in the vicinity of the critical point. Hence already the diagrams with one pion mode give bigger contributions than those with the σ. In the critical region the σ, as the critical mode, becomes massless. However, this regime is exceedingly small and a detailed analysis of its features will be presented elsewhere. This leaves us with the flow equation of h = h π , or rather that ofh, wherem 2 π (ρ) =V (ρ)/k 2 , and the subscript (qτ q)π indicates the projection on the pseudoscalar part of the Yukawa coupling. The diagrams contributing to Flow (3) (qτ q)π are shown in the third line of Figure 6. Note that, as discussed in Section III E, the hadronisation functionȦ explicitly enters the flow of the Yukawa coupling. This highlights the fact that dynamical hadronisation, in contrast to conventional bosonisation, stores the full dynamical information of the four-quark correlation in the constituent-quark-meson sector. The ratioh 2 /(2m 2 φ ) then is identical toλ q (without dynamical hadronisation) in the chirally symmetric phase, reflecting the smooth transition from fundamental to composite degrees of freedom in QCD.
However, instead of using the flow equation for Γ (3) qqπi [Φ 0 ], it is simpler to use the fact that the scalar part of the quark two-point function is proportional to Z q σ 0 h(ρ 0 ). Note that this only holds true within a dynamical hadronisation of the theory which keeps the formulation maximally symmetric: The only term that breaks explicitly chiral symmetry are the linear terms in σ and σ s . As shown in Appendix B this impliesĊ k ≡ 0 in (15). Then, the flow ofh can be deduced from that of the quark two-point function, see Figure 5, with where the subscript (qτ 0 q) indicates the projection on the scalar partqτ 0 q = (ūu +dd)/2 of the quark two-point function. More details can be found in Appendix B. Both, (114) and (115) provide the flow ofh. This is different from the flows of h π and h σ , which differ by a term proportional to h (ρ). In the present case we use (115) for ∂ th , as it has the simpler diagrammatic part.

D. Effective potential
The effective potential V k (ρ, L,L) in (39) receives direct contributions from the quark and meson loop in Figure 2. Via the quark loop of the flow equation for the gluon two-point function the gluon propagator also develops a φ-dependence. The latter φ-dependence then also propagates to the ghost loop. However, both dependences are negligible and are dropped in the following. This leaves us with a flow for the effective potential which only receives contributions from the quark and meson loop in (26), see also Figure 2. Owing to the coupling between the gluon background field A 0 and the quarks, the quark loop also induced a dependence on the Polyakov loop to the effective potential. The flow equation for ∂ t V k , (F1), is discussed in detail in Appendix F.
Here we only emphasise the important parts of Appendix F. In the present work we resort to a Taylor expansion of the effective potential. In general this can be done about a k-dependent expansion point, whereρ is the renormalised field (42). On the right hand side of (116) we have suppressed the dependences on L,L of the expansion coefficients: λ n,k = λ n,k (L,L) and κ k = κ k (L,L). Moreover, in (116) we have implicitly assumed its convergence. It has been shown in [70] that the most rapid convergence of this expansion is achieved if the expansion point κ is kept fixed, for applications see also [28,71,73]. Then the only k-dependence of the expansion point in (116) comes from the mesonic wave function Z φ .
In the present work we resort to an expansion about the flowing solution of the quantum equation of motion, Using the flowing expansion point (117) facilitates the access to the expectation value of the σ-field. The expansion is described in more detail in Section F. The price we pay for this simple access to the EoM of the σ field is that this expansion shows bad convergence properties for large chemical potential. This is but one of the reasons why our current analysis is limited to µ B /T 6. For general considerations concerning the stability of this expansion scheme see [70,178]. Respective k-independent expansion schemes, full effective potentials and a detailed stability analysis will be considered in future work.

V. RESULTS AND DISCUSSION
In this section we discuss a selection of numerical results obtained in the present setup. We concentrate on N f = 2 + 1 flavour QCD, and use the corresponding results for N f = 2 only for the comparison of the N f = 2 and N f = 2 + 1 flavour QCD phase structures. Confinement and chiral symmetry breaking is discussed in Section V A. Indications for an inhomogeneous phase at large chemical potential are shown in Section V B. The strong couplings and propagators at finite T and µ B are discussed Section V C. In Section V D we present phase structure of N f = 2 and N f = 2+1 flavour QCD at finite temperature and density.
In practice, we numerically solve the coupled set of flow equations discussed in the previous section together with the initial conditions discussed in Appendix E. As it has been discussed extensively in the literature, e.g. [27-29, 32, 61, 125, 179-181], the running couplings in the bound state sector are governed by an IR attractive fixed point in the chirally symmetric regime. As a consequence, if we initialise the flow in the perturbative regime, the initial conditions are fixed by the initial values of the strong coupling and the quark masses. As long as the initial meson masses are larger than the cutoff scale, our results are independent of the choice for the initial values in the meson sector, i.e. independent of any low-energy effective parameters. The low-energy sector emerges dynamically in our case, and the corresponding parameters cannot be tuned, but are predicted.
Note that the mutual coupling of the different n-point functions through their respective flow equations leads to a resummation of infinitely many diagrams upon integration of these equations. The renormalisation and regularisation provided by the fRG gives us access to the nonperturbative IR regime of QCD.
A. Chiral and confinement-deconfinement order parameters & signatures for the CEP Observables and order parameters are derived from fields with nonvanishing expectation values: The scalar condensatesσ EoM andσ s,EoM are directly related to dynamical chiral symmetry breaking, and the Polyakov loop and its conjugate, L,L, are related to confinement, as introduced in Section III F. Note however that these order parameters are not observables themselves.
In particular, the mesonic fields φ = (σ, π) are introduced as low energy effective fields that carry the dynamics of the scalar-pseudoscalar four-quark interaction. They carry the same quantum numbers as the related mesons but their direct physics interpretation in terms of onshell mesons has to be taken with a grain of salt. Nonetheless, the renormalised light chiral condensate and the reduced chiral condensate, ∆ l,R and ∆ l,s respectively, are given in terms of the unrenormalised expectation values σ EoM and σ s,EoM , for more details see Appendix A. In the present work we concentrate on the renormalised light chiral condensate ∆ l,R , with a normalisation −c σ /(2N R ) which is fixed to the lattice, see also (A3). The reduced condensate is also discussed in Appendix A, but carries a larger systematic error in the present approximation. In turn, the Polyakov loop is an order parameter for confinement in pure Yang-Mills theory, but does not carry a direct physics interpretation in QCD with dynamical quarks. Despite these deficiencies both,σ EoM and (L,L) EoM play a distinguished rôle for the low energy dynamics of QCD.
In Figure 7 we  breaking triggered by the scalar-pseudoscalar four-quark interaction. We see that both transitions are crossovers for small baryon chemical potential as is well-known due to explicit chiral and center symmetry breaking, see e.g. [182]. They are getting sharper and shift to smaller temperatures with increasing µ B . As L andL are related to the free energy of a single quark and antiquark respectively, they are different at finite µ B : Heuristically speaking, it is easier for an antiquark to propagate in a system with an excess of quarks, and henceL is larger than L at finite, positive µ B .
As discussed in particular in Section III E, the fourquark interaction plays a central rôle for chiral symmetry breaking. The physical picture is that, as the strong coupling increases in strength with decreasing energy scale, the scattering of quark-antiquark pairs, mediated by gluon exchange, becomes stronger and eventually resonant. If this resonance occurs in the channel which carries the quantum numbers of a scalar quark condensate, qq , chiral symmetry is spontaneously broken. Dynamical hadronisation allows us to accurately capture the dynamics of the resonant quark-antiquark scattering channels in both the symmetric and the spontaneously broken phase. As argued in Section IV C, the resulting effective four-quark interaction of the scalarpseudoscalar channel is the given by the ratioh 2 /(2m 2 π ). In Figure 8 we show it as a function of the RG scale k for several values of the temperature. Starting from large k, the coupling rises as expected from the increasingly strong quark-antiquark scattering. At T = 0 and around k ≈ 400 MeV, the increasingly steep rise indicates the resonance. The coupling turns into a pion-exchange interaction below this scale through dynamical hadronisation. Since the pion mass and the pion-quark Yukawa coupling are essentially constant in the chirally broken phase (cf. Figures 11 and 16 approximately constant as well. The thermal suppression of the strong coupling discussed below leads to a suppression of the four-quark interaction with increasing temperature. As a consequence, the resonance disappears and chiral symmetry is restored at sufficiently high temperature. This is discussed in more detail in [176]. In Figure 9 we show the effective four-quark interactionh 2 /(2m 2 π ) as a function of T for k = 0 for various chemical potentials. The insensitivity to µ B at vanishing temperature is a consequence of the Silver Blaze property discussed in Section IV A 1. For correlation functions below onset density, it implies that their dependence on µ B is given solely by a µ B -dependent shift of the frequencies of external legs that carry baryon number. This is discussed in detail in [96]. For the present case, see Appendix J. For increasing temperature, the four-quark interaction melts down earlier and more rapidly with increasing µ B . Consequently, the chiral transitions occurs also more rapidly and at lower temperatures.
The order parameter for chiral symmetry breaking used in the present work is the renormalised light chiral condensate ∆ l,R , see Appendix A, in particular (A7) and Figure 22. In Figure 10 (left panel) we show ∆ l,R as a function of temperature at various baryon chemical potentials. At µ B = 0 we compare our result to the continuum extrapolation of lattice results at the physical point from [183] and find very good agreement. The corresponding thermal susceptibility, ∂∆ l,R /∂T , is shown in the right plot of Figure 10. We use the peak position to define the pseudocritical temperature. We note that the pseudocritical temperature T c = 156 MeV at µ B = 0 is in quantitative agreement with the lattice prediction. Similarly, one can use the inflection points ofm l and the Polyakov loop in Figure 7  chemical potential. While these are not unique definitions of the crossover temperature, they match well with each other. For a more detailed comparison of chiral order parameters in the crossover regime see e.g. [70]. T,µ B =0 . As a result, z φm coincides with the pole mass defined in (88) and (98) in the vacuum. The in-medium behavior is that of the curvature masses defined in (43). For more details see Section IV A 1.
to the chiral condensate is clearly visible. Furthermore, sincem σ is directly related to the inverse correlation length of the system, it has a dip at the chiral crossover. The closer the system is to the CEP, the lighter the sigma becomes at the transition until it eventually dips below m π . This indicates that the system enters the critical region, as visible in Figure 11 for µ B = 630 MeV. At the CEP, the sigma is exactly massless since it is the critical mode of the chiral phase transition.
We also find that in the chirally symmetric phase m π and m σ become degenerate and grow rapidly with the increase of T . This entails that the mesonic degrees of freedom decouple quickly when the temperature of the system is above T c . In this regime the dynamics are governed by the fundamental degrees of freedom of QCD. This will be discussed in more detail below. We note that the inclusion of only the resonant scalar-pseudoscalar interaction channel in the two-flavour subspace of 2+1 flavours entails that U (1) A is maximally broken. Hence, the mesons in the light sector that receive a positive contribution to their mass from the axial anomaly, η and a 0 , are fully decoupled here. We defer a more complete study of chiral symmetry restoration, including the anomalous effects, to future work.

B. Inhomogeneous regime at large chemical potential
At finite chemical potential the situation changes significantly: The mesonic dispersion develops a minimum at nonvanishing momentum. This is signaled by Z φ (0) at k = 0 as defined in (85) and (90)   tive. Indeed this happens for baryon chemical potentials µ B 420 MeV for a temperature regime that widens with increasing chemical potential, see Figure 12. The left panel shows Z φ (0) as a function of temperature for different baryon chemical potentials, while the right panel shows 1/|Z φ (0)|, which highlights the negative regime bounded by the spikes. This property entails that the mesonic two-point function develops a negative slope at vanishing momentum implying a minimum at nonvanishing spatial momentum, p 2 min = 0. This signals a spatially modulated or inhomogeneous regime. Vanishing of the wave function renormalisation, as the coefficient of p 2 in the meson propagator, can be indicative of an instability towards the formation of an inhomogeneous quark bound state, such as a chiral density wave. Whether or not inhomogeneous condensation occurs cannot be answered on the basis of only the meson wave function renormalisation. Such an analysis requires a systematic study of the competing effects of potential homogeneous and inhomogeneous resonances in quark-antiquark scattering channels. This will be deferred to future work. The relation between a minimum of the dispersion relation of bound states and inhomogeneous phases has been discussed in detail in the past decade within low energy effective theories at and beyond mean field. For reviews see e.g. [184][185][186][187][188], for applications with the functional RG see e.g. [111,[189][190][191].
This regime is indicated in the phase diagram in Figure 21. The blue shaded area corresponds to the region with Z φ (0) < 0. The hatched red area denotes where this region overlaps with a sizable homogeneous chiral condensate. Note that the chiral phase boundary lies within this regime for baryon chemical potentials µ B 500 MeV. Hence, in particular in the hatched red region a competition between homogeneous and inhomogeneous quark-antiquark pairing is expected. This competing order effect can alter the phase structure. The hatched red region should therefore also be interpreted as a systematic error on our phase boundary. Accordingly, the present approximation should be upgraded by inhomogeneous four-quark interactions at larger µ B . Most notably, the CEP lies in this regime.
The connection between the CEP and inhomogeneous phases has been investigated within low energy effective theories, e.g. [189]. It has been shown that the relative position of the onset of the inhomogeneous regime and the chiral CEP is related to the σ mass gap, m σ , and the constituent light quark mass, m l . For m l = m σ /2 the two onsets agree on the mean field level. In this case the phases with restored chiral symmetry, the homogeneous phase and the inhomogeneous phase meet at a singe point, the Lifshitz point. For m l ≥ m σ /2 the location of the (potential) chiral CEP is at larger baryon chemical potential as that for the onset of the inhomogeneous regime. In the present case we have m l ≥ m σ /2 with m l /m σ = 347/510 ≈ 0.68 in the vacuum. However, it has been argued in [192] that strong IR fluctuations wipe out the Lifshitz point. This could even lead to the destruction of the CEP. We emphasise that a detailed investigation of this regime goes beyond the scope of the present work but will be discussed in a forthcoming publication. However, we emphasise again that our observation of a minimum at nonzero spatial momentum in the mesonic dispersion is a strong indication for very interesting physics in this regime.  [193], based on the tmfT configurations [194], for the tmfT simulation setup see also [195,196]. The lattice results depicted here are obtained for β = 2.1, a = 0.0646 fm and an approximate pion mass of mπ = 369(15) MeV leading to a pseudocritical temperature of T deconf = 193 (13)(2) MeV, for more details see [193].

C. Dynamics and sequential decoupling
The relevance of the different fields for the dynamics of the system at different energy scales is encoded in the strength of the vertices and the propagator dressings. For example, we have already seen, that the effective four quark couplingh 2 /(2m 2 π ) decays rapidly for large momentum or cutoff scales, see Figure 8 and Figure 9. Accordingly, the quark-self coupling gets irrelevant in the UV. Moreover, for increasing temperature and chemical potential the effective four quark coupling decreases further. In the present parameterisation with dynamical hadronisation this decoupling is solely triggered by the increase of the pion mass functionm 2 π , for the temperature dependence of the mesonic mass functions see Figure 11. In turn, the Yukawa interaction does not show a significant momentum scale and parameter dependence, see Figure 16 and discussion there. Similar observation apply to the gluonic sector. All avatars of the strong coupling show a rather strong scale dependence, as do the gluon and meson dressings. This will be discussed in the following. In combination, the different dynamics and scale dependence of the couplings, wave function renormalisations and mass functions allow us to provide a consistent picture of a sequential decoupling of gluon, quark and mesonic degrees of freedom with decreasing momentum and cutoff scales. This leads to the natural emergence of low energy effective theories (LEFTs) of QCD for scales below ∼ 1 GeV and is discussed at the end of the present Section.
We start our analysis of the dynamics and decoupling properties of the system with the gluon dressing. Here we also provide a further benchmark test of the present approach at vanishing baryon chemical potential: In Figure 13 we compare our results at finite T and µ B = 0 to lattice results on the magnetic gluon dressing function from [193] based on the tmfT configurations [194]; for this simulation setup see also [195,196]. Note that the lattice results are obtained for N f = 2 + 1 + 1 flavours, with β = 2.1, a = 0.0646 fm and an approximate pion mass of m π = 369(15) MeV leading to a critical temperature of T deconf = 193(13)(2) MeV. Hence, agreement with our results can only be expected at large temperatures and the corresponding spatial momentum scales. In this regime, the vacuum constituent quark masses are small against their thermal masses and the perturbative running of N f = 3 vs N f = 4 is not yet dominant. Furthermore, the effects of the additional quark flavour and the different masses are subleading. In order to facilitate this comparison, we rescaled our momentum scale by |p| → 1.07|p| as to match our results in the semiperturbative regime around p ≈ 2 GeV. This is necessary since, owing to the different mass scales, also the momentum scales of our results and the lattice results are inherently different. Indeed, we find good agreement for the temperatures, 191, 218, 254, 277 and 305 MeV, for 0.8 GeV |p| 3 GeV. For larger momenta, our result show the correct leveling-off expected from the logarithmic momentum dependence in the perturbative regime. For smaller momenta, both the different mass scales and the different infrared gauge fixing come into play. Note also that the lattice data are subject to lattice artefacts at large momenta, as is the case in the N f = 2 vacuum data, see Figure 25.
In Figure 14 we show the gluon dressing function 1/Z A (|p|) = 1/Z A,k=|p| at vanishing Matsubara frequency p 0 = 0 as a function of spatial momenta for different temperatures (left panel) and baryon chemical potentials (right panel). While the gluon dressing shows a mild dependence on the temperature, it is essentially independent on the baryon chemical potential. Note that the respective dependences are more significant in the gluon propagator G A = 1/(p 2 Z A (p)) shown in Figure 29 in Appendix H. However, we emphasise that the flows of the couplings and observables computed in the present work depend on the gluon dressing rather than on the gluon propagator.
The suppression of the thermal and chemical potential dependences discussed above are linked to the mass gap of the gluon propagator in the Landau gauge. This mass gap signals confinement, see [69,[128][129][130]164], and manifests itself through the nonmonotonicity of the curves depicted in Figure 14. The backbending of the gluon propagator at small momenta also signals positivity violation [197], and results in a negative spectral function at low frequencies [198]. With increasing temperature above T c , the chromoelectric part of the gluon is subject to increasingly strong Debye-screening. As discussed in Section IV A 2, we define the gluon dressing via the chro- momagnetic contribution. Hence, the thermal suppression in the left plot of Figure 14 is due to nonperturbative magnetic screening. For more technical details on the inmedium gluon propagator we refer to Section IV A 2 and Appendix H.
Overall, the temperature dependence of the gluon dressing is a subleading effect in the temperature region relevant for the phase transition. Below the critical temperature it is mostly affected in the infrared, where the mass gap suppresses gluon fluctuations anyway. While thermal corrections are also present in the pure gauge theory, density corrections are only triggered by the quark contributions discussed in Section IV A 2. We find that the gluon dressing function is almost in-sensitive to µ B . This is shown in the right plot of Figure 14, and also in Figure 29 for the propagators itself. For temperatures in the vicinity of the crossover, we find a slight enhancement of the gluon propagator with increasing µ B for momenta below ∼ 1 GeV. The reason is that the quark contribution to the gluon propagator, i.e. the vacuum polarisation, is suppressed at finite µ B , which can be read-off from the explicit expression in (H1). This counteracts the color screening of quarks. These findings support the qualitative or even semiquantitative approximations in functional applications (fRG and DSE), where either both or part of the thermal and density dependence of the gluon dressing is left out.  N1a) and (N1b), Z q,k and Z φ,k enter the system of flow equations only indirectly through the corresponding anomalous dimensions, (30). Consequently, the flows of the wave function renormalisations do not have to be integrated and their initial values are irrelevant. We have chosen Z k=Λ = 1 for convenience, see (41). One observes that there is a bump in Z q as a function of T during the phase transition, and it evolves into a sharp peak with the increase of the baryon chemical potential. This indicates that constituent quark fluctuations are enhanced in the vicinity of the phase transition at finite µ B . We also find that the meson wave function renormalisation decreases after chiral symmetry is restored. This effect drives the decoupling of mesons from the physical spectrum above T c , which is also visible in Figure 11 and discussed in the previous section.
Interestingly, the renormalised Yukawa coupling, (43), is relatively stable over all scales. Accordingly, the quarkmeson vertex running counterbalances the strong scale dependence of the meson wave function renormalisation. This is explicitly seen in Figure 16: On the left side we show the temperature and chemical potential dependence of the physical Yukawa coupling,h k=0 . On the right side we show its running with k at different temperatures. The Yukawa coupling is related to the expectation value of the renormalisedσ-field and the constituent mass of the light quarks viah k=0 = 2m q / σ . In the symmetric phase, the fast decrease ofh with increasing T triggers a more rapid melting ofm q as compared to σ . It is evident from the running of the Yukawa coupling in the right panel of Figure 16 that thermal corrections only set in if the RG scale is smaller than the temperature scale, k 2πT . For larger k and the 3d spatial regulators used here, see Appendix N, they are exponentially suppressed with exp(−k/(2πT )) and there is only the vacuum running. For a respective discussion of the thermal properties including also other regulators we refer to [199]. Furthermore, we have explicitly checked that our results are independent of the initial value of the Yukawa coupling: The Yukawa coupling is an irrelevant coupling that is generated by the flow. The memory of its initial condition is washed-out by an IR-attractive fixed point in the regime of small strong coupling. This is in line with the discussion at the beginning of this section.
This stability of the renormalised Yukawa coupling together with the rapid rise of the meson mass parameters m σ ,m π entails the decoupling of the scalar-pseudoscalar mesonic channel of the four-quark interaction in the chirally symmetric phase at large scales, discussed in Section V A.
For larger scales, k 1 GeV, the gluon exchange interactions are dominating the system, while for smaller scales, k 1 GeV, the gluons decouple. This is clearly seen in Figure 17, where we show the running of the different strong couplings: The quark-gluon and the threegluon couplings. We have distinguish between the quarkgluon coupling for u and d quarks, denoted by αl Al , and for the s quark, αs As . The quark-gluon couplings αl Al and αs As match for momentum scales k above ∼ 1 GeV. In turn, the strange quark coupling is slightly smaller at lower scales. This is expected since quark number conservation of the strong interactions implies that there is no flavour mixing in the quark-gluon interactions. Hence, only strange quarks contribute to the leading quantum corrections of αs As , resulting in a suppression relative to the leading corrections to αl Al . The quark-gluon couplings and the three-gluon coupling agree with each other, and are also consistent with perturbation theory, for scales above ∼ 5 GeV [28,32]. One also finds that for scales k 5 GeV, the light flavour coupling αl Al and the strange coupling αs As grows bigger than α A 3 , which is qualitatively consistent with the calculation in [32]. For scales k 1 − 2 GeV, αl Al , αs As and α A 3 deviate from each other pronouncedly, and all strong couplings are suppressed significantly. This behavior in the IR region is due to the running of the gluon dressing function shown in the left panel of Figure 14 and our definition of the strong couplings in (54) and (55). Since they effectively describe gluon exchange, the suppression of glue dynamics in the nonperturbative regime due to the dynamically generated gluon mass gap is transferred to the RG invariant strong couplings. Apparently, the three-gluon coupling is more suppressed than the quarkgluon vertex in the IR, since more gluons are attached to the former, which is verified in Figure 17. The dependence of the strong couplings on the temperature is consistent with our expectation inferred from the results of the gluon propagator. The strong couplings decrease with the increase of the temperature, which indicates the interaction between gluons and quarks gets weaker at high temperature. We have also investigated the dependence of the strong couplings on the baryon chemical potential, and find a small µ B -dependence in the µ B region of our interest.
We close this Section with a discussion of sequential decoupling and the natural emergence of low energy effective theories (LEFT) in the present fRG approach to QCD for low cutoff scales k 1 GeV: The flows of matter vertices and propagators (quarks and mesons) are driven by the tree-level four-point single field ex-change couplings with either quark or meson legs in the current approach. We emphasise that this is not due to the approximations used here but originates in the one loop completeness of the flow equation for the effective action, (26). Relevant examples are provided by the flows of the quark-gluon and Yukawa couplings depicted in Figure 6, and the flow of the four quark coupling depicted in Figure 4. These flows contain all single field exchange couplings considered in the present approach: The four-quark single gluon exchange coupling and the four-quark single meson exchange coupling. Note that the potential third single field exchange coupling, the two-quark-two-meson single quark exchange coupling simply comes from a reordering of the building blocks in the diagrams. It also originates in the fundamental building blocks in the QCD matter sector considered here, the scalar-pseudoscalar channel of the four-quark interaction. Accordingly, the two-quark-two-meson single quark exchange coupling simply measures the strength of the Yukawa coupling and has no physical interpretation. In particular, it does not entail a possible decoupling. This leaves us with with φ = σ, π and the dimensionless mass functionsm defined in (F2). Note also that these exchange couplings generate genuine four-field interactions in the flow, such as the four-quark and four-meson vertices. Consequently, the strength of these higher order couplings is tightly linked to that of the single field exchange couplings, the analysis of which suffices in the absense of additional resonant interaction channels besides the scalar-pseudoscalar one.
In Figure 18 we depicted the two fundamental single field exchange couplings defined in (119). The UV dominance of the gluonic exchange coupling for cutoff scales k 1 GeV is clearly visible. In turn, for cutoff scales k 1 GeV the mesonic exchange coupling takes over gradually with g 2 lAl =h 2 /(1+m 2 π ) at k ≈ 600 MeV, while the gluonic coupling decays more quickly. Note that the exchange couplings only provide the qualitative picture; the respective flow diagram also have different combinatorial factors from the different colour and flavour traces. For the four-quark flow a details comparison is done in Appendix L, and is summarised in Figure 30.
In the matter-dominated regime the best indicator for the further decoupling is given by the respective propagator gapping with φ = σ, π, q = l, s and them defined in (F2). These gapping functions are depicted in Figure 19, where also the full gluon dressing is displayed for the sake of comparison. Note that the latter does not only entail the gapping information for cutoff scales smaller than that of the peak postion of the dressing function, k k peak , but also the running of the gluon wave function for k k peak . Qualitatively, this carries the same information as that depicted in Figure 18: gluonic fluctuations decouple first. However, Figure 18 carries the full decoupling which also requires the coupling strengths of quark-gluon vertex and Yukawa vertex respectively.
In turn, after the decoupling of the gluonic sector the interactions between quark and mesons only depend on the Yukawa coupling, and the respective propagator gapping is sufficient for the decoupling information: the propagator gappings of quarks and mesons clearly entail that diagrams with quark propagators decouple first (first strange, than the light quark diagrams). The σ-field decouples at roughly the same scale as the light quarks and finally the contributions of pion diagrams tend towards zero below the pion mass scale.
In summary, the sequential decoupling is clearly seen. First the gluonic dynamics decouples from the matter sector, then the quark and σ exchange diagrams decouple and finally the pion decouples. For example, this already entails that in the vacuum chiral perturbation theory emerges naturally in this setup, for investigations of respective low energy parameters see e.g. [200,201]. Moreover, it allows us to also investigate the validity bounds of chiral perturbation theory. Note also that the ordering in the sequential decoupling is changes at finite density. In particular the σ mode decouples later in the vicinity of the CEP, where, as the critical mode, gets massless.
This leads us to an intriguing and simple picture: Initiating the QCD flow at a large perturbative cutoff scale k = Λ UV QCD , here Λ UV QCD = 20 GeV, the flow towards smaller scales leaves us with the dynamics of quarks and For comparison, we also show the gluon dressing function Z peak A /ZA with the solid blue line. Here, 1/ZA is normalised by its peak position at k peak , that is Z peak A = ZA(k peak ): for cutoff scales k k peak the k-dependence of ZA is dominated by the running of the gluon wave function renormalisation, for k k peak it shows the gluon gapping. In summary, the sequential decoupling of gluon, quark and meson dynamics towards the IR is evident here. mesons for k 1 GeV. In the present setup with dynamical hadronisation we flow into a Polyakov loop enhanced quark meson model. In turn, without dynamical hadronisation we arrive at a Polyakov loop enhanced NJL-type model. We also emphasise that while the offshell dynamics of gluons decouples, the gluon sector still leaves its imprint in terms of the quantum equations of motion for the gluonic background A 0 or L,L.
Hence, the underlying assumption of LEFTs, that high-energy degrees of freedom are integrated out, can be made explicit with the present setup. Moreover, this can also provide valuable reliability checks for LEFTs. In summary this provides an explanation and further reliability. The sequential decoupling demonstrated here adds further justification for the unreasonable effectiveness of LEFTs in QCD. For a selection of fRG work on QCD LEFTs see [30,34,202].
We emphasise that the emergence of QCD LEFTs in the infrared at scales k 1 GeV within the current approach allows for their systematic improvement in terms of QCD-assisted LEFTs, and constitutes important progress in the setup and the qualitative and quantitative reliability of low energy effective theories. Of course, QCD LEFTs in general rely on QCD information. For example, the parameters of QCD LEFTs used for QCD at finite temperature and density are typically fixed with vacuum observables such as the pion decay constant and meson masses as well as quark constituent masses. Further QCD input ranges from the enhancement of the models with a Polyakov loop potential that carries information about the confinement-deconfinement phase transition.
A step in the direction towards using dynamical or fluctuation information of QCD is done with the implementation of the QCD-running of external parameters such as the critical temperature parameter T 0 in the Polyakov loop potential as suggested in [203]. A first example for a QCD-assisted LEFTs is then provided by [88,89], where the phenomenological approach in [203] was confirmed and improved by results for the Polyakov loop potential from fRG-computations in Yang-Mills theory [69,129] and N f = 2 flavour QCD [26]. Another example is provided by [109], where the UV initial conditions for the QCD LEFT were computed from QCD flows.
The results and in particular the flows of the current work can be readily used as external input for QCD LEFTs, hence lifting them to QCD-assisted LEFTs. An interesting and relevant example under progress is provided by fluctuation observables such as cumulants of net-particle multiplicity distributions, e.g. [30,75]. A qualitative improvement of the computations there is given by utilising the quark-gluonic flows from the present work as external input. This provides quantitative reliability of the results at large density or small √ s relevant for the search of the critical end point. In Figure 20 we show our results on the phase structure in T -µ B plane for N f = 2 and N f = 2 + 1 flavours. In both cases, as discussed in Section V A, we define the pseudocritical temperature via the thermal susceptibility of the renormalised light chiral condensate, ∂∆ l,R /∂T . At vanishing baryon chemical potential, the respective chiral transition temperatures are The width of the transition is taken to be the full width at 80% of the maximum of this susceptibility. For N f = 2, the crossover phase boundary is given by the long-dashed red line and its width by the gray band. The corresponding CEP is denoted by the red star. For N f = 2 + 1, the short-dashed black line marks the phase boundary and the blue band its width. The black point shows our result for the CEP of N f = 2 + 1 QCD. At small chemical potential lattice results provide a benchmark test for the present fRG computations. This includes the curvature of the phase boundary. It is a sensitive measure for the correct relative strength of quantum, thermal and density fluctuations. The curvature coefficient κ of the phase boundary at vanishing baryon chemical potential is the quadratic expansion coefficient of T c (µ B ) around µ B = 0, i.e., with T c = T c (µ B = 0). We refer, e.g., to [7] for additional relevant discussions. As for the transition temperature, the curvature of the phase boundary of a crossover depends on its definition [70]. By using the renormalised light chiral condensate, we ensure comparability with the lattice results in [44,49]. To extract the curvature, we fit our numerical results for T c (µ B )/T c by the polynomial given in (122) for orders 2 and 4 in µ B /T c within the regions µ B /T ≤ 2 and 3, and read-off the corresponding quadratic coefficient κ. Our final result for the curvature with its numerical error is then given by the mean and the standard deviation of all these coefficients. The curvature for the two-flavour pseudocritical line is determined to be For the curvature of the 2+1 flavour phase boundary we find As has been emphasised before, our results are consistent with recent N f = 2 + 1 lattice results, e.g. κ = 0.0149 (21) in [44] and κ = 0.015(4) in [49]. Note also that This entails that the two dimensionless phase transition lines given by T c (µ)/T c (0) in (122) agree for small baryon chemical potentials. This even holds for a large µ B -range as can be already deduced from Figure 20.
reference curvature κ N f = 2 N f = 2 + 1 fRG: this work 0.0176(1) 0.0142 (2) lattice: [44] 0.0149 (21) lattice: [47] 0.0144 (26) lattice: [49] 0.015 (4) DSE: [37] 0.0238 DSE: [63] 0.038 lattice: [204] 0.0078 (39) lattice: [205] 0.0056(6) All results, including DSE, are summarised in Table II. For a comparison to other results and potential implications for experimental CEP searches, we also refer to the discussion in the introduction, Section I, and in particular Figure 1. At finite baryon chemical potential lattice simulations are hampered by the sign problem, resulting in increasing statistical and systematic errors with increasing µ B . The former are shown in Figure 1, the latter are unknown. The respective validity bounds are indicated by the dashed lines at µ B /T = 2 and µ B /T = 3 in Figure 20. For µ B /T close to the crossover line and larger than µ B /T = 3 also the approximations used in the present work start to lose reliability: In the presence of the nontrivial mesonic dispersions, the phase structure can be altered, cf. Section V B. Accordingly, results in this regime, including our result for the location of the CEP, have to be interpreted with care.
We find CEPs for N f = 2 and N f = 2 + 1 flavour QCD at and respectively. This amounts to rather large values of µ B CEP /T CEP of Our results for the CEP are extracted from the onset of critical scaling on the chiral phase boundary. As discussed in Section IV D and Appendix F, we expand the effective potential about its running minimum. The flow of the minimum is proportional to the correlation length 1/∂ 2 σVk EoM (cf. (A29b)), and hence becomes unstable in the critical region. We use this instability to determine where the system enters the critical region from below in µ B . To pin down the location of the CEP, we exploit the well established fact that, owing to the very small pion mass and strong mesonic fluctuations close to the CEP, the critical region is at most a few MeV wide in µ B ; see e.g. [79]. We emphasise that the critical properties are solely driven by the quark-meson fluctuations and hence QCD LEFT investigations are fully applicable here. Accordingly, to very high accuracy, the CEP is located at the point on the phase boundary where the running minimum becomes unstable.
We emphasise that the reduced reliability of our results at large µ B is not of conceptual origin, as opposed to the sign problem on the lattice. It is related to correlation functions that have been neglected in our truncation of the effective action. Improvements require the systematic inclusion of additional momentum-dependences of correlation functions. Moreover, a Fierz complete basis of the four-quark tensor structures should be taken into account, including potentially dynamical hadronisation of additional resonant channels. This is work in progress within the fQCD collaboration [207].
A physically intriguing source for part of the systematic error in the present work is related to the indications of an inhomogeneous regime we discussed in Section V B. To highlight the potential effect on the phase boundary, we show the region where the meson wave function renormalisation at vanishing spatial momentum is negative, Z φ (0) < 0, with the shaded blue area in Figure 21. Note that we only scanned the phase diagram for µ B < µ B CEP . Hence, the boundary of the shaded blue region at µ B = 635 MeV is of no physical significance. As discussed in Section V B, the possibility of inhomogeneous condensation leads to a competition between potential resonances in the homogeneous and inhomogeneous quark-antiquark interaction channels. Regarding the phase structure itself, this is most relevant when condensation occurs in the first place. The region where the mesonic dispersion with a minimum at nonvanishing spatial momentum overlaps with a sizable condensate, defined in Section V B and the caption of Figure 21, is given by the hatched red area in this figure. There, the phase structure could be altered significantly by possible inhomogeneities. To put his into perspective, we also added a selection of freeze-out points to Figure 21. It is very intriguing that the inhomogeneous regime could be probed in heavy-ion collisions at small beam energies of roughly √ s 6 GeV.

VI. CRITICAL END POINT AND SYSTEMATIC ERROR ESTIMATES IN FUNCTIONAL APPROACHES
Here we combine the up-to-date CEP results of functional approaches for a first estimate of the most likely region of the CEP. To that end we discuss the respective approximations and the potential effect of missing fluctuations. A full discussion of the respective systematic errors is deferred to future work as it goes far beyond the scope of the present work.
The results for the phase structure in the present fRG approach are summarised in Figure 1 and Figure 20. The respective DSE results and an exhaustive discussion can be found in the recent review [7], for original work see e.g. [35-37, 39, 63, 98, 132-135].
It goes without saying that the approximations applied in the different works differ vastly. However, all DSE works share a common feature with the present fRG work: The most general quark-gluon interactions are not fully taken into account and both, momentum dependencies as well as tensor structures are missing. Typically, this causes a lack of infrared strength of dynamical chiral symmetry breaking [27,32]. This is compensated for by an infrared enhancement similar in spirit to that described in Appendix E 2. Accordingly, the systematic error inherent to these approximations (apart from other error sources) relates to the strength of this infrared enhancement.
As also discussed in Appendix E 2, such a setup guarantees the value of the chiral condensate and the correct temperature dependence of the order parameters at vanishing density or chemical potential. However, it potentially lacks quantitative precision for the µ B -dependence, as the µ B -dependence of the infrared enhancement, that is the µ B -dependence of the missing tensor structures in the quark-gluon vertex, is not known. On the other hand, while the explicit µ B -dependence is only present in the quark-gluon correlations, the chemical potential dependence of the gluon propagator has been shown to be very small, see Figure 14. Accordingly, the respective systematic error of dropping this dependence that has been used in various works is presumably very small. Now we adopt these observations to the phase structure computations with the DSE. All of them exhibit a larger curvature κ DSE in comparison to the lattice results and the fRG results of the present work, which are compatible. The larger curvature at low density implies a stronger µ B -dependence of, e.g., order parameters, which suggests -in a linear error propagation-also a CEP at too small chemical potential and temperature. Indeed, the DSE curvatures and positions of the CEPs can be ordered accordingly, see Figure 1.
Let us now concentrate on [37]. In this work, as in the present one, the phase boundary is well described by only the leading contribution proportional to κ in (122). Assuming that all physical effects influence the µ B -dependence only linearly, we can rewrite the curva-ture term of the DSE phase boundary as follows: and we note that T c at vanishing chemical potential is practically the same in both cases. Using the value κ DSE = 0.0238 and µ B,DSE CEP = 488 MeV from [7,37], we are led to to be compared with µ B CEP = 635 MeV in the present work, see (126). This is promisingly close. However, we rush to add that such a 'linear' error analysis can be deceptive. For example, the analysis of the µ B -dependence of the baryonic contributions in [7,135] show that a smaller κ induced from baryonic fluctuations comes with a smaller µ B CEP in contradistinction to the discussion here, involving the direct µ B -dependence in the quarkgluon vertex. This calls for a more elaborate combined error analysis in functional approaches, and we hope to report on this in the near future. Still, our preliminary analysis of the systematic error in state-of-the-art functional approaches allows for a simple estimate of a region where the critical end point, or more precisely the onset of a new phase of matter, is most likely. For this combined estimate we rely on the DSE results from [37], as this computation has the most complete µ B -dependence of the gluon propagator and the quark-gluon vertex to date. This is also reflected by the fact that it features the smallest curvature κ DSE of all the DSE computations, as well as the largest ratio µ B CEP /T CEP . Now we correct the DSE curvature by hand with the factor κ fRG /κ DSE . As the above systematic error analysis for the DSE also suggests that the CEP is at too low chemical potential, the rescaled CEP from the DSE gives us a lower estimate on the CEP on the T c (µ B ) curve of fRG and lattice with µ B ≈ 480 MeV. Moreover, at µ B 500 MeV the potential inhomogeneous regime intersects the chiral transition region, see Figure 21, and the competition between homogeneous and inhomogeneous condensation could modify the phase structure there. Below these values we see no indications for a critical end point. In combination, this leads us to the estimate where it is indicated with ∼ that the upper and lowers bounds are not strict, as they carry additional systematic errors of our approximation that are difficult to estimate. Below the lower bound we see neither an indication for a CEP nor signals for the failure of the current approximation. This is also sustained by the investigations in [74,76,142], where the Fierz-complete basis of four-quark scattering vertex has been taken into account: For more details see Section V. The phase boundary globally agrees well with recent lattice results, and in particular the curvature of the phase boundary for small chemical potential, κ = 0.0142 (2), is consistent with recent lattice results. We find a critical end point at (TCEP, µB CEP ) = (107, 635) MeV. fRG: inhom: The blue area depicts the regime with a negative slope of the mesonic dispersion at vanishing spatial momentum, that is Z φ (0) < 0, see Figure 12. There, the the meson dispersion has a minimum at a nonvanishing spatial momentum. This is a strong indication for an inhomogeneous regime. We also show the minimum of Z φ (0) in T at a fixed value of µB with the blue dashed line. The red hatched area shows where the inhomogeneous regime overlaps with a sizable homogeneous chiral condensate. The latter is defined as the size of the (reduced) chiral condensate where the thermal susceptibility ∂∆ l,R /∂T at fixed µB reaches 80% of its peak height (from above), and larger. In this region, a competition between homogeneous and inhomogeneous condensation is expected and the phase structure could be altered significantly. Except of a diquark channel at very large chemical potential, µ B /T 7, only the scalar-pseudoscalar channel gives large contributions.
In turn, for chemical potentials above the regime in (130) either a CEP has already been observed or the standard chiral symmetry breaking pattern has broken down. For the upper bound on the CEP region we therefore assume that the system does not enter the critical region in the vicinity of the CEP we find here. A simple estimate of the phase boundary without a CEP at (T CEP , µ B CEP ) = (107, 635) MeV is then given by an ex-trapolation of the boundaries of the width of the chiral transition in Figure 20 from the crossover region towards larger µ B . The intersection of these boundaries then yields the r.h.s. of (130). To be more specific, we fit the upper and lower bound of the width of the (2+1)-flavour chiral phase boundary in Figure 20 by the polynomial in (122) up to fourth order for chemical potentials µ B ≤ 630 MeV. The two curves intersect at (103, 660) MeV.
We emphasise again, that these considerations are only the first step towards a combined functional analysis of the CEP. Moreover, (130) should be rather seen as the regime with either a CEP or interesting new phenomena such as a phase with an inhomogeneous condensate. As already mentioned in the beginning of this Section, in particular the present estimate is based on a linear error analysis in a nonlinear problem and hence has to be taken with a grain of salt. Despite the obvious deficiencies of the present analysis this combined functional approach at large chemical potential, while also utilising benchmark results from lattice simulations at smaller chemical potential will finally enable us to provide quantitative predictions for the large density regime including the potential CEP in the next few years.

VII. SUMMARY AND OUTLOOK
In this work we have studied the QCD phase structure at finite temperature and baryon chemical potential as it emerges from the fundamental dynamics of QCD. To this end, we extended the dynamical hadronisation technique within the functional renormalisation group approach to QCD towards the inclusion of in-medium effects. Our results therefore constitute a new state-of-the-art regarding the understanding of the QCD phase diagram based on the dynamics of quarks and gluons.
We have investigated various quantities, including their temperature and chemical potential dependence. In particular, we have computed order parameters for chiral symmetry breaking, e.g. the renormalised chiral condensate, see Figure 10, as well as the Polyakov loop expectation value, see Figure 7. In addition, we studied the constituent quark mass, the pion and σ-meson masses, the wave function renormalisations/propagators, the running strong couplings and the interactions of quarks and the mesonic low energy degrees of freedom. This allowed us to gain insights into the interplay of different degrees of freedom at different momentum-, temperature-and density-scales.
Most importantly, we have obtained the phase diagram of N f = 2 and N f = 2 + 1 flavour QCD at finite temperature and baryon chemical potential, see Section V D. Various aspects of our results, including comparisons to other recent works, are shown in Figures 1, 20 and 21. Our result for the curvature of the phase boundary for 2+1 quark flavours at small chemical potential is κ N f =2+1 = 0.0142 (2), which is fully compatible with the most recent lattice results [44,49]. In addition, we find a critical end point at (T CEP , µ B CEP ) = (107, 635) MeV. We compare our results for N f = 2 + 1 and N f = 2 in Figure 20. Owing to strange quark fluctuations, the (2+1)flavour phase boundary is systematically lower than the two-flavour boundary.
An intriguing aspect of our results is that we find strong indications for an inhomogeneous regime at large chemical potentials in the region µ B 420 MeV, shown in Figures 1 and 21. Since this regime intersects the phase boundary, it could have important implications for the QCD phase structure. For a detailed discussion see Section V B.
A preliminary analysis of the systematic error in our computation is done in Section VI. In combination with a similar analysis of the DSE as well as with the benchmark results from the lattice at small chemical potential, we put forward a first suggestion for the region of the critical end point, see (130). This combination of results from functional approaches at large chemical potential and lattice results at small chemical potential will finally allow us to pin down the location of CEP or, given potential inhomogeneous phases and their implications, the onset of a new state of matter at large density.
This calls for a systematic improvement of the approximations used here. Firstly, only quark-antiquark scattering in the σ-π channel has been taken into account for the dynamical hadronisation. Other channels, however, may play a sizable rôle in particular at large µ B , such as diquark and vector-meson channels, see, e.g., [74,76]. Therefore, the present approximation for the four-quark couplings is currently extended to a Fierz-complete basis, see [27,32,74,76]. In addition, we need to allow for inhomogeneous quark-antiquark scattering channels in order to supplement our present results in this direction. Secondly, we have employed a minimal extension to include strange quark dynamics in our calculation. In the current work we have been concentrating on specific light quark observables which are insensitive to the details of the strange sector. However, this is insufficient for the highly interesting investigations of strangeness at finite chemical potential and temperature; for an fRG point of view see [34,75]. For these aspects a full (2+1)flavour study, including the corresponding four-quark interactions and resonances, is required. Finally, we relied on external input for the momentum-dependent gluon and ghost propagators in vacuum, the in-medium screening of gluons in Yang-Mills theory and the gauge part of the gluon effective potential. Only the average momentum dependences encoded in the cutoff dependence have been considered. Full momentum dependences have so far been considered in vacuum QCD and finite temperature Yang-Mills theory, [27,28,[31][32][33]. Its extension to QCD at finite temperature and chemical potential is under way. This also includes a complete basis of tensor structures for in particular the quark-gluon vertex as well as the purely gluonic vertices. The selfconsistent determination of the gluon effective potential along the lines of [69] is also on our agenda.
Further interesting applications of the present work are manifold. In the following, we list selected active projects within the fQCD collaboration [207]. A detailed study of QCD thermodynamics and the equation of state is the most obvious application. This will provide input for phenomenological studies of heavy-ion collisions that is firmly rooted in QCD even at small beam energies. Since the chiral limit is easily accessible within the present approach, the magnetic equation of state is also being computed. Furthermore, fluctuations and correlations of conserved charges are evaluated. The present approach allows for both, the spectral reconstruction and the direct computation of real-time correlation functions in QCD along the lines of [104,198,[208][209][210][211][212][213][214][215][216][217][218][219], that are important for QCD-assisted transport, [220], and QCDassisted hydrodynamics, [221]. These applications play important rôles in experimental searches for the CEP, and more generally for our understanding of the QCD phase structure. We hope to report on these matters in the near future.

ACKNOWLEDGMENTS
We thank R. Alkofer, J. Braun, C. Fischer, J. Papavassiliou, R. Pisarski, J. Rodriguez-Quintero, B. J. Schaefer, S. Valgushev, L. von Smekal, N. Wink and S. Zafeiropoulos for discussions. We thank the authors of [144,145] for providing us with unpublished lattice data of the N f = 2 + 1 gluon propagator and a collaboration on related subjects, [146]. This work is done within the fQCD collaboration [207], and is supported by EMMI, the BMBF grant 05P18VHFCA, by the National Natural Science Foundation of China under Contract No. 11775041, and is part of and supported by the DFG Collaborative Research Centre SFB 1225 (ISOQUANT) as well as by the DFG under Germany's Excellence Strategy EXC -2181/1 -390900948 (the Heidelberg Excellence Cluster STRUCTURES). F.R. is supported by the DFG through grant RE 4174/1-1.

Appendix A: Chiral condensates
Here we explain of how to extract chiral condensates within the present approach. Results in the present approximations for the renormalised chiral condensate, (A3), and the reduced chiral condensate, (A4), are presented in Figure 22 and Figure 23 respectively. Naturally, these observables are derived from the thermodynamic potential density or free energy density defined in (67), evaluated on the EoM (6) with J = 0. This is a finite, renormalisation group invariant observable, as follows from its relation to the effective action Γ[Φ EoM ]. The chiral condensates are proportional to x q i q i with q i = u, d, s and vanishes in the chirally symmetric regime. Up to the normalisation it is given by the derivative of Γ w.r.t. the current quark mass m 0 qi , to wit where no sum over i is implied in (A1). The normalisation in (A1) includes the trivial volume factor T /V with the spatial volume V. The multiplication factor m 0 qi leading to the logarithmic current quark masses derivative removes possible ambiguities relative to the condensates from other formulations such as DSEs and the lattice. In the present work the masses of the light quarks u, d are the same, m u = m d = m l , and we define the light quark condensate ∆ l = ∆ u = ∆ d . Up to renormalisation terms the condensates ∆ qi read The necessity of the renormalisation of (A2) is typically circumvented by subtracting two condensates from each other. One convenient choice is to only consider the thermal and density part of the condensates, similar to the definition of the pressure . This leads us to the renormalised condensate ∆ qi,R with Indeed, ∆ qi,R = −m 0 qi ∂ qi p. In (A3) the normalisation N R renders ∆ qi,R dimensionless. N R is a convenient normalisation that is typically chosen to be one of the characteristic, well-determined, scales in the theory. A common choice is N R = m 4 π for the physical case. If also being interested in the chiral limit, N R = f 4 π is better suited.
Another common choice is the reduced condensate ∆ l,s , which is proportional to the weighted difference between light and strange quark condensate, ∆ l and ∆ s . It is normalised by its value in the vacuum and reads In the functional RG approach the necessity of renormalising (A2) is resolved by the very definition (A1) to use current quark mass derivatives of the the basic object in the approach, the finite effective action Γ[Φ EoM ]: We first use that ∆ l can be represented as a c σ -derivative, where we have used that c σ ∝ m 0 l carries the only dependence of Ω on the current quark mass m 0 l . This leads us to the finite chiral condensate which is simply the term in the effective action that explicitly breaks chiral symmetry. We emphasise that (A6) is an exact relation. Using (A6) in (A3) we arrive at This is readily computed within the present approach. In Figure 22 we show the results within the current approximation in comparison to the continuum extrapolated lattice result in [183], for results at finite baryon chemical potentials see Figure 10 in Section V A. The normalisation constant in (A3) and (A7) is chosen to match the scale in the lattice calculation.
In an approach with (dynamical) hadronisation also for the strange quark sector also the reduced condensate can be represented in terms of the meson field expectation values. Similarly to (A5) we get leading to With (A6) and (A9) we arrive at where we have used that see also (47) and Section III C for its derivation. It is implicitly understood that σ and σ s in (A10) are evaluated on the respective EoMs.
In the current work we utilise the simple approximation (49) for the respective strange condensate. Using (49) in the reduced condensate (A10) leads us to the simple expression Notably, the explicit dependence on the current strange quark mass m 0 s has dropped out. Equation (A12) also makes evident that it is even nonvanishing in the flavoursymmetric case with m 0 l = m 0 s even though both the numerator and the denominator in (A10) vanish.
Equation (A10) still contains explicitly the parameters of explicit chiral symmetry breaking in the initial action, c σ , c σs . Now we rewrite (A10) in terms of the full effective action Γ k=0 . To that end we utilise the equations of motion. For constant solutions σ EoM , σ s,EoM they read where V (1,0) k (ρ, ρ s ) = ∂ ρ V k (ρ, ρ s ) and V (0,1) k (ρ, ρ s ) = ∂ ρs V k (ρ, ρ s ). Inserting (A13) for k = 0 in (A10) yields where the fields ρ, ρ s are evaluated on the respective EoMs (A13) at finite T, µ q and in the vacuum. Equation (A14) depends on the unrenormalised pion mass m 2 π = V (1,0) (ρ EoM ), and the effective potentials in (A14) are evaluated at vanishing cutoff scale, V (ρ, ρ s ) = V k=0 (ρ, ρ s ). Equation (A14) only depends on the full effective potential at vanishing cutoff scale k = 0: the apparent dependence on the parameters of the initial effective action has been removed. In the chiral limit both the numerator and the denominator of the reduced condensate vanish as they should. A convenient reparameterisation of (A14) leads us to ∆ l,s = σ(T, µ q ) σ(0, 0) Its computation only requires the knowledge of the full effective potential V (ρ, ρ s ) at vanishing cutoff scales at the T, µ q -dependent EoMs σ EoM , σ s,EoM : the T, µ qdependent pion and masses. The flow equation for the effective potential and the relevance analysis in Section V C and Appendices K, L allows us also to provide a simple expression for V (ρ, ρ s ). The effective potential is generated from the σ, σ s dependences of the Yukawa terms. These terms have the full N f -symmetry which is unspoiled by the flow as neither c σ nor c σs are present. In the present approximation we only consider the diagonal terms in the full potential. This is consistent as long as the potential is mostly driven by the quark contribution: This approximation is well justified for larger temperatures and small baryon chemical potential, but may fail for larger chemical potential µ B and small temperatures. This is one of the reasons why we do not consider the reduced condensate as an order parameter for chiral symmetry breaking here, its reliability at large µ B requires a better approximation. In the present quark-dominated approximation the potential is simply a sum of the same potential V k for both ρ and ρ s separately. This yields which is a quantitative approximations for cutoff scales k 200 − 300 MeV. With (A16) we arrive at Equation (A17) can be evaluated within our current approximation, see Figure 23. Note that the reduced condensate, ∆ l,s , in contradistinction to the renormalised condensate, ∆ l,R , is sensitive to the correct relative treatment of strange and light quark sector. In the current work we have used a simple approximation for the strange quark sector as none of our observables is sensitive to the subleading differences between light quark and strange quark sector. Put differently, we may use the reduced condensate for fixing the strange current quark mass m 0 s or c σs , alternatively to the ratio of the decay constants f K /f π . Given the large value of the constituent strange quark mass in the Landau gauge, m 0 s 500 MeV, see e.g. [7], we expect m 0 s 150 MeV for this adjustment, for a detailed discussion see Section III C around (47). There we also referred to the N f = 2 + 1 flavour ratio from lattice simulations, f K /f π ≈ 27, which is obtained here for ∆m sl = 150 MeV , c σs = 97.5 GeV 3 . (A18) Note also that c σs and hence the ratio varies rapidly with ∆m sl from current quark mass ratios temperature-dependence of ∆ l,s at large temperature. This is clearly seen in Figure 23. Note that while this strange quark mass variation has no impact on the light quark and gluon correlation functions and observables, it of course is relevant for observables and correlations with strange quark content. In particular we see quantitative agreement of the reduced condensate for all temperatures studied here for current quark ratios close to the physical one.
We close this Appendix with a discussion of possible improvements of the computation of the condensates (or other observables) utilising the flow equation of the condensates instead of the their representation in terms of the effective action/potential at k = 0. Such an improvement is based on the fact that in a given approximation each flow-step generates terms (information) that is dropped when projecting on the approximation of the effective action at hand. Accordingly, a flow representation of a specific observables may keep this information.
Of course, without approximation both ways have to give the same results, which can be cast in form of an integrability condition, for more details see [56,118,222]. While being trivial for the full theory, (A19) is nontrivial within a given approximation. For the flow of the chiral condensate ∆ qi we use the flow equation for the effective action, Note that all the dynamical hadronisation terms in the flow vanish on the equation of motion within our choice as ∂ tφk EoM = 0. Seemingly, (A20) carries explicit and implicit dependences on the current quark mass. The implicit ones are that of the Φ EoM and that of the couplings. The latter dependences can be computed with the same recursive relations as has been put forward in [71] for the µ q -dependence. However, (A20) can be simplified with (A5), and we write for the light quark condensate where we have used that the diagrams depends on c σ only via the equations of motion Φ EoM . As in (A20) all the dynamical hadronisation terms in the flow vanish on the equation of motion within our choice as ∂ tφk EoM = 0. For the strange quark condensate we get In the current approximation this leads us to and For the explicit result we first concentrate on the light condensate. Within the approximation (A23) we get Equation (A26) leads us to Collecting all the ingredients this leads us to Note that the expression in the square bracket is nothing but ∂ σVk and m 2 σ is the second σ-derivative of V k , see (A27). Now we use that the flow of σ EoM,k follows straightforwardly from the cutoff dependence EoM which can be solved for the flow of σ EoM,k , Accordingly, (A28) is nothing but the flow of σ EoM,k up to a prefactor −c σ . This follows already from (A6), to wit which proves that (A19) holds nontrivially in the present approximation. Note also that at large cutoff scales the term in the second line of (A28) from the mesonic loop drops quickly proportional to σ EoM . The prefactor of the first term approaches the current quark mass cσ h k m 2 σ → 2m 0 l . Hence, in this limit (A28) reduces to the explicit m 0 l -derivatives in (A20). In turn, for smaller cutoff scales (A28) carries also the implicit ones in (A20).

Appendix B: Dynamical hadronisation
Here we carefully derive the hadronisation relations for the given hadronisation (15), which is recalled here, withê σ is a conveniently normalised vector in the σdirection,ê 2 σ = 1/ 2N f . Now we use the shift in the σdirection,Ċ k (Φ)ê σ , for absorbing the current quark mass of up-and down-quarks into the mesonic field. Here we allow a field dependence of the coefficientĊ k . Then, the scalar quark two-point function is proportional to σ and its symmetry part vanishes for all momenta p and cutoff scales k, This entails, that the flow of the symmetric part vanishes, Moreover, it also leads to vanishing flow diagrams, For the determination ofĊ k we evaluate the flow equation for the scalar part of the quark two point function, derived from (26) with the choice (B1). It reads (qτ 0 q)σ = Flow where the field-dependence of all correlation functions has been suppressed, and the subscript (qq) indicates the projection on the scalar part of the quark two-point functionqq = 2qτ 0 q. Upon symmetrising (B6) the flow of the quark two-point function, ∂ t Γ (qτ 0 q) and the diagrams, Flow (2) (qτ 0 q) , drop out due to (B4), (B5). This leaves us witḣ where we have used that (B3) leads us to Γ (qτ 0 q)σ [σ]. In the chiral limit the numerator in (B7) vanishes andĊ k ≡ 0 as expected. For finite current quark masses (B7) the symmetric sum Γ (1) This is the expected property following from the symmetry constraints on the σ-dependence we have imposed in (B3), see also the discussion in Section II A. With (B8) we are led to As (B9) depends on Γ (1) σ , we also study the respective flow in order to discuss the selfconsistence of the current expansion schemes as well as the optimisation of the approximation (39) used in the present work. The flow of Γ (1) σ at (q,q) = 0 reads and is that in the chiral limit: it does not depend on c σ , see Section 52. Accordingly the dynamical hadronisation put forward here carries the full chiral symmetry at each cutoff scale due to the shift in J σ . Explicit chiral symmetry breaking is entirely carried by the linear term in σ. Inserting (39) in (B9), we arrive at where m 2 π (ρ) = ∂ ρ V (ρ). In terms of dimensionless variables and after rescaling with the wave function renormalisations we are led to (115).
In summary we are left with the two key equations (B10) and (B11). In particular the latter one, (B11), relates the ρ-dependence of the pion mass function m 2 π (ρ) to that of the Yukawa coupling. This entails that the field-dependence of the Yukawa coupling is in one-to-one correspondence to that of the effective potential, see also [70]. Since we do not take into account the ρ-dependence of h k (ρ), the respective systematic error is minimised for an expansion point in the Taylor expansion in ρ, that minimises the higher Taylor coefficients of m 2 π (ρ): The Taylor coefficients in (B12) are minimised at the flowing minimum, the k-dependent EoM for κ k = σ 2 EoM,k /2. Note however that this expansion point does not lead to best convergence for the flow of the effective potential, (B10). The latter flow shows best convergence for a fixed expansion point close to the final minimum σ EoM = σ Eom,k=0 . Accordingly, for small cutoff scales a fixed expansion point fares better than the running one in terms of convergence. It is suggestive to switch of the running of the renormalised expansion point at see (117) in Section IV D. However, if the frozen expansion point is too far away from the physical configuration ρ EoM =ρ EoM,k=0 ≤ρ EoM,k peak , the latter may not be in the radius of convergence of the Taylor expansion. This is checked by stopping the flow of the minimum later and checking the stability of the results under this procedure. This subtlety can be avoided completely by a fielddependent Yukawa-coupling (and further field-dependent couplings such as the wave function renormalisations), see e.g. [70]. This will be studied in future work.
Appendix C: Strange quark and the minimal extension to N f = 2 + 1 In this work we adopt an minimal approach to extend the calculation of N f = 2 to that of N f = 2 + 1. First of all, the contribution of the strange quark loop to the gluon anomalous dimension is included, i.e., another term denoted by η s A,k is added in (100), which reads The in-medium gluon vacuum polarisation generated by the light quarks, ∆η l A , is given by the in-medium part of (H1) and (H2) withm q =m l , N f = 2 and the lightquark-gluon coupling, gl Al , defined in (106) and (M1). For the strange quark contribution, η s A , we use (H1) directly withm q =m s , defined in (45), N f = 1 and the the strange-quark-gluon coupling, gs As , defined below. Note that η s A not only contains the thermal contribution, as does ∆η q A , but also includes the vacuum contribution, since our input for the unquenched QCD calculation in vacuum, i.e. the first term on the r.h.s. of (C1), in only for two flavours [32].
Furthermore, the quark-gluon coupling for the strange quark is distinguished from that for the u-and d-quarks, and its flow equation is given by where the last term is defined in (M2). In comparison to (106), we note that in (C2) we have neglected the contribution of mesonic fluctuations to the running of gs As . This will be improved in our future work. All explicit expressions for the quantities discussed here are given in Appendix H and M.

Appendix D: Direct sum versus full back coupling
Here, we motivate the approximation we use for the gauge field propagator. The goal is to find a simple, yet quantitatively reliable scheme. This discussion is largely based on [28,131] and reiterated here for completeness. The most prominent feature in the gauge sector is confinement. In linear covariant gauges it manifests itself in the emergence of a gluon mass gap, cf. (32), together with an enhanced ghost propagator in the IR [223][224][225]. This implies that in order to capture the mass gap, the full momentum dependence of the gauge field propagators has to be resolved. This is a formidable task which requires sophisticated truncation schemes and considerable numerical effort, e.g. [27,31,33,[226][227][228][229]. Here, we motivate an approximation scheme which uses the full gauge field propagators from fist-principles computations of the pure gauge theory as input, while matter effects, i.e. quark and hadron fluctuations, and their backreaction on the gauge sector, are computed selfconsistently.
Within the truncation scheme described in Section III, the propagators enter the flow equations through the corresponding anomalous dimensions and mass parameters. The gluon anomalous dimension, for instance, is in general a function of the strong couplings and the mass gap, where we omitted the k-dependence on the couplings and masses for the sake of brevity. Within both, the fRG and DSE approaches this can always be split into a part where only diagrams of the pure gauge theory contribute, η glue (αc Ac , α A 3 , α A 4 ,m 2 A ), and the matter contribution, η quark A (αq Aq ,m 2 q ). We note that the couplings themselves are the ones of QCD, i.e. the couplings in the gauge sector also depend on the ones in the matter sector and vice versa. Hence, the splitting of the gluon anomalous dimension into gauge and matter part is purely formal. A first simplification can be made by using that the different gauge couplings only differ from each other in the low-energy regime where the mass gap is generated, i.e. where gluons decouple in QCD. Hence, we parametrise the anomalous dimension as where α s stands for α ccA , α A 3 or α A 4 . Note that the gluon mass gap is part of the gauge contribution. As discussed in Section III, the glue diagrams of the anomalous dimension in Yang-Mills theory and in QCD agree, and it can be parameterised in terms of the couplings α s and the mass gap m 2 A . Accordingly, the differences between the glue parts in QCD and Yang-Mills theory come from the differences in α s and m 2 A . We have where α s and m 2 A are the QCD couplings and mass gap respectively. This suggest a possible approximation scheme for the gluon propagator: One takes the Yang-Mills anomalous dimension η YM A in (D3) as an external input from a first-principles computation, computes α s , m 2 A and η quark A (αq Aq ,m 2 q ) as functions of the RG scale k and puts everything together according to (D2).
The advantage of this procedure is that the information about the full momentum dependence, and in particular the emergent mass gap, is encoded in the quantitatively reliable external input. The corrections from matter fluctuations as well as their feedback on the gauge contribution are computed 'locally', i.e. the full RG scale dependence of the correlation functions is resolved, but only one momentum configuration (typically vanishing momenta) is taken into account. A simple reduction of this procedure is to ignore the feedback from the matter sector onto the gauge sector and simply define This amounts to only adding the gluon vacuum polarisation to the input from YM.
To check the quality of both approximations we compare the resulting gluon propagators to results from lattice gauge theory [230,231] at vanishing temperature and chemical potential in Figure 24. In this figure we also show the benchmark fRG calculations with the state-ofthe-art truncation of pure Yang-Mills theory in vacuum, [31], (thin green dashed line), and N f = 2 flavour QCD [32] (thick cyan dashed line). The fRG results in the current approximation have been taken from [28,131]: The The open red circles are the corresponding lattice results [230].
The blue open squares are the lattice results for two-flavour QCD [231]. The solid blue thick line shows the fRG result using (D2) and the dotted black line shows the result based on (D4). Both results are taken from [28,131]. Furthermore, results from the fRG calculations with the state-of-the-art truncations at vacuum for the pure Yang-Mills theory [31] and N f = 2 unquenched QCD [32] are shown for comparison, which are denoted by the green and cyan dashed lines, respectively.
red line shows the input propagator from the pure gauge theory. The open circles are the corresponding lattice data in [230]. The open blue squares show the lattice results for two-flavour QCD [231]. The solid blue line shows our results within the approximation based on (D2). We used α s = αc Ac . Indeed, this approximation leads to a perfect agreement with the corresponding lattice results. To facilitate a comparison between p-and k-dependent quantities, we have identified k = p. Note that the momentum scale is arbitrary here. If the feedback to the gauge sector is ignored, i.e. the approximation in (D4) is used, we find the dotted black line in Figure 24. While it still is a good approximation to the full propagator, it is clearly outperformed by (D2).
This results motivates the use of a local approximation, i.e. resolving only the RG scale dependence of the correlation functions, for the effects of matter fluctuations, as long as the relevant nonlocal information, i.e. the genuine momentum dependence, is taken into account through the input. This is crucial for the present work, since we compute the thermal and density corrections to the gluon vacuum polarisation in a local approximation directly, while we use the gluon anomalous dimension of QCD at T = 0 as well as the thermal effects of YM as input, cf. (100).
such that all avatars of the strong coupling agree at the initial scale. The initial strong coupling α s (Λ) is constrained by the requirement of the RG consistency in (112), and we find that (112) is satisfied for The only other input at our perturbative initial scale are the current quark masses. We identify the up and down quark masses, m u = m d = m q , and fix their value by the pion mass in the vacuum to be m π = 138 MeV. Following our discussion in Section IV A 1, we emphasise that the curvature massm π of the pion has been shown to be in very good agreement with the pion pole mass m π,pol defined by Γ (2) ππ (p 2 = −m 2 π,pol ) = 0, see [95,215]. This can be traced back to the mild momentum dependence of the pion correlation function, [27,32,95] as well as its relatively small width. For the σ resonance this is less obvious and in particular the width is not expected to be small. Accordingly the curvature mass m σ may not be close to the position of the smallest scalar resonance. In the present work it is computed to be m σ = 485 MeV. As discussed in Section III, the light quark mass is absorbed in a shift in the σ-field and simply leads to a linear term in the potential, c σ σ, see Section III C. We arrive at Note thatc σ,k = c σ /Z 1/2 φ,k with c σ being independent of k. This leads to The strange current quark mass is fixed as m s,Λ = 120 MeV by the ratio f K /f π , see the discussion in Section III C and Table I. This amounts to c σs = 48.8 GeV 3 , with ∆m sl = 120 MeV . (E6) We have checked that a variation of the difference of the constituent quark masses 100 MeV ∆m sl 200 MeV does not influence our results for light quark and gluon correlation functions and observables.
In the perturbative regime all other couplings are irrelevant, so they are either vanishing or have irrelevant RG flows. This entails in particular that the full fourquark couplingλ q,eff (that on the EoM for φ) should be infinitesimally small, where we have used the trivial initial effective potential Equation (E7) fixes the ratio of the meson mass parameter and the Yukawa coupling. Dynamical hadronisation entails that results should only depend on the combination (E7), and not onm 2 φ,Λ andh Λ separately. Moreover, ifλ q,eff is small in comparison to its flow, λ q,eff /∂ tλq,eff ≈ 0, (E7) holds, and the results at k = 0 do not depend onλ q,eff . We have checked that this is indeed the case as well as the independence onm 2 φ,Λ and h Λ separately, see also [28]. In fact, this behaviour is guaranteed by an IR-attractive fixed point in the weakcoupling regime [125].
The above scenario holds for initial cutoff scales of Λ 10 GeV, as we are safely in the perturbative regime of QCD. Then the above initial conditions should be applicable with appropriately adjusted values for the RGconsistent strong coupling. We have investigated the dependence of our results on the value of Λ in this regime. We found that the Λ-dependence is negligible.

Infrared enhancement
In Section IV B 1 we have discussed a potential infrared enhancement of the quark-gluon coupling in order to compensate for the dynamics of the missing nonclassical tensor structures in the present approximation. We follow the approach proposed in Ref. [28]: The potentially missing interaction strength is compensated for by a phenomenological infrared enhancement of the vector tensor structure. This is done with the replacement, with the infrared enhancement function given by This function behaves as ς a,b (k) → 1 with k > b and ς a,b (k) → 1 + a with k < b. In this work b = 2 GeV and δ = 2 are chosen. The strength of the enhancement for N f = 2 + 1, a = 0.034, and for N f = 2, a = 0.008, is fixed by fitting the physical constituent quark mass. Note that this IR enhancement is only applied for the quarkgluon vertex, and the three gluon coupling g A 3 is not enhanced. This takes into account that the infrared enhancement takes into account missing tensor structures, which is relevant for the quark-gluon vertex, but not for the three-gluon vertex.

Scales
QCD in the chiral limit, that is with vanishing current quark masses m 0 l,s = 0, has no scale on the classical level, and all observables can only be measured in the dynamically created scales: The confinement scale, carried e.g. by Λ QCD or the string tension, and the chiral symmetry breaking scale, carried e.g. in the chiral condensate or f π,χ .
At the physical point with m 0 l,s = 0 all scales are measured in the above dynamical scales, for more details and the possible choice of observables see Section III C. In particular, this leads us to Naturally, the scales at the physical point may also be measured in units of the reduced condensate, Λ QCD or the string tension. However, the pion decay constant has a direct physical meaning and hence is preferable. These scales may also be set in comparison to lattice simulations, which has been done in [27,32] for N f = 2 flavour QCD. In the present work we took over the scales from [27,32], fixing the scales at Λ = 20 GeV in a direct extension to N f = 2 + 1. RG-consistency, see Section IV B 3, then determines the value of the N f = 2 + 1 flavour strong coupling α s,Λ at the initial scale, providing us with a prediction of the N f = 2 + 1 gluon propagator, see Figure 25. This prediction is in quantitative agreement with the continuum-extrapolated lattice data from [144][145][146] with a pion mass of m π = 139 MeV. This completes our scale setting for the physical point. All quantities are measured in the respective absolute units.
For the sake of completeness we also provide results for the pion decay constant. This computation is also interesting in terms of the predictive power as well as the limitations of the present approximation. The pion decay constant can be determined from the momentum dependence of the quark propagator dressings Z q (p) and M q (p), defined in (77). These dressings are linked to the Bethe-Salpether wave function of the pion, for reviews see e.g. [232][233][234], Equation (E12) is sensitive to the full momentum dependence of the quark mass function M q (p). This dependence is suppressed in the flows due to the presence of the regulator term in the propagators. The letter suppression  [32] is the input in the present work. It is in quantitative agreement with the respective lattice results. The lattice data here are taken from [231]. The N f = 2 + 1 flavour gluon dressing shown here is a genuine prediction of the present computation. It is shown to be in quantitative agreement with the respective lattice results [144][145][146].
is the argument for the momentum-independent approximation in the present work with m l,k = M q,k (p = 0). Note also that the results in [32] have been obtained for a fixed expansion point or background with ∂ t κ = 0 with a four-dimensional flat regulator. Both properties can be taken account in the present work by evaluating the mass function on a fixed background κ = 1/2σ EoM at vanishing cutoff k = 0 as well as rescaling the cutoff scale k → 3/4k. Moreover, we have to rescale the mass function M q,k → M q,k m l,0 /M q,0 in order to have the same infrared value. This leads us to Figure 26, which first of all confirms impressively the validity of our approximation scheme. Moreover it allows us to use the rescaled quark mass function M q,res (p) of [32] for an estimate of f π in the present work with (E12). The rescaled quark mass function reads and the result for the pion decay constant at the physical point is see also Table I. In the current approximation the N f = 2+1 flavour mass function has the same cutoff behaviour as the N f = 2 flavour one. This is related to the fact, that the strange quark feeds into the light quark mass function only indirectly over the gluon propagator and running coupling αl Al in the quark-gluon vertex. Rescaling the mass functions with their value at k = 0 makes this very apparent, see Figure 27. We conclude that the agreement of the rescaled M q,k (0) with the N f = 2 + 1 flavour light quark mass function m l,k works equally well as for the two-flavour case depicted in Figure 26. The respective plot is shown in Figure 28.
Consequently we can use the same estimate as in (E13) for the N f = 2 + 1 flavour case. This leads us to see also Table I. Finally we are interested in the pion decay constant in the chiral limit. The respective light quark mass function m l,χ at k = 0 is given by m l,χ = 319 MeV, leading to f π,χ ≈ 88.6 MeV (E16) In summary the scale setting in the present work is selfconsistent with both the confinement scales, here determined via the gluon dressing function, as well as the chiral symmetry breaking ones, here determined with the pion decay constant. This is further nontrivial evidence for the reliability of the present approximation for not too large chemical potential. are given in (N9) and (N10) in Appendix N, and the Polyakov loops L,L are not shown explicitly for brevity. The right hand side of (F1) also depends on the dimensionless renormalised quark and meson masses, i.e.,  [32] and the N f = 2 + 1 light quark mass function from present work in the fixed expansion.
Performing a Taylor expansion of the effective potential, one arrives at where the expansion point κ k can depend on the RG scale. Note that the expansion coefficients λ n,k should not be confused with the four-fermion coupling λ q in (39). It is more convenient to work with RG-invariant variables, which are denoted in this paper by quantities with a bar, except for the strong couplings in (54) and (55), for example,ρ =Z φ,k ρ,κ k =Z φ,k κ k , and λ n,k = λ n,k /Z n φ,k withZ φ defined in (92). The η Φi,k 's in (F1) are the anomalous dimensions for respective fields, which are given in (30). We reformulate the effective potential in terms of RG-invariant variables through V mat,k (ρ) =V mat,k (ρ), i.e., The expansion pointκ k can be chosen freely, provided that the physics of our concern is within the radius of convergence of the expansion [70]. Expanding about a fixed background is a convenient choice, see e.g. [28,70,71,73], where the unrenormalised expansion point is chosen to be independent of k. In this work we adopt the physical point expansion. In this expansionκ k is chosen to be the minimum of the effective potential at each scale k, i.e. it is given by the solution of the EoM, The flow equations of the expansion coefficients of the effective potentialλ n,k , which describe 2n-meson scattering, are ∂ tλn,k =∂ n ρ ∂ t ρV k (ρ) ρ=κ k + nη φ,kλn,k + ∂ tκk + η φ,kκk λ n+1,k .
In this work the maximal order of the Taylor expansion for the effective potential in (F4) is chosen to be N v = 5, and we have checked its convergence for the µ B -range studied here, 0 ≤ µ B 600 MeV. We employ the parameterisation put forward in [171] for the glue potential in (69). The advantage is that, besides the expectation value of the Polyakov loop and the pressure in Yang-Mills theory, quadratic fluctuations of the Polyakov loop are taken into account. For a detailed discussion we refer to [71]. The glue potential reads The temperature dependence of coefficients on the r.h.s. of (G1) enters through x(T ) = x 1 + x 2 /(t YM + 1) + x 3 /(t YM + 1) 2 1 + x 4 /(t YM + 1) + x 5 /(t YM + 1) 2 , The constants in G3 and G4 are taken from [171], and are also collected in Table III for convenience. t YM is related to the reduced temperature of Yang-Mills theory, and it has been found in [88,89,235] that unquenching effects in QCD are well captured through a linear rescaling of the reduced temperature of the pure gauge theory,   − F F (2,1) (m 2 q,k ,m 2 q,k ) + √ x cos 2 θ − cos θ with the fermionic regulator in (N1d). We have defined x = q 2 /k 2 , x = (q−p) 2 /k 2 . q and p are the loop and external momenta, respectively, and θ is the angle between them. In order to take into account the screening effect of the quark on the gluon propagator more efficiently, we choose the external momentum to be |p| = k in (H1). This is necessary for consistence, since we also evaluate the external input for the gluon anomalous dimension at |p| = k. The in-medium contribution, which we use for the light-quark contribution, is then The threshold functions in (H1) are given in (N20). Moreover, the second term on the r.h.s. of (100), viz. ∆η glue A , which is the thermal part of the pure glue contribution to the gluon anomalous dimension, has to be specified. Again, we build on the results in Yang-Mills theory in [33]. As discussed there, ∆η glue A is captured quantitatively through the inclusion of a thermal screening mass of the gluon, ∆m 2 scr (k, T ). In addition to the temperature, it has to depend on the RG scale since thermal effects are rapidly suppressed for k 2πT . Therefore, in the following we discuss the modification of the gluon anomalous dimension resulting from its thermal screening mass. Beginning with the inverse gluon propagator, one obtains a modified one with the screening mass, Differentiating both sides with respect to t leads tō whereη A = −∂ tZA /Z A is the anomalous dimension with the screening mass. Replacing Z A in (H5) by resorting to (H4), one is led tō Guided by the results of [33], we propose an ansatz for the screening mass as follows: Here c = 2 is adopted for N f = 2 + 1, which is consistent with the result in [33]. Furthermore, we choose n = 2 in (H7). We have also investigated the dependence of our results on the parameters in (H7), and find that variations of the parameters result in only a mild change of the pseudocritical temperature of the chiral phase transition.
In Section V C we show the momentum (or cutoff) dependence of the gluon dressing function 1/Z A (p) for different temperatures and baryon chemical potentials. In particular the dependence on µ B is very mild and only concerns the infrared. The in-medium behaviour is more pronounced in the gluon propagator depicted in Figure 29. However, we emphasise that the dressing functions are the relevant quantities for the flow equations in the present work.
where the ratioZ φ /Z φ (0) takes into account that the flow on the right hand side of (I1) has been defined with derivatives w.r.t. the renormalised fieldφ =Z φ φ. The term with the threshold function BB (2,2) arises from the mesonic loop of the flow equation for the meson propagator in Figure 5.
with x = q 2 /k 2 and x = (q − p) 2 /k 2 as same as in (H1). Note that in (I2) we have neglected the external momentum dependence of the meson loop for simplicity in the numerical calculations. This is a good approximation for both, small and large k, since meson loops are suppressed by the meson masses in the IR and decouple in the UV. Only in the phase transition region and in particular at the CEP, corrections to the mesonic part of (I2) could become relevant. treatment has been verified in [30] through a comparison with the fully frequency-dependent calculations.
In comparison to [71] we use the limit p 0,ex + i µ q → 0 instead of p 0,ex + i µ q → k for convenience. This modification allows us to compare our vacuum results more directly to that in [28], as there the flows are evaluated at p 0,ex = 0. We drop the suppression factor in the thermal part of the diagrams: p 0,ex → (πT ) − i µ q . This leads to simpler expressions, and the necessary suppression of thermal effects at small T /k is guaranteed by the thermal distribution functions. Both modifications have no sizable impact on the numerical results. simplify the flow in Eq. (L1) as well as numerical calculations, we approximate the expression above as q F · q F , which is exact for the 4d regulator. We have also investigated the reliability of this approximation through the comparison to the relevant four-quark flow in [28]. There a 4d regulator has been used, and we find that the difference is negligible.
In Figure 30 we compare (L1) and (L2) at T = 0 and µ B = 0 for the purpose of illustration. The flow divided by k 2 is shown in the left plot of Figure 30, and the ratio between the respective flow and their sum in the right plot. Note that the peak in the right plot corresponds to the position where the sum of (L1) and (L2) is vanishing, thus resulting in a divergent ratio there. This is because in the low k region meson exchange dominates, and this contribution changes sign at about k ≈ 75 MeV, as shown by the blue vertical line in both plots. One observes that gluon exchange, i.e., the diagrams in the first line in Figure 4, is the dominant process for k 300 MeV. Therefore, it is reasonable to expect that the mixed diagrams in the third line of Figure 4 play a subleading rôle for dynamical hadronisation. It is interesting to observe that gluon exchange becomes dominant already in (and below) the phase transition region at k ≈ 400 MeV, which can be read-off, e.g., from Figure 8. This gives an implicit upper bound on the validity of low-energy models, namely k 300 MeV. For first discussions in this direction, see [127,131]. Note that this bound is in quantitative agreement with the results in these references.
For the light and strange quarks we use the corresponding masses defined in (45) and set N f = 2 and N f = 1, respectively, above. Note that, as discussed in Appendix C, we only take (M2) into account for the strange-quarkgluon coupling.
We remark that the simple substitution n F (m 2 ; T, µ q ) → n F (m 2 ; T, µ q , L,L) in the flows is part of the current approximation. There are additional contributions from the nontrivial color trace in threshold functions that mix bosonic and fermionic contributions, see (N16) and (N19), which we neglect here. For L,L = 1 (N8) reduces to the Fermi-Dirac distribution of quarks, (N6). In turn, for L,L = 0, (N8) is the Fermi-Dirac distribution of baryons, that is (N6) with x → 3x. Note however, that this limit is not achieved for low temperatures as the Polyakov loop expectation value decays with L exp(2m q /T ) → ∞ for T → 0. Nonetheless, baryon number fluctuations show a baryonic regime for low temperatures as they should, for more details, see [71]. For the sake of brevity we drop the variables L and L in the quark distribution function.
In (N16) The threshold function L in the flow equation of the Yukawa coupling in (K1) is given by FB (1,2) (m 2 f , m 2 b ; T, µ q , p 0 ) Furthermore, we also need another class of threshold functions, appearing in the flow equation of the four-fermion coupling induced by the meson exchange in (L2). This leads us to define Again, all the functions FBB's in (L2) can be easily obtained from FBB (1,1,1) , whose analytic expression reads Finally, we would like to present the threshold functions in (H1), which is somewhat different from those mentioned above, since the external momentum p in the fol-lowing is nonvanishing, and the definition reads FF (n f 1 ,n f 2 ) (m 2 f1 , m 2 f2 ; T, µ q , p 0 , p) with G f given in (N4) and G f by and hereq = q/k andp = p/k. In the same way, one has and the explicit expression for FF (1,1) is given as follows FF (1,1) (m 2 f1 , m 2 f2 ) = k 3 4E q E q−p − 1 + n F (E q−p ; T, µ q ) + n F (E q ; T, −µ q ) with and the modified quark distribution function n F defined in (N8) with x = E q .