Fluctuating temperature and baryon chemical potential in heavy-ion collisions and the position of the critical end point in the effective QCD phase diagram

We use the linear sigma model with quarks to locate the critical end point in the effective QCD phase diagram accounting for fluctuations in temperature and quark chemical potential. For this purpose, we use the non-equilibrium formalism provided by the superstatistics framework. We compute the effective potential in the high- and low-temperature approximations up to sixth order and include the contribution of ring diagrams to account for plasma screening effects. We fix the model parameters from relations between the thermal sigma and pion masses imposing a first order phase transition at zero temperature and a finite critical value for the baryon chemical potential that we take of order of the nucleon mass. We find that the CEP displacement due to fluctuations in temperature and/or quark chemical potential is almost negligible.


I. INTRODUCTION
The study of the transition describing the phase change from nuclear to quark-gluon matter with increasing temperature (T ) and baryon chemical potential (µ B ), constitutes one of the most active fields of research of modern high-energy nuclear physics. The description of this transition is encoded in the so called QCD phase diagram [1] represented by a temperature vs. baryon chemical potential plane, where the different phases of strongly interacting matter can be identified. From the theoretical side, efforts to describe this diagram have been carried out employing a wide range of tools such as finite energy sum rules, Dyson-Schwinger equations, functional renormalization, holography and effective models . Lattice QCD (LQCD) calculations are also a useful technique [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40], although these cannot be used to access large values of µ B , given the severe sign problem [41].
On the other hand, relativistic heavy-ion collisions provide the experimental tool to access the properties of the QCD phase diagram. In recent years the STAR BES-I program has analyzed data from nuclear collisions in the energy range 200 GeV > √ s N N > 7.7 GeV [42], to explore deeper into the QCD phase diagram reaching nuclear matter at higher densities. Also, new experiments will soon enter into operation, providing data at lower collision energies with higher luminosity to better study the properties of baryon-rich matter [43].
It is well known that fluctuations play a relevant role for the analysis and interpretation of heavy-ion data [44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60]. In particular, when T and/or µ B are not uniform over the entire reaction volume, fluctuations from the average values can be accounted for using the so-called superstatistics scenario [61] where a non-extensive be-havior naturally emerges due to these fluctuations. Indeed, at the onset, the reaction starts off from regions the size of overlapping nucleon pairs. If overall thermalization is to be achieved, it seems natural to assume that these regions form subsystems from where thermalization spreads later over the entire reaction volume. In this scenario, the temperature and chemical potential between subsystems may not be the same. Thus, a superposition of two statistics needs to be considered: one in the usual Gibbs-Boltzmann sense for particles in each subsystem and another one for the probability to find particular values for T and µ B for different subsystem.
In this work, we study the QCD phase diagram implementing superstatistics with fluctuations both in T and µ B . We concentrate on the location of the critical end point (CEP) using the linear sigma model with quarks (LSMq) as an effective model to find the boundaries on the phase diagram from the chiral symmetry restoration/breaking point of view. We show that when conditions for a first order phase transition are enforced to happen for a critical value of µ B = µ c B at T = 0, the position of the CEP varies little as compared to an analysis where one only follows the evolution of the crossover transition from high T without requiring a first order phase transition at µ c B . The work is organized as follows. In Sec. II we review the superstatistics ideas when temperature and chemical potential fluctuate between different subsystems. We expand the generalized Boltzmann factor to first order in 1/N , where N is the number of subsystems that make up the whole system. In Sec. III we discuss the calculation of the effective potential in the LSMq, including screening effects together with the effective couplings computed at finite T and µ B . In order to explore a wide region in the phase diagram, the effective potential is computed analytically in two regimes; first at low temperature and high chemical potential, and then at high temperature and low chemical potential. In Sec. IV we spell out the conditions that give rise to the equations to find the values of the model coupling constants. In Sec. V we use these couplings to compute the critical T and µ B that define the transition curves and locate the CEP, for the cases of superstatistics with varying number of initial subsystems and comparing to the result in the thermodynamic limit. Finally, we summarize and conclude in Sec. VI. We reserve for the appendices the explicit computation of the vacuum stability conditions and the temperature and baryon chemical potential dependence of the effective coupling constants.

II. SUPERSTATISTICS
To include the superstatistics effects, we closely follow Ref. [62]. accounting also as an extra ingredient fluctuations in the ratio µ B /T . Recall that for a thermodynamic system which contains an space fluctuating intensive variable β, such as the inverse temperature or the product of the chemical potential and the inverse temperature, one may consider the full system as made up of subsystems where β is constant with a local Boltzmann factor e −βĤ , withĤ being the Hamiltonian. A generalized Boltzmann factor B(Ĥ) can be defined as the normal Boltzmann factor e −βĤ weighted by a distribution function f (β), Let us consider the two fluctuating intensive variables β = 1/T and η = µ/T . The local Boltzmann factor is e −β(Ĥ−µQ) , whereQ is the number operator and we have assumed that each subsystem is part of a Grand Canonical ensemble. The generalized Boltzmann factor may be defined in a similar fashion as in Eq. (1) with probability distribution F that depends on both variables, namely, F (β, η). Assuming that the two variables are statistically independent, then F (β, η) = f (β)g(η) and the generalized Boltzmann factor becomes To carry out the calculations, we consider that f (β) and g(η) correspond to a χ 2 distribution function F given by where Γ is the gamma function, N represents the number of subsystems in the full system, x corresponds either to T or η. The average of the random variable x is given by The χ 2 describes the distribution of the sum of N random variables X i , each of which are in turn distributed obeying a Gaussian distribution. Therefore, the χ 2 distribution is well-suited to describe random variables with positive definite values. To evaluate Eq. (2) with the distribution functions given by Eq. (3), we can use that whereÂ is any operator. The Taylor series expansion of Eq. (5) is Therefore, using thatQ commutes withĤ, Eq. (2) can be written as Expanding Eq. (7), up to first order in 1/N , we get the generalized Boltzmann factor Recall that the partition function is Z = T r[B(Ĥ,Q)], therefore, we explicitly get with where Ω and V ef f (β 0 , η 0 ) represent the volume and the effective potential, respectively. Changing the variable β 0 to T 0 , the partition function is explicitly obtained in terms of T 0 and η 0 , To obtain the effective potential with superstatistics corrections, recall that Thus The logarithm of the partition function in Eq. (11) is given by Substituting Eq. (14) into Eq. (13), we obtain the effective potential with corrections due to fluctuations in β and η Finally, using Eq. (10) into Eq. (15), we obtain the expression for the effective potential in terms of T 0 and η 0 We now proceed to compute the effective potential V ef f using the LSMq.

III. EFFECTIVE POTENTIAL
In order to explore the effective QCD phase diagram from the chiral symmetry prospective, we work with an effective model that takes into account spontaneous symmetry breaking/restoration at finite temperature and density; the LSMq. The dynamical degrees of freedom consist of the lightest quarks together with the pions and the sigma. The Lagrangian is given by where ψ is an SU(2) isospin doublet of u and d-quarks, π = (π 1 , π 2 , π 3 ) is an isospin triplet of pions and σ is an isospin singlet. λ is the boson's self-coupling, g is the fermion-boson coupling and a 2 > 0 is the mass parameter. To allow for an spontaneous breaking of symmetry, we let the σ field to develop a vacuum expectation value v which can later be taken as the order parameter of the theory. After the spontaneous symmetry breaking, the Lagrangian for the LSMq is given by The shift in the σ-field produces that the quarks, the sigma boson and the three pions acquire dynamical masses given by respectively. The tree-level potential is given by whose minimum is given by Since v 0 = 0, we notice that the symmetry is spontaneously broken. The coupling constants λ and g, as well as the mass parameter a are to be determined from physical conditions valid at the phase transition. The minimum of the effective potential v represents the order parameter that evolves when the system approaches chiral symmetry restoration at finite T and/or µ B and vanishes in the symmetric phase. In order to determine the transition conditions as function of temperature and quark chemical potential, we study the behavior of the effective potential.
The effective potential is computed beyond the mean field approximation. This means that we include, in addition to the tree-level contribution and the one-loop correction, both for bosons and fermions, also the ring diagrams contribution and the effective coupling constants. This correction accounts for the plasma screening effects [63]. Also, in order to consider a non-vanishing pion mass, we add to the Lagrangian an explicit symmetry breaking term and thus and take the vacuum pion mass as m π ≃ 139 MeV and the vacuum sigma mass m σ = 500 MeV. As a consequence of the explicit symmetry breaking, the effective potential at tree-level has a minimum at a value of v given by The one-loop contribution contains vacuum as well as matter terms. The vacuum piece can potentially distort the tree-level potential, shifting the position of the minimum and distorting its curvature. Since the properties of vacuum should be independent of the perturbative order of its description, we introduce vacuum stability countertemrs δa 2 and δλ so that the tree-level potential is now written as δa 2 and δλ are computed from requiring that the minimum and the curvature at the minimum of the sum of the vacuum piece coming from the one-loop effective potential and the tree-level potential in Eq. (25) do not change with respect to the values they had without loop corrections. The explicit calculations of the counterterms is given in Appendix A. These conditions yield In order to find the ring diagrams correction to the effective potential, one needs to compute the meson selfenergies, which serve as the temperature and density infrared regulators. Furthermore, an improved description of the properties of the phase diagram can be achieved when including temperature and density corrections to the couplings, which we carry up to one-loop order.
The effective potential up to the ring diagrams contribution, the boson and fermion self-energies, and the effective couplings, can be analytically computed in the low and the high-temperature expansions. In the former, one considers that the largest energy scale is provided by µ q , such that M/µ q , T /µ q ≪ 1, where M is any of the particle masses. In the latter, it is only necessary to consider that M/T ≪ 1, regardless of the relation between T and µ q . The explicit expression for the effective potential in the low-T expansion is given by where for the matter contribution to the one-loop effective potential, we have included terms up to O(T ) 6 . We use the expansion technique described in Ref. [64] for the fermion case. Also, we have adopted the MS regularization scheme using µ c e 1/2 as the renormalization scale, with µ c the critical quark chemical potential at T = 0 in the QCD phase diagram. We have made tests scaling the renormalization scale µ c e 1/2 → 3µ c e 1/2 and found that the results vary only slightly.
In the high-T expansion, the effective potential is given where for the matter contribution we have included terms up to O(M ) 6 , using the expansion technique described in Ref. [63]. The expression for the boson self-energies, Π b , in the low-and high-temperature approximations are given by and respectively. The temperature and density corrections to the couplings, accounting for the modification of the intensity of the interaction around the phase transition region, are explicitly computed in Appendices B and C, and given by and where for the fermion mass m f in the low temperature approximation of g ef f LT we have used the fermion selfenergy given by These effective couplings enter the effective potential through the boson and fermion masses which are now written as as well as through the counterterms δa 2 and δλ. Before proceeding to explore the properties of the phase transition including the superstatistics effects, we first need to determine the free parameters of the LSMq which are appropriate to the conditions of finite temperature and density around the transition lines.  FIG. 1. QCD phase diagram obtained from the low-and hightemperature approximations with µc = 300 MeV and λ = 1.5, and the corresponding g = 1.77 and a = 102 MeV. The CEP is located within the full circles on each curve obtained using the low temperature expansion.

IV. DETERMINATION OF THE FREE PARAMETERS
The effective potential contains three free parameters that need to be fixed, namely the couplings λ and g, and the mass parameter a. These parameters are to be determined using physical input valid at the phase transition. For this purpose, we enforce that at T = 0 and for a critical value of the quark chemical potential µ c = µ B c /3 [65], with µ B c the critical baryon chemical potential, the effective potential describes a first order phase transition. In analogy with the Hagedorn's limiting temperature concept [66] extended to finite µ B , we take µ B ≃ m B , where m B ≃ 1 GeV is the typical value of the baryon mass. Recall that for a first order phase transition, the effective potential develops two degenerate minima. The value of v ≡ v * = 0 at the minimum becomes a new quantity that needs to be also determined. However, as we proceed to show, the conditions to describe a first order phase transition provide only three equations to determine the four unknowns. Our strategy to find the solutions consists of finding the model parameters when varying one of them and choose λ as the parameter to vary. First, in order to determine a, we use the relation between the σ and π dynamical masses at T = 0 and µ q = µ c . This involves Eqs. (20), including the matter corrections coming from the self-energies, Eq. (29). Thus, we have  The remaining two equations for the two unknowns, v * and g, are given by Equations (38) describe the conditions for the effective potential to show degenerate minima at T = 0 for µ q = µ c , one of them at v = 0 and the second one at v = v * . The appropriate expression to be used for the self-energy corresponds to the case of the low temperature regime. Although the above described procedure fixes the model parameters for conditions valid at the putative first order phase transition, we also require consistency of the description with the crossover transition at µ B = 0 and finite T . From LQCD [67] calculations, it is wellknown that this transition happens at T c 0 ≃ 155 MeV for 2+1 light flavors and at T c 0 ≃ 170 MeV for 2 light flavors. In order to get solutions of Eqs. (37) and (38) satisfying the conditions at both ends of the transition curve with T c (µ = 0) < 200 MeV we find that λ is required to lie between the very restricted range 1.4 < λ < 1.5. The lowest limit for these λ values corresponds to a choice for µ c = 300 MeV, whereas the upper limit corresponds to µ c = 290 MeV.

V. QCD PHASE DIAGRAM
We now, proceed to use the effective potential we have determined, into the supersatistics formalism. Following the procedure developed in Sec. II, we substitute either Eq. (27) or Eq. (28) into Eq. (10) for the low and hightemperature descriptions, respectively. In this manner, we can obtain the partition function using Eq. (11) and from its logarithm, we write the effective potential, including the superstatistics corrections, as To explore the QCD phase diagram, we identify the transition lines where chiral symmetry is restored. The procedure consists of finding the transition temperatures and chemical potentials using first the high temperature approximation for the effective potential. We start from µ q = 0 and stop when the ratio of the critical chemical potential and temperature is ≈ 0.8 from where we start using the low-temperature approximation for the effective potential to continue finding the critical curve.
From this procedure, we find either second or first order transitions. For the former, the vacuum expectation value continuously moves from a finite value to zero and the phase transition occurs when this vanishes. For the latter, the phase transition is identified when the two developed minima become degenerate. Figure 1 shows the effective QCD phase diagram using λ = 1.5 and µ c = 300 MeV with a = 102 MeV and g = 1.77. In the region where the HT approximation is valid, we find only second order phase transitions (our proxy for the crossover phase transition) and these are shown as the curve with 0 ≤ µ q < 150 MeV. In contrast, both second and first order phase transitions are found deep in the phase diagram where the LT approximation is valid. This happens for 260 MeV < µ ≤ µ c . The CEP location, where the phase transition changes order, is indicated with a full circle at µ q ≈ 291 MeV and T ≈ 40 MeV. Figure 1, also shows the variation on the transition line and the CEP due to the superstatistics effects in temperature (S T ), in µ q /T (S η ) and in both parameters at the same time (S T η ), with N = 100 and N = 300 for all combinations. We notice that as N increases, the corresponding CEP moves towards the corresponding CEP in the thermodynamic limit, labeled as N = ∞. The changes due to any of the superstatistics effects are quite small and for all practical purposes they are neglegible. Figure 2, shows the effective QCD phase diagram using other set of allowed parameters, λ = 1.6 and µ c = 290 MeV with, a = 93.4 MeV and g = 1.85, the general features of the phase diagram in this case do not change with respect to the previous case. However, it is relevant to mention that the CEP is now located at µ q ≈ 271 MeV and T ≈ 51 MeV. Nevertheless, the superstatistics effects are also negligible. Finally, we notice that the systematics of the CEP displacement for the two explored sets of parameters is the same.

VI. SUMMARY AND CONCLUSIONS
In this work we have used the LSMq to locate the CEP in the effective QCD phase diagram taking into account fluctuations in the temperature and the quark chemical potential. For this purpose, we have implemented the superstatistics scenario assuming that both fluctuations in the temperature and quark chemical potential are described according to a χ 2 distribution. We computed the superstatistics effective potential up to order 1/N . To numerically find the effects, N has been taken as about half the number of nucleons in the collision of heavy ions.
In order to locate the CEP we have used the LSMq, computing the effective potential in the high-and lowtemperature approximations up to sixth order, including the ring diagram contribution which accounts for the plasma screening effects. We fix the model parameters, namely the coupling constants and the mass parameter, imposing a first order phase transition at zero temperature and a finite critical quark chemical potential, using the Hagedorn limit temperature concept applied to finite baryon density. We fix this critical baryonic chemical potential µ c B ∼ 1 GeV. We find that the CEP displacement due to fluctuations in temperature and chemical potential are of order ∼ 0.1 MeV and are much smaller compared with the case when only fluctuations in temperature are considered, and the model parameters are fixed in the high-temperature effective potential [62]. As noted also in Ref. [62], N can be associated with the number of participants in the heavyion collision, the specific heat and the smallest of the mass numbers of the colliding nuclei. Thus, in order to give a more accurate estimation of how much the CEP is displaced under appropriate experimental conditions, these need to be included for the estimation of the parameter N . values, even after including the vacuum pieces stemming from the one-loop corrections. These conditions are 1 2v where V vac is the one-loop vacuum piece of the effective potential. The solution for the counterterms δa 2 and δλ is given by (A3)

B. Effective coupling constants (λ)
We now compute the one-loop correction to the coupling λ, including thermal effects in the high temperature and low temperature regimes. Figure 3 shows the Feynman diagrams that contribute to this correction. Columns (a), (b), (c), (d), (e) and (f) contribute to the correction to the σ 4 , (π 0 ) 4 , (π + ) 2 (π − ) 2 , σ 2 π + π − , (π 0 ) 2 π + π − and σ 2 (π 0 ) 2 terms of the interaction Lagrangian in Eq. (17), respectively. Since each of these corrections lead to the same result, we concentrate on the diagrams in column (a). Each of the three diagrams involves two propagators of the same boson. For the first two diagrams the intermediate bosons are neutral and for the third one the intermediate bosons are charged. Therefore the expression for the diagrams can be ob-tained from where P i ≡ (ω, p) is the total incoming four-momentum, K ≡ (ω n , k), ω n = 2nπT is the bosonic Matsubara frequency and D is the boson propagators defined as Considering the permutation factors and the contribution from the s, t and u-channels, the correction to the self-coupling λ to one-loop order is given by we have to compute Eq. (B1), explicitly The sum in Eq. (B4) is calculated using Ref. [68]. We obtain where Calculating in the infrared limit (E 1 = E 2 ≡ E), we get First, we study the high temperature (HT) case. We focus on the matter term of Eq. (B7), notice that in the second of Eqs. (B8) we retain only the dominant term. We define y ≡ m/T and x ≡ k/T and get The integrals over x can be expressed in terms of the well known functions [69] h n (y) = 1 Γ(n) ∞ 0 dx x n−1 which satisfy the differential equations Therefore, Using the high temperature expansions for h 1 (y) and f 1 (y) [69] h 1 (y) = π 2y + 1 2 ln y 4π and keeping the leading terms, we get Using Eq. (B14) into Eq. (B3), we finally obtain Effective boson self-coupling λ eff , in the lowtemperature approximation, as a function of T .
We now calculate the low temperature (LT) case. We focus on the matter term of Eq. (B7) Note that in the second equality of Eq. (B16), there is a term that is zero when T → 0, if we make the change of variable w = √ k 2 + m 2 , then dk = wdw √ w 2 −m 2 . Therefore, we get as T → 0, then e −w/T << 1 and we can use the geometric series in the form thus, using Eq. (B18) into Eq. (B17), we obtain (B19) Notice that the integral of (B19) gives where in the last equality we have expanded in a Taylor series for T << 1. Therefore we get and thus, Finally, using Eq. (B22) into Eq. (B3), we get C. Effective coupling constants (g) We now turn to the calculation for the one-loop correction of the coupling g. Figure 6 shows the Feynman di- agrams that contribute to this correction. Columns (a), (b) and (c) contribute to the correction to the quarkσ, quark-π 0 and quark-π ± terms of the interaction Lagrangian of Eq. (17), respectively. Since each of these corrections lead to the same result, we concentrate on the first diagrams in column (a). The expression for this diagram is written as where D(K) is given by Eq. (B2) and S(K) is the fermion propagator therefore, considering all diagrams, the correction of the self-coupling g to one-loop order is given by In order to proceed, we have to compute G(P i ; m 2 ). Using trace algebra and applying the Dirac equation in Eq.(C1), we get In order to apply the imaginary-time formalism of finite temperature field theory, we make the substitutions p 0 → iω m + µ for fermion part, k 0 → iω n for boson part, where ω n = 2nπT , and therefore, we obtain where Both sums in Eq. (C6) are calculated using Ref. [68]. We obtain To study the high (HT)-and low-temperature (LT) regimes, we analyze Eq. (C8) and Eq. (C9) in each regime. First, we compute the high temperature regime. For the case of I HT , the dominant terms behave as ∼ T 2 . This behavior is obtained for the case when s 1 = −s 2 in (C8). Considering only the matter piece, we get where have we used 1 + f (−E) = −f (E) andf (−E 1 − µ) = 1 −f (E 1 + µ). Taking E 1 = E 2 , we obtain  Therefore, I HT does not contribute to terms of order T 2 . Now we compute II HT . Once again, terms that behave like ∼ T 2 are obtained when s = −s 1 = −s 2 in Eq, (C9). Therefore we get, where have we used 1 + f (−E) = −f (E) andf (−E 1 − µ) = 1 −f (E 1 + µ), p 0 1 = iω 1 − µ and p 0 2 = iω 2 − µ. Taking E 1 = E 2 , we obtain Taking p 0 1 = 0 and p 0 2 = 0, we get Therefore, using Eq. (C14) and Eq. (C11) in Eq. (C6), we obtain In the high-temperature case, we can approximate k 2 + m 2 π ∼ k in the denominator of the previous integral. Then, we have an analytical integral of the kind