Thermodynamics of a gas of hadrons with attractive and repulsive interaction within S-matrix formalism

We report the effect of including repulsive interactions on various thermodynamic observables calculated using a S-matrix based Hadron Resonance Gas (HRG) model to already available corresponding results with only attractive interactions \cite{Dash:2018can}. The attractive part of the interaction is calculated by parameterizing the two body phase shifts using K-matrix formalism while the repulsive part is included by fitting to experimental phase shifts which carry the information about the nature of interaction. We find that the bulk thermodynamic variables for a gas of hadrons such as energy density, pressure, entropy density, speed of sound and specific heat are suppressed by the inclusion of repulsive interactions and more sensitive for second and higher order correlation and fluctuation, particularly for the observable $\chi^2_Q$, $\chi^2_B-\chi^4_B$ and $C_{BS}$. We find a good agreement for the $C_{BS}$ and lattice QCD simulations. Additionally, assuming that the value of interacting pressure versus temperature for a gas of hadrons calculated in S-matrix formalism is same as that from a van der Waals HRG (VDHRG) model, we have quantified the attractive and repulsive interactions in our model in terms of attractive and repulsive parameters used in the VDHRG model. The values of parameters thus obtained are comparable to that which reproduces the properties of nuclear matter in the ground state.


I. INTRODUCTION
One of the primary goals of relativistic heavy ion collision is the study of QCD (Quantum Chromo Dynamics) phase diagram [2]. There are at least two phases in the phase diagram, one where the degree of freedom are quarks and gluons called the Quark Gluon Plasma (QGP) phase and other where the degrees of freedom are hadronic. An approach to study the properties of hadronic phase formed by hadronization of the QGP is through a statistical model of a gas of hadrons called the hadron resonance gas model (HRG) [3]. The hadron resonance gas (HRG) [2,[4][5][6][7][8][9][10][11][12][13][14] models have successfully described the hadron multiplicities produced in relativistic nuclear collisions over a wide range of center of mass energy. The main result of such an investigation was the observation of rise in the extracted chemical freeze-out temperature values from lower energies to almost a constant value of temperature T = 160 − 165 MeV at higher energies, supplemented with the decrease of the baryon chemical potential (µ B ) with increasing energy [15]. The saturation of temperature supports the Hagedron's limiting temperature hypothesis [16], suggesting the possibility of a phase boundary. Similarly, theoretical investigation of QCD on lattice (LQCD) at vanishing µ B indeed predicts a sharp increase of thermodynamical quantities near deconfinement temperature T c [17][18][19][20][21][22][23][24]. The HRG model is also successful in describing LQCD data related to the bulk properties of hadronic matter in thermal and chemical equilibrium below T c [18,20,21,[23][24][25].
The phenomenal success of the ideal HRG (IDHRG) model in predicting the hadronic yields can be attributed to a theorem by Dashen and Ma [26] which states that the partition function of an interacting hadronic gas, can be decomposed into a free and an interacting part. Considering that only resonances contribute to the interacting part, it can be shown that in a narrow resonance width approximation, the net effect of the interacting part is equivalent to considering all such hadronic resonances as free particles. However, relaxing the above assumptions by including overlapping resonances and resonances of finite widths, it has been seen that the variation of thermodynamic variables with temperature changes substantially [27][28][29][30][31]. Further, it can be argued that such interaction contribute only to the attractive part of partition function and the inclusion of a repulsive part could partially negate the effect of the attractive part. For example, in Refs. [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50] the authors have used an excluded volume approach which only had the repulsive part whereas Refs [51][52][53][54][55][56][57][58] considered an van der Waals' (VDW) type of interaction, which has both the attractive and repulsive part and a comparison of the calculated thermodynamic pressure from both the approaches shows the feature as discussed above.
In the present work, we extend the K-matrix formalism for HRG model that we had previously developed with attractive interactions among the hadrons to include the effect of repulsive interactions between the hadrons in the S-matrix formalism. In Ref. [1], we used K-matrix formalism to calculate the phase shifts of the resonance spectral function in contrast to the popular Breit-Wigner parametrization. It has been argued previously that the K-matrix formalism preserves the unitarity of the scattering matrix (S-matrix) and neatly handles multiple resonances [1,30,59]. However, the formalism fails to handle any repulsive channel in the scattering matrix. Therefore in this work we include the repulsive part by fitting to experimental phase shifts that encodes the information about the nature of interaction. We use the phase shifts data form Scattering Analysis Interactive Database (SAID) partial wave analysis for nucleon-nucleon (N N ), pion-nucleon (πN ) and kaon-nucleon (KN ) interaction in their respective isospin channels [60][61][62]. Additionally, we have also included the repulsive isotensor channel in the pion-pion (ππ) scattering, as have been pointed in many earlier works [28,63].
After constructing the interacting hadron resonance gas model with both attractive and repulsive interactions using phase shift information for various hadronic interactions we calculate the various thermodynamic observables like pressure, energy density, entropy density, interaction measure, specific heat, speed of sound and susceptibilities. The temperature dependence of these observables are then compared with corresponding results from Lattice QCD, IDHRG and HRG with attractive interactions using K-matrix formalism.
The paper is organized in the following manner. In the next section we discuss the formalism used to introduce repulsive interactions to our HRG model developed earlier using Kmatrix approach with attractive interactions [1]. In Sec.III we discuss the results from the new interacting HRG model with both attractive and repulsive interactions among the hadrons. The temperature dependence of our results are compared to those from LQCD and IDHRG (with different hadron spectrum). Finally in Sec. IV we summarize our findings.

II. FORMALISM
The equation of state for an interacting gas of hadrons can be computed by using the method of virial expansion. Specifically, the pressure of such a gas can be written as [64], where ξ = (m/2βπ) 3/2 e βµ and the inverse temperature, chemical potential, mass are denoted by β, µ, m. The term J i takes into account the interaction between groups of i hadrons and which are given as, etc., where U 12 is the interaction energy. Differentiating Eq.
(1) with respect to µ, we obtain the expression for number density i.e., Eliminating ξ to the first order from Eq.(1) and (3) gives us the ideal equation of state P = nT , where T is the temperature. For a relativistic non-interacting quantum gas the expression for the pressure is given in [28]. In the present work the contribution to the non-interacting part comes from all the stable hadrons. In the second order we retrieve the van der Waals' equation of state P = nT (1+nB(T )), where B(T ) = −J 2 /2 is called the second virial coefficient. In this work while calculating the virial coefficients we will be using the S-matrix approach to statistical mechanics, which has also been used previously in Refs. [28,[65][66][67] to study the thermodynamics of interacting hadrons.
In the S-matrix formalism, the second virial coefficient is related to the scattering amplitude or alternatively to the scattering phase shifts δ I l for a given spin l and isospin I channel. The correction to the ideal pressure i.e., the second term in the expansion of Eq.(1) called P int is given as where the terms z i , g I,l and ε stand for the fugacity, the spinisospin degeneracy factor and the total center of mass energy respectively. The function K 2 (x) stands for the modified Bessel function of second kind and the term M is the invariant mass of the interacting hadron pair at threshold. Additionally, there is a sum over all possible spin-isospin channels and the prime over the summation sign denotes that for given l the sum over I is restricted to values consistent with statistics. From the resulting Eq.(4), it can be seen the second virial coefficient gives positive (attractive) or negative (repulsive) contribution depending on whether the derivative of phase shifts are positive or negative. The phase shifts are obtained from experiments or from theoretical calculations. In the present work, we determine the attractive phase shifts using the K-matrix formalism which takes the masses and widths of resonances from the PDG (Particle Data Group) [68] as input. Since the K-matrix formalism is not applicable for handling the repulsive phase shifts, these are obtained by fitting to experimental data. We would like to note here that since we do not have the information of masses and widths of resonances (mentioned in PDG) that decay into a pair of nucleons, we extract phase shifts in such situation by fitting to experimental data.

A. K-matrix Formalism
A theoretical way of calculating phase shifts is to use the K-matrix formalism. The K-matrix formalism preserves the unitarity of S-matrix and neatly handles multiple resonances [59]. In addition to that, widths of the resonances are handled naturally in the above formalism. In contrast to the notion of ideal HRG is only valid for narrow resonances and not for broad resonances but the K-matrix formalism can be applied quite generally. Similarly for overlapping resonances the K-matrix gives a more accurate description of the phase shifts than the Breit-Wigner parametrization. In Ref. [30] the K-matrix formalism was used to study an interacting gas of hadrons and it was extended further in [1].
The resonances, contributing to the process ab → R → cd, appear as a sum of poles in the K-matrix, where a, b and c, d are hadrons and the sum on R runs over the number of resonances with mass m R . The sum is restricted to the addition of resonances for a given spin l and isospin I. The residue functions are given by where √ s is the center of mass energy and Γ R→ab ( √ s) is the energy dependent decay width.
Furthermore, once one computes the K-matrix by providing the relevant masses and widths of resonances, the phase shift can obtained using the relation:

B. Experimental Phase shifts
As mentioned earlier for repulsive interactions and for interactions where the information about m R and Γ R are not available, the K-matrix formalism is not applicable and we resort to extraction of phase shifts from experimental data. In our extraction of repulsive (πN , KN ) and nucleon-nucleon (N N ) interaction phase shifts we use the data from the SM16 partial wave analysis [60]. For the repulsive isotensor channel δ 2 0 in the π − π scattering we use the data from Ref. [69]. However, the S-matrix formalism elucidated here is only applicable for elastic scattering and the inelastic part that enters into the analysis by fitting to experimental data has to be removed. To get around this problem, we make an estimate of the contribution coming from the inelastic part by first defining a generic l dependent scattering amplitude f l ( √ s), where η l is the inelastic parameter and the elastic crosssection, and the inelastic cross-section are given by where q is center of mass momentum. The total cross section σ is the sum of Eq.(9) and Eq.(10). We can approximate the contribution to the elastic part of the phase shift δ el by the following expression where δ is the total phase shift that is obtained from fit to experimental data [60][61][62].

N-N interactions
For the nucleon-nucleon (N N ) interaction we have included the phase shifts for l ≤ 7 in both I = 0 and I = 1 isospin channels. Combinations of l, S and I are chosen so that the total wave function for N N interaction is asymmetric as dictated by Pauli's principle. We have restricted the range of energies up to the pion (π) production threshold. Beyond this threshold, the contribution from the inelastic channels become dominant and the present formalism fails to disentangle the contribution from the elastic and inelastic part. However, below this threshold where the contribution from inelasticities are sub-dominant, we can extract the contribution from the elastic part using the approximation Eq. (11). In order to use Eq.(11) we need a parametrization of the cross-section as function of energy which in the present study are used from Ref. [70].
where p lab is the laboratory momentum. Similarly the elastic cross-section σ el can be parametrized as By comparing Eq.(12) and Eq. (13), we can see that that contribution from inelastic processes is small below p lab < 0.8 GeV and increases further with energy. In Fig. 1 we have plotted the experimental NN phase shifts from the SAID partial-wave analysis [60] as a function of center of mass energy ( √ s). Dominant contribution comes from lower l values for e.g the 1 S 0 phase shift which peaks at lower √ s and then falls sharply or the rapidly falling and largely repulsive 3 S 1 phase shift. An interesting case to observe are the triplet Pwaves which can have J = 0, 1, 2 corresponding to phase shifts 3 P 0 , 3 P 1 , 3 P 2 . The behavior of the phase shifts are quite different in the above three channels, from zero crossing to purely repulsive and purely attractive as seen in Fig.1. This could be accounted because of the spin-orbit coupling which splits them in to the triplet states having different behavior depending on the sign and strength of the coupling. However, most of the phase shifts become negative at higher √ s energies signifying the hard core nature of NN interaction.

π-N interactions
For pion-nucleon (πN ) interaction we have included only those phase shifts [61] which are purely repulsive and the attractive part are from the K-matrix parametrization. Here, we where q is the center of mass momentum. In the range of momenta 0.5 GeV< p lab < 1.5 GeV, the inelastic channel πN → ππN is the most dominant whose cross section can be parametrized as The dominant repulsive contribution in the πN interaction comes from the S 31 (l 2I,2J ) phase shift corresponding to ∆(1620) resonance. We would like to stress here that in our previous study of resonances in Ref [1], using K-matrix formalism, resonances such as ∆(1620), ∆(1910), N (1720) etc. were included in the attractive part of the S-matrix. However, a comparison to experimental phase shifts has rendered that, this is not the case. Thus, such resonances are included in the repulsive part by fitting to experimental phase shifts.

K-N interactions
For the KN interaction, the dominant repulsive contribution comes from the S 11 (l I,2J ) phase shift containg the Σ(1660) resonance. Here, as in case of πN resonances like etc. were considered attractive in [1], but here we include them in the repulsive part as their experimental phase shift have confirmed [62]. The cross sections are parametrized from Ref. [71] σ

π-π interactions
For the pion-pion (ππ) interaction we have included the dominantly repulsive phase shift from Ref. [69], in the isotensor channel δ 2 0 , as have many previous studies indicated [28,63]. This phase shift is known to cancel the isoscalar channel δ 0 0 containing the broad f 0 (500) (σ meson). The relevant energies have been restricted to pion production threshold.

III. RESULT
In Fig.2, we show the temperature variation at zero chemical potential for various thermodynamic quantities such as scaled pressure, energy density, entropy density, speed of sound and the specific heat capacity at constant volume. Results of attractive K-matrix (KM) based HRG model from Ref. [1] are compared with the total contribution (Total) from both attractive and repulsive channels done in the present work. It should be noted that both KM and Total contain noninteracting part as well. We observe that the effect of repulsive interactions cancels some of the contribution from attractive channels, there by slightly lowering the net result for Total relative to KM for the observables studied here. A comparison with ideal HRG, that considers all the confirmed hadrons and resonances consisting of up, down, and strange flavor valence quarks listed in the PDG 2016 Review [68], shows that effect of repulsion bring the results closer to the ideal limit. However, it is worth mentioning here that the agreement of ideal HRG with the LQCD data is because of the increase in the number of degeneracies and not due to some inherent interaction that is naturally present in the system revealed within the S-matrix formalism. On the whole, we conclude that the effect of repulsive channels suppress the bulk variables studied here, compared to K-matrix (KM) approach and are shown in Fig. 2.  [68]. Results are compared with lattice QCD data of Refs. [18] (WB), [20] (HotQCD) and [24] (Lattice).
The effect of repulsive interactions are most prominent when we calculate second order diagonal and off diagonal susceptibilities. Results for χ 2 B , χ 2 Q , χ 11 BS and χ 11 BQ (B, Q stand for baryon and electric charge respectively, definition of susceptibilities can be found in Ref. [1]) agree better with the LQCD data, in the case when both attraction and repulsion is taken into account than in the K-matrix formalism. The effect of repulsion is mostly visible in the baryonic sector because of the inclusion of dominantly repulsive N N interactions as discussed previously. We have also checked for the remaining second order diagonal and the off diagonal susceptibilities, and the difference between Total and the K-matrix formalism for χ 2 S is small and for χ 11 QS is negligible. The contribution from repulsive interaction can be explored further by considering certain combinations of diagonal and off diagonal susceptibilities which are identically zero for non-interacting or ideal HRG but not for non-interacting gas of quarks and gluons [72]. The quantity χ 2 B − χ 4 B = 0, for a hadron gas which has baryon number ±1, but for noninteracting QGP for which χ 2 B − χ 4 B > 0, since all quarks carry a baryon number of ±1/3. However for an interact-ing gas, the inclusion of N N interaction which carries a net baryon number ±2 might give us a non zero result. This is shown in Fig.4(a), and the result shows that sign of χ 2 B −χ 4 B > 0 again indicating the hard core nature of N N interaction. In Ref. [55,67] the same increasing trend of χ 2 B − χ 4 B with temperature was also found using repulsive mean field in a multi-component hadron gas and excluded volume approach. Our results using the S-matrix formalism validate the previous results. Other observables like v 1 = χ 31 [72] are trivially zero in our analysis since we do not have the information about interactions (phase shifts) among baryons which have |B| > 1 and |S| = 1 or vice-versa.
The correlation between the strangeness S and baryon number B is a sensitive probe of the relevant microscopic degrees of freedom. The quantity C BS [73] defined as C BS = −3χ 11 BS /χ 2 S is one such observable. For a gas of non-interacting QGP, C BS = 1 but for a gas of kaons and anti-kaons where a light quark is always correlated with its strange partner (kaons) or vice versa (anti-kaons) C BS < 1.  with strange quark (anti-quark) and hence have C BS > 1. Therefore, for large baryon chemical potential, C BS could be larger than unity in a hadron gas. Moreover, significant difference between LQCD and ideal HRG has been reported previously [74]. It has been argued that such discrepancy can be cured by allowing additional strange hadrons which have not been discovered but are predicted in various quark models [75,76]. Fig. 4(b) shows that the difference between LQCD and ideal HRG can accounted by including interaction with-out invoking any additional hadrons.
Finally, we match the second virial coefficient obtained using S-matrix formalism with the virial coefficient B(T ) of a Van der Waals gas and extract the VDW parameters a and b. For a VDW gas the coefficient B(T ) is given as [64] where b = 16πr 3 /3, where r is the hard core radius and a is a positive constant denoting attraction. Thus, the interacting pressure P V DW int for a VDW equation of state is given as Matching Eq.(4) and Eq. (18), we extract the values of a and b. Fig.(5) shows the P int /n 2 int calculated using Eq. 4 as a function of temperature. This study indicates that the simple (constant) parametrization of the VDW parameters is not correct in a realistic situation, where both the attractive parameter a and the repulsive parameter b could in general be temperature dependent. This fact also supports models [28] where a temperature dependent radius was used. However assuming the VDW parameters are temperature independent, a straight line fit to the results in Fig. 5 with a functional form of −bT + a is carried out to extract the VDW parameters. The values of the VDW parameters are a = 0.54 ± 0.014 GeV fm 3 and the hard core radius r = 0.56 ± 0.005 fm. This is in qualitative agreement with a previous work [52], where VDW interaction were used in the HRG model with a = 0.329 GeV fm 3 and r = 0.59 fm by fixed to describe the ground state properties of nuclear matter. However the value of a and r are smaller from those extracted by fitting VDWHRG model results to LQCD data at zero chemical potential [58].

IV. SUMMARY
To summarize we have included repulsive interaction by fitting to experimental phase shifts which carry the information about the nature of interaction. The attractive part of the interaction are also included and was calculated by parameterizing the two body phase shifts using K-matrix formalism [1] which is known to preserve the unitarity of S-matrix. Since experimental phase shits for attractive part of the interactions was available for the N N scattering, those are used in the calculations. Thermodynamic quantities like pressure, energy density, trace anomaly, specific heat and speed of sound etc. were calculated using the S-matrix formalism. The results indicate that the effect of repulsive channels is to suppress the bulk variables studied here.
The most prominent effect of repulsive interactions are seen when we calculate the second and higher order fluctuations and correlation. Particularly the two most interesting observations which resulted from the current work are as follows. We find that the observable χ 2 B − χ 4 B > 0 for the interacting HRG model discussed in this work, contrary to the expectation χ 2 B − χ 4 B = 0 for an uncorrelated gas of hadrons like of IDHRG model. Similarly, for the observable C BS which is a sensitive probe of the relevant microscopic degrees of freedom of a system, the HRG model in the present formalism very well describes the LQCD data. Also seen from Fig. 4(b) that IDHRG model with additional strange hadrons which has not yet been discovered agrees with the LQCD data at a similar level [74]. The difference in physics interpretation is the following: while the IDHRG model with additional strange hadrons attributes the matching with LQCD data relative to normal IDHRG model, to the increase in hadronic degrees of freedom for the system, our results in contrast attributes the matching to be due to interactions among the constitutes that is captured naturally through the formalism used in the current work.
Finally we have tried to quantify the attractive and repulsive interactions in our model in terms of the VDWHRG attractive and repulsive parameters a and r, respectively. In doing so we assume that the parameter values do not change with temperature and the interactive pressure values are same in the two models at a given temperature. The extracted values of a and r are of the similar order as reported for a VDWHRG model that describes the LQCD data with its parameters fixed to those that describe the ground state properties of nuclear matter [52]. It may be noted that our results as shown in Fig.  5 indicates a and r could be temperature dependent. We will end by saying that as an outlook it would be interesting to calculate various transport co-efficients in a S-matrix based HRG model and compare to other different type of HRG model and corresponding Lattice QCD results.
ACKNOWLEDGEMENT BM acknowledges financial support from J C Bose National Fellowship of DST, Government of India. SS and AD acknowledge financial support from DAE, Government of India.