Baryon number fluctuations in the 2+1 flavor low energy effective theory

High order cumulants of the baryon number distribution are calculated in a 2+1 flavor low energy effective theory. Quantum fluctuations are encoded through the functional renormalization group approach. The chiral and deconfinement phase transitions are investigated at finite temperature and baryon chemical potential. The equation of state for the QCD matter and the kurtosis of the baryon number distribution are calculated in the low energy effective theory, and one finds the results are consistent with those from lattice QCD.


I. INTRODUCTION
Studying the phase structure of the hot and/or dense QCD matter has attracted lots of attentions in recent years. In the QCD phase diagram spanned by the temperature and the chemical potential, the critical end point (CEP), which separates the first-order phase transition at high chemical potential from the continuous crossover at high temperature, plays a central role due to its uniqueness and significance. Unfortunately, location, and even existence, of the CEP is hitherto unclear. Lots of efforts, however, have been made to be aimed at unravelling the mysterious veil. In the experiments significant progress has been made in the Beam Energy Scan (BES) Program at the Relativistic Heavy Ion Collider (RHIC) [1][2][3], see a review in [4].
In the meantime theoretical studies of the phase structure of dense QCD matter, as well as confrontation of theoretical calculations with experimental data, are currently ongoing research frontier. In recent year, the QCD equation of state (EoS), fluctuations and correlations of conserved charges, etc. have been computed and studied systematically in lattice QCD at finite chemical potential, see e.g. [5][6][7][8][9][10][11][12]. In particular, freeze-out parameters in heavy ion collisions are already accessible through comparison between theoretical calculations and experiments [5,10,11]. Functional continuum field approaches besides lattice QCD also provide us with a wealth of knowledge about the QCD phase structure [13][14][15] and the properties of QCD matter [16,17], especially in the regime which can not be reached by lattice simulations.
Given the importance of the baryon number fluctuations in the phenomenology of CEP searching in the experiments [41][42][43][44], in this work we would like to extend relevant studies in Refs. [20,21,41] from two to 2+1 flavors, based on the 2+1 flavor Polyakov-loop improved quark-meson (PQM) effective model, see, e.g. [45][46][47][48][49][50][51][52] for more details about the model, and Refs. [53][54][55] for other relevant 2+1 flavor low energy effective models of the same class. Note that including the strangeness is nontrivial, since it introduce the degrees of freedom of not only the strange quarks but also open-strange mesons, such as kaon mesons.
From another point of view, the 2+1 flavor low energy effective model is the hadronic sector of the 2+1 flavor rebosonized QCD, and is readily embedded into the glue sector within the FRG approach, see [26,27] for the two flavor rebosonized QCD. Therefore, the 2+1 flavor low energy effective model investigated in this work is also aimed to facilitate for the construction of the 2+1 flavor rebosonized QCD in the near future. Therefore, we would like to describe the theoretical framework in more detail. The QCD thermodynamics, such as the pressure and the EoS at finite temperature and baryon chemical potential, will be investigated, and then we will compare our calculated results of the baryon number fluctuations with those from lattice simulations. This paper is organized as follows. In Sec. II we describe the 2+1 flavor low energy effective model. The FRG and the flow equation for the effective potential are given in Sec. III. We will discuss the thermodynamics and the baryon number fluctuations in Sec. IV. Numerical results are presented in Sec. V, and then a summary and an outlook are presented in Sec. VI. Furthermore, more technical details about the flow equation of the effective potential, the glue potential, and the numerical setup are given in Appendix A, Appendix B, and Appendix C, respectively.

II. 2+1 FLAVOR LOW ENERGY EFFECTIVE MODEL
In this work we adopt the 2+1 flavor Polyakov-loop improved quark-meson model [48][49][50][51], and quantum fluctuations of different scales, as well as the thermal and density fluctuations, are encoded successively through the evolution of the renormalization group (RG) scale dependent effective action Γ k , which in the Euclidean formalism reads with Φ denoting all the field dependence, where we have used a shorthand notation x = 1/T 0 dx 0 d 3 x with the temperature T . The quark chemical potential µ is related to the baryon chemical potential through µ B = 3µ. In this work only the baryon chemical potential is taken into account, and interested readers are referred to, e.g., [51] for relevant discussions about the influences of the strangeness chemical potential on the QCD thermodynamics and phase structure. Note also that the local potential approximation (LPA) to the effective action is implicitly implied in Eq. (1), viz. only the effective po-tentialŨ k (Σ) is dependent on the RG scale k. For more discussions about truncations beyond LPA, taking into account, for instance, the nontrivial dispersion relations and the running Yukawa coupling, etc., see e.g. [50].
The temporal gluon background field A 0 in Eq. (1) is related to the color confinement and its phase transition, which could be formulated as the Polyakov loops for convenience, to wit, with where P on the r.h.s. denotes the path ordering. The Polyakov loop, from the viewpoint of statistics, can be regarded as the order parameter of the Z(3) symmetry for the deconfinement phase transition. The dynamics of the Polyakov loop is governed by the glue potential, also called as the Polyakov loop potential V glue (L,L) in Eq. (1). Usually the glue potential and its dependence on the external parameters, such as the temperature, can be parameterized by employing first-principle QCD calculations, for instance the lattice computation [56] and the FRG [57]. Then, the parameterized glue potentials are applied in low energy effective theories, for more relevant discussions, see, e.g. the review article [58] and reference therein.
In the effective action Eq. (1), quarks are coupled to the scalar and pseudoscalar meson nonets via a chirally symmetric Yukawa term with where T a are the generators of the flavor U (N f ) group, which for the superscript a = 1, ..., 8 can be represented by the Gell-Mann matrices, i.e., T a = λ a /2, and T 0 = The kinetic term for the mesons in Eq. (1) is formulated in terms of the adjoint representation of U (N f ), which reads In Eq. (1)Ũ k (Σ) is the meson effective potential, which consists of several parts as follow with U k (ρ 1 ,ρ 2 ) on the r.h.s. being an arbitrary func- Therefore, U k (ρ 1 ,ρ 2 ) has the maximal symmetry of flavors, and the invariants ρ 1 andρ 2 are defined as ξ in Eq. (6) is the Kobayashi-Maskawa-'t Hooft determinant which reads and the relevant term breaks the U A (1) symmetry, which stems from the axial anomaly due to quantum fluctuations of QCD. The coefficient c A is scale independent. The last two terms linear in the sigma fields in Eq. (6) break the chiral symmetry explicitly, which results directly in mass acquirement for the Goldstone bosons, such as pions and kaons, etc. In Eq. (6) we have employed the light-strange basis, which is related to the singlet-octet one through the relation as follows with φ denoting scalar and pseudoscalar mesons collectively. Apparently, the strength of explicit breaking of the chiral symmetry, in another word, how massive the Goldstone bosons are, is connected to the magnitude of the coefficients j L and j S in Eq. (6), which in this work are regarded as the RG scale independent parameters to be determined in the following. The constituent quark masses for the light and strange quarks are given by respectively, and the pion and kaon decay constants read [59] whereσ L andσ S denotes the expected values of σ L and σ S fields. The meson mass squares are obtained by diagonalizing the Hessian matrix of the effective potential, which reads The Hessian matrix is block diagonal in the scalar and pseudoscalar channels, and in each block the only nonvanishing nondiagonal element is H 80 or H 08 , which resulting in the mixing of mesons between the octet and singlet. We will not go into the details in this work, and for more relevant discussions as well as explicit expressions for the meson masses, see, e.g. [50].

III. QUANTUM FLUCTUATIONS WITHIN FRG
In the RG scale dependent effective action Γ k in Eq. (1), quantum fluctuations of wavelength 1/k are suppressed, and the scale k is in fact an infrared (IR) cutoff scale. Therefore, by lowering the IR cutoff scale toward the IR limit, i.e., k → 0, one arrives at a full quantum effective action Γ k→0 . As a consequence, one also has to specify the initial ultraviolet (UV) scale where the evolution begins. An ideal choice is a scale deep in the perturbative regime, for instance the mass scale of the Z boson, which, however, entails the embedding of the glue dynamics. The dynamics of glue part, including the gluon and ghost in the Landau gauge for example, is indispensable to the evolution of Γ k , when k is above ∼ 1 GeV, thus the effective action in Eq. (1) is insufficient at high scale. For more discussions about the first-principle FRG QCD calculations and recent progresses whereof, see e.g. [17,[26][27][28][29][30]60]. It has been known that, with the evolution of the RG scale k from the UV to IR regime, the glue dynamics decouples from other degrees of freedom when k is reduced to about 1 GeV, since a mass gap of the gluon field develops therein [27][28][29][30]. Therefore, one can start the evolution of flow equations at a scale Λ ∼ 1 GeV, safely neglecting the glue part, and the quantities at the initial scale are the parameters of the low energy effective model, which are needed to be determined.
The flow equation for the effective action in Eq. (1), i.e., the Wetterich equation [61], is given by with the RG time being t = ln(k/Λ), where Λ is the initial evolution scale as mentioned above, and it is also called as the UV cutoff. The two terms on the r.h.s. of Eq. (14) corresponds to contributions from the quarks and mesons, respectively, and in our case there are three flavor quarks, meson nonets in the scalar and pseudoscalar channels. In Eq. (14) G k 's are the propagators of quarks and mosons, and the IR regulators R k 's suppress quantum fluctuations of wavelengths 1/k. With the truncation of LPA, the Wetterich equation in Eq. (14) is equivalent to the flow equation for the effective potential, since all the k-dependence is encoded in such potential, which reads with the dimensionless mass squarem 2 i,k ≡ m 2 i,k /k 2 and the threshold functions l where all the scalar and pseudoscalar mesons in the nonets are shown explicitly, and the numbers in front of l (B) 0 's denote their respective degeneracy of isospin. The dependence of the threshold functions on the temperature and chemical potential is also shown in Eq. (15). Note that the baryon chemical potential does not enter into l (B) 0 's, since the mesons do not carry the baryon number. However, this is not the case for other chemical potentials, such as those for the electric charge and strangeness, which can be carried by a meson, for more relevant discussions, see, e.g. [51]. In this work we employ the Taylor expansion around the physical point to solve the flow equation for the effective potential in Eq. (15), which is presented in detail in Appendix A.

IV. THERMODYNAMICS AND THE BARYON NUMBER FLUCTUATIONS
The thermodynamical potential density is connected to the effective action Γ k in Eq. (1), more exactly its IR limit with k = 0, through the relation as follows whereΦ are the physical values of fields, i.e., the solutions of their equations of motion, and V is the volume of the system. Upon inserting Eq. (1) into Eq. (16) and considering the expected values of all the fields, one arrives at where it is assumed that Ω has been normalized to be vanishing at vacuum as Eq. (16). In order to investigate the dependence of our calculated results on the parameterization scheme for the glue potential, in this work we employ two different glue potentials, which are commonly used in literatures, i.e., the polynomial potential in Eq. (B1) and that with the Haar measure in Eq. (B5), respectively. More detailed descriptions and discussions about these two potentials are deferred to Appendix B.
With the thermodynamical potential in Eq. (16) in hands, one can obtain other equilibrium thermodynamical quantities, such as the pressure and the entropy density: The equation of state for the QCD matter is well describe by the interaction measure, i.e., the trace anomaly which with the energy density . The n-th order cumulants of the net baryon number N B distributions are given by is the probability distribution of N B , and N B its mean value. Therefore, theoretical calculations of the cumulants are feasible, if the probability distribution P (N B ) is obtained, that is what has been done in our former work [62]. However, one can also chooses another equivalent approach, which is more commonly used, viz. computing the n-th order derivative of the pressure w.r.t. the baryon chemical po- tential, which reads They are related to the cumulants of the baryon number distribution through relations, for instance up to the fourth order, as follow In the experimental measurements, the mean value M , variance σ 2 , skewness S, and the kurtosis κ of the net proton or baryon number distribution are usually used, which are connected to the generalized susceptibilities in Eq. (22) by

V. NUMERICAL RESULTS
In this section we would like to give our numerical results, but before that, the initial conditions of the flow equation for the effective potential in Eq. (15) as well as other parameters, e.g. those in Eq. (6), have to be specified. They are fixed by fitting hadronic observables at vacuum, and relevant discussions are presented in Appendix C in detail, and we also discuss how to reduce the influence of the UV cutoff Λ on observables at large T or µ therein.
In Fig. 1 we show the meson and quark masses as functions of the RG scale k at vacuum, and as functions of the temperature. When we refer to results at finite temperature or chemical potential, the RG scale k = 0 is assumed, which is still applicable in the following. Comparing the left and right panels in Fig. 1, one finds that the dependent behaviors of the masses on k and the temperature are similar, which are reasonable and self-consistent. With the increase of k or T , the dynamically broken chiral symmetry is restored, and the chiral partners, such as π-σ, a 0 -η , K-κ, become degenerate with each other, and furthermore, the constituent quark masses decrease pronouncedly, especially the light quarks. Note, however, that if we compare the results obtained here in the low energy effective model with those in the QCD calculations, including the quantum fluctuations of the glue sector, e.g. in [27,63], the distinction is remarkable, and one would find that the mesons in the QCD calculations decouples much more quickly than the results here, once the scale or the temperature is located in the chiral symmetrical phase. And also the quark masses approach their bare masses more slowly in the effective model than those in the QCD. This feature of the effective model will affect the thermodynamics to be discussed in the following.
In Fig. 2 we show the pion and kaon decay constants as functions of T with µ B = 0. As shown in Eq. (12), f π and f K are linked to the expected values of the σ L and σ S fields, and therefore serve as order parameters for the QCD chiral phase transition, or more exactly the continuous crossover. Alternatively, one can also use ρ 1 defined in Eq. (7) as the order parameter, which is reduced to after expected values for all mesonic fields have been inserted, and note that only those of σ fields are nonvanishing. The pseudocritical temperature extracted from the peak of |∂ρ 1 /∂T | at µ B = 0 is T χ c = 194 MeV, where the superscript χ denotes the chiral crossover, which is intended to be distinguished from the pseudocritical temperature for the deconfinement phase transition T d c . It is found in our calculations that T d c = 177 MeV, where T d c is obtained from the peak of the derivative of the Polyakov loop w.r.t. the temperature. Note that T χ c in this work is larger than that of lattice QCD by HotQCD Collaboration 154 ± 9 MeV [6] and Wuppertal-Budapest Collaboration 156 ± 9 MeV [9,64]. We think this remarkable difference is attributed to several reasons as follow. Firstly, the absolute scale of the effective model is inherently different from that of QCD, and the former is larger. A larger critical temperature is also found in  [7,8] and by Wuppertal-Budapest Collaboration [10] are also presented for comparison.
another effective model calculation with some different setup [51]. Secondly, in the calculations we have used the glue potential V glue-Haar in Eq. (B5) with two parameters T glue c = 270 MeV and α = 0.52, and although the value of T c is allowed, but a bit larger. Finally, the σ-meson mass is fixed to be m σ = 510 MeV in our calculations, as discussed in detail in Appendix C. This value is located in the experimental measured mass region of f 0 (500), i.e., 400 − 550 MeV [65], but almost touches the allowed maximal value. However, because of the numer- ical instability for the Taylor expansion of the effective potential around the physical point in Eq. (A10), see Appendix A for more relevant discussions, it is difficult to decrease m σ further. A smaller σ-meson mass will results in a smaller T χ c . We will try to overcome this numerical instability and present relevant results elsewhere.
Although there is an absolute scale difference between the effective model and the lattice simulation, it will not hamper the comparisons of results obtained from the two calculations, if the relative scale is used, such as the temperature in unit of the critical one, as will be discussed in what follows. We show the pressure and the trace anomaly in Fig. 3 and the kurtosis of the baryon number distribution, i.e., the ratio χ B 4 /χ B 2 , in Fig. 4. Two different glue potentials in Appendix B are used in the calculations. Our calculated results are also compared with relevant lattice results by HotQCD Collaboration [6][7][8] and Wuppertal-Budapest Collaboration [9,10]. In order to eliminate the difference of absolute scale between the lattice QCD and the effective model as discussed above, we rescale the temperature by their respective pseudocritical temperature. The pseudocritical temperature T c for the lattice calculations is chosen to be central value, i.e., 154 MeV from 154±9 MeV by HotQCD Collaboration [6], and 156 MeV from 156±9 MeV by Wuppertal-Budapest Collaboration [9,64], and the errors are not taken into account in this work. Note that, in the low energy effective model, two pseudocritical temperatures T χ c and T d c are defined, which are related to the chiral and deconfinement phase transitions, respectively. In Fig. 3 we choose T c = T χ c for the effective model. One finds that our calculated pressure and trace anomaly agree with the lattice results. This is not a surprise, since two parameters in the glue potential are employed to fit the pressure and the trace anomaly, see Appendix B for more detailed discussions. Therefore, it is more valuable to compare our calculated baryon number fluctuations with the lattice simulations as shown in Fig. 4. However, there is still an intricacy needed to be fixed before the comparison. It has been known that the ratio χ B 4 /χ B 2 is linked to the degrees of freedom [20,54,66,67], and therefore this ratio is more sensitive to the deconfinement phase transition. Thus in Fig. 4 we relax T c of the effective model to be an appropriate value between T χ c and T d c , and we find T c = 185 MeV for the Haar potential which gives an optimal agreement. This is nontrivial, and it would be clearer if one looks at the red dashed curve calculated with the polynomial glue potential. The best choice for the polynomial potential is T c = T χ c . Nevertheless, the agreement with the lattice results for the polynomial potential is not as good as that for the Haar potential, specifically at high temperature.
In this work, we also perform calculations at finite baryon chemical potential. The phase diagram of the 2+1 flavor low energy effective model in the plane of T and µ B is shown in Fig. 5. We use the color to indicate the value of ∂ρ 1 (T, µ B )/∂T with the chiral order parameter ρ 1 given in Eq. (24). The two solid lines are determined by 90% peak height of this derivative at each value of µ B . This two lines approach toward each other with the increase of the baryon chemical potential, and crosses at the critical end point. The location of the CEP is found to be (T CEP = 84 MeV , µ B CEP = 840 MeV). The method of Taylor expansion of the effective potential results in numerical instability around the first-order phase transition, so we do not show results at very high baryon chemical potential. The ideal approach to reveal the global properties of the effective potential is to lattice the potential in the field space, which therefore can be employed to study the first-order phase transition. Interested readers are referred to, e.g., [68].
In Fig. 6 we investigate the dependence of the pressure and trace anomaly on the baryon chemical potential. One finds that the pressure scaled by T 4 as a function of the temperature moves up globally with the increase of µ B , and the height of the trace anomaly increases as well. Note that our calculated results are obtained with a glue potential without direct dependence on the baryon chemical potential. When the influence of the chemical potential on the glue potential is taken into account, for instance, through the µ-modification of the parameter T glue c in the glue potential [51], It is reasonable to expect that the behavior of the trace anomaly would change considerably, cf. Fig. 8. Unfortunately, a conclusion about the dependence of the glue potential on the chemical potential has not yet been arrived at, but it might be inferred through a detailed comparison of the trace anomaly at finite µ B between the low energy effective theory and the lattice simulations, which is though very interesting but beyond the scope of this work, and we will report it elsewhere in the future.
In Fig. 7 we show as functions of the temperature at several values of the baryon chemical potential. χ B 1 /χ B 2 and χ B 3 /χ B 2 are linearly dependent on µ B when µ B is not far from zero, so they are sensitive to the chemical potential in this region, which is also verified in our calculations. Therefore, is usually employed to extract the freeze-out chemical potential [8,10]. Furthermore, one finds that with the increase of the baryon chemical potential, the kurtosis of the baryon number distribution, i.e., the ratio χ B 4 /χ B 2 develops a minus value as well as a peak, which is more significant when the critical end point is approached. But the region of minus χ B 4 /χ B 2 shrinks and vanishes at the CEP [47].

VI. SUMMARY AND OUTLOOK
We have studied the QCD phase transition, thermodynamics, and the baryon number fluctuations in the 2+1 flavor low energy effective theory. In the calculations quantum fluctuations are included through the functional renormalization group approach. The flow equation for the effective potential is solved by Taylor-expanding it around the physical point.
The QCD phase transition, as well as the behaviors of the masses of mesons and quarks, the pion and kaon decay constants during the phase transition, has been investigated in the 2+1 flavor effective model at finite temperature and chemical potential.
The equation of state for the QCD matter, including the pressure and the trace anomaly, and the kurtosis of the baryon number distribution are calculated and compared with lattice results. We find the agreement between the low energy effective model and the lattice QCD is acceptable, especially for the calculations with the glue potential which takes the Polyakov loop fluctuations into account. We also calculate the EoS and baryon number fluctuations up to the four order at finite baryon chemical potential.
It should be noted that this work is our first calculation of high order cumulants of the baryon number distribution in the 2+1 flavor low energy effective model within the FRG approach. Lots of things have to be done in the future, for instance, going beyond the LPA truncation employed in this work; including the degree of freedom of gluons and extending the low energy effective theory to the 2+1 flavor rebosonized QCD. Furthermore, in our calculations only the baryon chemical potential are included, two other chemical potentials, i.e., the electric charge and strangeness chemical potentials, are also needed to be taken into account, when we calculate their fluctuations and correlations, or the strange neutrality [51] and the ratio Z/A are demanded, where Z and A are electric and mass number for a nucleus. Therefore, more sophisticated calculations and comparison with lattice QCD and experimental data are especially desirable.

ACKNOWLEDGMENTS
We thank Jan M. Pawlowski and Fabian Rennecke for valuable discussions, and Heng-Tong Ding for providing us with lattice data. We are grateful to Bernd-Jochen Schaefer and Heng-Tong Ding for reading this manuscript and giving us a number of valuable comments and suggestions. We also thank the members of the fQCD collaboration [69] for work on related projects. The work was supported by the National Natural Science Foundation of China under Contracts Nos. 11775041. The flow equation of the effective potential in Eq. (15) is obtained with the flat or Litim 3d IR regulators [70,71], as follow for the quark and meson fields, respectively, where the shape functions read with the Heaviside step function Θ(x). The threshold functions in Eq. (15) are given by where the bosonic distribution function reads and the Polyakov-loop modified fermionic distribution function n F (m 2 ; T, µ, L,L) = 1 + 2L e x/T + L e 2x/T 1 + 3L e x/T + 3L e 2x/T + e 3x/T , with In this work, we employ the method of Taylor expansion to solve the flow equation for the effective potential in Eq. (15), which reads with N ≥ m + 2n being the maximal order of the expansion, and N = 5 is adopted in this work, which is found to suffice for the convergence of calculations. κ 1,k and κ 2,k are the expansion points for ρ 1 andρ 2 , respectively. The expansion points usually can be chosen with a freedom given being in the convergence regime, in which the physical results are independent of the choice of the expansion points. In literatures there are two commonly employed choices. One is the fixed point expansion with the expansion points, i.e., κ 1,k and κ 2,k here, independent of the RG scale k [35,50,51]. The prominent advantage of the fixed point expansion lies in its excellent numerical stability. The other is the physical point expansion, see, e.g. [62], in which κ 1,k and κ 2,k are the physical points for every value of k, and therefore they are k dependent. As the convergence is concerned, the physical point expansion is superior to the fixed point expansion, while the former loses the excellent numerical stability. In this work we employ the approach of physical point expansion. Differentiating both sides of Eq. (A10) w.r.t. the RG time t, one arrives at which yields Specially, for the expansion coefficients λ 10,k and λ 01,k , one has In the meantime, in order to find the evolution equations for the expansion points, the stationarity condition is implemented as follows which leads us to the following equations which read and It follows from Eq. (7) and Eq. (8) that one arrives at which result in Equations (A13), (A14), (A16), (A17), (A20), (A21) constitute a closed set of linear equations, which can be solved straightforwardly to obtain the flows equations for κ 1,k , κ 2,k , λ 10,k , and λ 01,k . Since the final expressions for these flows are lengthy and the set of linear equations can be quite easily solved, we will not present them here. In this work two different parameterizations for the glue potential, i.e., the Polyakov-loop potential V glue (L,L) in Eq. (1), are employed. One is the polynomial potential, which is widely used in literatures, see e.g. [72], and given bȳ withV glue-poly = V glue /T 4 , and the temperature dependence is only encoded in the coefficient b 2 which reads where t ≡ (T − T c )/T c is the reduced temperature. Note that the Polyakov-loop potential is parameterized by employing the thermodynamics and the behaviors of Polyakov loop in the Yang-Mills (YM) theory, thus T c is the critical temperature for the deconfinement phase transition in the YM theory. Interestingly, it is found that the YM Polyakov-loop potential is also applicable in QCD [73], and the unquenching effect is quantitatively taken into account, given the reduced temperature is appropriately rescaled, for instance with where α is the linear scale factor and T glue c is the pseudocritical temperature for the deconfinement phase transition in QCD. It is found that α 0.57 in the case of flavor N f = 2 [73]. T glue c can also be estimated by the dependence of Λ QCD on N f via the QCD RG running [45], and T glue c = 208 MeV for N f = 2 and 187 MeV for N f = 2 + 1. In this work we will relax the restriction on the values of parameters α and T glue c , and study the dependence of the QCD thermodynamics on them systematically.
In recent year, another parameterization for the Polyakov-loop potential is proposed [56], in which fluctuations of the Polyakov loop are also taken into account besides usually employed quantities. Furthermore, this parameterization employs the SU (N c ) Haar measure, which solves the problem, as the polynomial potential in Eq. (B1) has, that the Polyakov loop exceeds unity at high temperature. We denote this potential as V glue-Haar which reads where coefficients in Eq. (B5) are hatted with a bar in order to be distinguished from those in Eq. (B1). Different from the polynomial potential, all the coefficients in Eq. (B5) are dependent on the temperature, and their dependence is given by for x ∈ {ā,c,d}, and Values of all these coefficients as well as those in the Eqs. (B1) and (B2) are collected in Tab. I. In Fig. 8 and Fig. 9 we have investigated the dependence of the thermodynamics on the two parameters T glue c and α in the glue potential, respectively. One finds that T glue c shifts the curves of the dimensionless pressure and trace anomaly along the T direction, and affects the height of the trace anomaly as well, while variation of α only result in the change of height of the trace anomaly. Therefore, we would like to employ these properties to

Appendix C: Numerical setup
In an ideal case, we wish to start the evolution of the FRG flows at an initial scale Λ well above the scale of interests, such as ∼ 2πT related to the Matsubara gap of the temperature. In a low energy effective theory this, however, is restricted by the lack of the glue dynamics as we discuss in Sec. III, since the gluon mass gap disappears rapidly when the RG scale is above ∼ 1 GeV [26][27][28][29]. Therefore, the value of Λ is limited to a small window around ∼ 1 GeV. In this work without loss of the generality, we choose Λ = 1 GeV.
At the initial scale Λ, the effective potential in Eq. (15) is classical, and it is reasonable to assume that its irrele- vant terms, i.e. those with couplings of dimension minus, in the Taylor expansion in Eq. (A10) are vanishing, which leads us to U k=Λ (ρ 1 ,ρ 2 ) = λ 10,Λ ρ 1 + λ 20,Λ 2 ρ 2 1 + λ 01,Λρ2 , where the three coefficients together with those in Eq. Note that the σ-meson mass is chosen to be m σ = 510 MeV in this work; however, in the experiments the σmeson related scalar meson f 0 (500) with I G (J P C ) = 0 + (0 ++ ) is located in a broad mass region of 400 − 550 MeV [65]. Therefore, it is necessary to investigate the influence of the σ-meson mass on the QCD thermodynamics, and the relevant results are presented in Fig. 10 through Fig. 12 for the EoS, baryon number fluctuations, and the kurtosis of the baryon number distribution, respectively. As we discuss in Sec. V, due to the numerical instability for the Taylor expansion of the effective potential around the physical point, we can not decrease the value of m σ to be below 500 MeV. With the variation of m σ in the interval 510 − 600 MeV, we only find a mild dependence on the σ-meson mass.
As we have mentioned above, the initial scale is chosen to be Λ = 1 GeV. Making a crude estimate from the Matsubara gap as follows where T Λ is the temperature, above which the UV cutoff effect becomes significant and can not be neglected. It follows from Eq. (C2) that T Λ ∼ 160 MeV, which is around the critical temperature of QCD phase transitions. Therefore, we have to reduce the influence of the UV cutoff on observables of interest, especially at high temperature. One way to diminish the UV cutoff effect is to introduce an appropriate modification for the effective action at the initial scale, i.e., where the effects of external parameters, such as the temperature and chemical potential, can be taken into account to some degree at the second term on the r.h.s. of Eq. (C3). For instance, one can integrate the flow of the effective action from infinity down to Λ, while with contributions of the vacuum subtracted, to wit, In the low energy effective model, mesons decouple and are irrelevant in the scale above the UV cutoff, thus ∆Γ k=Λ in Eq. (C4) receives contributions only from free quarks. Therefore, it follows from Eq. (C4) that with the Polyakov-loop modified fermionic distribution function n F given in Eq. (A8), andm 2 l,k = m 2 l,k /k 2 and m 2 s,k = m 2 s,k /k 2 . When the RG scale k is above Λ, the quark masses are approximated as the values of those at k = Λ in this work, viz. m l,k>Λ = m l,k=Λ and m s,k>Λ = m s,k=Λ , for values of m l,k=Λ and m s,k=Λ , see Fig. 1. Recently, Braun et al has proposed the concept of RG consistency, which is employed to analyze the cutoff effects, and is aimed to enhance the predictive power of low energy effective theories [74].