Chiral condensate and the equation of state at nonzero baryon density from the hadron resonance gas model with a repulsive mean field

We study the QCD equation of state and the chiral condensate using the hadron resonance gas model with repulsive mean-field interactions. We find that the repulsive interactions improve the agreement with the lattice results on the derivatives of the pressure with respect to the baryon chemical potential up to eighth order. From the temperature dependence of the chiral condensate we estimate the crossover temperature as a function of baryon chemical potential, $T_{pc}(\mu_B)$. We find that the chiral crossover line starts to deviate significantly from the chemical freeze-out line already for $\mu_B>400$ MeV. Furthermore, we find that the chiral pseudocritical line can be parametrized as $T_{pc}(\mu_B)/T_{pc}(0)=1-\kappa_2 (\mu_B/T_{pc} (0))^2-\kappa_4 (\mu_B/T_{pc} (0))^4$ with $\kappa_2=0.0150(2)$ and $\kappa_4=3.1(6) \times 10^{-5}$, which are in agreement with lattice QCD results for small values of $\mu_B$. For the first time we find a tiny but non-zero value of $\kappa_4$ in our study.


I. INTRODUCTION
Understanding QCD matter at nonzero temperatures and baryon densities presents several challenges on the theory side.Apart from asymptotically large temperatures and baryon densities when perturbation theory can make predictions for thermodynamics, the problem is nonperturbative.For zero baryon density, lattice QCD methods have provided continuum extrapolated results for physical quark masses for many thermodynamic quantities.However, due to the infamous sign problem, the lattice Monte Carlo techniques based on the method of important sampling break down at finite densities.Attempts to extend the lattice QCD calculations to nonzero baryon densities through Taylor expansion in baryon chemical potential, µ B , or the imaginary chemical potential approach has allowed us to understand QCD thermodynamics µ B < 2.5 T [1][2][3][4], which can be extended in a certain temperature regime to µ B /T ≈ 3.5 [5].
Below the chiral crossover temperature, QCD thermodynamics can be fairly well understood using the hadron resonance gas (HRG) model.The main idea of this model is that the interactions between hadrons can be taken into account through hadronic resonances and the interacting gas of hadrons can be replaced by a gas of noninteracting hadrons and hadronic resonances treated as stable particles.This approach can be justified within the framework of relativistic virial expansion [6,7] and the comparison between the HRG model results and the lattice data [8][9][10][11][12][13].Indeed, the most recent lattice QCD calculations of the QCD pressure and the second order fluctuations and correlations of conserved charges agree well with the the HRG model results [2,[14][15][16][17][18].For some quantities like baryon-strangeness or baryon-charm correlations, it is important to include additional hadron states that are not listed by the Particle Data Group (PDG) but are predicted by lattice QCD and quark models [18][19][20][21].For strangeness-baryon number correlation, this idea is further supported by the calculation within the relativistic virial expansion with state-of-the-art phase shift analysis [22].
The temperature dependence of the chiral condensate has been studied within the HRG model [23][24][25][26][27], and attempts to use the HRG model to estimate the chiral crossover temperature at small baryon density have been made recently [27].There is a widespread expectation that the chiral crossover may become a true phase transition at large baryon density, namely that the crossover will turn into a first-order phase transition at µ B = µ CEP B corresponding to the critical end-point (CEP).Lattice QCD results strongly disfavor µ CEP B < 400 MeV for a range of temperatures around 0.85 T pc , while functional renormalization group studies estimate µ CEP B to be even larger, around 635 MeV [28] at T = 107 MeV.Therefore, it would be interesting to estimate the chiral crossover temperature at larger values of baryon density using the HRG model framework.
The agreement between lattice QCD and the HRG model results no longer holds for higher-order fluctuations of conserved charges [2][3][4]29].This may highlight the importance of the inclusion of non-resonant attractive interactions as well as repulsive interactions among the baryons, which could become very important at higher baryon densities.The most common approach to include repulsive interactions within the HRG model is through the inclusion of a hard-core repulsive excluded volume for hadrons [30][31][32][33][34][35][36][37].Recently these models were refined to include both repulsive and attractive nonresonant interactions via a van der Waals-like potential [38][39][40][41][42]. Another approach to include the repulsive interactions is through the mean-field approximation.Here the strength of the repulsive interaction is proportional to the baryon density [43][44][45].It has been shown that this model can explain the deviations from the HRG model seen in the lattice QCD calculations of higher-order derivatives of the QCD pressure with respect to the baryon chemical potential [46].Following this work the use of the HRG model with repulsive mean-field interactions was explored to explain the lattice results on the QCD equation of state [47,48] and fluctuations and correlations of conserved charges [49,50].
The aim of this paper is to study the temperature and the µ B dependence of the chiral condensate, and using this to estimate the chiral crossover temperature as a function of µ B .We will also revisit the description of the higher-order baryon number fluctuations and the QCD equation of state within the HRG model with repulsive mean-field interactions in light of the newest lattice QCD results.The paper is organized as follows: In the next section, we review the HRG model with repulsive interactions.In that section we also present the comparison of this model with the state-of-the-art lattice QCD results on fluctuations of conserved charges.In Sec.III we discuss the calculation of the chiral condensate as a function of temperature and baryon chemical potential.In Sec.IV we present our results, i.e., the estimate of the chiral crossover temperature as a function of baryon chemical potential, and the curvature of the pseudocritical line.In Sec.V, we discuss the temperature and baryon density dependence of the speed of sound for the strangeness neutral system.Finally, our conclusions are presented in Sec.VI.In the Appendices, we give some technical details of the calculations.

II. INCLUDING REPULSIVE BARYON INTERACTIONS WITHIN THE HRG MODEL
It is well known from nucleon-nucleon scattering experiments that there are repulsive interactions between nucleons at short distances (or equivalently at high energies).From the lattice results we know that this is also true for other baryons [51].These interactions will become important at large baryon density.Since we are interested in extending the HRG model to high baryon density the effect of the repulsive baryon-baryon interactions has to be included.The repulsive baryon-baryon interactions could play a role also at zero or small baryon densities as the temperature increases toward T pc , since the abundance of baryons and antibaryons increases.But from previous analysis we know that the corresponding effect is small for the pressure and energy density [47].However, higher-order derivatives of the pressure with respect to the baryon chemical potential are sensitive to the effect of the repulsive interactions even at zero baryon density [46].This is because these derivatives carry information about the physics at high baryon density even when evaluated at zero baryon chemical potential.Comparison of HRG calculations with the lattice results on higher-order derivatives of the pressure with respect to baryon chemical potential allows one to constrain the contribution of the repulsive interactions.As already mentioned in the introduction in this work we use the repulsive mean-field to include the effect of repulsive baryon-baryon interactions in the HRG model.
The pressure for an interacting ensemble of baryons and anti-baryons with densities n b and nb respectively can be represented within the repulsive mean-field model as [48,49] Here the effective chemical potential for the ith species is defined as µ ef f = B i µ B − Kn b{ b} , with B i being the baryon number.B i = +1 for baryons, and B i = −1 for antibaryons respectively, and β = 1/T , where T is the temperature of this system.The number densities for baryons and anti-baryons can be solved self-consistently from the following pair of equations: The total pressure within the HRG model is the sum of the contributions from this interacting ensemble of the (anti)baryons and the noninteracting ensemble of mesons.We include the quark model (QM) predicted states [52,53] in addition to the Particle Data Group (PDG) listed hadrons up to mass 3 GeV following [54,55].The PDG list was made considering all the states, even the ones with large uncertainties, and to avoid double counting the QM states are replaced by the experimentally determined states if they have the same mass and quantum numbers.We denote this HRG model with the extended list of particles as the QMHRG model.For µ B = 0 we can use the Boltzmann approximation for baryons and antibaryons due to their large masses.Therefore, for the pressure and the number density we could write the following set of equations [46]:  [57] results and the circles corresponding to Wuppertal-Budapest (WB) results [2].The gray and yellow bands represent the fluctuation data in the continuum limit from HotQCD [57] and WB [58] Collaborations.
This approximation greatly simplifies the calculations of the baryon number susceptibilities, which are defined as The explicit expressions of χ B n are presented in Appendix A. In this work, we have included the partial pressures of the QMHRG states with repulsive interactions among the (anti) baryons at the mean-field level.Typically in the earlier works the mean-field coefficient has been chosen to be K = 450 MeV fm 3 which is equivalent to 56.25 GeV −2 [46,56] from phenomenological considerations.Here we suggest a new way of estimating this constant.Our approach is motivated by the fact that we have high precision lattice QCD data for χ B 2 and χ B 4 extrapolated to the continuum limit [4,18,29].Therefore, we adjust the value of K to reproduce these lattice QCD results.In Fig. 1 we show the comparison of our results for the second-and fourth-order baryon number fluctuations with the corresponding lattice QCD results.The lattice QCD results on χ B 2 results disagree with the QMHRG model results for T > 150 MeV, while the lattice QCD results for χ B 4 start to disagree with QMHRG model results already for T > 140 MeV, when no repulsive interactions are taken into account.By including repulsive mean-field interactions within the QMHRG model with K = 33 GeV −2 , we obtain a good agreement with the lattice QCD results for both χ B 2 and χ B 4 up to T ≈ 155 MeV.In addition there are lattice QCD results for N τ = 8 and N τ = 12 from the HotQCD Collaboration [57] and N τ = 12 results from the Wuppertal-Budapest Collaboration [2] for spatial volume V 1/3 T = 4.Here N τ is the temporal lattice extent that is related to the temperature and the lattice spacing a, as N τ = 1/(aT ).Very recently continuum extrapolated lattice QCD results on for sixth-and eighth-order baryon number fluctuations were released for small spatial volume V 1/3 T = 2 [58].
In the lower panels of Fig. 1 we show the comparison of the QMHRG model results with the corresponding lattice QCD results for χ B 6 and χ B 8 .We see from the figure that the QMHRG model results with repulsive mean-field and K = 33 GeV −2 agree well with the lattice QCD results for χ B 6 and χ B 8 within errors.The agreement with Wuppertal-Budapest lattice results is especially remarkable as these results have relatively small errors.These lattice QCD results clearly disagree with the QMHRG model without repulsive interactions; see Fig. 1 (lower panels).We thus conclude that, in order to apply the QMHRG model at large values of baryon number density to describe QCD, repulsive baryon-baryon interactions have to be taken into account.Here we also note that based on the lattice QCD estimates of the charm quark pressure [59,60] we do not expect significant contribution from quark degrees of freedom to thermodynamics quantities below T pc .Therefore, thermodynamic estimates based on purely hadronic models will also agree with the lattice data up to T pc .
In the above discussion, we assumed that both baryons and baryon resonances contribute equally to the mean field.This is certainly a very simplistic assumption.One may expect that not all baryon resonances contribute to the mean field at the same level because of their short lifetimes.As an extreme assumption, we can assume that only ground state baryons contribute to the repulsive mean field.As shown in Appendix A, we can describe the lattice results on χ B 2 , χ B 4 , χ B 6 , and χ B 8 also in this case if the parameter K is increased from K = 33 GeV −2 to K = 100 GeV −2 .We also calculated the QCD pressure and energy density using the QMHRG model with repulsive mean-field interaction.We find that repulsive mean-field interaction has a significant effect on the pressure and energy density for µ B /T ≤ 2, and its inclusion improves the agreement with the lattice QCD data.This is shown in the Appendix B. We will return to the calculations of the equations of state within the QMHRG model with repulsive interactions in Sec.IV D for a strangeness neutral system.

III. CHIRAL CONDENSATE IN THE MEAN-FIELD QMHRG MODEL
The light quark condensate at nonzero temperature and density is defined as where m l is the light quark mass and P is the pressure of the thermodynamic medium described by QCD.We work with two degenerate light quarks, m u = m d = m l .In the HRG model, the derivatives of the pressure with respect to the light quark mass can be written as the mass derivative of the interacting baryonic gas pressure and the ideal mesonic gas pressure.Using Eqs. 1 and 2 the derivative of the interacting baryonic pressure in Eq. 1 can be further written as Here f b{ b} i is the Fermi-Dirac distribution function corresponding to the baryons and antibaryons with the modified chemical potential µ ef f = B i µ i − Kn b{ b} .In the limit K = 0 one obtains the ideal gas result of Ref. [27].The calculations of the mass derivative of the mesonic pressure are the same as in Ref. [27].
The nontrivial input needed for the calculation of the chiral condensate is the dependence of the hadron and resonance masses on the light quark mass m l .In this work, we have estimated these mass derivatives following Ref.[27].The uncertainties in the derivatives of the hadron masses with respect to m l have been estimated in Ref. [27], and we propagate these uncertainties when evaluating the temperature and µ B dependence of the chiral condensate.As in our previous work, we use the renormalized chiral condensate defined by the HotQCD Collaboration [61] as where the parameter r 1 is derived from the static quark potential [62] and d = r 4 1 m s (lim m l →0 ⟨ ψψ⟩ l,0 ) R .In the chiral limit, the light quark condensate has only a multiplicative renormalization factor.The superscript R denotes the renormalized quantity.Taking into account the fact that (lim m l →0 ⟨ ψψ⟩ l,0 ) R = 2Σ and using the values of the low energy constant of SU (2) chiral perturbation theory (χPT), Σ 1/3 = 272(5) MeV and m s = 92.2(1.0)MeV in the MS scheme at µ = 2 GeV from the FLAG 2022 review for the 2+1 flavor case [63], as well as r 1 = 0.3106 fm [64], we obtain the value of d = 0.022791.
We estimate the chiral crossover temperature, T pc as the temperature where the renormalized chiral condensate drops to half of its vacuum value.From the lattice QCD calculations at µ B = 0 we know that this definition of the chiral crossover temperature agrees well with the usual definition of the chiral crossover temperature, defined as the peak position of the chiral susceptibility [61].In our previous work, we showed that, using this definition of T pc and the temperature dependence of ∆ R l within the QMHRG model, one obtains a value of T pc at µ B = 0 as well as for small values of µ B that is only a few MeV larger than the lattice results [27].We expect that the above criterion to estimate T pc should work at larger values of µ B as long as we are far away from a true phase transition.Therefore, in this work, we assume that the chiral transition remains a crossover transition also for large values of µ B and estimate the chiral crossover temperature using the QMHRG model and the above considerations.This is reasonable since there is no definite indication from lattice QCD for a true phase transition at large µ B , and based on symmetry arguments the transition could be a crossover even for very large values of µ B [65].

IV. RESULTS ON THE CHIRAL CROSSOVER LINE
A. Renormalized chiral condensate in the mean-field formalism: One of the goals of this study is to study the dependence of the renormalized chiral condensate on temperature and baryon chemical potentials for a large range of µ B , including values not accessible in lattice QCD calculations.With the already estimated value of the mean-field parameter K = 33 GeV −2 , we have evaluated ∆ l R as a function of µ B for different temperatures, shown in Fig. 2. We have also shown the ∆ l R calculated in the noninteracting QMHRG model for comparison.In the figure, the uncertainty in the values of renormalized chiral condensate obtained in the QMHRG model with repulsive mean-field interactions due to the uncertainties of the quark mass derivatives is represented as the width of the lines.For the noninteracting QMHRG model, we do not show the corresponding uncertainties for better visibility.For temperatures T < 100 MeV, we observe a prominent plateau in the value of ∆ l R for a range of baryon chemical potentials, the extent of which reduces with increasing temperature.Furthermore, we observe that the renormalized condensate falls faster with increasing µ B from the plateau region for the noninteracting QMHRG model compared to the QMHRG model with repulsive mean-field interactions.Thus the repulsive mean-field pushes the transition point toward larger values of µ B .As stated above we use the condition of the chiral condensate dropping to half of its vacuum values to estimate the crossover point.Thus using the values of µ B where ∆ R l is at half of its vacuum value for each temperature we obtain the estimate of crossover line as a function of baryon chemical potential, T pc (µ B ).
We are also interested in estimating the chiral crossover temperature as a function of the net baryon density, n B .Therefore, we have also studied the dependence of ∆ R l on n B .The right panel of Fig. 2 shows the change in ∆ l R at different temperatures as a function of net baryon number density n B scaled by the nuclear saturation density n 0 = 0.16 fm −3 .At a fixed temperature, ∆ l R decreases as the net baryon density increases because baryons significantly contribute to the reduction of the chiral condensate.From this figure, we observe that the density corresponding to the chiral crossover transition for T = 100 MeV slightly exceeds that for T = 25 MeV, for example.This seemingly contradicts the usual expectation that n B decreases as T increases along the chiral crossover line.This is an artifact of the hadron resonance gas model at very low temperatures, as discussed in Sec.IV C.

B. Curvature of the crossover line in the repulsive mean-field QMHRG model
As discussed in the previous subsection, the temperature at which the observable ∆ l R falls to half its zero-temperature value is used as an estimate of the pseudocritical temperature T pc [27] for each value of µ B .The extracted pseudo-critical temperature as a function of µ B , T pc (µ B ), is shown as a function of µ B in the left panel of Fig. 3 using QMHRG model for K = 0 and K = 33 GeV −2 .The uncertainty in the value of ∆ R l results in an uncertainty in the value of the crossover temperature, which, however, is quite small, about the size of the symbols in Fig. 3 or even smaller.In the same figure, the vertical axis is scaled by the estimated pseudocritical temperatures at zero baryon density, which are T pc (µ B = 0) = 161.2(1.7)MeV for the ideal case [27] and T pc (µ B = 0) = 161.5 (1.6) MeV for the repulsive QMHRG model with the mean-field coefficient K = 33 GeV −2 .Our analysis is restricted to µ B ≤ 750 MeV for reasons explained below.As one can see from the figure, the effects of repulsive baryon interactions start to have an impact on the value of T pc for µ B > 300 MeV, pushing the chiral crossover temperature toward larger values and leading to a better agreement with the lattice QCD results on the curvature of the pseudocritical line.On the other hand for µ B < 300 MeV the effects of the repulsive mean-fields are small, and good agreement for the curvature of the pseudocritical line is obtained using QMHRG model with or without the repulsive mean fields included.
In Fig. 3 we also show the recent parametrization of the chemical freeze-out line in heavy ion collisions from Ref. [67].The freeze-out line is normalized by the freeze-out temperature T f = 158 MeV at µ B = 0 [67].For µ B < 300 MeV the freeze-out line agrees well with the pseudocritical line.However, we observe significant differences between the pseudocritical line and the chemical freeze-out line for µ B > 400 MeV, which become larger as the value of the baryon chemical potential increases.Namely, the chiral crossover temperature is always larger than the chemical freeze-out temperature.For µ B = 400 MeV, the chiral crossover temperature is more than 8 MeV larger than the freeze-out temperature while at µ B = 750 MeV it is more than 20 MeV larger than the freeze-out temperature.Thus, the systems produced in heavy ion collisions in the fixed target mode of the beam energy scan (BES) program at the Relativistic Heavy Ion Collider (RHIC) as well as in heavy ion collisions at GSI are expected to undergo a significant nonequilibrium chemical evolution in the hadronic phase unlike the matter produced in heavy ion collisions at higher center-of-mass energies, like the heavy ion collisions at the Large Hadron Collider (LHC) and at RHIC in the collider mode.
We next performed a fit to the obtained values of T pc (µ B ) with the following widely used ansatz: It turns out that a good fit can be obtained by this ansatz for µ B ≤ 750 MeV which results in the following values of the curvature coefficients: κ 2 = 0.0150(2) and κ 4 = 3.1(6) × 10 −5 .The value of κ 2 is consistent within errors with the two most recent continuum estimates from lattice QCD studies with physical masses, κ 2 = 0.012(4) [66] and κ 2 = 0.0153(18) [68], as well as with the earlier lattice QCD estimates [69][70][71].For the first time, we could estimate a nontrivial value of κ 4 which is distinctly different from zero given the estimated errors.Lattice QCD calculations report a κ 4 which is compatible with zero within errors [66,68].The above estimate for the leading curvature coefficient is somewhat smaller than the one obtained from the relativistic nuclear mean-field model calculation, while κ 4 has the opposite sign [72].Allowing for a nonzero value of κ 6 does not improve the quality of the fit, and only leads to much larger errors in κ 2 and κ 4 .Therefore, it is fair to say that, within our accuracy, κ 6 is consistent with zero.

C. Net-baryon density, energy density, and particle composition along the pseudocritical line
The right panel of Fig. 3 shows the estimated values of net-baryon density normalized by the nuclear saturation density n 0 as function of µ B along the pseudocritical line in the T − µ B plane.The baryon density increases monotonically with µ B as expected, but for sufficiently large µ B it reaches a maximum and then decreases with a further increase in µ B .This behavior is clearly unphysical and therefore, the applicability of the HRG model is restricted for large values of µ B .The maximum in n B of 1.53(13)n 0 is reached at µ B = 700 MeV when K = 0 and at µ B = 750 MeV for K = 33 GeV −2 .Thus the inclusion of repulsive mean-field interaction in the QMHRG model calculation helps to extend the model to slightly larger values of µ B .However, our QMHRG model with a single repulsive mean field is too simplistic and is not valid for very large values of µ B .The large value of n B obtained at relatively small values of µ B is due to the contribution from baryon resonances.Since the chiral crossover temperature decreases with increasing µ B , the contribution from these states will eventually decrease and the large value of the net-baryon density near the crossover point should come from nucleons and possibly their attractive interactions.But the simple HRG model cannot describe this switchover.Thus for temperatures close to the chiral crossover temperatures, the QMHRG model is only applicable for µ B < 750 MeV.
The nonmonotonic variation of n B with µ B is also observed in the Nambu-Jona-Lasinio model [73], where the mass-gap equation determines the phase transition line, resulting in a qualitatively similar trend.A similar maximum of the net-baryon density as a function of µ B was earlier observed in a HRG motivated study along the freeze-out boundary in Ref. [74], which was determined from the hadron yields at different collision energies.
It is interesting to study the variation of energy density along the pseudocritical line, which is shown in the left panel of Fig. 4. One can see that the energy density along the pseudocritical line slightly decreases with increasing µ B .The relative contribution to the energy density coming from the baryon sector increases till µ B = 650 MeV as the baryon chemical potential increases, resulting in a shift from a meson-dominated to a baryon-dominated scenario, beyond µ B = 350 MeV.A comparable phenomenon can be observed in heavy-ion collisions, where a transition from meson to baryon dominance during freeze-out takes place at lower energies [75], where the net-baryon density is larger.
Furthermore, we would like to understand the relative abundances of different baryon species with increasing µ B .In the right panel of Fig. 4  relative contribution of the nucleons also increases.The relative contributions of the lowest lying strange baryons and of the lowest ∆ resonances do not change significantly with increasing µ B , either.While the net contribution from higher lying resonances decreases with increasing µ B , this contribution remains significant up to µ = 750 MeV.Even for the largest value of µ B = 750 MeV, where our HRG model is applicable, nucleons contribute less than 50% to the total net-baryon density.Thus the contribution of the baryon resonances remains important even at large values of µ B , and this contribution is responsible for the relatively large values of n B .At the same time, the treatment of the repulsive interactions for baryon resonances is quite simplistic.In future, a more refined QMHRG mean-field model, which treats the repulsive interactions differently in the nucleon, strange baryon, and baryon resonance sectors will be needed to extend the reach of the QMHRG model.

D. Strangeness neutrality and the chiral pseudo-critical line
The system created in heavy ion collisions has zero netstrangeness, since the incoming nuclei do not contain strange particles.Therefore, it is important to estimate the pseudocritical line for the conditions realized in heavy-ion experiments, namely for zero net strangeness, n S = 0.At nonzero values of µ B , one needs to tune the strangeness chemical potential µ S to achieve zero netstrangeness.We performed this tuning to estimate the crossover temperature and the equation of state within our QMHRG model.The values of µ S needed to achieve strangeness neutrality are shown in Appendix C.
In Fig. 5 (left) we show the crossover temperature in the case of n S = 0 as a function of the baryon chemical potential, while in Fig. 5 (right) we show the crossover temperature as a function of the net-baryon number.For a comparison, we also show our results for T pc for the unconstrained, µ S = 0 case in the same figure.We observe that imposing strangeness neutrality pushes the crossover temperature to a slightly larger values compared to the case of µ S = 0 discussed in the previous subsections, whether or not one includes the effect of the repulsive mean field.The effect of imposing strangeness neutrality on the crossover temperature is considerably smaller than the effect of the repulsive interactions.Interestingly enough, the effects of the repulsive interactions are not visible on T pc when plotted as a function of the net-baryon density.This can seen from Fig. 5 (right).The clear difference between the pseudocritical line and the freeze-out line is also apparent for the zero net-strangeness case, cf.Fig. 5, especially when the results are shown in the T − n B plane.
The maximum net-baryon number density that can be reached within our model is slightly reduced compared to the unconstrained, µ S = 0 case to about 1.4n 0 compared to 1.53n 0 as obtained in the previous subsection.This is due to the fact that a positive nonzero value of µ S required by the condition n S = 0 reduces the contribution of strange baryons to the partition function.
K=33, µ S =0 Freeze-out line    In heavy-ion experiments, the chemical freeze-out surface is typically defined as a surface of constant entropy per baryon density, denoted as s/n B [76].Hence it is interesting to study the evolution of the QCD matter formed in heavy-ion collisions along the lines of constant s/n B .The evolution of this matter is subject to several constraints, including strangeness neutrality (n S = 0) and a fixed ratio of net electric charge (n Q ) to net-baryon number (n B ) density due to charge conservation.We study a particular case of strangeness-neutral system with isospin symmetry which is realized when n Q /n B = 0.5, corresponding to a nonzero µ B and µ S but µ Q = 0.It has been shown that the differences in bulk thermodynamic observables that arise from deviations in n Q /n B from the isospin-symmetric value of 0.5 are minimal [29].We have examined the isospin-symmetric scenario for a range of entropy per baryon density values, s/n B = 15, 20, 30, 50, 100, and 400.These values effectively cover the range of center-of-mass energies explored in the RHIC BES program in the collider mode [29].The left panel of Fig. 6 illustrates the isentropic trajectories for different values of s/n B and compares to the lattice QCD results available up to µ B /T = 2.5 from Ref. [29].We have also included the pseudocritical line estimated by us in the same plot for comparison.Our estimates obtained within the mean-field repulsive QMHRG model show a good agreement with the lattice results at temperatures below the pseudo-critical temperatures.For completeness, we have also included the isentropic trajectories for s/n B = 15 and s/n B = 20 calculated within the repulsive QMHRG model.These trajectories are beyond the reach of lattice QCD calculations and hold significance for lower collision energies of the RHIC BES program in the fixed target mode [77].
In the right panel of Fig. 6, we have shown the speed of sound calculated within the mean-field repulsive QMHRG model calculated along these constant s/n B trajectories.The speed of sound is defined as [29] Here, all the derivatives are taken along the lines of constant n S = 0, n B /n Q = 0.5 for different values of s/n B .We show the speed of sound for our QMHRG model and compare it with the latest lattice QCD results from Ref. [29].The c 2 s calculated within the mean-field QMHRG model decreases as the temperature is increased towards T pc for s/n B > 30.The QMHRG results with repulsive mean field for c 2 s are in agreement with the lattice QCD results for T < 150 MeV.Beyond this temperature, the lattice QCD results show an increasing trend with the temperature, which cannot be reproduced within a QMHRG model, even when the effect of the repulsive interactions are included.We also observe that in the temperature range 100 < T < 150 MeV the speed of sound decreases with increasing temperature both in the lattice QCD calculations and in the QMHRG model.This leads to a minimum of the speed of sound in the range of temperatures 140-150 MeV.From the inset of Fig. 6 (right) we also observe that the speed of sound has a maximum in the low temperature region for s/n B ≥ 30.
The behavior of the speed of sound for s/n B = 20 and s/n B = 15, currently inaccessible through lattice QCD techniques, is somewhat different.From Fig. 6 (right) we observe that for s/n B = 20 the speed of sound has a very shallow minimum around T = 145 MeV, while for s/n B = 15 it is always monotonically increasing with increasing temperatures.This implies that the maximum in c 2 s , evident from the plot at low temperatures T < 75 MeV and small baryon densities, should disappear when baryon densities are increased.We thus do not observe a softening of the equation of state near the chiral crossover transition at large values of baryon densities corresponding to s/n B ≈ 15 within this mean-field treatment of repulsive baryon interactions.

VI. SUMMARY
In this paper, we have studied the QCD equation of state and the chiral crossover at non-zero net baryon density using the QMHRG model with repulsive mean-field interaction.We have shown that this model can describe the lattice results on the fluctuations of net baryon number up to the eighth order in the vicinity of the chiral crossover temperature if the parameter characterizing the strength of the mean-field repulsive interactions is properly chosen.This is contrary to the usual QMHRG which can describe the higher-order net baryon number fluctuations only for considerably lower temperatures.
We extended our previous calculation of the chiral crossover temperature within the QMHRG model limited to the region of small baryon chemical potential [27] to much larger values of the baryon chemical potential.We showed that for µ B > 400 MeV there is a significant effect of the repulsive mean field on the value of the chiral crossover temperature.We found that the µ B dependence of the chiral crossover temperature can be very well parametrized by the following form: 4 with κ 2 = 0.0150(2) and κ 4 = 3.1(6) × 10 −5 .The value of the leading curvature coefficient, κ 2 , agrees well with the lattice QCD results within errors.The smallness of κ 4 is consistent with the upper bounds from lattice QCD, but for the first time we found that the corresponding value is significantly different from zero.
To realize the relevance of this study in the context of heavy-ion collision experiments, we have estimated the chiral crossover line imposing the strangeness-neutrality condition.The separation between the freeze-out curve and the pseudocritical line increases toward high density, which implies a longer-lived interacting hadronic phase at lower collision energies.Furthermore, we have used the repulsive mean-field QMHRG model to calculate the speed of sound along various isentropic ( constant s/n B ) lines pertinent to heavy-ion collision experiments, constrained by the strangeness neutrality and n Q /n B ratio.In the hadron phase, there is good agreement between our calculation and the latest available data from lattice QCD.We have also predicted the isentropic trajectory and the speed of sound within this repulsive mean-field model for s/n B = 15, where no lattice QCD data are available.We see that starting from this value of s/n B the previously seen characteristic softening of the QCD equation of state near the chiral crossover regions disappears.
We pointed out that the QMHRG model with a single repulsive mean field for all baryons only works for µ B < 750 MeV.For larger values of the baryon chemical potential, a more refined mean field approach is needed, with different mean fields for different baryon species.
Taking the first-order derivative with respect to βµ B , we get The above equation ensures that the net baryon number is given by the difference of the number of baryons and the number of antibaryons, and vanishes for µ B = 0. Using the above expressions, it is easy to extend the calculations to higher-order derivatives with respect to µ B /T .For completeness here we list all the derivatives up to the eighth order. 13 .
It is obvious from the above equations that the odd fluctuations of net baryon number vanish for µ B = 0. Using the above expressions we calculated the fluctuations of the net baryon density up to the eighth order at µ B = 0 and compared them to the available lattice results for different values of the parameter K.The value K = 33 GeV −2 leads to a good agreement with the lattice QCD results as one sees from Fig. 7.Here it is assumed that all baryons and baryon resonances contribute to the mean field.If we assume that only ground-state baryons contribute to the mean field then we do not get such a good agreement with the lattice QCD results for the same value of K.However, using a larger value of K, namely K = 100 GeV −2 a reasonably good agreement with the lattice QCD results can be obtained; cf.Fig. 7.

Appendix B: Calculation of the equation of state within the mean-field HRG model
The energy density is defined from the pressure as Taking the derivative of the number density with β we find represent the ideal QMHRG model, whereas the K = 33 data are calculated with mean-field repulsive interactions among all the baryons with a mean-field coefficient of K = 33 GeV −2 .Similarly, MF in OD represents the case where we have included repulsive mean-field interaction only among the baryons in octets and decuplets.Lattice QCD results for Nτ = 8 and Nτ = 12 are shown in red and green colors respectively, where the squares and circles correspond to data from the HotQCD [57] and Wuppertal-Budapest (WB) [2] Collaborations.The gray and yellow bands represent the fluctuation data in the continuum limit from the HotQCD [57] and WB [58] Collaborations.
Rearranging the above equation we get Putting the derivative terms together in Eq.B1, one can derive the final expression for the energy density, In the ideal limit of K = 0 the above formula reduces to the ideal HRG expression for energy density.In Fig. 8  FIG. 8.The pressure P/T 4 (upper row) and energy density ϵ/T 4 (lower row) normalized by T 4 calculated for ideal (red) and repulsive mean-field (blue) QMHRG models.Lattice results (black points) are from Ref. [1].The results are displayed for µB/T = 0, 1, and 2, in the left, middle, and right panels respectively.
show the pressure and energy density in units of T 4 in QMHRG and compare them with the lattice QCD results from Ref. [1].The results for both the ideal case and repulsive mean field case with K = 33 GeV −2 at µ B /T = 0, 1, and 2 are displayed.For all baryon chemical potentials, the QMHRG estimations align well with lattice results up to 150 MeV.At high µ B /T , the effect of repulsive mean field is more pronounced due to the increased thermal abundance of baryons at high density; cf.Fig. 8.

33 FIG. 1 .
FIG.1.The second (top left), fourth (top right), sixth (bottom left), and eighth (bottom right) order baryon number fluctuations results compared between ideal QMHRG (red line) and mean-field QMHRG (blue line) models and lattice QCD results with physical quark masses.Lattice results for Nτ = 12 and Nτ = 8 are shown as red and green points respectively with the squares representing HotQCD[57] results and the circles corresponding to Wuppertal-Budapest (WB) results[2].The gray and yellow bands represent the fluctuation data in the continuum limit from HotQCD[57] and WB[58] Collaborations.

FIG. 2 .
FIG. 2. Left:The µB dependence of the renormalized chiral condensate at different temperatures T = 25-150 MeV.The bands are results from the QMHRG model with repulsive mean-field interaction and the lines with points are the results from the ideal QMHRG model.Right: The variation of the ∆ l R with net-baryon number density normalized by the nuclear saturation density, n0 = 0.16 fm −3 , for the same range of temperatures.The widths of the bands represent the uncertainty of the chiral condensate due to the uncertainties in the quark mass derivatives of the hadron masses; see text for details.

7 T
FIG. 3. Left:The pseudocritical line calculated from the QMHRG model is compared with the lattice QCD results from the HotQCD Collaboration (light blue) from Ref.[66].The ideal gas (K = 0) results are shown as a line and results including the mean-field repulsion (K = 33) are shown as open symbols.These are compared with the freeze-out line (magenta) using the parametrization in Ref.[67].The gray line represents the fit performed by us to obtain the pseudocritical line; see text for details.Right: The net-baryon density in units of n0 = 0.16 fm −3 is shown as function of µB along the pseudocritical line.The horizontal gray band represents the variance in the nB/n0, which is due to the variation of Tpc arising from the uncertainties in the sigma terms.We also show the variation of the net baryon density along the freeze-out line.The color scheme is the same as in the left panel.

FIG. 4 .
FIG.4.Left: Relative contributions to the energy density from mesons (blue) and baryons (red) along the pseudo-critical line compared to the total energy density (gray) calculated in the mean-field QMHRG model.Right: The relative number densities of various baryon species shown along the pseudo-critical line.

FIG. 5 .
FIG. 5.The chiral crossover temperature at finite density µB normalized by its corresponding value at µB = 0 is shown as function of µB (left) and as a function of nB (right) for µS = 0 and under strangeness neutrality condition, nS = 0. Ideal QMHRG model results are shown in blue.Coral and green represents mean-field QMHRG results for nS = 0 and µS = 0 respectively. 0

FIG. 6 .
FIG. 6. Left: The isentropes for s/nB = 15, 20, 30, 50, 100, and 400 respectively are shown in the T − µB plane.Right: The speed of sound as a function of the temperature is shown along these isentropes.The bands are the lattice QCD data from Ref. [29] and the points are the QMHRG model calculations with repulsive mean field.The isentropes and the speed of sound in the low temperature region 60 < T < 150 MeV are shown in the insets of the left and the right figures respectively.The black line represents the pseudocritical line evaluated within our mean-field QMHRG model.
FIG.7.Top row: second-order (left) and fourth-order (right) baryon number fluctuations; bottom row: sixth-order (left) and eighth-order (right) baryon number fluctuations, shown as a function of temperature and µB = 0.The K = 0 results represent the ideal QMHRG model, whereas the K = 33 data are calculated with mean-field repulsive interactions among all the baryons with a mean-field coefficient of K = 33 GeV −2 .Similarly, MF in OD represents the case where we have included repulsive mean-field interaction only among the baryons in octets and decuplets.Lattice QCD results for Nτ = 8 and Nτ = 12 are shown in red and green colors respectively, where the squares and circles correspond to data from the HotQCD[57] and Wuppertal-Budapest (WB)[2] Collaborations.The gray and yellow bands represent the fluctuation data in the continuum limit from the HotQCD[57] and WB[58] Collaborations.