Chemical Freeze-out Parameters via a Non-perturbative QCD Approach

By analyzing the calculated baryon number susceptibility ratios ${\chi_{1}^{B}}/{\chi_{2}^{B}}$ and ${\chi_{3}^{B}}/{\chi_{1}^{B}}$ in two-flavor system via the Dyson-Schwinger equation approach of QCD, we determine the chemical freeze-out temperature and baryon chemical potential in cases of both thermodynamic limit and finite size. We calculate the center-of-mass energy dependence of the ${\chi_{4}^{B}}/{\chi_{2}^{B}}\, (\kappa \sigma^{2})$ at the freeze-out line and find an excellent agreement with experimental data in $\sqrt{S_{NN}^{}} \geq 19.6\,$GeV region when taking into account the finite size effect. Our calculations indicate that the $\kappa \sigma^{2}$ exhibits a non-monotonic behavior in lower collision energy region.


I. INTRODUCTION
Phase transitions of strong interaction matter have been explored for more than forty years since the research may reveal the nature of the early universe matter evolution [1][2][3]. The transitions include chiral phase transition (from dynamical chiral symmetry to dynamical chiral symmetry breaking) which generates more than 98% of the mass of visible matter and the confinement transition (hadronization) which slaves the quarks and gluons to hadrons. They are driven by the temperature (T ) and the baryon density (ρ B ) or chemical potential (µ B ). Since the strong interaction can be well described by Quantum Chromodynamics (QCD), the above mentioned phase transitions are usually referred to as QCD phase transitions. Moreover, many calculations (see, e.g., Refs. [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]) have shown that the chiral phase transition at low chemical potential is a crossover at physical quark mass. Theoretical calculations (see, e.g., Refs. [3,6,[8][9][10][11][12][13][14][15][16][17][18][19]) also indicate that the chiral phase transition at high chemical potential is first order. Therefore, there would exist a critical end-point (CEP) in the T -µ B plane at which the first order phase transition turns to crossover. The position of the CEP or even its existence becomes thus one of the most significant topic in both theories and experiments. Besides the efforts in theories, the Beam Energy Scan (BES) program at RHIC, the FAIR at GSI and the NICA at DUBNA all take the search of the CEP as their investigation focus (see, e.g., Refs. [20][21][22]) and some meaningful information has been provided by the RHIC experiments [23][24][25][26].
In experiments, one can measure only the states after the hadronization but not the phase transition directly, and thus the chemical freeze-out line which is defined as the set of states ceasing the inelastic collision of the * qwertylou@pku.edu.cn † Present address: Institute of Theoretical Physics, Chinese Academy of Science, Beijing 100081, China ‡ Corresponding author: yxliu@pku.edu.cn newly formed hadrons plays the essential role. Especially, as the chemical freeze-out line approaches to the CEP, nonmonotonic behavior of conserved charge fluctuations could be observed [9,[27][28][29][30][31][32]. The freeze-out temperature and chemical potential have then been studied in statistical hadronization model (SHM) [33][34][35][36][37], hadron resonance gas (HRG) model [38,39], lattice QCD simulations [40][41][42][43] and other models [44,45]. In fact, the matter system generated in relativistic heavy ion collision (RHIC) experiment has a finite size and cools in a finite time [30,31,46,47]. The finite size and finite time prevent the correlation length ξ from diverging near the CEP, and smoothen the fluctuations [46]. Model calculations have shown that the finite size influences both the phase diagram and the thermodynamical properties drastically [48][49][50][51][52][53], the surface of the system may also play the role [54][55][56][57][58][59]. The effects of the finite size and the surface on the chemical freeze-out parameters will then complement the information for searching the CEP in experiments. However, different models give contradictory results. It is therefore imperative to investigate the finite size and the surface effects on the chemical freeze-out parameters with sophisticated QCD approaches.
It has been known that Dyson-Schwinger equations (DSE), a nonperturbative approach of QCD [60][61][62][63][64][65][66][67][68], have been successful in describing QCD phase transitions (see, e.g., Refs. [8-13, 64, 69-74]) and hadron properties (For recent reviews, see Refs. [64,65]). We then, in this paper, take the DSE approach to investigate the chemical freeze-out parameters with the finite size and surface effects being taken into account. We calculate the baryon number susceptibilities in two light-flavor quark system. By comparing the obtained baryon number susceptibility ratios χ B 1 /χ B 2 and χ B 3 /χ B 1 with the experimental data of the net-proton distribution cumulant ratios C 1 /C 2 and C 3 /C 1 at different collision energies, we determine the freeze-out parameters. We observe that with the finite size and surface effects being included, the calculated collision energy dependence of the χ B 4 /χ B 2 agrees with the experimental data excellently, and the calculated κσ 2 shows a nonmonotonic behavior in lower collision energy region. Moreover, we propose that hyper-order cumulant ratio such as χ B 6 /χ B 2 also shows a nonmonotonic dependence on the collision energy.
The remainder of this paper is organized as follows: In Sec. II, we describe briefly the Dyson-Schwinger equation approach and its relation to the chemical freeze-out parameters. In Sec. III, we calculate the freeze-out parameters by adopting the DSE approach as well as the experimental data. In Sec. IV, we give the phase diagram in the vicinity of CEP and reveal the effect of finite size and surface. In Sec. V, we give a summary and discussion.

A. Dyson-Schwinger Equation Approach
The Dyson-Schwinger equations are inifinite number of coupled equations. In this paper, we focus on the DSE for quark propagator S(ω j , p). The corresponding equation is: (1) where Z 1 , Z 2 and Z 4 are renormalization constants. m 0 is the current quark mass.ω j = ω j + iµ q , with µ q being the quark chemical potential, and ω j = (2j + 1)πT the Matsubara frequency for quarks. Σ(ω j , p) is the selfenergy of quark, and reads: where D µν is the dressed-gluon propagator, Γ ν is the dressed quark-gluon vertex, Ω jl = ω j − ω l is the Matsubara frequency for gluon. In principle, the gluon propagator D µν and the quarkgluon vertex Γ ν should be solved by corresponding DSEs, which depend on higher order correlation functions. And a truncation must be applied in order for numerical solution. In this paper, for the quark-gluon vertex, we adopt at first stage the rainbow approximation for the vertex Γ ν ( p,ω m , q,ω l ; T, µ q ) = γ ν .
The gluon propagator has the general form where P T,L µν are the transverse and longitudinal projection operators, respectively: where k = Ω nl , k . D T and D L are the effective interactions and can be represented using models. Note that the coupling constant g and the renormalization constant Z 1 has been absorbed into the effective interaction.
In this paper, we adopt the infrared constant model (QC model) [9,13,68], which reads: where D and ω are the parameters of the model. α pQCD is the ultraviolet perturbation term and reads: where The gluon screening mass m g is also considered in the longitudinal part of gluon model [13,75]: whose value is determined by leading-order perturbative QCD [76]: In order to fix the renormalization constant Z 2 and Z 4 , we need to specify the renormalization condition. In this paper, the renormalization condition is: where ζ is the renormalization point. We choose ζ = 19 GeV, m 0 = 3.4 MeV, D = 1.024 GeV 2 and ω = 0.5 GeV as in Ref. [68]. The quark propagator can be decomposed according to its Lorentz structure. At finite temperature, the decomposition is: There should be, in principle, a fourth term in this decomposition. However, its contribution to order parameter is extremely small and is usually omitted in practical calculations [61,77]. The mass function of the quark propagator can then be defined as: In vacuum, the physical solution to the quark DSE has a non-zero mass function even if the current quark mass is zero. At high temperature, the mass function gradually approaches the current quark mass. Therefore, the mass function at zero momentum, M (ω 2 0 , 0), is often used as the order parameter of the QCD crossover.

B. Number Density and susceptibilities
Experimental observations indicate that the yields of pion and proton are much larger than that of kaon [36], we can then simplify the matter generated in RHIC experiments as that including mainly two light flavor quarks. In the system of u and d quarks, baryon number density n B and electric charge density n Q can be fixed with quark number density n u,d as: From Eq. (12) we notice that the u and d quarks are in exact isospin symmetry, if only the baryon number is considered. In this sense, both the u quark and d quark hold the same quark chemical potential µ q = µ B /3, and the quark number density n q = 3n B . In view of statistical physics, the quark number density can be determined as where Z 2 is the quark wave-function renormalization constant, N c = 3 the color number, and N f = 2 the flavor number. Notice that the flavor number here represents the flavor degeneracy, and is different from the one we used in the ultraviolet perturbation term in the gluon model. In Eq. (14), the summation runs over an infinite number of Matsubara frequencies. In practice, we can only carry out the calculation with a finite number of Matsubara frequency, and the summation must have a cut-off, N . However, the convergence of summation Eq. (14) is slow, especially at low temperature. On the other hand, we observe that the scalar functions A, B and C defined in Eq. (10) converge quickly to the free-quark-propagator scalar functions A = C = 1, B = m 0 , as the Matsubara frequency grows large. For free propagator, the distribution function is: where E = p 2 + m 2 0 . Therefore, the contribution of missing Matsubara frequencies that exceed the cut-off N in Eq. (14), can be approximated using the distribution function of free prop-agator: Using this technique, the calculated quark number density has better convergence at moderate and high temperature, whose chemical potential dependence is shown in Fig. 1. After calculating the quark number density, we can calculate the susceptibility by taking derivatives. The kth order baryon number density susceptibility (fluctuation) is obtained as where β = 1/T and k = 2, 3, 4, · · · . The susceptibilities are related to the moments of the multiplicity distributions of the corresponding conserved charges as where M , σ 2 , S and κ are the mean, the variance, the skewness and the kurtosis of the multiplicity distribution, respectively. By comparing the theoretical net-baryon number fluctuations in terms of temperature and chemical potential with the experimental data one can determine the freeze-out parameters [43].

C. Finite Size and Surface Effect
The system created in RHIC exists in finite size, rather than thermodynamical limit. To determine the freezeout parameters in experiment one has to take the finite size and the surface effects into account. Assuming the system as a cube of size L, and adopting the anti-periodic condition, the momentum of a fermion should be p j = (2j + 1)π/L. In principle, one should sum over discrete momentum values. For simplicity, the finite size effect is roughly incorporated by a non-zero momentum cutoff |p| min = π/L [51,52]. It corresponds to an infrared momentum cut-off in Eqs. (1) and (13). This can also be understood by the quantum uncertainty principle, that L −1 is the energy scale for a system of finite size L. It is remarkable that such an L is not exactly the same as the size of fireball, but an effective scale that the ingredients of the quark matter can interact.
We also incorporate the effect of the surface through the multiple reflection expansion (MRE) approximation. In the MRE approximation, the thermodynamical quantities of a droplet composed of quarks can be derived from a density of states in the form [54-56, 78, 79] dN dp = 6 where V is the volume of the droplet, S = 4πL 2 and C = 8πL are the area, the extrinsic curvature of the surface of the droplet, respectively. The f s and f c are the contributions to the density of states from the surface and curvature, given explicitly as [50,78]: with p being the momentum and M the constituent quark mass. Rigorously, one should solve the coupled equations for the constituent mass M [50]. For simplicity we set M in Eq. (19) to be M = Re M (ω 2 0 , 0) from the mass function defined in Eq. (11), calculated from Eq. (1) with a finite system size L.
The modified density of states is then [50,79]: After taking into account the finite-size as well as the surface effect, the momentum integration in Sec. II A should be converted as follows:

III. FREEZE-OUT PARAMETERS
We have carried out calculations with L = ∞ (thermodynamical limit) and various finite values of L. The calculations manifest that the fluctuations (skewness, kurtosis, etc.) in the T -µ B plane behave qualitatively the same as those given in Ref. [9], respectively. The obtained µ B dependence of the baryon number susceptibility ratios χ B 1 /χ B 2 and χ B 3 /χ B 1 in case of L = 2.2 fm at several values of temperature are shown in Fig. 2. It is evident that our results agree with the lattice QCD results [43] qualitatively very well. In order to extract the freezeout parameters, we plot the experimental values of the cumulant ratios C 1 /C 2 = M/σ 2 and C 3 /C 1 = Sσ 3 /M of net-proton multiplicity distributions in central collisions [26] as horizontal lines.
By fitting our calculated χ .5, 7.7 GeV given in Ref. [26]. The stars label our assigned freeze-out points.    Table I. It appears that our theoretical results in the thermodynamical limit (L = ∞) do not fit the experimental values well at low collision energy, whereas the deviations are smaller if the finite size parameter L changes. We illustrate the presently calculated relation between the baryon chemical potential µ f B and the center-of-mass energy of the collision, S N N , and the comparison with those given in lattice QCD simulations (e.g., Ref. [43]) and model calculations (e.g., Refs. [34,39]) in Fig. 3. We see from (color online) Comparison of presently obtained S N N dependence of the baryon chemical potential in cases of L = ∞ and L = 2.2 fm with those given in lattice QCD simulation [43], HRG model [39] and the parameterized one in SHM model [34].
Only the freeze-out points with small deviation between theory and experiment in Table I are used for fitting the freeze-out conditions, and the obtained best-fitted parameters are listed in Table II. The fitted µ f B ( S N N ) curve is also displayed in Fig. 3.
With the parametrization, one can predict the freezeout parameters (µ f B , T f ) of the system generated in any collision energy. For example, with L = 2.2 fm, S N N = 5.8 GeV, 7.7 GeV, 11.5 GeV, 14.

IV. PHASE DIAGRAM AND FURTHER PREDICTION
With the quark propagator obtained by solving the DSE, we can get the temperature and chemical potential dependence of the quark condensate and the quark dynamical mass, which are commonly regarded as appropriate order parameters of chiral phase transition. Taking the chiral susceptibility criterion [8,75] we determine the upper and lower boundary of the chiral phase crossover region by the full width at half maxima (FWHM) of the susceptibility. The chiral susceptibility is defined as: where the constituent quark mass is defined as the real part of M (ω 2 0 , 0) in Eq. (11), the same as in Sec. II C. The obtained crossover regions in cases of L = ∞, 2.5 fm and 2.2 fm are shown as the shadowed regions in Fig. 4. With the chiral susceptibility criterion [8,75] or the fluctuation criterion [9], we determine the boundaries of the first order transition region and the location of the CEP. The obtained results in the two cases are displayed in Fig. 4. We illustrate also the presently obtained chemical freeze-out lines in these cases in Fig. 4. The figure manifests that the chemical freeze-out happens in the obtained chiral crossover region. Quantitatively, the freeze-out temperature T f = 155.7 MeV is higher than the T c = 144.0 MeV in the thermodynamical limit (i.e., with L = ∞) at µ B = 0, but the deviation is smaller in case of L = 2.2 fm where T f = 128.4 MeV, T c = 121.6 MeV at µ B = 0. We also notice that the finite size effect shifts the location of the CEP to higher baryon chemical potential and lower temperature drastically: in case of L = ∞, 3.5 fm, 2.8 fm, 2.5 fm and 2.2 fm, (µ CEP B , T CEP ) = (296,132), (360,121), (426,111), (510,101), (711,77.5) MeV, respectively. This is consistent with the behavior given in phenomenological model calculations [48,51]. Compared with the location of CEP calculated from functional renormalization group (FRG) approach [80,81], it indicates that the inclusion of the finite size effect is reasonable for 2.2 fm ≤ L ≤ 2.8 fm. (color online) Calculated collision energy √ SNN dependence of κσ 2 = χ B 4 /χ B 2 at the freeze-out line. The black circles are the experimental values [26], the gray triangles stand for our results in case of infinite volume, the blue, green and red points denote our results in the case of L = 2.8 fm, 2.5 fm and 2.2 fm respectively. The shadowed region(s) displays the numerical uncertainties.
It is known that the κσ 2 = χ 4 /χ 2 is a direct observable in experiment and may demonstrate the property of the states around the CEP well. We calculate χ B 4 /χ B 2 in the T -µ B plane and pick out the value along the freezeout line to get the S N N dependence of χ B 4 /χ B 2 . The obtained results in case of thermodynamical limit, finite size of L = 2.8 fm, 2.5 fm and 2.2 fm are depicted in Fig. 5. In case of L = 2.2 fm, the lowest collision energy calculated is 5.8 GeV. It is apparent that, without considering the finite size effect, our calculated χ B 4 /χ B 2 decreases more rapidly than the experimental data as the S N N descends. With the finite size effect being taken into account, we can reproduce the experimental data excellently. In the lower collision energy region, the κσ 2 exhibits a nonmonotonic behavior, whose minimum is reached at S N N ≈ 10 GeV, and then increases drastically as S N N further decreases.
Since the calculated kurtosis fit the experimental data best in case of L = 2.2 fm, we further calculate the S N N dependence of χ B 6 /χ B 2 from our fitted freeze-out line under this finite size parameter, depicted in Fig. 6, whose numerical uncertainties are much larger than the results of κσ 2 though.  [82] for 0 -10% and 30 -40% centralities respectively, and the red points denote our results in case of 2.2 fm. The shadowed region(s) displays the numerical uncertainties.

At
√ S N N = 200 GeV, we obtain a negative value of χ B 6 /χ B 2 = −0.97 ± 0.36, which is qualitatively consistent with the experimental results [82][83][84] and lattice calculation [85]. For √ S N N 10 GeV, our result shows a nonmonotonic behavior with a shallow minimum, in agreement with the FRG calculation [86]. We also show that χ B 6 /χ B 2 may have a large maximum along with a sharp and deep minimum when √ S N N 10 GeV. As a result, we predict a complex nonmonotonic behavior of χ B 6 /χ B 2 as a function of collision energy.

V. SUMMARY
In summary, we have calculated in this work the baryon number susceptibilities in a two-flavor quark system via the DSE approach of QCD in case of not only thermodynamic limit but also finite size. By comparing the calculated ratios χ B 1 /χ B 2 and χ B 3 /χ B 1 with the experimental data of the net-proton multiplicity distribution in BES at RHIC, we obtained the temperature and the baryon chemical potential at the chemical freeze-out states. We calculated also the collision energy dependence of the κσ 2 at the freeze-out line and observed an excellent agreement with experimental data when taking into account the finite size effect. It shows that the finite size effect is significant in studying the QCD phase transitions with RHICs. The obtained collision energy S N N dependence of the κσ 2 exhibits a nonmonotonic behavior in lower collision energy region. We also predict that the collision energy dependence of hyper-order cumulant ratios such as χ B 6 /χ B 2 may also be nonmonotonic.