Chiral condensate from a hadron resonance gas model

In this work we address the question of how well the chiral crossover transition can be understood in terms of a noninteracting hadron resonance gas model. Using the latest results on the variation of hadron masses as a function of the pion mass from lattice quantum chromodynamics, we study the temperature dependence of the renormalized chiral condensate in 2+1 flavor QCD. Furthermore, we suggest a better criterion to estimate of the pseudocritical temperature, which gives $T_c = 161.2 \pm 1.7$ MeV, which is much improved compared to all the earlier results within the hadron resonance gas model or chiral perturbation theory. For the curvature of the pseudocritical line we find $\kappa_2 = 0.0203(7)$, which is in very good agreement with continuum extrapolated lattice results.


I. INTRODUCTION
The breaking of chiral symmetry is one of the key features of quantum chromodynamics (QCD) underlying much of our understanding of hadron physics. The breaking of this symmetry is signaled by the nonzero value of the quark (chiral) condensate, ψ ψ = 0. It is well known that at high temperature the chiral symmetry is restored and this restoration happens as an analytic crossover at a temperature T c = 156.5 ± 1.5 MeV [1]. The chiral condensate rapidly decreases at and above T c . Below the crossover temperature, it is expected that the thermodynamics of QCD matter can be understood in terms of hadronic degrees of freedom. There is a wide belief that the description of QCD matter as a gas of hadrons is valid up to the crossover temperature, T c . This belief is the basis of our current understanding of the freeze-out condition in heavy-ion collisions [2][3][4][5][6][7][8]. These considerations also implicitly assume that for T < T c the decrease in the chiral condensate can be understood in terms of a hadron gas. Therefore, an important question to clarify is to what extent we can understand the decrease of the chiral condensate for T < T c in terms of hadronic degrees of freedom. This question is also important for our understanding of the QCD phase diagram at large values of baryon density. Recent studies suggest that lattice QCD calculations based on Taylor expansion in chemical potential or analytic continuation may only work for baryon chemical potential µ B < 2.5 T , [9][10][11][12] though another study claims that lattice QCD calculations may be useful for baryon chemical potential as high as 3.5 T [13]. In any case, the study of the temperature dependence of the chiral condensate at large values of µ B is of great interest and can be performed within the framework of a hadron gas.
The hadron resonance gas (HRG) model turned out to be very successful in describing QCD thermodynamics below the crossover temperature. The main idea of the HRG model is that the interacting gas of hadrons can be replaced by a noninteracting gas of hadron and hadron resonances, i.e., all interactions between hadrons can be approximated by the creation of additional hadronic resonances. This approach can be justified in terms of relativistic virial expansion, where interactions between hadrons are expressed in terms of phase shifts [14]. Using experimentally known phase shifts within this approach, it has been shown that the nonresonant part of the phase shifts largely cancel out in the QCD pressure, and indeed the interacting part of the pressure can be well described by the contribution of resonances treated as stable hadrons [15].
There has been an extensive effort to test the validity of the HRG model by comparison to the lattice QCD results. Early works in this direction were presented in Refs. [16][17][18][19] and confirmed the validity of HRG approach. Recent continuum extrapolated lattice QCD results on the equation of state show that HRG works reasonably well below T c if all known resonances from the Particle Data Group (PDG) are included [20][21][22]. A precision lattice study of the second derivatives of the pressure with respect to the chemical potentials again confirmed the validity of the HRG model [23]. However for the quantities involving strangeness, for example, strangeness-baryon-number correlations, including additional resonances not listed by the PDG but predicted by QCD (or QCD inspired models) turned out to be necessary [23,24]. The need for additional resonances to describe the baryon-strangeness correlations was also confirmed using relativistic virial expansion with state-of-the-art partial wave analysis of the phase shifts [25]. Furthermore, the baryon-charm correlation calculated on the lattice requires taking into account resonances not listed by the PDG [26]. For higher-order derivatives of the pressure with respect to the chemical potential the HRG model does not work as well. It has been argued that repulsive baryon-baryon interactions play an important role in this case and the contribution of these repulsive interactions increases with the order of the derivative. Including these interactions in some phenomenological way improves the agreement between the lattice results and HRG model [27][28][29][30].
The temperature dependence of the chiral condensate has also been studied in the HRG model [20,[31][32][33][34] and within the chiral perturbation theory [35,36]. These calculations show a significant change in the chiral condensate in the vicinity of the crossover temperature. The pseudocritical temperature extracted from the point where the chiral condensate of a pion gas within the next-to-next-to-leading order (NNLO) chiral perturbation theory falls to zero is about 250 MeV, which is lowered to about 190 MeV when additional hadrons are included [35]. However, the uncertainty in the HRG calculations of the chiral condensate is difficult to quantify. This is due to the fact that, in order to obtain the chiral condensate in the HRG model one has to know the precise dependence of the hadron masses on the light quark mass. Except for the few lowest-lying hadron states this dependence is not well known. The main goal of this paper is to revisit the temperature dependence of the chiral condensate in the HRG model and quantify the uncertainties related to the poorly known quark mass dependence of hadron masses. We will show that with systematic improvements in the procedure we can make a much more precise estimate of the T c within the HRG model.
The paper is organized as follows: In the next section, we discuss the general formalism and the relevant observables. In particular, we focus on two renormalized definitions of subtracted chiral condensate, which we calculate within the HRG model. We show how much contribution to these observables comes from different ground state mesons and baryons as well as their higher excited states in Sec. III. We find that the dominant contribution to the renormalized chiral contribution at small baryon densities comes from the mesons, in particular, the ground state pseudoscalar and vector mesons. We have extracted the pseudocritical temperature T c = 161.2 ± 1.7 MeV and a curvature of the pseudocritical line κ 2 = 0.0203 (7) corresponding to the chiral crossover transition which not only is in very good agreement with the latest lattice QCD result but can be obtained with considerably less computational effort. Furthermore, our estimates of the pseudocritical temperature are much improved compared to the corresponding calculations within the chiral perturbation theory. Though pseudoscalars dominantly contribute to the chiral condensate, extending this formalism to the chiral limit has its own limitations. We discuss this in detail, showing how we could improve from the previous works in Ref. [35]. We conclude by discussing the important implications of our work, in particular, the possible extension of these calculations to understand the QCD phase diagram at finite baryon densities where lattice QCD calculations are severely limited due to the infamous sign problem; see e.g. Ref. [37] for a recent review.

A. Renormalized definitions of chiral condensate
The total pressure due to a noninteracting ensemble of several species of hadron labeled by α, enclosed within a volume V and at an ambient temperature T , is given as, Here, +(−) corresponds to baryons (mesons) respectively. The degeneracy factor and the energy of each particle α are denoted by g α and E α , respectively. Furthermore, E α = p 2 + M 2 α with p, M α being the magnitude of the particle momentum and its mass respectively. We consider here a grand canonical ensemble with each species of hadron denoted by a chemical potential µ α = B α µ B + Q α µ Q + S α µ S , which depends on the quantum numbers B α , Q α , and S α denoting the conserved quantum numbers baryon number, electric charge, and strangeness respectively. Within this model the interactions between hadrons are taken care of by inclusion of hadron resonances as free stable particles [14,15], i.e., the sum in Eq. [1] also includes resonances. The light quark condensate, also called the chiral condensate, at nonzero temperature is defined as where m l is the light quark mass. We consider the case of two degenerate light quarks, m u = m d = m l . The zero temperature light quark condensate ψ ψ l,0 can be thought of as the derivative of the vacuum pressure with respect to m l , and it contains all the multiplicative as well as additive divergences (for m l = 0 case). To have a finite observable which is also renormalization group invariant we consider the following combination where m s is the strange quark mass. A definition of the renormalized chiral condensate that has been discussed in the literature [20] constructed out of the quantity defined in the Eq. [3] but with a different normalization factor, is Another, possibly a more natural way to define a dimensionless chiral condensate at nonzero quark mass [38] is where the parameter r 1 is derived from the static quark potential [39] and d = r 4 1 m s (lim m l →0 ψ ψ l,0 ) R . This definition takes advantage of the fact that, in the chiral limit, the light quark condensate has only a multiplicative renormalization. The superscript R denotes the renormalized quantity. Taking into account 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 [40], as well as r 1 = 0.3106 fm [41], we obtain the value of d = 0.022791.
In our calculations we have included all hadron resonances predicted from the quark models [42,43] in addition to those mentioned in the Particle Data Group 2016 lists, as tabulated in Ref. [44]. The HRG model constructed out of this augmented hadron list is called the QMHRG, and it turned out to be successful in describing the temperature dependence of cumulants of different conserved charges obtained from lattice QCD for T 150 MeV [44][45][46][47]. Recently, the HotQCD Collaboration constructed a HRG model with a similar list of states for their recent comparative study [23] on the second-order fluctuations of net baryon number, strangeness, and electric charge and their correlations. We have also verified that the HRG model used in this work can reproduce the results from Ref. [23] with a very good precision.
It is clear that one of the most important ingredients needed for the calculation of the chiral condensate is the dependence of the hadron masses, including those of the resonances, on the light quark mass m l . Although this variation is well understood only for a few primary hadrons, more specifically the ground states, we will discuss in this work the detailed systematics for the resonance states. In the following subsections, we will separately discuss the contributions of the pseudo-Goldstone modes, i.e., those of pions and kaons, and of all other heavier hadrons and resonances to the chiral condensate.

B. Contributions from pions and kaons
We first consider the contribution of pions and kaons to the chiral condensate. The contribution of pions and kaons to m s (∂P/∂m l ) can be written as Here n B is the Bose-Einstein distribution. Extensive lattice studies of the quark mass dependence of the pseudoscalar meson masses found excellent agreement with SU(2) χPT, whereas the agreement with the SU(3) χPT is not so great. For more comprehensive reviews, we refer to Refs. [48][49][50][51][52][53]. Evidence from lattice studies suggests that an extended version of the SU(2) χPT, considering the strange quarks as heavier degrees of freedom [48,54] can be useful in understanding the quark mass dependence of the kaon mass and decay constant. We will henceforth refer to the two-flavor χPT results for the pion mass from Ref. [55] and its extended version from Ref. [48] for the kaon mass. The dependence of the pion and kaon masses on light and strange quark masses m l , m s described in terms of SU(2) low energy constants F, B,l 3 , F π , are given as The χPT relations for the pseudoscalar mesons are for the two-degenerate flavor case m u = m d = m l , where most lattice QCD calculations exist in the continuum limit. The bag constant B K (m s ) and the NLO low energy constants λ 1,2 (m s ) in the kaon mass can be expressed in terms of the SU(3) low energy constants L 4 , L 5 , L 6 , L 8 , B 0 , F 0 as [48], We have set the scale appearing in the above formulas Λ χ = 770 MeV, i.e., at the ρ meson mass. We have taken the values of the different χPT low energy constants from the latest FLAG review [40]. In particular we use the value Σ 1/3 = 272(5) MeV for 2+1 flavor QCD. The other SU(2) low energy constant F is obtained from the relation F π /F = 1.062 (7), where F π = 92.2(1) MeV. The various SU(3) low energy constants that appear in the the expression of the kaon mass squared for 2 + 1 flavor QCD are Σ K . We also took the latest value of m l = 3.381(40) MeV from the same review. We then calculated the derivative of the pseudoscalar meson masses in Eq. (7) with respect to m l . These derivatives were evaluated at the physical point as well as for the vanishing light quark mass. The latter are needed for the subsequent discussions in Sec. II D. Lattice QCD studies show that, at the physical point, m l ∂M π /∂m l = M π /2 to a very good approximation [56]. This is equivalent to the assumption that close to the physical point, M 2 π = 2Bm l , and hence the error on the m l derivative of M π can be estimated by propagating the errors on B and m l as quoted in the FLAG review.

C. Contributions from heavier hadrons
To obtain the contribution from heavier hadrons to the chiral condensate we can write with n α being the Bose-Einstein distribution for α ≡ mesons and the Fermi-Dirac distribution for α ≡ baryons. Here we introduce the so-called σ terms: For the ground state baryons, the σ terms are known from the fits of the low energy data of masses within the chiral effective theory and the lattice. These are summarized in Table I. For other hadrons we will use the lattice results on the pion mass dependence of the hadron mass for constant strange quark mass set to its physical value, to evaluate M 2 π ∂Mα ∂M 2 π . A summary of these evaluations is given in Table II and the corresponding details are given in Appendix A. For the ground state pseudoscalar meson η and vector mesons ρ(770) and K * (892), we have collected the lattice data for the variation of these masses as a function of the pion mass from Refs. [57][58][59] and have extracted the corresponding σ terms, details of which are discussed in Appendix A.
Apart from the vector and pseudoscalar mesons, we have categorized the other meson-resonances in three groups according to their quantum numbers as shown in Table II. Contribution to the chiral condensate from mesons with the same quantum numbers as ρ (like a, b, and higher excited states of ρ) have been included considering their σ values similar to ρ(770), whereas the contribution from higher excited states of K and K * have been included considering their σ terms to be the same as K * (892). The σ terms for the isoscalar mesons and their excited states in our calculations are kept identical to those of the corresponding ground states mesons (η, ω, φ). We have also included the contribution of the η to the renormalized chiral condensate separately, considering its σ term, details of which are again mentioned in Appendix A. The contribution of the scalar meson f 0 (500) is not considered, as it was observed earlier that the pole due to attractive interactions between them in the I = 0 sector are exactly canceled by the repulsive interactions [60]. I. Sigma terms of baryon octet and decuplet ground states used in this work taken from Ref. [61]. (1)(2) 25(1)(1) 15(1)(1) 29(9)(3) 18(6)(2) 10(3)(2) 5(1)(1) Following the same procedure as for the mesons, we have carefully categorized the higher excited states for the baryons as well. For the nucleon and ∆ resonances we have checked that their σ terms, extracted from the lattice data [62] are within 1σ with the corresponding ground state values. Excited strange baryons with strangeness content S = 1, 2, 3 have been included by setting their corresponding σ term same as the ground state baryons Σ * , Ξ * and Ω − respectively. We end this section by summarizing the σ terms used in our calculations for the ground states and higher excitations in the various quantum number channels for baryons and mesons in Table II. TABLE II. A summary of the σ terms for various meson and baryon channels used in this work.

Mesons
Meson NLO χPT Refs. [48,56] Octet (N, Λ, Σ, Ξ) Table I Similar to ρ(770) Excited states of N and ∆ Similar to σN and σ∆ Resonances |S| = 0, (K, K * ) Similar to K * Excited states labeled as Λ, Σ, Σ * Similar to Σ * Resonances |S| = 0, |I| = 0, (φ, η, f ) Similar to ground states (η, η , ω, φ) Excited states labeled as Ξ, Ξ * , Ω Similar to Ξ * and Ω D. The condensate in the chiral limit In QCD with two light quark flavors, the crossover transition goes over to a real phase transition in the chiral limit. This happens because the critical fluctuations in the free energy start dominating over the regular (analytic in the light quark mass) terms towards the chiral limit. Within a HRG model, the information about critical fluctuations may be entirely missing in the absence of additional inputs about the thermal widths, particularly of the pseudo-Goldstone partners π, f 0 . Unlike in QCD, the ψ ψ T within the HRG model will not be zero at the critical temperature T 0 c . In order to calculate the chiral condensate we need to calculate the derivatives of the square of the hadron masses with respect to the light quark mass, but now in the limit m l → 0. These derivatives for the pion and kaon mass squared can be calculated using Eq. (7), which for the m l → 0 case read: From the latest FLAG review [40], the value of B = Σ/F 2 = 2.76 GeV, where F and Σ are the pion decay constant and light quark condensate respectively measured in the chiral limit of SU (2) χPT. Since pions are the Goldstone modes of QCD with two massless flavors, they will have zero mass. The kaon and η meson masses decrease only slightly in the chiral limit as their masses are determined primarily by the strange quark mass. We then calculate the chiral condensate defined in Eq. (2) in the limit m l → 0 and normalize by its zero temperature value. Except for the pions and kaons the remaining hadrons and resonances have been included in the definition of renormalized chiral condensate by rearranging Eq. (3) such that it could be generalized easily to the chiral limit. This is done through the relation Here we have used the same definition of the σ term as introduced for the baryons in Eq. (10) and hence we had to normalize the right-hand side of Eq. (12) by the square of the physical pion mass. The derivative term lim m l →0 ∂M 2 π /∂m l = 2B and the sum in Eq. (12) is over all the hadron species except for the pions and kaons. In the Eq. (12), we have used the masses of the ground states of the baryon octet and decuplet which were extracted in the chiral limit from Ref. [61]. For the ρ meson, the fit to its mass data as a function of M 2 π , as shown in Appendix A, yielded M ρ = 690(18) MeV in the chiral limit which is about 10% lower than its physical mass. However, we could not find any data for the masses of most other resonances in the chiral limit, but instead took their physical masses. To estimate the systematic uncertainty due to this approximation, we varied the mass of these resonances by 10% but found a ≈ 1% change in the value of the temperature at which the chiral condensate vanishes. We thus anticipate that the effect of the mass modification of these resonances will be tiny in the chiral limit.

E. Estimating the errors in our results
The major source of error in the calculation of the renormalized chiral condensate comes from the errors in estimating the σ terms. In this work, we have improved upon the errors considerably by estimating the σ terms from the lattice QCD data for the ρ, K * and the η, η mesons and the excited states of pions. For the baryon sector, the σ terms for the octet and decuplet ground states have been extracted to a very good precision [61]. However, large uncertainty exists about the values of the σ terms of the excited baryon states. We show in Appendix A how extraction of the σ term of the excited nucleon state is prone to larger errors due to uncertainty in calculating the masses of these states. Although we have assumed that the higher mass resonances have a similar value of σ term as their ground states hadrons, various studies hint that their σ terms may deviate from those of their corresponding ground states [63,64]. To reliably account for this uncertainty, we have taken the relative errors in the σ terms of all the excited states of the baryons to be 50%. We have found that this variation does not result in a significant change in the renormalized chiral condensate. Due to their relatively heavier masses, the condensate is mainly dominated by the ground state hadrons in the temperature region till 155 MeV. Although the resonances have a significant cumulative contribution to the condensate, the uncertainties in their σ terms effect the estimates of T c by less than 1%.
The largest sources of error are still due to the uncertainties in the derivatives of the pion and kaon mass with respect to m l , which we have estimated within the χPT. For example, the relative error for the m l derivative of the pion mass is almost 6%, which arises mainly due to the uncertainty in the low energy constant B. For kaons, the error is even higher ≈ 23%, when we use Eq. (7). This is because the low energy constants in SU (3) χPT are not very precisely determined, and the majority, ≈ 18%, of the error arises from the uncertainty in Σ 1/3 0 MeV and F 0 that appear explicitly in the m l derivative of kaon mass. The errors due to these two SU(3) χPT low energy constants also contribute to the error that comes from the term λ 1 (m s ) + λ 2 (m s ), which accounts for the remaining 5% contribution to the error on the m l derivative of the kaon mass squared. As the ground state pseudoscalar mesons dominantly contribute in the chiral condensate at T ≈ 156 MeV, these errors also dominate the total error of the renormalized chiral condensate and hence the estimation of T c .

III. RESULTS
A. Comparing the HRG model results for the renormalized chiral condensate with lattice QCD We compare our HRG model calculations with the continuum extrapolated lattice QCD results for ψ ψ R , from Ref. [20]. The lattice results were calculated for 2+1 flavor QCD with staggered discretization for quarks with one-link stout improvement and physical quark masses such that m s /m l = 28.15. To compare our HRG model calculation to the lattice data we have accordingly set the ratio of quark masses to the same value in the model as well.
We show the comparison of the renormalized chiral condensate as defined in Eq. (4) calculated in lattice QCD [20] and the HRG model in Fig. 1 (left). This quantity shows a smooth analytic rise as a function of temperature, as is evident from the continuum extrapolated lattice results from Ref. [20], gradually saturating at T ≈ 200 MeV. Since there is no real phase transition for physical quark masses, we do not observe any sharp jump at the chiral crossover transition temperature T c . Our results calculated within the noninteracting HRG model agree well with the lattice results below T = 140 MeV. Above this temperature, the HRG model results show a rapid rise with temperature, a trend similarly visible within the continuum extrapolated lattice data, however, the HRG estimates are lower than the lattice QCD results. As will be discussed later, while for T > 140 MeV pion and kaons provide the dominant contribution to the renormalized chiral condensate, the relative contributions of the vector mesons and baryons start to increase near the chiral crossover, thus influencing the location of the corresponding pseudocritical temperature. The ψ ψ R calculated from HRG shows a diverging behavior at high temperatures T ≈ 170 MeV, which is similar to the singularity in the Hagedorn spectrum [65], due to the density of states increasing exponentially with temperature. This is primarily driven due to a large number of excited states of different spin-parity channels in the baryon sector. There is no sign of saturation of ψ ψ R in the HRG model at high temperatures unlike on the lattice.
We have also compared our results with the calculations from NLO chiral perturbation theory in the two-loop approximation in the presence of a pion heat-bath taken from Refs. [66,67]. The χPT results agree with our HRG π ψ ψ l,T − ψ ψ l,0 calculated within QMHRG is compared to the continuum extrapolated lattice data from Ref. [20]. Right: The relative contribution to the subtracted light quark condensate, −ms ψ ψ l,T − ψ ψ l,0 /T 4 , due to different meson and baryon channels is shown and the resultant total contribution within QMHRG (orange band) is compared to the lattice data from Ref. [20].
results at T 100 MeV, beyond which the latter increases more strongly with temperature leading to a larger difference between the two.
A similar comparison of lattice QCD data with the HRG model predictions was shown in Ref. [20]. However, the uncertainty in the HRG predictions have been reduced by an order of magnitude in our results. This is possible now because of a more precise determination of σ terms for the ground state baryons [61]. Furthermore, with a significant improvement in the calculation of ρ meson mass on the lattice as a function of the pion mass [68], we could extract the corresponding σ term more precisely. Similar improvement of the σ terms of the pseudoscalar isosinglet meson excited states has also led to a more accurate HRG model estimate of the renormalized chiral condensate.
To check the robustness of our predictions we have also calculated a different definition of the renormalized chiral condensate −m s ψ ψ l,T − ψ ψ l,0 /T 4 , shown in the right panel of Fig. 1. The temperature dependence of this quantity is similar to that of the trace anomaly calculated in Ref. [22], and it is also useful in order to understand the relative importance of different hadron species contributing to the chiral condensate. As one can observe from Fig. 1, the heavier hadrons become important for T > 150 MeV. These are the states which are responsible for the increase of this observable as a function of temperature, since the contribution from the pseudoscalar states monotonically decreases beyond T ≈ 100 MeV. We can also see that the contributions of the higher-lying baryon and meson resonances are important for T > 150 MeV and become comparable in magnitude to the individual contributions due to octet and decuplet baryons and vector mesons. The quark mass dependence of these excited hadron resonances are not very well constrained by lattice QCD calculations and the errors on the corresponding σ terms are comparatively large. Hence the error band for this observable also increases with temperature. Based on the discussion in Appendix A we assign a generous 50% relative error for the σ terms of these higher-lying meson and baryon resonances. In Appendix B we further scrutinize the relative contribution of these resonances to the chiral observables.
In the left panel of Fig. 2, we show the temperature dependence of another definition of the renormalized chiral condensate, ∆ l R as calculated from the HRG model and compared with the lattice data. This quantity is related to the above definition of the chiral condensate through ∆ l R = d−(m s /m l )(r 1 m π ) 4 ψ ψ R . The study of ∆ l R is motivated by our desire to extract the pseudocritical temperature corresponding to the chiral crossover transition, within our HRG model calculations. A comparison with the lattice QCD data for ∆ l R leads to similar conclusions as earlier, namely for T ≤ 140 MeV the difference between the lattice and HRG model result is small, while for 140 < T < 160 MeV the HRG model results for ∆ l R drop slower than the lattice QCD results. At still higher temperatures the lattice results for ∆ l R flattens off signaling a change in the degrees of freedom, while the HRG model data keep decreasing. Let us now discuss the technique of determination of the pseudocritical temperature. In lattice calculations, the pseudocritical temperature is defined as an inflection point of the renormalized chiral condensate or the peak of the chiral susceptibility as a function of temperature. Due to the breaking of the exact chiral symmetry due to a finite quark mass, the value of the pseudo-critical temperature will depend on the quantity used to define it, including FIG. 2. Left: ∆ l R , as defined in Ref. [38], ∆ l R = d+msr 4 1 ψ ψ l,T − ψ ψ l,0 is compared to the lattice data taken from Ref. [20].
We have used d = 0.022791 and (r1mπ) 4 = 2.23 × 10 −3 following [38]. The dotted line denotes the half value of ∆ l R compared to its magnitude at the lowest temperatures, which is used as a criterion to determine Tc. Right: ∆ l R has been plotted for different values of baryon chemical potentials µB/T ≤ 2. In both the figures an estimate of Tc is obtained from where the magnitude of ∆ l R falls to half its zero-temperature value.
the normalization of the corresponding quantity [69]. In a recent lattice QCD study with highly improved staggered quark discretization, it was found that the spread of the pseudocritical temperature obtained from chiral observables is surprisingly small, ≈ 1.5 MeV. A somewhat larger spread was obtained in another study with stout-improved staggered fermions [20]. Here the spread was about 8 MeV when results from the chiral susceptibility were taken into account. The main difference is due to a different normalization of the chiral susceptibility suggested in Ref. [20].
The HRG model results on the renormalized chiral susceptibility, on the other hand, do not have an inflection point. This is obvious from Fig. 2. Therefore, one cannot define a pseudocritical temperature from the HRG model results alone. However, as mentioned in the Introduction, the breaking of the chiral symmetry plays a crucial role in hadron physics. If at some temperature the renormalized chiral condensate becomes small compared to its vacuum value the hadronic matter should gradually melt to a new state. We do not know a priori how small the chiral condensate should be for this to happen because of the inherent nonperturbative nature. Lattice QCD calculations show that around the pseudocritical temperature the value of the chiral condensate ∆ R l is around half of its vacuum value [38,70]. Motivated from this insight from QCD, we therefore set the same criterion to estimate T c within the HRG model calculations. With this procedure we estimate the pseudocritical temperature T c = 161.2 ± 1.7 MeV, which is only 2σ away from the recent high precision lattice QCD results on T c in the continuum limit [1,71]. This is not completely surprising as the temperature dependence of the renormalized chiral condensate in the HRG model and in lattice QCD is similar, and the obtained estimate of T c is not completely independent of the lattice determination. However, this analysis represents a marked improvement over earlier estimates T c obtained from the HRG model [31][32][33] or from χPT [35,36], where the central values are more than 10σ away from the latest lattice result.
It is more interesting to extend the above analysis to nonzero baryon chemical potential and see how the transition temperature defined within the HRG model depends on baryon density. For that purpose, we calculate the temperature dependence of ∆ l R for different values the baryon chemical potential 0 ≤ µ B /T ≤ 2. The corresponding results are shown in Fig. 2 (right). At finite netbaryon densities the contribution from baryon states and resonances is expected to become more decisive, which results in the ∆ l R decreasing faster with increasing temperature compared to the µ B = 0 case, as can be seen from Fig. 2 (right). In turn, this leads to a decrease in the pseudocritical temperature with increasing baryon density as expected. However, with the increase in netbaryon density, repulsive interactions between baryons will be more probable and hence this simple noninteracting HRG model will cease to effectively describe the QCD medium. The effect of repulsive interactions on the pseudocritical temperature will be discussed in detail in our forthcoming publication [72].
B. Curvature of the chiral crossover line from the HRG model As the next step we will investigate the dependence of the chiral crossover line as function of the baryon chemical potential µ B within the HRG framework. As discussed in the previous sub-section we expect that the chiral crossover temperature will decrease with increasing µ B . We define the pseudocritical temperature at nonzero µ B again, as the temperature where the renormalized chiral condensate drops by factor of 2, since we do not expect qualitative change in the chiral crossover unless µ B is too large. At moderately large values of µ B the pseudocritical temperature can be written as The lattice calculations of the curvature κ 2 and its higher order corrections κ 4 are quite challenging [73]. For several years κ 2 measured using different lattice discretizations and systematics, did not agree with each other [74,75], which was recently understood to be due to the fact that the continuum limit needs to be taken carefully [76]. Presently all continuum extrapolated lattice measurements of κ 2 from different groups using Taylor expansion [1] or from imaginary chemical potential techniques [71,76] agree quite well with each other. However, the extraction of κ 4 remains a challenge since it involves delicate cancellation between noisy operators, due to which the signal-to-noise ratio in this observable is extremely low [1,71].
In this context, we will study how well we can calculate the curvature terms κ 2 and κ 4 within the HRG model. For that, we first extract the pseudocritical temperature as a function of the baryon chemical potential µ B . We have then fitted the T -µ B crossover curve with the ansatz in Eq. (13) for values of µ B /T < 1. The resultant fits yielded a κ 2 = 0.0203±0.0007 within the HRG model which also is in a good agreement with the κ 2 = 0.0150(35) extracted from the continuum extrapolated (subtracted) chiral condensate calculated in lattice QCD [73]. However, our extracted value of κ 4 = −3(2) × 10 −4 is again quite noisy, consistent with the findings from lattice QCD [1,71]. At small baryon densities, the curvature of the pseudocritical line is mainly determined by κ 2 , which makes the extraction of κ 4 so difficult. To extract higher-order curvature coefficients one needs to extend our analysis to µ B /T > 2, where it is anticipated that this simple noninteracting HRG approximation may break down due to repulsive baryon interactions. We summarize our main findings on the curvature coefficients in the left panel of Fig. 3. The ratio of the chiral condensate at a finite temperature to its zero temperature value, compared between our calculation within QMHRG (orange band) and three-loop χPT results up to O(T 6 ) from Ref. [35] which are augmented by the contribution of heavier hadrons with the spectrum that has been used in our work (shown as the gray band).
C. Temperature dependence of the chiral condensate in the massless limit The chiral transition in the limit when m l → 0 becomes a second-order phase transition with O(4) or a U (2) L ×U (2) R critical behavior depending on whether or not the anomalous U A (1) part of the chiral symmetry remains broken substantially [77][78][79][80][81][82][83]. A high precision lattice QCD calculation of the chiral phase transition temperature with highly improved staggered quarks gives a T 0 c = 132 +3 −6 MeV [84]. More recently there is a lattice QCD calculation using Wilson fermions [85] reporting a T 0 c = 134 +6 −4 MeV which in good agreement with the earlier estimate. The T 0 c estimated in 2 + 1 flavor QCD within the functional renormalization group method yields a slightly higher value of T 0 c ≈ 142 MeV [86] but is reasonably close to the lattice calculation. Even though pions will contribute dominantly to the condensate in the chiral limit, extracting the T 0 c is not possible within the conventional QMHRG models which do not have the information of the precise critical universality class [87,88]. It was recently shown within linear sigma models with O(4) universality or in U (3) χPT, that it is necessary to include finite temperature corrections to the f 0 (500) spectral function to accurately account for the temperature dependence of the scalar susceptibility near the transition [89].
The earliest effort towards extracting the chiral transition temperature from ψ ψ T = 0 within a purely hadronic model, the three-loop χPT at finite temperature, gave a T 0 c ≈ 190 MeV [35]. This large disagreement with the latest lattice estimates of T 0 c can be understood from the fact that the chiral condensate due to an interacting pion gas calculated at three loops i.e., O(T 6 ), is reduced to half by T = 150 MeV [35] thus limiting the validity of the calculations beyond this temperature. Nonetheless the authors in Ref. [35] extended their calculations to higher temperatures, also including the effects of a dilute gas of heavier hadrons that the pions increasingly encounter. The presence of heavier hadrons reduces the chiral condensate to zero at a lower temperature, ≈ 170 MeV, compared to a pure self-interacting pion gas. Since in this work we have a substantially improved the procedure for inclusion of heavier hadrons and resonances, we included the effect of these states on the chiral condensate in addition to the three-loop pion gas contribution. The corresponding ratio of chiral condensates ψ ψ T / ψ ψ 0 is shown as a gray band in Fig. 3. We observe that the chiral condensate reduces to zero at a temperature 162 MeV, which is about 8 MeV lower compared to the earlier observed T 0 c in Ref. [35]. If we extend our earlier calculation of the chiral condensate within QMHRG model but now in the chiral limit, following the procedure outlined in Sec. II D, we observe it to fall to zero at a comparatively higher temperature ≈ 168 MeV. This can be understood from the fact that only the attractive resonant interactions between pions are included in our QMHRG model. With this analysis, we conclude that extending the existing χPT calculation beyond three loops and carefully including the contributions of heavier hadrons may reduce the estimated T 0 c , but this estimate will be still far from the lattice QCD result.

IV. SUMMARY AND OUTLOOK
Motivated from several pieces of evidence of a reasonably good representation of QCD thermodynamics at temperatures T 0.9 T c by a noninteracting gas of hadrons and resonances, we revisit its validity for chiral observables in this study. We provide an updated analysis of the renormalized chiral condensate mainly using extensive lattice data available for the masses of several hadrons and resonances as a function of the pion mass, and some recent updates of the so-called σ terms for ground state baryons. Furthermore, we include the effects of additional resonances, mainly strange baryon resonances, which are not yet measured experimentally but predicted from quark model and lattice studies. With the systematic inclusion of these effects we observe that a noninteracting HRG description for the renormalized chiral condensate is adequate to describe QCD data from first-principles lattice studies, up to T 140 MeV. However, using the fact that the lattice studies observe a drop in the value of renormalized chiral condensate to half its zero temperature value at the chiral crossover transition, we made an estimate of the pseudocritical temperature T c within this simple noninteracting HRG model. The extracted T c = 161.2 ± 1.7 MeV is in good agreement with the latest continuum extrapolated lattice QCD result, although this determination requires insights from lattice studies.
Encouraged by this effort we further estimate the curvature of the critical line at vanishingly small baryon density, yielding κ 2 = 0.0203 (7), which is in excellent agreement with the lattice results, which is highly nontrivial. Furthermore, with the more precise knowledge of the hadron σ terms, we could achieve about 8 MeV reduction in the earlier estimation of the chiral transition temperature in a system comprising a self-interacting pion gas at three loops and a gas of noninteracting heavier hadrons [35]. However, from our study we conclude that extracting the chiral transition temperature in this model would require further inclusion of higher order interactions between pions. Nonetheless our study has a deeper implication, since it allows us to carry forward our calculations within the HRG model with physical masses at higher baryon densities µ B /T > 3.5, which is currently out of bounds of the lattice calculations, and calculations are in progress in this direction [72]. As discussed in the main paper, the m l derivative of the masses of the ground state pions and kaons, could be very well described within the χPT framework. However, the σ terms of other mesons are not known very precisely. This is true for the σ terms for some ground state mesons like the vector mesons and η, η mesons. The σ terms for excited baryon states are also not very well known. This fact contributes to the uncertainty in the estimation of the chiral condensate. In this section we will evaluate the σ terms with as much precision as possible using the most recent values of the hadron masses available from extensive lattice QCD studies. We have collected the most precise values of the hadron masses as available in the literature, as a function of the pion mass along the lines of constant and physical strange quark mass. We then fit these data to the ansatz M α (M π ) = m exp α + b α · M 2 π and extract its slope b α corresponding to each hadron species labeled by α. Since we need the σ term at the physical point, we have taken the product of the slope b α extracted after performing the fit with the square of the average physical pion mass M phys π = 138 MeV. We first compile the dependence of the mass of the ground states of the isoscalar J P C = 0 −+ and vector J P C = 1 −− mesons, on the pion mass. The η meson mass data were taken from Ref. [56], and the results for η and ρ meson masses were compiled from Refs. [68,90], respectively. From the fit shown in Fig. 4, we extract a σ = 10 ± 7 MeV for the η-meson, and a similar value of the σ term for the η meson i.e., σ = 9.5 ± 4.5 MeV. For the ρ mesons we could extract σ = 27 ± 4 MeV with very good precision, which has enabled us to better constrain the T c from the HRG model. Using the data from Ref. [57], the fit of the mass of K * as a function of M 2 π , which is shown in the left panel of Fig. 5. From the fit, we extract a σ K * = 10 ± 1 MeV. We also determined the σ terms for the ω and φ mesons using the lattice data from Ref. [64], where the mass spectrum was calculated in 2 + 1 flavor QCD for two pion mass values M π = 391, 524 MeV respectively again along the lines of constant m s set to its physical value. We used the same procedure as above. The results are shown in Table. III. Ref. [56]. The middle and the right panels show the analogous fits for the η , ρ mesons with the data taken from Refs. [90] and [68] respectively. All fit results are in units of MeV.
We next extract the σ terms of the excited meson states. For the isotriplet excited states, we could again use the lattice data from Ref. [64]. Since there are only two data points in the fit, in order to estimate the goodness of the fit we have taken the experimentally determined masses as the intercept. This is a reasonable guess since we expect that the reduction of the mass of these excited states in the chiral limit will be very small, within ≈ 5-10% of its physical value and not as dramatic as the pseudo-Goldstone pion modes. We also performed the fit without constraining the intercept in order to estimate the systematic uncertainty in the results of the σ term. The results of the fits are shown in the middle panel of Fig. 5. We find that the σ terms corresponding to the first excited state of the pion are σ = 26(3), 16 MeV for the one-and two-term fits respectively. We see that the different estimates agree with each other if we assume the error on the constrained fit is 50% of its central value. We have done a similar study for the first excited state of the η meson using the lattice data from Ref. [64]. The results are shown in the right panel of Fig. 5. We again observe that the σ term obtained from the unconstrained fit lies within a variation of 50% of the value of σ = 15(3) MeV obtained from the constrained fit.
In our fit estimates we could not observe any significant difference of the σ term of the excited states of mesons in different spin, parity sectors from their ground state values within error. For example the excited states of both strange and nonstrange isoscalars have σ = 10(3), 33(5) MeV, also consistent with the values corresponding to their ground states. Hence we have taken the σ terms of the excited states to be the same as the ground state mesons in different spin, parity sectors but with a 50% error to take into account the systematic uncertainties. The data has been taken from Ref. [57]. The middle and the right panel shows the fit to the mass of the excited states of the isoscalar π and η-mesons as a function of the square of the pion mass, using the data from Ref. [64]. All fit results are in units of MeV.
In the baryon sector, we could find the data for the mass of the ground and excited states as a function of the pion mass only for the nucleon (octet) both for 2 + 1 [91] and 2 + 1 + 1 flavor [92] QCD, which are shown in Fig. 6. For the 2 + 1 + 1 flavor data which correspond to pion masses of 313, 225.9, 135.6 MeV respectively a fit to the ground state yields σ = 36(2) MeV, whereas for the excited state the data are available only for the first and last values of the pion mass stated above, which gives a σ = 68 (27) MeV. For the 2 + 1 flavor data, a fit to the positive parity nucleon mass data and its first excited state gives σ = 27 ± 2 and 36 ± 6 MeV, respectively, which are quite close to each other. For the 2 + 1 flavor fits we assumed that the ground state mass of the nucleon is reduced by 70 MeV in the chiral limit and the difference of its mass with the first excited state remains similar, ≈ 770 MeV, as seen in Ref. [61]. A more sophisticated fit strategy discussed in Ref. [61] gives a ground state octet σ ≈ 44 MeV, which we use in our paper. We nevertheless performed these simple fits to understand how much the σ terms of the excited state deviate from the ground state of the baryon octet. We find that the central value varies at most by 1.5 times from the ground state estimates but agree with each other within errors. We summarize our findings of our fits in the Table III. We have also roughly estimated the σ terms for the higher spin excited states of the baryon octet and decuplet using the data from [62]. Here again we have the mass data for the excited states for two different pion masses 396 and 524 MeV, respectively. For the spin 1 2 , 3 2 , 5 2 , and 7 2 first excited nucleon states we have σ = 37, 51, 51, and 51 MeV, respectively by performing a two-term unconstrained fit. For the first excited state of ∆ we have σ = 40, 39, and 38 MeV, for spins 1 2 , 3 2 , and 5 2 , respectively. Since there are only two data points in each fit, we could not estimate the errors of the σ terms but they all lie within 1σ error of the ground state estimates [61]. Henceforth in this work the σ terms for these excited states of baryons are also taken to be the same as their ground states. Similarly to our prescription for the excited states of mesons we have considered a 50% variation of the central values of these σ terms when estimating the error on the chiral condensate. We discuss here the relative contributions of the mesons and baryons with different quantum numbers to the renormalized chiral condensate as a function of temperature. The relative contributions of different hadrons are shown in Fig. 7. Within the meson sector, we find that the largest relative contribution to the renormalized chiral condensate near the crossover transition temperature comes from the ground state pions and kaons followed by the vector meson. The combined relative contribution from the excited states of both strange and nonstrange mesons is only 1/3 of the total ground state contributions. The overall meson contributions to the renormalized chiral condensate is almost 75% at a temperature T = 160 MeV. The remaining contributions come from the baryon sector, dominantly from the ground state octet and the excited states, primarily those of the nucleons.