Baryon number probability distribution at finite temperature

The probability distribution of the net baryon number is investigated within the functional renormalization group approach. We find that the Roberge-Weiss periodicity related to the Z(3) symmetry of the gluon fields results directly in that, states of the net baryon number $N_B=N\pm 1/3$ with $N\in\mathbb{Z}$ are prohibited, and only those of $N_B=N$ are possible. By employing the probability distribution of the net baryon number, we also compute the cumulants of the baryon number distribution, which are found to be well consistent with those obtained from the generalized susceptibilities. A question about the relation between the color confinement and the probability distribution of net baryon number is put forward.


I. INTRODUCTION
Studies of the QCD phase structure have attracted lots of attentions in recent years. The Beam Energy Scan (BES) program at Relativistic Heavy Ion Collider (RHIC) is aiming at locating the critical end point (CEP) of the QCD [1], which separates the first-order phase transition at high densities from the continuous crossover at high temperature in the phase diagram [2]. Interesting and promising results have been arrived at in the experiments, in particular a nonmonotonic behavior of the kurtosis of the net proton number distribution with the variation of the collision energy has been observed [3,4].
In this work we will investigate the probability distribution of the net baryon number, which is intimately relevant to the fluctuations of the net proton number, and therefore plays a significant role in the experiments of searching for the CEP of QCD [1]. We will study the influence of the glue dynamics on the probability distribution, and an interesting interrelation between the Roberge-Weiss periodicity [34] due to the Z(3) symmetry of the gluon fields and the resulting probability distribution will be revealed. Furthermore, different from the * wjfu@dlut.edu.cn cumulants of the baryon number distribution, the probability distribution of the baryon number calculated in this work, can be directly used as input for some transport simulations in heavy-ion collisions, see e.g., [35,36], through which one can assess whether the critical behavior related to the chiral symmetry can survive through the subsequent evolution due to some noncritical effects, for instance, volume corrections, detector acceptance cut, etc. [37].
This paper is organized as follows. In Sec. II we give a brief introduction about the formalism of the probability distribution. The low energy effective theory and the FRG approach are presented in Sec. III. In Sec. IV the cumulants of the net baryon number distribution are discussed. We present our numerical results and discussions in Sec. V. Then a summary and conclusion is given in Sec. VI.

II. THE BARYON NUMBER PROBABILITY DISTRIBUTION
We follow the formalism presented in, e.g., Refs. [26,27], and consider an thermodynamical system with a volume V , temperature T and a baryon chemical potential µ B , which is homogeneous and well equilibrated, the net baryon number N B probability distribution reads where Z(T, V, N B ) and Z(T, V, µ B ) are the canonical and grand canonical partition functions, respectively. Note that, in contradistinction to the ordinary case where net fermion number should be integer, the net baryon number N B here could be multiples of the baryon number arXiv:1805.12025v1 [hep-ph] 30 May 2018 of a quark, i.e., N B = N q /3 with N q ∈ Z being the net quark number. This is necessary, because if one would like to describe the baryon number probability distribution during the QCD phase transition within one theoretical framework, where the degrees of freedom transform, both quarks and hadrons should be taken into account. Therefore, one arrives at the normalization for the probability distribution as follows Immediately, upon inserting Eq. (1), it follows that with the fugacity λ = exp(µ q /T ) and the quark chemical potential µ q = µ B /3. Obviously, the canonical partition function is just the expansion coefficient of the Laurent series of the grand canonical partition function as a function of fugacity, which leads us to obtain with C indicating a closed contour with the origin point in it in the complex λ plane. Choosing the unit circle for the contour C, i.e., λ = e iθ , we are led to Therefore, the canonical partition function can be deduced from the grand canonical one Z(T, V, µ q ) with an imaginary chemical potential, i.e, µ q = iθT . Thus, the task of calculating the net baryon number probability distribution in Eq. (1) is converted to compute the thermodynamic potential density Ω with an imaginary chemical potential, which is related to the grand canonical partition function as follows

III. THE LOW ENERGY EFFECTIVE THEORY WITHIN THE FRG APPROACH
As we have shown in the section above, the thermodynamical potential for a system with an imaginary chemical potential is indispensable to the investigation of the net baryon number probability distribution in Eq. (1). In this work, we will focus on two ingredients which affect the baryon number distribution: one is the chiral critical behavior, which is relevant to the chiral phase transition; and the other is the confinement information, that is encoded in the glue dynamics and related to the deconfinement phase transition. It should be emphasized that the glue dynamics plays a significant role in the determination of the probability distribution, as we will show in what follows. This is consistent with the fact that the kurtosis of the baryon number distribution is pronouncedly influenced by the glue dynamics as well, see e.g. [29,33]. In this work we employ the two-flavor Polyakov-quark-meson (PQM) low energy effective model within the FRG approach, for more details see, e.g., [29][30][31]. As a nonperturbative continuum field approach, FRG encodes quantum fluctuations of different wavelengths successively, through the running of RG scale k from the UV to IR regime, and the scaledependent effective action for the low energy effective model is given by are the mesonic fields, and the effective potential Z q,k and Z φ,k are the anomalous dimensions for the quark and meson, respectively; h k is the Yukawa coupling; the term −cσ in Eq. (7) breaks the chiral symmetry explicitly. Moreover, in order to take into account the deconfinement information, we include the temporal component of the gluon background field A 0 in the effective action. It can also be transformed into the formalism of the traced Polyakov loop L and its conjugateL, i.e., with where P on the r.h.s. is the path-ordering. V glue in the effective action in Eq. (7) is the glue potential which will be specified in the following. The evolution equation for the scale-dependent effective action in Eq. (7), i.e., the Wetterich equation [38], reads with t = ln(k/Λ) and Λ is the initial evolution scale, i.e., the UV cutoff scale. G k 's are the propagators for the quark and meson fields; R k 's are the IR regulators that suppress specific quantum fluctuations whose wavelengths are larger than 1/k. In this work, the local potential approximation (LPA), i.e., Z q,k = Z φ,k = 1 and ∂ t h k = 0, is adopted. Then following from the Wetterich equation, one immediately obtains the flow equation for the effective potential in Eq. (7), i.e., where the dimensionless masses for the meson and quark readm The bosonic distribution function is given by When the chemical potential is imaginary, such as µ q = iθT as discussed in Sec. II, the Polyakov loops L and L in Eq. (8) are complex conjugate to each other exactly. Note, however, that a difficulty will arise when the chemical potential is real valued, for example, when we calculate the grand canonical partition function in the denominator in Eq. (1) with a real-valued µ B . Inserting µ q = iθT , L = |L|e iφ ,L = |L|e −iφ , one obtains the fermionic distribution functions as follow andn with E q,k = k 1 +m 2 q,k . The thermodynamic potential density is given by where the expectation value of the meson field is determined through its equation of motion (EoM), and the normalization Ω = 0 at vacuum is assumed. Since the glue potential V glue has Z(3) symmetry, the symmetry of the thermodynamic potential in Eq. (16) with an imaginary chemical potential has been well known, see e.g., [34,58]. We introduce modified Polyakov loops as follow with φ = φ+θ. Substituting the modified Polyakov loops above into Eq. (16), one immediately recognizes that Ω is invariant with the replacement θ → θ + 2π 3 , which is also known as the Roberge-Weiss periodicity. The periodicity of the thermodynamic potential leads to as well. Furthermore, it has been found that L (θ) is not continuous at θ = (2n + 1)π/3 with n ∈ Z, when the temperature is above some critical value, which is also called as the Roberge-Weiss phase transition [34]. the Roberge-Weiss periodicity of the thermodynamic potential, leads directly to the same periodicity of the grand canonical partition function with an imaginary chemical potential µ q = iθT in Eq. (6). Then we can rewrite Eq. (5) as follows where the expression in the parentheses in the integrand has an interesting property, i.e., it is nonvanishing, only when the remainder of N q over 3 is zero. Therefore, the probability to find a system in a state with non-integer baryon numbers is zero, which is apparently a manifestation of the color confinement.

IV. CUMULANTS OF THE BARYON NUMBER DISTRIBUTION
With the probability distribution in Eq. (1) in hand, one can immediately obtain the statistical average for a quantity, denoted here symbolically as O(N B ), through the equation as follows and the n-th order cumulant of the baryon number dis- In fact, the cumulants can also be calculated in another commonly employed method, which resorts to derivatives of the thermodynamic potential with respect to the real-valued baryon chemical potential, i.e., the generalized susceptibilities, see e.g., [29]. The relevant definition reads The relations between χ n and the cumulants are given by taking the quadratic and quartic orders for instance.
In the following, we will compare results of the cumulants calculated in these two different approaches. The comparison is nontrivial, especially when the glue dynamics and, thus the information of confinement, is encoded in the calculation. This is because when the baryon chemical potential is real valued, the Polyakov loop L and its conjugateL are ill-defined, since they are not complex conjugate to each other. One method commonly employed to bypass this difficulty is to ignore the phase of L andL, and treat them as two different real quantities, which are determined through their respective equations of motion.

V. NUMERICAL RESULTS
We employ the Taylor expansion around the physical point to solve the flow equation for the effective potential in Eq. (11), which reads Here κ k is the solution of its EoM for every RG scale k, which fulfills Inserting Eq. (23) into Eq. (11) and constraining the expansion point with Eq. (24), one obtains the flow equation for κ k , i.e., and those for the expansion coefficients λ n,k 's, which reads and where ∂ n ρ denotes the n-th order partial derivative with respect to ρ, and ρ the differentiation with ρ fixed. In our calculations, the maximal order of the Taylor expansion N in Eq. (23) is chosen to be N = 5, which has guaranteed the required convergence very well, for more discussions about the convergence, see e.g. [23] .
We evolve these flow equations from an initial UV scale, which is chosen to be Λ = 700 MeV in this work. The effective potential at the initial scale is chiral symmetric, except the explicit breaking induced by the −cσ term. Hence, the effective potential at k = Λ can be approximated as Parameters in the effective potential above, to wit, λ Λ and ν Λ , as well as the Yukawa coupling h and the explicit chiral symmetry breaking related coefficient c, are fixed by fitting hadronic observables: the π-and σmesons masses m π = 135.9 MeV and m σ = 502.2 MeV, the π-meson decay constant f π = 92.0 MeV, and the constituent quark mass m q = 299.1 MeV. Their values are given by λ Λ = 5.7, ν Λ = 0.23 GeV 2 , h = 6.5, c = 1.7 × 10 −3 GeV 3 , respectively.
It is left to specify the glue potential in Eqs. (7) and (16). In this work, we adopt the glue potential parameterized firstly by Lo et al. in Ref. [59], which has also been used in Ref. [31]. This potential is characteristic of its ability to describe the quadratic fluctuations of the Polyakov loop, which reads with where M H is the Haar measure of the SU (3) group in color space. The parameterization of the coefficient a(T ) is given by  and similar parameterizations also apply to c(T ) and d(T ), but for b(T ) it reads All the constants required for this parameterized glue potential can be found in [31,59]. We collect them once more in Tab. I for the convenience and completeness. t YM in Eqs. (32) and (33) are the reduced temperature for the Yang-Mills gauge theory at finite temperature. It has been found that the QCD glue potential, which encodes the backreaction of the matter on the glue sector, can be remarkably well parameterized through the Yang-Mills one [20,22,44], just by making the following replacement: where the scaling constant α = 0.57 is found for the two-flavor QCD, and we adopt the critical temperature T glue c = 250 MeV for the glue potential. In Fig. 1 we show the dependence of the modified Polyakov loop L in Eq. (17) on the θ with µ q = iθT , which is calculated in the PQM within the FRG approach at several values of the temperature. The magnitude and phase of L are presented in different panels. One can see that the Roberge-Weiss periodicity is quite pronounced in all panels of Fig. 1. We also investigate the dependence of the temperature, and confirm that when the temperature is above a critical value, that is found to be T RW c = 208 MeV in our calculations, a discontinuity appears in the curve of sin(φ ), as the solid line shows in the bottom panel of Fig. 1. A recent lattice simulation by the Wuppertal-Budapest Group also find a similar value of T RW c [14]. Figure 2 shows our calculated probability distributions of the net baryon number, whose formula is given in Eq. (1). We have discussed a lot about the canonical partition function in the numerator of Eq. (1) in former sections. The grand canonical partition function in the denominator can be easily obtained from its thermodynamic potential with a real-valued baryon chemical potential, which has been studied heavily in literatures, and we refer interested readers to such as Refs. [29,30] etc. In Fig. 2 we only consider the case of vanishing baryon chemical potential. Therefore, the probability distributions are symmetric with respect to the zero point in N B . The four subplots in Fig. 2 correspond to different values of temperature, and for every temperature, we choose two volumes. One can observe that with the increase of the volume or temperature, the profile of distribution becomes wider or "fatter", which can be understood reasonably. We find that when the chemical potential is vanishing, the pseudo-critical temperature for the deconfinement phase transition is T Poly c = 178 MeV, and that for the chiral phase transition T χ c = 194 MeV in our calculations, which are identified by the peaks of ∂L/∂T and ∂ρ/∂T , respectively. When the temperature is low, such as T = 110 MeV as shown in the top-left panel of Fig. 2, on one hand, the chiral symmetry is broken dynamically and the constituent quark mass is large, on the other hand, quarks are confined inside baryons. Consequently, the possible net baryon number N B is restricted to a narrow region as the results of T = 110 MeV show. The situation, however, is drastically changed with the increase of the temperature, especially when the temperature is above the T Poly c and T χ c . One can see that it is much easier to excite state with large |N B | in the plots of T = 200 MeV as well as T = 250 MeV.
In Fig. 3 we compare the probability distribution calculated in the PQM and that in the quark-meson (QM) model, in order to illustrate the significance of the glue dynamics in the studies of baryon number distributions. As we discussed above, the Roberge-Weiss periodicity encoded in the glue dynamics entails that only integer N B is possible, as the blue circles show. On the contrary, loss of the periodicity, which is the case in the QM model, leads to the loss of the restriction as well. In another word, states of N B = N ± 1/3 with N ∈ Z can be excited as equally as those of N B = N , as the solid dots show in Fig. 3.
By employing the probability distribution of the net baryon number obtained in this work, we calculate the quadratic and quartic cumulants of the baryon number through Eq. (20), and relevant results are presented in  n 's, we label results relevant to the temperature below and above T RW c with circles and squares, respectively. Furthermore, we also show the result of χ B n , calculated directly from the norder derivative of the thermodynamic potential with respect to the real-valued baryon chemical potential, i.e., the generalized susceptibilities in Eq. (21), which is presented in Fig. 4 with the solid line. One can see that the two method agree with each other remarkably well. An more interesting observable is the kurtosis of the baryon number distribution, given by R B 42 = χ B 4 /χ B 2 , since it is closely related to the degree of freedom in a system. In Fig. 5 the kurtosis R B 42 computed in these two methods are presented as well. As one can see, R B 42 approaches to 1 in the low temperature regime, which indicates that it is the hadronic degree of freedom in the hadronic phase, which is in sharp contrast to the quark dominated system, where R B 42 would have been 1/9 in the low temperature limit, for more detailed discussions, see e.g., [29]. Once more, we find the calculated R B 42 in these two different approaches agrees with each other very well.
Note that the agreement between these two methods is nontrivial. This is because, in the approach of the generalized susceptibilities, the thermodynamic potential with a real-valued baryon chemical potential is required. The Polyakov loop, however, is ill-defined in the SU(3) gauge theory when µ B is real, due to the notorious sign problem. Therefore, in the actual calculations the Polyakov loop L and its conjugateL are treated as two independents quantities, and their phases are ignored. Thus, one can even say that the generalized susceptibility is an approximate approach. On the contrary, when we employ the probability distribution to calculate the cumulants of the net baryon number, only an imaginary chemical potential is needed, especially when the probability distribution is symmetric with respect to N B = 0. Thus, the calculations in the approach of the probability distribution are exact. Another noteworthy phenomenology in Fig. 4 and Fig. 5 is that both the circles and squares coincide with the solid line very well, which indicates that the Roberge-Weiss phase transition, i.e., the discontinuity observed in the Polyakov loop when the temperature is above T RW c , does not affect the observables, e.g., χ B 2 , χ B 4 , R B 42 etc.

VI. SUMMARY AND DISCUSSIONS
In this work we investigate the probability distribution of the net baryon number in the low energy effective model within the FRG approach. Emphases are put on the influence of the glue dynamics, in particular the Roberge-Weiss periodicity, on the probability distribution of the baryon number. We find that the Roberge-Weiss periodicity directly results in that states of N B = N ± 1/3 with N ∈ Z are prohibited, and only those of N B = N are possible, which is an indication of the color confinement from another viewpoint.
By employing the probability distribution of the net baryon number obtained in our calculations, we compute the quadratic and quartic fluctuations of the baryon number, and the kurtosis of the baryon number distribution. The obtained results are compared with those from the derivatives of the thermodynamic potential with respect to the real-valued baryon chemical potential, i.e., the generalized susceptibilities. We find that these two different approaches yield consistent results.
The probability distribution of the net baryon number obtained in this work can be used as input for some transport simulations in heavy-ion collision experiments, see e.g., [36], so the critical behavior near the chiral phase transition can be combined with noncritical effects, such as the volume corrections, detector acceptance cut, resonance decays, etc., see e.g. [37] for more discussions. All these effects are needed to be identified, before the signal of the QCD critical end point is pinned down at the BES program of RHIC [1]. Related work is in progress.
There is a question or problem raised in our studies. As we have discussed above, the probability distribution of the net baryon number is only possible when N B = N with N ∈ Z, because of the Roberge-Weiss periodicity. It is reasonable in the hadronic phase at low temperature. The periodicity, however, is still there in the high temperature regime as shown in Fig. 1; hence states of N B = N ± 1/3 are still prohibited at high temperature, which seems that it is not consistent with the picture of deconfined quarks. Another more amazing finding is that, even with the hadronic probability distribution of the net baryon number in the high temperature limit, the predicted cumulants of the baryon number distributions in Fig. 4 and Fig. 5 are consistent with those obtained from the approach of generalized susceptibilities, which is exotic. More studies are needed to answer these questions in the future.