Hyper-order baryon number fluctuations at finite temperature and density

Fluctuations of conserved charges are sensitive to the QCD phase transition and a possible critical endpoint in the phase diagram at finite density. In this work, we compute the baryon number fluctuations up to tenth order at finite temperature and density. This is done in a QCD-assisted effective theory that accurately captures the quantum- and in-medium effects of QCD at low energies. A direct computation at finite density allows us to assess the applicability of expansions around vanishing density. By using different freeze-out scenarios in heavy-ion collisions, we translate these results into baryon number fluctuations as a function of collision energy. We show that a non-monotonic energy dependence of baryon number fluctuations can arise in the non-critical crossover region of the phase diagram. Our results compare well with recent experimental measurements of the kurtosis and the sixth-order cumulant of the net-proton distribution from the STAR collaboration. They indicate that the experimentally observed non-monotonic energy dependence of fourth-order net-proton fluctuations is highly non-trivial. It could be an experimental signature of an increasingly sharp chiral crossover and may indicate a QCD critical point. The physics implications and necessary upgrades of our analysis are discussed in detail.


I. INTRODUCTION
Some of the most challenging questions of heavy-ion physics are related to the transition from the early, non-equilibrium, state of quarks and gluons to the final hadronic states after chemical freeze out, which is observed in experiments. Unravelling this dynamics necessitates a thorough grasp on the physics in the QCD phase structure close to the confinement-deconfinement and chiral transitions. This regime is strongly correlated with highly non-trivial dynamics. Understanding this part of the phase structure, including the location and dynamics of a potential critical end point (CEP), plays a pivotal role in understanding phases of strongly interacting nuclear matter under extreme conditions. For works on the phase structure of QCD, covering experiment and theory see, e.g., [1][2][3][4][5][6][7][8][9][10][11][12][13], where theory covers first principles functional approaches and lattice simulations.
Fluctuations of conserved charges are very sensitive to the physics of the strongly correlated regime that governs the transition from the quark-gluon plasma (QGP) to the hadronic phase. They provide detailed information on the underlying dynamics. This includes, but is not limited to, possible experimental signatures of a CEP [4]. It has for example been proposed in [14][15][16] that non-monotonic variations of conserved charge fluctuations as functions of the beam energy can arise from critical physics in the vicinity of a CEP. During the last decade, significant fluctuation measurements have been performed in the first phase of the Beam Energy Scan (BES-I) program at the Relativistic Heavy Ion Collider (RHIC), involving various cumulants of net-proton, netcharge and net-kaon multiplicity distributions [17][18][19][20][21]. Remarkably, very recently the STAR collaboration has reported the first evidence of a non-monotonic variation in the kurtosis (multiplied by the variance) of the netproton number distribution as a function of the collision energy with 3.1 σ significance for central collisions [22]. The measurements have been extended to the sixth-order cumulants of net-proton and net-charge distributions, for preliminary results from STAR see [23,24].
Recent first-principle QCD calculations at finite temperature and density, within both the functional renormalisation group (fRG) and Dyson-Schwinger equations (DSE), show that the transition from the QGP to the hadronic phase is a crossover which becomes sharper with increasing baryon chemical potential, µ B , for µ B /T 4 [7,8,[25][26][27][28]. Beyond this region, a CEP might occur, but quantitative reliability of the theory computations cannot be guaranteed within the present approximations [7,8,25,28]. In addition, critical physics may only be observable in a very small region around the CEP, see e.g. [29]. Since the available RHIC data is limited to µ B /T 3, it is important to understand how conserved charge fluctuations are affected by the increasingly sharp crossover away from a regime with critical scaling.
To address these open questions related to the physics of strong correlations in the QCD phase diagram, we study in detail the T -and µ B -dependence of net-baryon number fluctuations in the range µ B /T 3. We present results for fluctuations up to tenth order (where we refer to everything above fourth order as hyper-order ), including comparisons to available results from RHIC [22][23][24] and predictions for the beam energy dependence of fluctuations where no experimental results are available yet.
To facilitate the comparison between theory and experiment, T and µ B can be mapped onto the beam energy per nucleon, √ s NN , via phenomenological freeze-out curves [30]. While these curves are expected to be close to the QCD crossover at large beam energies (corresponding to small µ B ) [31], they may move away from the transition region at low energies (i.e. larger µ B ) [32]. This can affect the beam energy dependence of particle number fluctuations, and requires a detailed understanding of the physics also outside the critical region. Aside from their phenomenological relevance, netbaryon number fluctuations at finite µ B can also be used to assess the reliability of extrapolations of thermodynamic quantities to finite µ B based on a Taylor expansion at µ B = 0. Such a strategy is commonly used in lattice QCD simulations, where a sign problem prevents direct simulations at finite µ B , see e.g. [9,[33][34][35][36][37][38]. By comparing the results of direct computations at finite µ B to the ones obtained from extrapolations at µ B = 0, we study the range of validity of a Taylor expansion at a given order self-consistently. Understanding the limitations of such an extrapolation is also relevant for phenomenologically constructed equations of state, as, e.g., in [39], where the non-critical physics at finite µ B crucially rely on this extrapolation.
All this is addressed within a QCD-assisted low-energy effective field theory (LEFT) which is described in detail in the next section. We use first-principles QCDresults on the T -dependence of the kurtosis and the µ Bdependence of the chiral phase boundary to map the inmedium scales of the LEFT onto QCD. This improves the reliability of our predictions, in particular at finite µ B . Non-perturbative quantum-, thermal-and density fluctuations are taken into account with the functional renormalisation group (fRG). This work therefore is a significant upgrade of previous work in [40][41][42], where net-baryon number fluctuations up to fourth order have been studied. The present QCD-assisted LEFT approach has various advantages. Most importantly, it is directly embedded in QCD as the relevant low-energy degrees of freedom emerge dynamically from systematically integrating-out the fast partonic modes of QCD [8,[43][44][45][46]. In addition, this approach allows us to capture both critical and non-critical effects in the QCD phase diagram. This entails in particular that our results agree with the results of lattice QCD at small µ B and show the correct universal behaviour of QCD in the vicinity of the CEP, i.e. 3d Ising universality.
Concerning the existence and location of the latter we add that the first-principles results in [7,8,[25][26][27][28] include a CEP in a region of 450 MeV µ B 650 MeV and therefore outside the regime of quantitative reliability of these computations. This suggests that the experimental detection of a CEP requires explorations of the high-µ B in the region with µ B /T c 4. Moreover, the direct experimental measurement of the CEP may be very challenging as it requires very high statistics, and predictions that signal critical dynamics can be further complicated by non-equilibrium effects. In the present work we shall therefore also outline how the location of a CEP could be constrained based on data in the crossover region, with- out the necessity of observing critical scaling. On the theoretical side this asks for first principles QCD studies for µ B /T 4. In turn, a first experimental step towards this goal is the solidification of experimental observation of the non-monotonic energy dependence of fourth-order net-proton fluctuations. Both is safely beyond the scope of the present work. This paper is organised as follows: In Section II we give a brief introduction to the fRG-approach to QCD and low-energy effective theories, including their mutual relationship. Thermodynamics and the hyper-order baryon number fluctuations are discussed in Section III. In Section IV, we first introduce a systematic scale-matching procedure between QCD and the low-energy effective theory. We then present our numerical results and compare them to lattice QCD simulations and experimental measurements. A summary with conclusions is given in Section V. Technical details regarding the flow equations are presented in the appendices.

II. QCD AND EMERGENT LOW-ENERGY EFFECTIVE THEORIES
At low momentum scales the quark-gluon dynamics of QCD successively decouple due to the QCD mass gap and spontaneous chiral symmetry breaking. This decoupling also applies to most dynamical (hadronic) low energy degrees of freedom at even lower energies, finally leaving us with dynamical pions and hence with chiral perturbation theory. Indeed, this successive decoupling is at the root of the success of chiral effective field theory.
The functional renormalisation group approach to QCD with its successive integrating-out of momentum modes is ideally suited to follow and study this decoupling. Diagrammatically, this is already seen within the flow equation for the QCD effective action, depicted in Figure 1. The different lines stand for the full nonperturbative propagators of gluons, ghosts, quarks and emergent low-energy degrees of freedom (hadrons in our case), where the loop momentum q is restricted by the infrared cutoff scale k, q 2 k 2 . In this setup, emergent bound states can be incorporated systematically by dynamical hadronisation [47][48][49][50]. For quantitative QCD applications in the vacuum see [43][44][45][46], for further conceptual developments and the application to the QCD phase structure important for the present work see [8].
The decoupling is apparent in this framework as the propagators carry the mass gaps m gap of gluons and quarks and for cutoff scales k m gap of a given field the respective loop tends towards zero.
More importantly, in this way the emergent low-energy effective theory is naturally embedded in QCD, and its ultraviolet parameters (at Λ 1 GeV) as well as further input may be directly computed from QCD, leading to QCD-assisted low-energy effective theories. We emphasise that this procedure does not lead to a unique LEFT. The dynamical degrees of freedom of QCDassisted LEFTs at Λ 1 GeV depend on the dynamical hadronisation procedure applied within QCD-flows. This setup and the QCD-embedding entail that, provided the relevant quantum, thermal and density fluctuations of low energy QCD are taken into account in the QCDassisted LEFT at hand, all QCD-assisted LEFTs encode the same physics, namely that of low energy QCD.
This leads to an equivalence relation of QCD-assisted Polyakov-loop-enhanced NJL-type LEFTs (PNJL), Polyakov-loop-enhanced QM LEFTs (PQM) and variations including higher meson multiplets and/or diquarks and baryons. We emphasise again that this equivalence relation only holds if low energy quantum, thermal and density fluctuations are taken into account. For more details see in particular [8], and the recent review [51]. Most prominently this embedding has been used for determining the temperature-dependence of the Polyakov loop potential, see [52,53]. This setup was then applied to the computation of fluctuations in [40-42, 54, 55].
In summary this entails, that for sufficiently small momenta k, temperatures T , and also density or quark chemical potential µ q , the gluon (and ghost) loop in Figure 1 decouple from the dynamics, and only provide a non-trivial glue background at finite temperature and chemical potential. The latter is taken into account with the Polyakov loop potential discussed in detail below.
In the present work, we build upon previous investigations of the skewness and kurtosis of baryon number distributions [40][41][42], and baryon-strangeness correlations [56,57] within QCD-assisted LEFTs with the fRG. The present LEFT is an upgrade of those used in the works above, and includes the quantum, thermal and density dynamics of quarks, pions and the sigma mode in a Polyakov loop background. It is a QCD-assisted PQM. As argued above, for low enough chemical potential, this model sufficiently close to QCD, and leads to results that are independent of the LEFT at hand.

A. 2-flavour setup
For the physics of fluctuations we are interested in scales below approximately 1 GeV. We restrict ourselves to k 700 MeV and temperatures and quark chemical potentials T, µ q 200 MeV. In this regime the only relevant quarks are the light quarks q = (u, d) and the strange quark s. The latter, while changing the momentum-scale running of the correlation functions, has subleading effects on the form of the fluctuations. Hence, the effect of the momentum-scale running induced by strange fluctuations will be mimicked here by an appropriate scale-matching detailed in Section II B.
We also include the lowest lying hadronic resonances, the pion π = (π ± , π 0 ), and, for symmetry reasons, the scalar resonance σ as effective low energy degrees of freedom. Within QCD flows these fields are emergent low energy degrees of freedom at cutoff scales k 1 GeV, that are taken care of with dynamical hadronisation in e.g. [8]. At the present low energy scales k ≤ 700 MeV, they are fully dynamical, and hence are part of the effective action at the initial cutoff scale. The other members of the lowest lying multiplet as well as further hadronic resonances produce rather subleading contributions to the offshell dynamics and hence are dropped. The mesonic fields are stored in an O(4) scalar field φ = (σ, π) with the corresponding chiral invariant ρ = φ 2 /2.
Quantum, thermal and density fluctuations with scales k Λ = 700 MeV are taken into account within the fRG, whose dynamics is now reduced to the last two loops in Figure 1. The respective effective action of QCD in the low energy regime is approximated by . We assume isospin symmetry and the corresponding chemical potential flavour-matrix is given by µ = diag(µ q , µ q ) = 1 3 diag(µ B , µ B ). Z q,k and Z φ,k are the wave function renormalisations for the light quarks and the meson respectively. Further running couplings considered are the Yukawa coupling h k , the scattering between quarks and mesons, as well as the effective potential V k (ρ, A 0 ), that describes the multi-scattering of mesons in the non-trivial glue background present at finite temperature and chemical potential.
The flow equation for the effective action Equation (1), and that for V k , h k , Z φ,q is described in Appendix A and Appendix B. The initial condition for Γ k at the initial cutoff scale k = 700 MeV is described in Appendix C.
The potential V k (ρ, A 0 ) has contributions V glue,k (A 0 ) from offshell glue fluctuations (first two diagrams in Figure 1), and contributions V mat,k (ρ, A 0 ) from the quark loop (third diagram in Figure 1). This leads us to The first contribution is typically reformulated in terms of the Polyakov loop L(A 0 ), while the latter is directly computed from the present low energy flow. This allows us to trade the A 0 -dependence for that of the traced Polyakov loop, L(A 0 ),L(A 0 ), see Appendix D, Equation (D2), leading us to to the final form of our potential, More details can be found in Appendix D.
In conclusion, the QCD-assisted LEFT described above and used in the present work, is a PQM-type model, e.g. [40-42, 52-59, 66, 67, 69-78]. Quantum, thermal and density fluctuations below Λ = 700 MeV are taken account with the functional renormalisation group, and the setup is well-embedded in functional QCD. As argued above, within the present, and analogous, elaborate approximation, the respective results (for fluctuation observables) for all QCD-assisted LEFTs match those of QCD for sufficiently low density. Therefore, we will refer to this model from now on as generic QCD-assisted LEFT.
The current QCD-assisted LEFT setup enables us to compute thermodynamic observables and in particular hyper-order baryon number fluctuations. However, as already briefly discussed in Section II A, we have dropped the dynamics of the strange quark. While we expect sub-dominant effects on hyper-fluctuations, the s-quark influences the momentum running of the correlations in the ultraviolet. Importantly, in [8] it has been observed on the basis of genuine N f = 2 and N f = 2 + 1 flavour computations in QCD, that the latter effect is well approximated by a respective universal scale-matching of the 2-flavour results even in QCD. Such a scale-matching has already led to a quantitative agreement of thermodynamics and kurtosis within the current LEFT setup with lattice results, see [40][41][42]. Thus, the present scale-matching entails that we use information on the T -and µ B -dependence of welldetermined quantities in QCD. This leads to an improved reliability of our results of finite T and µ B , as in-medium effects in the QCD-assisted LEFT are directly connected to in-medium effects in QCD.

2-to 2+1-flavour scale-matching in QCD
Given its relevance for the predictive power of the present LEFT within a QCD scale-matching procedure we briefly recall the respective results in [8]: There, the phase boundaries of 2-and 2+1-flavour QCD have been computed within the fRG approach. These results allow us to evaluate the reliability of even linear scalematchings of temperatures and chemical potentials in 2and 2+1-flavour QCD introduced by With such a linear scale-matching the scaling factors c T , c µ B can be determined by evaluating the relations at a specific temperature and chemical potential. For the thermal scale-matching we naturally take (T, µ B ) = (T c , 0), the crossover temperature at vanishing chemical potential. In [8] the crossover temperatures have been determined with thermal susceptibilities of the renormalised light chiral condensate. Then, the linear rescaling of the 2-flavour chiral crossover temperature to the 2+1-flavour crossover temperature is done with For the matching of the chemical potentials we use the curvature κ of the phase boundary at vanishing µ B = 0, Adjusting the 2-flavour curvature −κ µ 2 B /T 2 c to the 2+1flavour one leads us to the relation The value c QCD µ B ≈ 1 entails that the change in the curvature coefficient κ is balanced by that of the temperature.
The fourth-order expansion coefficient λ is found to be very small in both functional, [8,27,28] as well as lattice computations, [10,79]. Moreover, the results for the phase boundary at finite chemical potential in [7,8,27,28] reveal that the phase boundary is still described well by the leading order expansion with µ 2 B -terms. We estimate, that this prediction is quantitatively reliable within µ B /T 4, using results from [7,8,25,27,28]. This covers the regime studied in the present work.
Applying the two scale-matching relations in Equation (4) with the coefficients Equation (5) and Equation (7) to the 2-and 2+1-flavour data of the QCD phase boundary in [8] leads us to Figure 2. In conclusion, this impressive agreement provides non-trivial support for the scale-matching procedure in QCD.

2-to 2+1-flavour scale-matching in LEFTs
The convincing quantitative accuracy of the linear scale-matching analysis presented for QCD in the last section also sustains its use in the LEFT within the present work. Note however, that we cannot simply Analogously to QCD we choose the chiral crossover temperature at vanishing chemical potential, (T, µ B ) = (T c , 0) for fixing the scale factor c T . Moreover, in the present work we are interested in fluctuations of conserved charges. Hence, instead of the renormalised condensate we use the kurtosis of baryon number fluctuations, or rather R B 42 = χ B 4 /χ B 2 , for the definition see Equation (13) and Equation (14) with Equation (16), Equation (18). This leads us to the following determination of c T : While the temperature-dependence of R B 42 is a prediction of the LEFT, its absolute temperature has to be adjusted. This is done by minimising the χ 2 of the difference between the lattice result and the LEFTprediction as a function of the rescaled absolute temperature c T T c , leading us to The respective result for R B 42 is shown in Figure 3 in comparison to the lattice result from [38]. The two curves match quantitatively supporting the predictive power of the LEFT.  (14), using the 2+1-flavour lattice results of [38]. This leads to c T = 1.247 (12) in Equation (5). The T /Tc-dependence of R B 42 is a prediction of the QCD-assisted LEFT, and agrees quantitatively with the lattice results.
Having adjusted the temperature with results from the WB-collaboration, [80], we use κ = 0.0153(18) from [10] for internal consistency. Note that the results presented here do only change marginally if using one κ in the range κ = (0.0142 − 0.0153). Within the current LEFT we obtain κ LEFT = 0.0193. In comparison, κ LEFT is larger than the 2-flavour QCD result in [8] with κ = 0.0179 (8). This reflects the lack of glue-dynamics in the LEFT. We use this in the relation Equation (7) instead of κ (N f =2) , and arrive at with the LEFT-c T from Equation (8).
In summary, as our first step towards a quantitative prediction for hyper-order baryon number fluctuations, in this work we do not deal with the strange quark as a dynamical degree of freedom for the moment, but rather take into account its effect on the modification of the momentum-scale running via an appropriate scalematching as shown in Equation (4). The validity of the scale-matching relations between 2-and 2+1-flavour QCD has been well verified in this section by means of results from the first-principle functional QCD in [8]. These relations were applied to the present QCD-assisted LEFT. The scale-matching was done with two observables relevant for the fluctuation physics studied here: 42 as a function of T and the curvature of the phase boundary κ, both at vanishing chemical potential. This led us to the coefficients Equation (8) and Equation (9) in Equation (4). The errors in these coefficients determine the errors of our results in Section IV.

III. THERMODYNAMICS AND HYPER-ORDER BARYON NUMBER FLUCTUATIONS
The thermodynamic potential in the LEFT at finite temperature and baryon chemical potential is readily obtained from the effective action in Equation (1), or rather from its integrated flow: we evaluate the effective action on the solution of the quantum equations of motion (EoMs). In the present work we consider only homoge- while the quark fields vanish on the EoMs, q,q = 0. We also note that the assumption of homogeneous solutions has to be taken with a grain of salt for larger chemical potentials with µ B /T 4, see [8]. Such a scenario has been investigated in LEFTs, see e.g. the review [81] and references therein.
With these preparations we are led to the grand potential Ω[T, µ B ] = V k=0 (ρ, L,L), the effective potential, evaluated at vanishing cutoff scale k = 0. It reads where the gluonic background field A 0 in Equation (2) has been reformulated in terms of the Polyakov loop L and its complex conjugateL. As mentioned before, the matter sector of the effective potential is integrated out towards the IR limit k = 0, for details see Appendix B. In turn, the glue sector is independent of k, see Appendix D.
The pressure of the system follows directly from the thermodynamic potential, The generalised susceptibilities of the baryon number χ B n are defined through the n-th order derivatives of the pressure w.r.t. the baryon chemical potential, to wit, To remove the explicit volume dependence, it is advantageous to consider the ratio between the n-and m-th order susceptibilities, defined by, The generalised susceptibilities are related to various cumulants of the baryon number distribution, which can be measured in heavy-ion collision experiments through the cumulants of its proxy, i.e., the net proton distribution, see, e.g., [4] for details. For the lowest four orders we get, with · · · denoting the ensemble average and δN B = N B − N B . Thus the mean value of the net baryon number of the system is given by , and the kurtosis , respectively. In this work the emphasis is, however, put on the baryon number fluctuations of orders higher than the fourth, i.e. χ B n>4 , which are named hyper-order baryon number fluctuations. As the low-order ones, the hyperorder susceptibilities are also connected to their respective cumulants, and their relations, taking the fifth through eighth ones for instance, are given as follows, Different aspects of hyper-order fluctuations have been studied in mean-field approximations in the past, see e.g. [66,67,82]. However, due to the decisive role of nonperturbative quantum fluctuations for these quantities, as functions of the temperature at vanishing baryon chemical potential (µB = 0). Results from the QCD-assisted LEFT are compared with lattice results from the HotQCD collaboration [9,36,37] and the Wuppertal-Budapest collaboration (WB) [38]. The inset in the plot of R B 82 shows its zoomed-out view. Our results agree quantitatively with the WB-results , and are qualitatively compatible with the HotQCD results. We also compare to the predictions of a hadron resonance gas (HRG) [30], which predicts only a very mild increase of R B n2 from unity with increasing T .
a treatment beyond mean-field, as in the present work, is necessary for their accurate description. A first study in this direction, discussing hyper-order fluctuations up to χ 8 within a PQM model with the fRG at small µ B /T can be found in [60].

IV. NUMERICAL RESULTS AND DISCUSSIONS
In this section we present and discuss our numerical results for hyper-order fluctuations on the freeze-out curve. At vanishing chemical potential the lower orders are compared to results from lattice calculations. We then discuss the implications of our predictions for the hyperorder baryon number fluctuations for decreasing collision energies (increasing chemical potential) for heavy-ion collision experiments.  We start our discussion of the numerical results in our QCD-assisted low-energy effective theory with benchmark results at vanishing chemical potential, µ B = 0. We have already seen in Section II B that the fourth order fluctuations R B 42 , Equation (14), agrees quantitatively with the respective lattice result, see Figure 3 and Figure 4, left panel. We emphasise again that the thermal dependence of R B 42 is a prediction of the present LEFT. Now we also compare the temperature dependence of the hyper-fluctuations R B 62 and R B 82 with the corresponding lattice results in the middle and right panels of Figure 4, respectively. We depict both our numerical results and lattice results from the HotQCD collaboration, [9,36,37], and the Wuppertal-Budapest collaboration, [38]. Note that lattice results in Figure 4 by the Wuppertal-Budapest collaboration in [38], and R B 62 and R B 82 by the HotQCD collaboration in [36] are not continuum extrapolated.
With the increase of the order of fluctuations, the uncertainties of lattice results increase significantly. Moreover, the eighth-order fluctuations R B 82 obtained by the two collaborations show a significant quantitative difference, although their form is qualitatively consistent with each other.
The hyper-order baryon number fluctuations computed in the current setup are in qualitative agreement with both lattice results. However, our results single out the lattice results of the Wuppertal-Budapest collaboration, with which we observe quantitative agreement. This situation is very reminiscent of the pressure prediction in [53]: similarly to the current situation with lattice predictions of hyper-fluctuations, the pressure predictions of the lattice collaborations had not converged yet. A less advanced version of the current QCD-assisted LEFT framework then predicted the correct pressure re- sult. We have also computed the hyper-order fluctuations within a simple hadron resonance gas model [30]. Essentially, they are all constant with only a very minor monotonic increase with T for T 140 MeV, starting from unity at T = 0. This is in quantitative agreement with our findings at low temperatures. In summary, the current setup passes all benchmark tests quantitatively and provides the full temperature-dependence of hyperfluctuations. We have also computed even higher order baryon number fluctuations. In Figure 5 we show our result for the temperature-dependence of the tenth order ratio R B 10 2 = χ B 10 /χ B 2 at vanishing chemical potential, µ B = 0. So far no lattice results for the tenth-order fluctuation are available, and the dependence of R B 10 2 on the temperature in Figure 5 is a prediction by the current QCD-assisted LEFT and awaits confirmation by other calculations, e.g., lattice QCD, in the future.

B. Hyper-order baryon number fluctuations at finite density
With successfully passing the benchmark tests, we proceed with the results for baryon number fluctuations at finite chemical potential. This will allow us to finally compare the theoretical predictions on the freeze-out curve with the experimental measurements in Section IV D.
Equally relevant is the self-consistent evaluation of the reliability of a Taylor expansion in baryon-chemical potential that underlies the extension of lattice results at vanishing chemical potential to µ B = 0. This is particularly important for predictions of the location of the critical end point based on such an expansion. Here we can investigate the reliability range of the Taylor expansion around µ B = 0 by comparison to our direct computation at finite µ B .
First we investigate the temperature-dependence of the baryon number fluctuations for different chemical potentials. This also allows us to discuss the reliability bounds of the current LEFT-setup for increasing chemical potential. In Figure 6 we show the temperature-dependence of the ratios R B 42 , R B 62 and R B 82 for chemical potentials µ B = 0, 100, 160, 200, 300, 400 MeV. First, we note that at small temperatures the thermodynamic properties of the QCD medium are well described by a dilute gas of hadrons, where the net-baryon number follows a Skellam distribution. Thus, all ratios approach unity at vanishing temperature. At very large temperature the system is governed by asymptotically free quarks, where fluctuations approach to the trivial Stefan-Boltzmann limit, and R B n2 goes to zero for all n > 4 at large T . Consequently, the non-trivial behaviour of the fluctuations between these two limiting cases shown in Figure 6 is directly related to the crossover from the hadronic-to the quark-gluon regime of QCD. The magnitude, but also the error of the fluctuations grow with increasing chemical potential. Both effects are more pronounced for higher order fluctuations. The increase in magnitude is directly linked to the sharpening of the chiral crossover with increasing chemical potential, cf. Figure 2. We expect that the current LEFT-setup is gradually loosing its predictive power for fluctuations on the freeze-out curve due to the rapid increase of the computational error for higher-order fluctuations at large baryon chemical potential, e.g., R B 82 with µ B 200 MeV. All results of the subsequent investigations have to be evaluated with this estimate on our systematic error.
For the evaluation of the reliability regime of the Taylor expansion about vanishing chemical potential we consider the Taylor expansion of the pressure in Equation (12) in powers ofμ B ≡ µ B /T aroundμ B = 0. This leads us to with the expansion coefficients χ B 2i (0) = χ B 2i (µ B = 0), the hyper-order fluctuations of the baryon charge. In Equation (23) we have suppressed the temperaturedependence of all functions for the sake of readability. Truncating the Taylor expansion in Equation (23) (24). Our results from the QCD-assisted LEFT are compared to those from lattice QCD by the HotQCD collaboration [9] and the Wuppertal-Budapest collaboration [38]. We show the comparison to HotQCD in the inlays, as these results deviate considerably from both ours and the WB results.
tain the expanded baryon number fluctuations, In Figure 7 we show the ratios R B 42 = χ B 4 /χ B 2 and R B 62 = χ B 6 /χ B 2 based on the Taylor expansion for two fixed temperatures: T = 155 MeV (close to the crossover temperature T c at µ B = 0) and T = 160 MeV (slightly above T c ). As an input we use χ B 2i (0) (i = 1, 2, 3, 4) from the current setup as well as from the lattice (HotQCD collaboration [9] and Wuppertal-Budapest collaboration [38]), depicted already in For more details about effects of these constraints on the fluctuations and correlations of conserved charges, see the relevant discussions in, e.g., [9,36] in lattice QCD and [56,57] in fRG.
Since we are not restricted by a sign problem within the fRG approach, the χ B n (µ B )'s in Equation (13) can also be computed directly for the current QCD-assisted LEFT without resorting to a Taylor expansion. By comparing this to the results of the Taylor expansion, we can study its range of validity. The results are presented in the left panel of Figure 8, again for T = 155 MeV and T = 160 MeV and with the Taylor expansion up to eights order inμ B .
We observe that the result for R B 42 from the Taylor expansion in Equation (24)  We emphasise that this is not an artefact of our model, but rather a generic, physical feature of these fluctuation observables. It reflects the increasingly pronounced non-monotonic temperature dependence of R B n2 due to long-range correlations in the crossover region, as seen in Figure 6. In particular, R B n2 develops distinctive areas around the crossover at larger µ B where its value varies significantly, even including sign changes. By following a line of fixed T close to T c and increasing µ B in the phase diagram, these areas are crossed eventually, leading to the oscillatory behaviour seen in Figure 8. This is also evident in the right plot of Figure 11, where we show the magnitude of R B 42 in the phase diagram. Since these strong fluctuations only occur at larger µ B , the resulting characteristic qualitative features cannot be captured by a (low-order) Taylor expansion at µ B = 0; it is bound to fail at the onset of the oscillatory behaviour, i.e. around µ B /T 1.
It is also interesting to evaluate to what extend a higher-order Taylor expansion can improve its reliability. We therefore include our prediction for R B 10 2 from Figure 5, hence extending the Taylor expansion in Equation (24) to the tenth order. The resulting comparison is shown in right panel of Figure 8 (25) might be increased a bit when terms of orders higher than the tenth one are included in the Taylor expansion. However, as we have discussed above, the full results of the baryon number fluctuations show a non-trivial oscillatory behaviour, which is generic, and stems from the fact that the system crosses over different phases with the increase of µ B , cf. the right panel of Figure 11. For a discussion of the radius of convergence based on mean-field theory we refer to [66]. In general, it is given by the distance to the nearest singularity of the equation of state in the complex µ B /T plane. Hence, possible candidates for this singularity are the CEP, the Roberge-Weiss endpoint at imaginary µ B [83], or the Yang-Lee edge singularity in the complex plane [84]. Evidently, the distance to the CEP is far larger than the radius estimated here. In turn, the closest endpoint at imaginary chemical potential is at |µ B /T | ≤ π. For physical quark masses it is most probably close to |µ B /T | = π. For a discussion within QCD-flows see [85], for lattice results see e.g. [86]. Hence, a particularly intriguing option is the Yang-Lee edge singularity, for a discussion see e.g. [87,88]. How-ever, while the location of the edge singularity has been determined for critical O(N ) theories [89], it is still unknown for QCD.
This interpretation also implies that the results from the Taylor expansion fail to agree even qualitatively with the correct µ B /T -dependence for µ B /T c [µ B /T ] Max (T ), see Figure 8. Interestingly, [µ B /T ] Max (T ) seems to grow for smaller temperatures. Whether or not this holds true requires a more systematic study, which will be considered elsewhere. In conclusion, the extrapolation of fluctuations of conserved charges in the vicinity of the chiral crossover within a Taylor expansion loses its predictive power for µ B 200 MeV, at least at tenth order.
In Figure 9 we show our full results for the temperature-dependence of R B 31 , R B 51 , and R B 71 with different values of baryon chemical potential. A further relevant odd fluctuation observable is R B 32 , depicted in Figure 10. Its experimental analogue, the proton number fluctuation R p 32 has been already measured in Au+Au central (0-5%) collisions at STAR, a comparison will be presented and discussed in Section IV E.

C. Determination of the freeze-out curve
The quantitatively successful benchmark tests analysed in Section IV A, and the evaluation of baryon number fluctuations at finite chemical potential in Section IV B allow us to discuss our main goal: the comparison of theoretical predictions on the baryon number fluctuations with experimental measurements.
All these different effects and experimental restrictions give rise to non-critical contributions to fluctuation observables in experiments, and pinning down their contributions plays a pivotal role in identifying the critical signals in the BES experiment. Additionally, due to critical slowing down, non-equilibrium effects become important in the vicinity of the CEP [111], which necessitates a theoretical description of the dynamics of critical fluctuations. For more details about recent progress on the dynamics of critical fluctuations in QCD, see [112] and references therein. We emphasise, however, that the present results of QCD-assisted LEFT model are well outside the critical region. Therefore they are not subject to critical scaling in the vicinity of the CEP.
In this work we will not take into account the noncritical and dynamical effects discussed above. Instead, we assume that the measured cumulants of the netproton multiplicity distribution at a given collision energy are in one-to-one correspondence to the calculated fluctuations in Equation (13) with single values for T and µ B (with other collision parameters e.g., the centrality and rapidity range fixed). Then, it is suggestive to identify the values of T and µ B with the ones when the chemical freeze-out occurs, viz. T CF and µ B CF . Such an approach for the comparison is usually employed in fluctuation studies of equilibrium QCD matter within functional methods or lattice simulations, see e.g. [9,26,41,42,62].
We adopt the freeze-out temperatures and baryon chemical potentials from [90] and from the STAR experiment [91], which are shown in the left panel of Figure 11 by the blue pentagons and red circles, respectively. They are both obtained from the analysis of hadron yields in the statistical hadron resonance gas model, see the aforementioned references for more details. The freeze-out data in [90] has also been parametrised as functions of the collision energy as follows with a = 1307.5 MeV, and with T (0) CF = 158.4 MeV. This parametrisation is depicted with the blue dashed line in the left panel of Figure 11. We use the same parametrisation functions in Equation (26) to fit the freeze-out data in STAR experiment, i.e., the red circle points. For this fit we invoke two procedures, called STAR Fit I & II in the following: For the first one, STAR Fit I, we simply take all 7 data points. The corresponding freeze-out curve is depicted by the red solid line in the left panel of Figure 11.
For the construction of the second one, STAR Fit II, we shall argue that some of the data points are potentially flawed, or rather await a physics explanation, and should be dropped accordingly in a fit based on Equation (26). Accordingly, we drop the first two data points at small chemical potential as well as the last one at the largest available chemical potential µ B ∼ 400. From general considerations we do not expect the freeze-out curve to rise with increasing chemical potential. Moreover, the physically motivated fit formula does not describe signchanges of the curvature of the freeze-out curve. For a respective discussion and possible explanation for the only apparent rise see [113]. The last data point is also not well-described by the fitting procedure described here. This may indicate the onset of a regime with different physics/phases. In this case, Equation (26) would not be an appropriate fit function. It may also indicate the onset of a regime of rapidly worsening systematics. In this case more points are needed in this regime. The freeze-out line of STAR Fit II is depicted by the green dotted line in the left panel of Figure 11. In comparison to STAR Fit I, STAR Fit II is located at slightly lower temperatures, which is more pronounced when µ B 200 MeV. In the right panel of Figure 11, we show the baryon number fluctuation R B 42 in the T − µ B plane. It can be observed that a narrow blue band, indicating the regime of negative R B 42 , develops around the crossover starting at µ B ∼ 250 MeV. The freeze-out curve STAR Fit II is approaching towards the boundary of the blue region firstly at small µ B , and then deviates a bit from it at large µ B . We emphasise that the large chemical potential region, and in particular asymptotically large µ B 500 MeV, is beyond of the reliability bound of the current computation, µ B /T = 4, see Figure 2. For a detailed discussion see Section II B. The determination of the freeze-out curve completes our setup, which enables us to compute hyper-order baryon number fluctuations R B nm along the freeze-out line within the QCD-assisted LEFT. These results are then used to compare with the experimental measurements of cumulants R p nm of the net-proton distribution from STAR experiment.
Before we discuss the numerical results, we also emphasise once more, that it follows from the analysis of Section IV B, that the simple extrapolation with the Taylor expansion about µ B = 0 lacks predictive power for µ B 250 MeV, that is √ s NN 15 GeV, see Equation (25). Moreover, it even lacks predictive power for the qualitative behaviour.
In the left panel of Figure 12 we show the √ s NNor chemical potential dependence of the baryon number fluctuations R B 42 , R B 62 , and R B 82 for the freeze-out lines from Andronic et al. [90] and STAR Fit I. The freeze-out line from Andronic et al. is obtained from an interpolation of the freeze-out data, the grey squares in Figure 11.
In the right panel of Figure 12 we show the same observables for the freeze-out line of STAR Fit II. As discussed in Section IV C, we have singled out the results for this freeze-out curve as the best-informed computation.
In both panels of Figure 12 we also show the experimental measurement of cumulants of the net-proton distributions in the beam energy scan experiments from the STAR collaboration. The fourth-order fluctuations, R p 42 , of the net-proton multiplicity distributions are measured in Au+Au collisions with centrality 0-5%, transverse momentum range 0.4 < p T (GeV/c) < 2.0, and rapidity |y| < 0.5, cf. [22] for more details. Moreover, results for the sixth-order cumulant of the net-proton distribution, R p 62 , are also presented in the middle plot of Figure 12 out curves considered are compatible with the respective experimental measurement of the κσ 2 of net-proton distributions in 0-5% central Au+Au collisions. In particular, the theoretical results feature a non-monotonic √ s NN -dependence: R B 42 first decreases with decreasing beam energy and then increases. The details of this behaviour, in particular how pronounced it is, is highly sensitive to the precise location of the freeze-out. For example, the increase at small √ s NN is larger for smaller freeze-out temperatures. Thus, the weak increase for STAR Fit 1 originates in the slightly larger freeze-out temperature of this freeze-out fit. This shows that even small variations in the freeze-out temperature have a substantial effect on the fluctuations in this regime. The underlying reason is that the freeze-out happens in or close to the crossover region, where the fluctuations vary significantly, see Figure 6. Importantly, this regime cannot be accessed within the extrapolation of the Taylor expansion at least within the current order. This entails that extrapolations based on a Taylor expansion are bound to fail to describe the data in this regime reliably. Consequently this calls for qualitatively improved direct theoretical computations at small beamenergies. This is work in progress and we hope to report on the respective results soon.
Our results for the hyper-order fluctuations R GeV. Another interesting property of the current LEFT setting is that the non-monotonic behaviour of our results for R B n2 at large µ B in Figure 12 does not arise from critical physics: in the LEFT used here, the CEP is at significantly larger µ B 700 MeV. Moreover, it is well established that the critical region is only very small. It is already small within mean-field approximations of lowenergy effective theories, and additionally shrinks considerably if quantum, thermal and density equilibrium fluctuations are taken into account, see [29]. Moreover, this does not change if transport processes are taken into account, see [114].
In the present LEFT the increasing trend at large µ B region originates from two effects: First, fluctuations are enhanced since the chiral crossover becomes sharper with increasing µ B . This leads to a stronger non-monotonic behaviour of R B 42 as a function of T , see Figure 6. This sharpening is also present in the vicinity of a CEP. Second, the freeze-out temperature is shifted away from the pseudo-critical temperature towards small  [90] (grey) and the STAR experiment [91] (red). Right panels: The freeze-out curve, STAR Fit II, is obtained from the freeze-out parameters of the STAR experiment [91]. The theoretical error bands show a highly correlated error, and should be interpreted as a family of curves with the same qualitative behaviour as the central curve. For more explanations see Section IV C with Figure 11. STAR data: R p 42 (top) is the kurtosis of the net-proton distributions measured in Au+Au central (0-5%) collisions [22]. R p
beam-energies, thereby probing different regimes of the cumulants. However, it should be noted that the uncertainty of our results increases significantly in the low energy region. These uncertainties include an estimate for the systematic error of the QCD-assisted LEFTapproach. This systematic error stems from the uncertainty in the matching of the in-medium scales of the LEFT and QCD, encoded in the coefficients c T and c µ in Equation (4). Moreover, for larger chemical potential the current LEFT lacks the back-reaction of the µ B -dependence of the glue-dynamics. While inherently small, it might still play a rôle. Furthermore, it is found that R B 42 can also be suppressed in the regime of low collision energies due to the effect of global baryon number conservation, cf. [107,108,115] which will be included in our future work.

E. Search for the CEP
The analysis in the previous sections entails that in the present QCD-assisted LEFT the non-monotonic behaviour of baryon-number fluctuations is triggered by the sharpening of the chiral crossover. This is highly nontrivial, since it is evident, e.g., from Figure 6 that neither a system only in the hadronic-nor only in the QGP phase could produce the beam-energy dependencies shown in FIG. 13. Baryon number fluctuation R B 32 as a function of the collision energy in comparison to STAR-data for R p 32 (0-5%) centrality [22]. Left panel: the freeze-out points are those from Andronic et al. [90] (grey) and the STAR experiment [91] (red). Right panel: The freeze-out curve, STAR Fit II, is obtained from the freeze-out parameters of the STAR experiment [91]. The theoretical error bands show a highly correlated error, and should be interpreted as a family of curves with the same qualitative behaviour as the central curve. For more explanations see Section IV C with Figure 11. In conclusion, the agreement of R B 42 , computed in the QCD-assisted LEFT and the measured beam-energy dependence of R p 42 shows, that the latter could be a signature for the presence of a sharpening crossover between these two phases. Whether or not it also signals the onset of the critical region will be subject of a future improved study. In the context of this latter study we also emphasise that the non-universal properties of the LEFT such as the existence and location of the CEP may not quantitatively agree with QCD as the latter regime lies outside the LEFT-regime with quantitative reliability. Still, the present LEFT probably has the same qualitative nonuniversal properties at large chemical potential, and it certainly has the same universal ones.
A non-monotonic energy dependence for the fluctuations is a highly relevant experimental observation, since this behaviour has been proposed as an experimental signature of a CEP [14,16]. The present analysis based on QCD-assisted LEFT model demonstrates that the nonmonotonic behaviour of fluctuations can serve as an indication of a CEP, but is not necessarily a smoking gun signature for it. The latter requires the extraction of critical scaling, or similar definite signatures such as the detection of a first order regime for large µ B , etc. . Still, the non-monotonic behaviour observed in both theory and experiments is a clear signature for interesting strongly correlated physics, whose uncovering requires joint and intensified effort of both, theory and experiment. Of course, whether or not these properties carry over completely to QCD remains to be seen.
Note also that the non-monotonic regime is far away from that covered by a simple extrapolation of the Taylor expansion at µ B = 0. It might be covered by a resummation of the latter, which can already be investigated within the present QCD-assisted LEFT. Constraints on such a resummation should also make use of odd hyper-order fluctuations at finite chemical potential, that are readily computed in the present setup: A prominent and relevant example is R 32 , already measured in the STAR experiment. In Figure 13 we show our predictions for R B 32 computed in the current QCDassisted LEFT on the different freeze-out curves defined in Section IV C, see in particular Figure 11. In this section it also has been argued, that our best-informed freeze-out curve is given by STAR Fit II. The respective results are shown in the right panel of Figure 13 in comparison with the STAR data for R p 32 (0-5% centrality). Indeed, these results show the best compatibility with the experimental data. Moreover, within the respective systematic and statistical errors the theoretical results with the freeze-out curve STAR Fit II and the experimental data agree down to collision energies of √ s NN ≈ 14.5 GeV or µ B ≈ 250 MeV.
Interestingly, below √ s NN ≈ 14.5 GeV the experimental data show a plateau, which is not present in the theoretical prediction. While this is in the large chemical potential regime, in which the LEFT gradually looses its predictive power, also the respective functional first principles QCD computation in [7,8,27,28], based on a grand potential, do not show any sign of new physics in this regime. This suggests that for √ s NN 14.5 GeV at least one of the implicit assumptions underlying the identification of R B nm with R p nm within a grand canonical ensemble with variable baryon charge (density) for given beam energies breaks down. As discussed before, this asks for a re-assessment of the identification The freeze-out curve, STAR Fit II, is obtained from the freezeout parameters of the STAR experiment [91]. The theoretical error bands show a highly correlated error, and should be interpreted as a family of curves with the same qualitative behaviour as the central curve. For more explanations see Section IV C with Figure 11. of baryon and proton number fluctuations, finite volume effects and finite volume fluctuations, the determination of the freeze-out curve for smaller collision energies, the evaluation of non-equilibrium effects such as transport, and finally the use of the grand potential in the theory computations. While highly relevant and interesting, this goes far beyond the scope of the present work and we defer any further investigation to future work.
The above example of R 32 demonstrates very impressively, that the odd (hyper-order) fluctuation observables encode highly relevant physics information which may be difficult or even impossible to extract from the even or-ders. As a first step in this direction, finally aiming at a resummation of the µ B -expansion that allows us to go beyond the validity regime of the Taylor expansion, we also have computed the fluctuation observables R 31 , R 51 , R 71 on the freeze-out curve STAR Fit II in Figure 14. An experimental confirmation of the respective predictions at least for the lower orders would be highly desirable.
The discussion in this section leaves us with the highly exciting possibility of unravelling the location and properties of a potential CEP within a combined experimenttheory analysis: First principle QCD at finite density should provide us with a prediction for the location of the CEP in terms of hyper-order fluctuations, Loc CEP (R nm ). This would allow us to use the experimental data on hyper-order fluctuation observables R p nm as input. We emphasise that this prediction does not necessitate the observation of critical behaviour in the R nm , but utilises the details of the non-monotonicity of the R nm .
In summary, such an analysis does explicitly not rely on the universal property of critical scaling measured in the R nm . Indeed, it uses the non-universal properties of the R nm to predict the non-universal location of the CEP, and hence is far more robust. We hope to report on this matter in the near future.

V. CONCLUSIONS
In this work we have computed baryon number fluctuations up to tenth order with a QCD-assisted lowenergy effective theory. This LEFT incorporates quantum, thermal and density fluctuations from momentum scales less than 700 MeV within the functional renormalisation group approach, and is embedded in QCD, for details see Section II. The quantitative predictability has been benchmarked with a comparison of baryon number fluctuations at µ B = 0 up to the eighth order from the lattice, see Section IV A. Our results are in quantitative agreement with that from the Wuppertal-Budapest collaboration, and are compatible with that of the HotQCD collaboration, as shown in Figure 4.
Our direct computation at finite µ B , presented in Section IV B, has allowed us to assess the range of validity of the Taylor expansion of the free energy of QCD around µ B = 0. Such an expansion is commonly used to extrapolate lattice results to finite density. We have shown that the expansion up to tenth order in µ B /T is only valid for µ B /T 1.5 in the chiral crossover regime, see Figure 8. Beyond this range, the Taylor expansion, at least to this order, fails to even capture the qualitative behaviour of the fourth-and sixth-order baryon number fluctuations. Thus, results for fluctuations at the freezeout curve based on a Taylor expansion around µ B = 0 should be interpreted with great caution for µ B /T 1.5, as relevant physical effects might not be captured by this extrapolation.
The main goal of the current work was the computation of baryon number fluctuations and in particular hyper-order fluctuations along the freeze-out curve at collision energies √ s 7.7 GeV. The respective results are discussed in Section IV D. They have been compared to experimental data of net-proton number cumulants from STAR for different estimates for the freeze-out curve, see Figure 12. Our result for the kurtosis, R B 42 , is in good agreement with the experimental data for collision energies √ s 7.7 GeV. In particular the increasing trend at lower beam energies √ s 19.6 GeV is captured well.
This non-monotonicity is also present in the hyper-order fluctuations R B 62 , R B 82 . We also note that a comprehensive comparison for the higher order cumulants is not possible due to the lack of experimental data. Accordingly, our results in Figure 12 are predictions that await experimental verification.
We have also investigated the twofold origin of the non-monotonicity for √ s 19.6 GeV in the present LEFT: First, for increasing chemical potential the chiral crossover gets sharper. Secondly, for smaller beam energies the freeze-out temperature may move away from the pseudo-critical temperatures. In the current setup both phenomena happen far away from a potential critical end point in the LEFT located at √ s CEP 3 GeV. The latter regime is also safely outside the reliability regime of the current setup, which gradually looses reliability for In summary, the non-monotonicities of hyper-order fluctuations, observed both in experiment and theory, are important signatures for interesting physics in the border regime between quark-gluon plasma and the hadron phase. This of course can include a potential CEP, and in any case deserves further investigation from both exper-iment and theory. In particular, we envisage that experimental data of fluctuation observables and their dependence on collision energy allow us to constrain the onset regime of this strongly correlated physics/CEP. Importantly, such a prediction does not rely on the observation of critical scaling in the hyper-order fluctuations, but is far more robust, for more details see Section IV E. We hope to report on this in the near future. The work is also supported by EMMI, and the BMBF grant 05P18VHFCA. It is part of and supported by the DFG Collaborative Research Centre SFB 1225 (ISOQUANT) and the DFG under Germany's Excellence Strategy EXC -2181/1 -390900948 (the Heidelberg Excellence Cluster STRUCTURES).

Appendix A: The fRG-approach to QCD & LEFTs
The functional renormalisation group or flow equation for QCD provides the evolution of its effective action Γ k with an infrared cutoff scale k. Here we use the setup with dynamical hadronisation, [8,[47][48][49][50]116]. The formulation used here has been developed in [8,[43][44][45][46]. Its current form has been described and further developed in [8], and for further details we refer to this work. The flow equation of the QCD effective action reads, In Equation (A1), the Φ = (A, c,c, q,q, φ) is a superfield that comprises all fields. This also includes hadronic (composite) low energy degrees of freedom introduced by dynamical hadronisation. The G's and R's are the propagators and regulators of the different fields, respectively. Diagrammatically it is depicted in Figure 1. For more works on QCD-flows at finite temperature and density see [8, 43-46, 85, 116-121], for reviews on QCD and LEFTs for QCD see [49,51,[122][123][124][125][126][127][128].
For scales k 1 GeV, the gluon decouples from the system due to its confinement-related mass gap. For these momentum scales, the (off-shell) dynamics of QCD is dominated by quarks and the emergent composite hadronic degrees of freedom. In particular, the lowest lying meson multiplet, and specifically the π meson is driving the dynamics. The pion is the pseudo-Goldstone boson of strong chiral symmetry breaking, and hence is the lightest hadron with a mass ∼ 140 MeV in the vacuum.
Consequently, in this regime with k 1 GeV, the flow equation of the QCD effective action in Equation (A1) is reduced to, where R q,k and R φ,k are the regulators for the quark and meson fields, respectively. The full propagators in Equation (A2) read, with Γ (2) In this work we employ 3d-flat or Litim regulators [129][130][131], with where Θ(x) denotes the Heaviside step function. Inserting the effective action Equation (1) into the flow equation Equation (A2), we arrive at where the threshold functions l (B/F,4) 0 as well as other threshold functions used in the following can be found in e.g., [8,55]. The dimensionless renormalised quark and meson masses read, The anomalous dimensions for the quark and meson fields in Equation (A6) are defined as respectively. Accordingly, the flow equation for the mesonic anomalous dimension is obtained from the (spatial) momentum derivative w.r.t. p 2 of the pion twopoint function, to wit, (A10) The approximation Equation (1) to the effective action together with Equation (A10) are based on two approximations: Firstly, in Equation (1) we have dropped the field-dependence of Z φ , which would lead to different Z π and Z σ . In Equation (A10) we have identified Z φ = Z π , and hence also Z σ = Z π . This is motivated by the fact that the meson dynamics are only dominant in the broken regime where the three pions are far lighter than the single sigma mode, which quickly decouples. Hence, the three pions drive the dynamics. Furthermore, in Equation (1) we do not distinguish between spatial and temporal components of Z φ . For finite temperature and density, the Euclidean O(4) rotation symmetry is broken, as the heat bath of density singles out a rest frame. This entails, that η φ,k splits into η ⊥ φ,k and η φ,k , the components transverse and longitudinal to the heat bath/density. We have used the approximation η φ,k = η ⊥ φ,k as we have three spatial directions. The influence of the splitting of η φ,k on the thermodynamics and baryon number fluctuations has been investigated in detail e.g. in [55]. There it has been found that the impact is small, supporting the reliability of the present approximation.
Similarly, the quark anomalous dimension is obtained by projecting the relevant flow onto the vector channel of the 1PI quark-anti-quark correlation function, In Equation (A11), the spatial momentum is set to zero, p = 0 as in the mesonic case: vanishing momentum is most relevant to the flow of effective potential in Equation (A6). Note, that the lowest fermionic Matsubara frequency is non-vanishing. We use p 0,ex = 0, its value is further described in Appendix B, based on [40,42]. As is implicit in Equation (A11), the flow of the quark two-point function is complex-valued at non-vanishing chemical potential. This originates in the Silver-Blaze property of QCD at T = 0. For quark correlation functions this entails that they are functions of p 0 − iµ q already before the onset of the baryon density, for a discussion in the present fRG-approach see [40,42,132]. In turn, the couplings are still real (i.e. real functions of the complex variable p 0 −iµ q ) in particular below the density onset. Hence, couplings (i.e. expansion coefficients in a Taylor expansion in momenta) are real. This is readily seen in a resummation of the external frequency of the quark propagator [42]. Without resummation they are obtained from a projection on the real part of the flow, see Equation (A11).
This projection is also used for the Yukawa coupling. Within the present approximation, the flow equation of the (real) Yukawa coupling is given by, The explicit expressions for the meson and quark anomalous dimensions, as well as the flow of the Yukawa coupling can be found in Appendix B.
In the threshold function FB's we have employed p 0,ex = πT for the finite temperature sector and p 0,ex = πT exp{−k/(πT )} for the vacuum sector. This choice guarantees a consistent temperature dependence for all k, which is particularly relevant for the thermodynamics in the low temperature region [40]. This can be resolved by means of a full frequency summation of the quark external leg [42], and the present procedure mimics this physical behaviour. The flow of the Yukawa coupling in Equation (A12) reads, (1,1) (m 2 q,k ,m 2 σ,k , η q,k , η φ,k ; T, µ, p 0,ex ) − (N 2 f − 1)L (1,1) (m 2 q,k ,m 2 π,k , η q,k , η φ,k ; T, µ, p 0,ex ) . The parameters are fixed with the pion pole mass m π,pol , the mass of the sigma resonance, mσ, and the pion decay constant, fπ ≈σEoM, the expectation value of the sigma-field. The input parameters are those of the initial effective potential VΛ: the meson self-coupling λΛ and the meson mass parameter νΛ.
The pion mass is tuned by cσ, the parameter of explicit chiral symmetry breaking. Finally, the constituent quark mass is fixed via the initial value for the Yukawa coupling, hΛ.
Explicit expressions of all the threshold functions mentioned above, such as BB, F's, FB's, and L can be found in e.g., [8,55]. In summary, the flow equations Equation (A6), Equation (B3), Equation (B5), Equation (B8), supplemented with Equation (B6) and Equation (B7), constitute a closed set of ordinary differential equations, which is evolved from the UV cutoff k = Λ to the IR limit k = 0.

Appendix C: Initial conditions
To solve the flow equation, we need to specify initial conditions. To this end, we choose initial values at a scale k = Λ such that known observables of QCD in the vacuum at k = 0, such as the pion mass and decay constant, are reproduced. The effective potential at the UV cutoff reads, We initialise the flows at Λ = 700 MeV. The input parameters of the QCD-assisted LEFT are given in Table I.

Appendix D: Glue potential
The dynamics of the glue sector in QCD is partly imprinted in the glue potential V glue,k (A 0 ), see Equation (2). This has been discussed in Section II. This potential is only needed for the determination of the expectation value of the Polyakov loop. Its inherent glue correlation functions are gapped and their backcoupling is suppressed for k 1 GeV. Accordingly, we can simply drop the scale dependence of the glue potential for the present purposes. This leads us to, where P on the r.h.s. stands for path ordering.
In this work we adopt the parametrisation of the glue potential in [136], which reads Both, the parametrisation of glue potential in Equation (D4), as well as determination of relevant parameters in Table II, is based on lattice results of SU(3) Yang-Mills theory at finite temperature. This potential does not only reproduce the lattice expectation value of the Polyakov loop and the pressure, but also the correct quadratic fluctuations of the Polyakov loop, [136]. These fluctuations, and higher ones, are important for the fluctuation observables discussed here [40,42]. The coefficients in Equation (D4) are temperature-dependent, x(T ) = x 1 + x 2 /(t r + 1) + x 3 /(t r + 1) 2 1 + x 4 /(t r + 1) + x 5 /(t r + 1) 2 , with x = a, c, d, and In Equation (D6) and Equation (D7) we have used the reduced temperature t r = (T − T c )/T c . The parameter values are taken from [136], and are collected in Table II for convenience.
The parameters in Table II are that of the glue potential in Yang-Mills theory. It has been argued and shown in [52,53,137] that unquenching effects in QCD are well captured by a linear rescaling of the reduced temperature in the regime about T c , very similar to the rescaling discussed in Section II B 2. This leads us to,