Covariant interacting Hadron-Resonance Gas model

The Hadron-Resonance Gas (HRG) approach - used to model hadronic matter at small baryon potentials $\mu_B$ and finite temperature $T$ - is extended to finite and large chemical potentials by introducing interactions between baryons in line with relativistic mean-field theory defining an interacting HRG (IHRG). Using lattice data for $\mu_B=0$ as well as information on the nuclear equation of state at $T=0$ we constrain the attractive and repulsive interactions of the IHRG such that it reproduces the lattice equation of state at $\mu_B=0$ and the nuclear equation of state at $T=0$ and finite $\mu_B$. The formulated covariant approach is thermodynamically consistent and allows us to provide further information on the phase boundary between hadronic and partonic phases of strongly interacting matter by assuming constant thermodynamic potentials.


I. INTRODUCTION
The phase diagram of matter is one of the most important subjects in physics since it also has important implications on chemistry and biology. Accordingly, the phase diagram of strongly interacting matter is a topic of utmost interest since decades and substantial experimental and theoretical efforts have been invested to shed light on this issue. The early universe went through different phases at practically vanishing baryon chemical potential µ B when expanding to its present size. Relativistic and ultra-relativistic heavy-ion collisions nowadays offer the unique possibility to study some of these phases, in particular a quark-gluon plasma (QGP) phase and its phase boundary to the hadronic one. Lattice Quantum-Chromo-Dynamics (lQCD) calculations suggest that at vanishing baryon chemical potential (µ B = 0) there is a crossover phase transition from hadronic to partonic degrees of freedom [1][2][3][4][5][6][7] for the deconfinement phase transition as well as for the restoration of chiral symmetry. However, at some finite baryon chemical potential the crossover might turn to a first-order phase transition implying a critical endpoint in the QCD phase diagram [8]. Since lattice calculations so far suffer from the fermion-sign problem, no first-principles information on the phase boundary can be extracted from lQCD at large µ B , whereas at low µ B Taylor expansions of the thermodynamic potential (in powers of µ B /T ) provide an alternative solution as demonstrated in Refs. [9,10]. Accordingly, heavy-ion reactions at Relativistic Heavy-Ion Collider (RHIC) and Large Hadron Collider (LHC) energies show no evidence for a phase transition and experimental studies at much lower bombarding energies are needed to explore the high µ B region of the phase diagram. To this end new facilities -such as the Facility for Antiproton and Ion Research (FAIR) and the Nuclotronbased Ion Collider fAcility (NICA) -have been planned and are presently under construction. * Thorsten.Steinert@theo.physik.uni-giessen. de The QCD phase diagram has two regions that are relatively well known, while the rest is more or less unknown. These two regions are the temperature axis at µ B = 0, which can be studied by lQCD calculations, and the region of T = 0, which is described by nuclear physics. Usually one constructs models that can only be applied to one of the two cases, but we intend to set up an approach that is able to describe both regions simultaneously and to allow for more stringent extrapolations on unknown areas of the phase diagram. A common model used in the first regime is the Hadron-Resonance Gas (HRG) model that treats hadrons as a gas of non-interacting particles. This model works at vanishing chemical potential but fails for the description of nuclear matter due to the lack of repulsive and attractive interactions. The latter are included in relativistic mean-field theories (RMFTs) whose interactions are based on meson-exchange potentials. While these models can describe infinite nuclear matter with the right properties of the binding energy, they fail for the QCD equation of state at vanishing chemical potential µ B . We will combine both approaches in the following to set up a model that is consistent with lQCD (µ B ≈ 0, T > 0) and the nuclear equation of state (T ≈ 0, µ B > 0) while using only hadronic degrees of freedom. Of course, nothing can be predicted for the partonic phase at high temperature or chemical potential, but we aim at defining a phase boundary between the hadronic and partonic phases by approaching it from the hadronic side.
It is commonly expected that the HRG gives a good description of hadronic 'QCD matter' at moderate chemical potentials due to the success of the statistical hadronization model in describing particle ratios from relativistic heavy-ion collisions. This model assumes that the medium -created in a heavy-ion collision -equilibrates (to some extent). The system will then continue to interact until it becomes too dilute and freezes out. The particle yields then are fixed by the temperature and the chemical potential at the chemical freeze out (apart from the volume, which drops out in particle ratios). This simple model can describe the particle abundances for various collision energies [11][12][13][14] from Alternating Gradient Synchrotron (AGS) up to LHC energies. It becomes even more precise if one assumes additional corrections for non-equilibrium effects [15][16][17][18]. The model is applicable also for p + p [19] and even e + e − collisions [20] and was also applied to the production of hypernuclei in Ref. [21]. Furthermore, the HRG equation of state and susceptibilities were compared to state of the art lQCD calculations in Refs. [22][23][24][25][26][27][28][29] and it was found that the ideal HRG leads to a satisfying description of the thermodynamics for temperatures below T ≈ 170 MeV. The quality of the description is improved if one includes an exponential increasing mass spectrum [30], as predicted by Hagedorn, and repulsive interactions [31][32][33][34][35][36].
On the other hand, the main application for relativistic mean-field theories is the calculation of ground-state properties of infinite nuclear matter and finite nuclei. The model works well for spherical and deformed nuclei [37,38], and can also be used to investigate neutron star properties [39,40], in particular the mass-radius relation. Furthermore, the RMFT provides the basis for nonequilibrium transport approaches when applied to heavyion collisions at lower energies or merely the hadronic phase [41,42].
This work is organized as follows: In Sec. II we recapitulate the equation of state from the HRG at vanishing chemical potential while in Sec. III we recall the framework of the relativistic mean-field model and discuss the nuclear equation of state at vanishing temperature. In Sec. IV we will combine both approaches in a thermodynamic consistent manner and define the covariant interacting hadron resonance gas (IHRG) approach, which will allow us to define a phase boundary in the (T, µ B ) plane from the hadronic side. This study is summarized in Sec. V.

II. REMINDER OF THE HADRON-RESONANCE GAS (HRG) MODEL
The most frequently used model for the thermodynamics of hadrons at finite temperature T and baryon chemical potential µ B is the HRG. The approach is based on the work of Dashen, Ma and Bernstein [43], who found that one can describe the thermodynamics of a system of particles, which interact through resonant scatterings, by including the resonances as stable particles in the partition sum. This is always possible if the spectral widths of the resonances are small compared to the temperature γ ≪ T . As an example we quote the interacting pion gas, which is thermodynamically equivalent to a free gas of pions and ρ-mesons [44,45]. The HRG generalizes this approach to all possible hadrons such that the thermodynamic potential for the hadronic system Ω HRG is given by the sum over all stable hadrons and all known hadronic resonances, without mutual interactions, In Eq. (2) n B/F denotes the Bose/Fermi distribution while d i is the particle degeneracy. The approach incorporates attractive interactions (for the dynamical formation of resonances) but discards repulsive interactions that describe the short-range repulsion between the hadrons. However, the effects from short-range repulsion can be introduced by assuming a finite volume of the particles [46][47][48][49] which is excluded in the thermodynamic analysis and thus leads to an increase of the pressure P (for fixed particle number, cf. Eq. (3)). The model presented in Refs. [46,47] assumes the same volume for each particle such that the excluded volume is proportional to the total particle density. The approach in Refs. [48,49] assumes the excluded volume for each particle to be proportional to its energy, which leads to a limiting energy density similar to the Hagedorn model (with a maximum temperature [50]). A non-relativistic version of this approach is the van der Waals model with the pressure P and internal energy U given by, where b characterizes the excluded volume and a the strength of the attractive interaction. The repulsive interactions, which are incorporated in excluded volume models, are quite different from repulsive interactions originating from vector mesons exchange as, for example, in the Walecka model [51]. In the latter case the strength of the vector repulsion is proportional to the net particle density and thus vanishes for µ B = 0. On the other hand, the excluded volume repulsion is proportional to the total particle density (or pressure) and is finite also at vanishing chemical potential. The only parameters in (1) are the masses of the hadrons, which are usually taken as the vacuum masses. In principle, the model is parameter free, but one has to decide on the amount of 'non-interacting particles' to include. The most fundamental hadrons are the spin-1/2 baryons and 0 − mesons as well as the spin-3/2 baryons and the 1 − mesons as resonances. All of these hadrons are important for the dynamics of heavy-ion collisions and have to be included to achieve a reasonable description of finite temperature QCD in the confined hadronic phase. In the context of chiral symmetry restoration at high temperature and/or baryon density other hadrons also play an important role, i.e. the chiral partners of opposite parity, e.g. the a 1 -meson and the N (1440) and N (1535) baryon or the scalar 0 + mesons. Among the 0 + mesons the most prominent is the f 0 (500) or σ-meson with a mass of 400 − 550 MeV [52], which is now established as a particle and contained in the latest version of the Particle Data Group (PDG) [53]. The other mesons in the scalar 0 + -multiplet are the f 0 (980) and The dimensionless pressure P/T 4 for a HRG with varying amounts of hadrons as a function of the temperature T . The dotted magenta line considers only pions, kaons and η's, the dashed blue line is the basic HRG-0 and the full red line shows the pressure for all hadrons listed by the particle data group [60] with a mass below 2.6 GeV (HRG-1). The lQCD results are taken from the Wuppertal-Budapest collaboration from Ref. [24].
the a 0 (980), which are both listed by the PDG, as well as the strange κ(720) meson. However, there is no consensus whether to include the σ and the κ in HRG models or not since calculations for the thermodynamic potential of an interacting pion gas in terms of experimental phase shifts show that the attractive pressure contribution from the scalar σ-mesons gets exactly canceled by the repulsive isotensor channel [45,54]; the same happens also for the κ-mesons. However, in these calculations the vacuum phase shifts have been employed. On the other hand, effective hadron field theories suggest that the σ-meson is much stronger affected by finite temperature or chemical potential effects than the pions and ρ-mesons [55][56][57][58]; accordingly vacuum phase shifts should no longer be valid for temperatures above T ≈ 100 MeV. Arguments in favor of the scalar mesons come from the statistical hadronization model in Ref. [59], where the σ-meson was explicitly included to improve the description of the K + /π + ratio that was observed experimentally in central Au+Au (Pb+Pb) collisions. Furthermore, it was shown in Ref. [52] that a neglect of the 0 + -mesons in the thermodynamic partition sum is inconsistent with respect to causality and unitarity. Although the 0 − , 0 + and 1 − mesons, the spin-1/2 and spin-3/2 baryons together with e.g. the a 1 , N (1440) and N (1535) are the most basic hadronic degrees of freedom required to describe the hadronic medium, they do not produce the required pressure at higher temperatures (above ∼ 150 MeV) to describe the lQCD equation of state. It is therefore necessary to include additional hadrons and a standard choice is to incorporate all hadrons listed by the Particle Data Group [53,60] with a mass below a certain threshold M max . For illus-tration we show in Fig. 1 the dimensionless pressures P/T 4 = −Ω HRG /T 4 at vanishing chemical potential calculated with different numbers of particles in comparison to the recent lQCD data from the Wuppertal-Budapest collaboration [24]. We denote the HRG with the most basic hadrons (quoted above) as HRG-0 and show also the result when including all hadrons from the PDG [60] with a mass below 2.6 GeV (HRG-1). A compact list of the hadrons, but without the σ and the κ-meson can be found in Ref. [61]. We see that at low temperatures the system is dominated by the lightest mesons, i.e. pions, kaons and the η-meson, which describe the equation of state up to T ≈ 120 MeV. In general the equation of state is meson dominated and baryons contribute only above T ≈ 140 MeV. The basic HRG-0 describes the lQCD data only up to T ≈ 150 MeV and underestimates the pressure at higher temperatures. The additional hadrons in HRG-1 provide the necessary pressure to describe the equation of state up to T ≈ 180 MeV within the error bars of the lQCD result. However, the speed of sound defined by with the energy density turns out to fail the lQCD results. This is demonstrated in Fig. 2 where the speed of sound (squared) c 2 s -for the different HRG versions -is compared with the lQCD data taken from the same simulation as the pressure [24]. Whereas the lQCD data show a clear minimum for T ≈ 150 MeV, the HRG versions do not show a substantial increase for higher temperature regardless of the hadron content. Any inclusion of further (heavier) resonances leads to an overall decreasing speed of sound for T > 150 MeV.
Extending the HRG to even heavier resonances with masses beyond 2.6 GeV has almost no effect on the pressure because additional hadronic degrees of freedom have only small effects on the equation of state. However, a comparison between the HRG and lQCD data shows that the experimentally established hadrons are not sufficient to describe strangeness fluctuations [62], which indicates that a full description of QCD (above about temperatures of 150 MeV) in terms of hadronic degrees of freedom requires an even stronger interaction.

III. RELATIVISTIC MEAN-FIELD THEORY
As mentioned above, an increase of the pressure can be achieved by an excluded volume in the van der Waals model that mimics the effects of short-range repulsive The speed of sound squared c 2 s for a non-interacting hadron gas as a function of temperature. The meaning of the lines is the same is in Fig. 1. The lQCD results are taken from the Wuppertal-Budapest collaboration from Ref. [24].
forces. The latter are naturally included in covariant effective theories with baryons and mesons by the massive vector mesons. We recall that the limit of T = 0 and finite chemical potential µ B is known as the "infinite nuclear matter" limit, a scenario found in the interior of all larger atomic nuclei, where the nuclear density is almost constant. The binding energy of nuclear matter E B /A = E/ρ B − m N has a minimum of E B /A ≈ −16 MeV at normal nuclear density ρ 0 ≈ 0.16 fm −3 to reproduce the stable nuclear matter inside finite atomic nuclei. For T = 0 the Fermi-distribution function becomes while the distribution function for anti-fermions and the Bose-distribution functions for mesons vanish. In Eq. (6) the momentum p F is the Fermi-momentum which specifies the largest occupied momentum state; at T = 0 it fixes the baryon density via and it is therefore common to describe nuclear matter properties in terms of the density ρ B . In (7) a degeneracy factor of 4 has been introduced (for protons and neutrons with two spin projections). Note that without anti-baryons the net-baryon density and the total baryon density are the same; this changes for T = 0.
Since only baryons with a mass larger than the chemical potential µ B can populate the system at T = 0, the conventional HRG can not describe the nuclear equation of state and in particular the minimum at ρ 0 . Baryons other than nucleons can appear only at very large (energy) densities that most likely have to be described by partonic degrees of freedom. The HRG -without inter-actions -reduces essentially to a gas of nucleons and thus fails not only at high T and µ B = 0, but also for T = 0.
The interactions between nucleons conventionally are assumed to be mediated by meson exchange and are described by relativistic mean-field theories [51,[63][64][65]. One usually includes isoscalar interactions -mediated by the scalar σ-meson and the vector ω-meson -and isospindependent interactions that are mediated by the ρ-meson and the δ-meson. The σ-meson describes the attractive part of the nucleon-nucleon interaction while the ωmeson is responsible for the short-range repulsion. The ρ and the δ-meson are important for asymmetric nuclear matter and neutron star physics, but give no contribution in isospin symmetric matter. Since this is approximately the case for the hot and dense medium created in heavyion collisions, we can neglect them in the following. The Lagrangian of relativistic mean-field theory then reads, where the functions U (σ) and O(ω µ ω µ ) describe selfinteractions of the σ-and the ω-fields that are introduced to incorporate non-linear density dependences as arising from Dirac-Brueckner calculations. The couplings Γ σN and Γ ωN are not constants but can be functions of the nucleon field [66][67][68] in order to better comply with results from Dirac-Brueckner calculations, too. To preserve Lorentz invariance the couplings have to be Lorentz scalars. The easiest way to ensure this is to write them as a function of a density Γ(Ψ, Ψ) = Γ(ρ 0 ), which is a Lorentz scalar itself. Two physical reasonable choices arê ρ 0 =ΨΨ andρ 0 =Ψu µ γ µ Ψ, where u µ is the four-velocity with u µ u µ = 1. The first one is denoted by scalar density dependence (SDD) and will lead to a dependence on the scalar density ρ s , the second one is called vector density dependence (VDD) and will lead to a dependence on the baryon density ρ B . It has been shown, that the application of the VDD gives better results when applied to finite nuclei in Refs. [66,67]. Furthermore, the density dependence of the couplings has the advantage that one can parametrize a realistic nucleon-nucleon interaction -as obtained from Dirac-Brueckner (DB) calculationswith less numerical effort [67,68], which allows us to apply DB calculations also to finite systems. We point out that only with density-dependent couplings we can reproduce the nuclear equation of state and the lQCD equation of state simultaneously in the same approach. Note, that for the thermodynamics one can describe the non-linear density dependence of the DB-interactions either with density-dependent couplings or with non-trivial mesonic selfinteractions.
Since the Lagrangian (8) is too complicated to be solved on the many-body level, we will use the mean-field approximation, where the meson fields are no longer in-dependent degrees of freedom but determined by their expectation values. When evaluating the equations of motion one ends up with the following two coupled equations [66], which have to be solved simultaneously: with Here d is the degeneracy of the fermion field, which in case of nucleons is d = 4. The density in the couplings is now the normal ordered expectation value of the density ρ 0 = :ρ 0 : . The distribution functions n F depend on ω * p = p 2 + m * 2 with the effective mass and on the effective chemical potential which both get affected by the interactions with the mesons. The mass m * gets modified by the scalar selfenergy Σ s that originates from the interactions with the σ-meson and the chemical potential gets modified by the vector selfenergy Σ 0 that originates from the interactions with the ω-meson. The selfenergies are split into a conventional Σ (0) and a rearrangement selfinteraction Σ (r) , which arises from the density dependence of the couplings. Their actual forms depend on the choice of ρ 0 . If the couplings are independent from the fields and just constants, the rearrangement selfenergies vanish, but otherwise they are mandatory for thermodynamic consistency (cf. Ref. [66]). In this work we will employ constant scalar couplings and VDD vector couplings. The rearrangement selfenergies for this case read The pressure and the energy density of the densitydependent relativistic mean-field model are given by where P 0 and E 0 are the pressure and energy density for a non-interaction particle evaluated for the effective quantities µ * and m * . The model is thermodynamic consistent as long as the selfconsistent equations of motion (12) and (13) are fulfilled. The binding energy per nucleon (for T = 0) as a function of ρ B has been the subject of extensive studies for decades and is not so well known for densities above about 3 ρ 0 . Note, however, that for a baryon density of 3 ρ 0 one obtains an energy density ∼ 0.5 GeV/fm 3 , which corresponds to the critical energy density in case of µ B = 0. In this work we will use the binding energy per nucleon from Ref. [69], which is consistent with microscopic Dirac-Brueckner calculations and the experimentally known momentum dependence of the nucleonnucleus optical potential (full black line in Fig. 3). We recall that the covariant approach in Ref. [69] is energymomentum conserving and most importantly also thermodynamically consistent. The conventional HRG result is shown by the red dash-dotted line in Fig. 3 and shows no binding as pointed out before. A suitable parametrization of relativistic mean-field theory, which reproduces the binding energy for densities up to 0.6 fm −3 , is given by: with the meson selfenergies given by The couplings here are taken as independent from the nucleon fields and the rearrangement selfenergies thus vanish. The minimum of this binding energy is We show in Fig. 4 the pressure from the equation of state from Ref. [69] (black solid line) together with other parametrizations for relativistic mean-field models taken from Refs. [70] (NL1 and NL3), [71] (TM1) and [72] (MTEC). The sets TM1 and MTEC are both stiff enough to allow for neutron stars with two solar masses (when extended by the isovector ρ exchange) as shown in Ref. [40].
We close the discussion of the nuclear equation of state with a short remark on the thermodynamic potential. So far we have shown the equation of state at T = 0 as a function of the density ρ B and not of the chemical potential µ B since it is more convenient to describe the system in terms of the Fermi momentum p F (ρ B ∼ p 3 F ). It is therefore natural to express all thermodynamic quantities as a function of ρ B , which implies, that we are not working in a grand-canonical ensemble but in a canonical ensemble. In this case the thermodynamic potential is the free energy, F = U − T S, which is equal to the internal energy U for T = 0. The pressure of the system in this special case follows from the relation, and is no longer proportional to the thermodynamic potential (as in case of the grand-canonical ensemble). As one can see from Eq. (24) this implies also negative pressures if the binding energy per nucleon decreases with ρ B , i.e. for all densities below the saturation density ρ 0 (cf. Fig. 4). The system for ρ B < ρ 0 is unstable, since the grand-canonical thermodynamic potential is larger than the vacuum. Another consequence of the canonical nature of the nuclear equation of state is the uncertainty in the chemical potential, which is no longer a natural variable. The system is described by the effective chemical potential µ * (17), i.e. the real chemical potential gets shifted by the repulsive interaction. It might therefore be impossible to define the nuclear equation of state as a function of the chemical potential in a unique manner.

IV. INTERACTING HADRON-RESONANCE GAS (IHRG)
In this Sec. we introduce interactions in the covariant hadronic model by constraining the Lagrange density at vanishing chemical potential (µ B = 0) by the most recent equation of state from the Wuppertal-Budapest collaboration [24] and at vanishing temperature (T = 0) by the nuclear equation of state from Ref. [69]. At finite temperature mesons will appear and interact with other mesons and baryons through resonant scatterings, which we describe in terms of the HRG by including several important resonances as non-interacting particles. We restrict the particles here to the most basic hadrons summarized by the list HRG-0 in Sec. II. We, furthermore, incorporate meson-exchange interactions in terms of relativistic mean-field models, which introduces additional attractive interactions as mediated by the σ-meson, which will account for the missing higher resonances in the HRG. The σ and the ω-meson appear in this model also as non-interacting particles, which may seem as a double counting, but the non-interacting contribution plays the role of an s-channel resonant-scattering amplitude. This channel -which is neglected in the mean-field limit -is missing in the meson-exchange model and the "exchange particles" appear only in the t channel. The nucleons appear in both channels, but in their case it is important to omit the non-interacting contribution in the HRG to avoid a true double counting. The thermodynamic potential of the interacting Hadron-Resonance Gas (IHRG) then is defined by the sum of the relativistic mean-field model and the regular HRG without nucleons, All other thermodynamical quantities, which follow from Ω IHRG by differentiation, then are also just the sum of the HRG and RMF contributions, but without the noninteracting nucleons. For the description of the nuclear equation of state it is not necessary to extend the model towards more interacting particles, since additional baryons will only appear for very large densities ρ B ≥ 2 − 3 ρ 0 [73][74][75]. This changes at finite temperature where other baryons populate the system and start to interact via meson exchange. However, the interaction has to be tuned in such a way, that it does not change the nuclear equation of state. The inclusion of additional interacting baryons such as hyperons and ∆'s is a frequently discussed question in the context of neutron-star physics, see Ref. [75] and references therein. We will use the findings from this field to extend the IHRG to include further interacting baryons while the mesons are kept non-interacting. Especially important in this context -and also for the description of heavy-ion collisions -are the ∆-resonances, which we describe by the Lagrangian [73,74,76], which is added to the Lagrangian of the relativistic meanfield theory (8). The couplings Γ ∆ may depend on an arbitrary Lorentz scalar or stay constant. The spinor Ψ ν ∆ is not a Dirac spinor but a Rarita-Schwinger spinor with 4 × 4 components that describes a spin-3/2 particle [77], however, the mean-field limit of the theory behaves just like Dirac spinors [73,74]. The selfconsistent equations (12) and (13) become Here ρ ∆ s and ρ ∆ B are the scalar and the particle density for non-interacting ∆-baryons. They depend on the effective mass m * ∆ and the effective chemical potential µ * ∆ that are defined by the selfenergies of the ∆'s, The actual form of the rearrangement selfenergies follows from the density-dependence of the couplings Γ σ∆ and Γ ω∆ . The pressure and the energy density of the system -without the HRG contribution -are given by and The entropy and the particle density are simply given by the non-interacting expressions, but with the respective effective quantities µ * ∆ , m * ∆ , The approach is thermodynamically consistent if the selfconsistent equations (27) and (28) are fulfilled. The thermodynamic potential of the IHRG with interacting nucleons and ∆'s is where Ω RMF is now the relativistic mean-field theory with nucleons and ∆-baryons. This extension introduces two additional couplings Γ σ∆ and Γ ω∆ . As for nucleons these couplings are not fixed by theory, but one can impose several constraints. The introduction of additional particles such as ∆'s or hyperons can create a second minimum in the binding energy [73,74], but since there are no ∆'s in the ground state of nuclear matter, this minimum can only describe a metastable state. Furthermore, any contribution from the ∆'s has to vanish at saturation density. There is also some guidance from finite density sum-rules, which show that the scalar selfenergy of the ∆'s is larger and the vector selfenergy smaller than the corresponding values for the nucleon selfenergies [78]. In Ref. [76] all these conditions are used to constrain the model in case of constant couplings. These findings are summarized by: (36) A simple choice for the couplings -in line with the relations (36) -are the conditions Γ σ∆ /Γ σN = m ∆ /m N and Γ ω∆ = Γ ωN . They are based on the argument that the ω-meson has a real quark-antiquark structure and the σmeson does not [74]. This choice leads to a fixed ratio of the effective masses m * ∆ /m * N = m ∆ /m N and equal chemical potentials for both baryons, µ * ∆ = µ * N . We will employ this choice whenever we treat the ∆'s as interacting particles.
The generalization to even more interacting baryons is straight forward. Spin-1/2 and spin-3/2 particles behave equally in the mean-field limit. We fix the scalar couplings by the ratio of the bare masses and keep the vector couplings identical, The selfconsistent equations (27) and (28) in their generalized form become, and the pressure reads The sum runs over all baryons that we include as interacting particles. These baryons have to be omitted in the HRG contribution in case of the IHRG. We will only discuss the cases of interacting nucleons as well as interacting nucleons and ∆'s in the following. Another reasonable choice are all baryons in the spin-1/2 octet. However, the results are similar to the case of interacting nucleons and ∆'s since the masses of the ∆, the Σ, Λ and Ξ are all in the vicinity of m ≈ 1200 MeV. The ∆'s are 16-times degenerated, the Σ's, Λ's and Ξ's have in total a degeneracy of 12, thus both cases are fairly similar, but the ∆'s are more important in low-energy heavy-ion collisions due to the reactions π + N ↔ ∆.
We will now fix the parameters of the IHRG. The right-hand side of the selfconsistent equation (39) is proportional to the net-baryon densities of the interacting baryons, which have to vanish for µ = 0. Since symmetries demand that O(ω) is an even function, the left-hand side of the equation vanishes for ω = 0. This fixes ω = 0 for µ = 0 and the repulsive interaction contributes only at finite chemical potential. As noted before this is different from repulsive interactions introduced through excluded volume effects that contribute even at vanishing chemical potential [46][47][48][49]. We can therefore fix the scalar interaction solely with the lQCD equation of state at µ B = 0 and then tune the repulsive interaction to reproduce the nuclear equation of state at T = 0 and µ B = 0.
We use the following strategy to define the scalar interaction. We subtract the non-interacting HRG from the lQCD equation of state and define in this way the contribution from the attractive mesonic interactions. We keep the scalar coupling as a constant Γ σN , and obtain from Eq. (37) a constant ratio for the effective masses m * X /m * N = m X /m N . The entropy density of the interacting model for a given temperature T is then a function of only the effective nucleon mass m * N , We demand that the interacting entropy density s Int is equal to the missing entropy density (in the HRG) to reproduce the lQCD result. This determines the effective mass m * N (T, µ B = 0). With m * N (T ) fixed we can easily calculate the scalar densities and use the selfconsistent equation (38) to determine ∂U/∂σ(T ) as a function of temperature. The value of the σ-field as a function of temperature follows from the effective mass relation σ = (m N − m * N )/Γ σN . We can, furthermore, fit ∂U/∂σ as a function of σ and by integration define the σ-selfinteraction. The polynomial ansatz for the selfinteraction (Eq. (22)) is able to reproduce the interaction for both nucleons as well as nucleons and ∆'s. The value of the scalar coupling Γ σN is arbitrary, since σ has no physical meaning, only Γ σN σ = m N − m * N . If one rewrites the selfconsistent equations (12) and (38) in terms of m * instead of σ, one finds for the polynomial ansatz of U (σ) that the equation is determined by the ratios m σ /Γ σN , B/Γ 3 σN and C/Γ 4 σN . We fix Γ σN by setting the σ-mass to its physical value m σ ≈ 550 MeV. The parameters in Tab. I give a good representation of the scalar selfinteraction for the temperature range between T ≈ 130 MeV and T ≈ 160 MeV and define the attractive interaction of the IHRG at µ B = 0. We denote the set for interacting nucleons by 'NLDD1' and for nucleons and ∆'s by 'NLDD2'. When comparing these parameters with those for other relativistic mean-field models one notices the large quartic coefficient C. Whereas the scalar selfinteraction in conventional mean-field models is just a small correction, it gives the dominant contribution in our approach. Another difference to conventional meanfield models is the much larger scalar density ρ s probed by the approach, since ρ s increases with ρ s ∼ T 3 for low µ B . This may lead to an unphysical phase transition if ∂U/∂σ is not strictly increasing monotonically, i.e. if the cubic or quartic coefficients B or C are negative. Note that the lQCD equation of state has no real phase transition and just a crossover, such that any model employed to describe lQCD cannot show a critical behavior for µ B = 0. Both parametrizations, NLDD1 and NLDD2, have negative cubic interactions, but the dominant quartic interaction is positive such that the model is regular. In the same way -as a non-monotonic behavior in ∂U/∂σ can introduce a phase transition -also a non-monotonic behavior in ∂O/∂ω potentially introduces one.
We show in Fig. 5 the ratio of the effective masses (for the nucleon and ∆) to the vacuum masses as a function of the temperature for vanishing chemical potential. The additional interactions do not show up for temperatures below T ≈ 100 MeV, such that the effective masses stay at their vacuum values before decreasing with increasing temperature. We note that a smaller effective mass The lQCD results are taken from the Wuppertal-Budapest collaboration [24].
results in a larger σ-field and more interaction strength compared to the non-interacting case. The effective mass for the parameter set NLDD1 decreases more rapidly than for NLDD2, because the whole additional interaction strength has to come from the nucleons alone while for NLDD2 also the ∆-contribution is included. We show the corresponding scaled entropy densities s/T 3 in Fig. 6 and compare them to the non-interacting HRG with the same degrees of freedom, HRG-0, and the lQCD entropy density from Ref. [24] used to determine the attractive interaction. At small temperatures the interacting models are similar to the non-interacting HRG since the interactions give no contribution. The additional interaction becomes visible for T ≈ 125 MeV. Up to temperatures of T ≈ 155 MeV both interacting models (by design) give the same result and describe the lQCD data within the error bars, however, differ at higher temperatures where one expects partonic matter anyhow. The IHRG with ∆'s increases too fast for higher temperatures and exceeds the lQCD entropy while the model with only nucleons reproduces the entropy density even up to T = 200 MeV. This is surprising since we fixed the interaction only for smaller temperatures. Nevertheless, the IHRG can not (and should not) describe the dynamics for temperatures beyond T ≈ 160 MeV, since it does not use the proper degrees of freedom; however, both models work well in the region where we expect a hadronic medium.
We now compare the equation of state from the two parametrizations with the lattice data from Ref. [24] in Fig. 7 (for NLDD1) and Fig. 8 (for NLDD2) scaled by powers of the temperature. We find an excellent agreement between lQCD and the model NLDD1, which de- scribes the whole equation of state within the error bars of the data, even at temperatures T > T c where lQCD becomes more reliable and the error bars shrink. The thermodynamic consistency of the approach ensures that we get the correct behavior in the pressure and the energy density once the entropy density is fixed. The first differences will appear only for T > 200 MeV. The interaction The dotted magenta line line shows c 2 s for the non-interacting HRG with all the hadrons in the IHRG (HRG-0) and the dashed blue line for all hadrons listed by the particle data group [60] with a mass below 2.6 GeV (HRG-1). The lQCD results are taken from the Wuppertal-Budapest collaboration [24].
measure I/T 4 has its maximum around this temperature and will then start to decrease in lQCD, but increases further in the hadronic model. The IHRG equation of state then exceeds the lQCD data. As mentioned earlier the parameter sets NLDD1 and NLDD2 give the same results for temperatures below T = 155 MeV and the model reproduces the lQCD data within the error bars for all temperatures below T ≈ 160 MeV, where we expect a dominantly hadronic system (cf. Fig. 8). At larger temperatures the results from NLDD2 differ substantially from the lQCD results and rise too fast for T > T c as already seen in Fig. 6. One might improve the description at larger temperatures by using a different parametrization for the scalar selfinteraction and employ a larger temperature interval, but this is of no relevance for the present study where we have fixed the interaction only up to T = 160 MeV. The excellent results from the parameter set NLDD1 at larger temperatures come out as a surprise.
We have seen for the case of the conventional HRG that a reasonable reproduction of the equation of state does not imply a correct behavior in the speed of sound, cf. Fig. 2. We now compare the speed of sound squared for the IHRG as a function of the temperature in Fig.  9 and show again the corresponding results from the non-interacting HRG (HRG-0) and also the HRG with hadrons up to a mass of 2.6 GeV (HRG-1). We find that only the parameter set NLDD1 describes the data properly. It reproduces the minimum at T ≈ 150 MeV and is within the error bars up to T = 180 MeV. Nevertheless, it benefits from the huge error bars at low temperatures. The parameter set NLDD2 can only describe the data up to T ≈ 150 MeV; it has also a minimum in c 2 s , but at a too high temperature, which is also too deep. On the other hand HRG-0 is completely off the data and has also the wrong T -dependence. The other version, HRG-1, with much more particles gives a better description but only up to T ≈ 155 MeV. From there on it has also the wrong slope. All models except for NLDD1 fail to describe the rise in the speed of sound at T ≈ 150 MeV. For this it is necessary that the models reproduce the equation of state up to the inflection points of the scaled equation of state. Thus one will always find a decreasing speed of sound in T if the scaled pressure P/T 4 has an increasing slope, which is a problem in most hadronic models.
With the scalar interaction defined at µ B = 0 we can now discuss the repulsive interaction in addition. We fix it in the same way as the scalar interaction using the nuclear equation of state at T = 0 as input. In this limit the HRG contribution of the IHRG vanishes and it reduces to a normal (density-dependent) relativistic meanfield model. The scalar interaction defines already the effective masses as a function of the density at T = 0 and therefore the selfinteraction U (σ), the scalar selfenergy Σ s and the non-interacting part E 0 of the energy density of the relativistic mean-field model (20). The remaining contributions to the energy density depend on the repulsive interaction that we determine as follows: We omit the selfinteractions in the ω-field and keep only the mass term as in Eq. (23), however, we describe the repulsive interaction with a density-dependent vector coupling, which depends on the net-baryon density, Γ ωN (ρ B ). Note that it is important for the consistency of the model that ρ B contains only the interacting baryons and not the whole baryon density of the IHRG; however, this is naturally the case for the nuclear equation of state. The coupling can not depend on the scalar density, since this would lead to a scalar rearrangement selfenergy that alters the effective mass and therefore the equation of state at finite temperatures and vanishing chemical potentials. The selfconsistent equations for the ω-field (13), (28) and (39) simplify to m 2 ω ω = Γ ωN (ρ B )ρ B . We use this form of the ω-selfinteraction to rewrite the equation for the energy density (20) as, The right-hand side of the equation depends on the repulsive interaction but the left-hand side is determined by the scalar interaction and the equation of state that one wishes to reproduce. Since both are already defined, we can use Eq. (42) to determine the vector coupling Γ ωN . Note that we describe the equation of state at T = 0 as a canonical system, where the density is a natural variable (instead of the chemical potential) and the chemical potential follows as a derivative of the thermodynamic potential with respect to density. Therefore, Eq. (42) defines the coupling directly as a function of the density ρ B , where the repulsive interaction is actually determined by the ratio Γ ωN /m ω . This is similar to the scalar interaction where we fixed m σ to its physical value. We do the same here and fix m ω = 783 MeV. With Eq. (43) one can employ any possible nuclear equation of state. Even if it is possible to directly use the function in the numerical calculations, it is convenient to have a parametrized form. We will use the ansatz This function is similar to the ansatz employed in Refs. [67,68] to fit the density-dependence of Dirac-Brueckner calculations, but with two important differences: The model in Refs. [67,68] was only applied to nuclear matter and finite nuclei, i.e. to T = 0, but we want to employ our model also at finite temperatures and vanishing chemical potential. It is therefore mandatory that the coupling is an even function of the density, which is guaranteed by the absolute values of the density that lead to Γ ωN (ρ B ) = Γ ωN (−ρ B ), but the function has to be continuously differentiable at ρ B = 0. We ensure this by taking the same linear coefficient in the numerator and denominator of Eq. (44). To compensate for the lost coefficient in the linear term we extend the polynomials in the function to third order. The coefficients for the fit are summarized in Tab. II and we show the couplings as a function of ρ B in Fig. 10. It might not be possible to see by eye, but both functions fulfill the condition Γ ′ ωN (ρ B = 0) = 0. The differences in the couplings are not due to the appearance of ∆'s but follow from the different scalar interactions. Since the parameter set NLDD1 has a stronger scalar, i.e. attractive interaction, it needs a stronger repulsive interaction to balance it.
We can now calculate the nuclear equation of state of the IHRG to check the quality of the fit (44). The binding energy per nucleon E B /A is shown in Fig. 11 in Even if their mass gets reduced due to the scalar interaction, they are still too heavy to give a thermodynamic contribution. The same holds also if one includes even more interacting hadrons, since they would additionally weaken the scalar interaction. The nuclear equation of state of the IHRG is therefore predominantly determined by the nucleons. Fig. 11 shows some odd behaviors at very small densities: NLDD1 has a very small peak and NLDD2 exhibits a second minimum at ρ B = 0.012 fm −3 that is higher than the global minimum at saturation density implying that it describes just a metastable state. These deviations can appear also in conventional relativistic mean-field models, but at even smaller densities. They occur due to the singular behavior of the binding energy per nucleon E B /A = E/ρ B − m N , if the energy density differs from the form E ≈ m N · ρ B at small densities. In case of the IHRG they result from the condition Γ ′ ωN (ρ B = 0) = 0 that constrains the functional form of Γ ωN (ρ B ). However, the effects are small and negligible if one considers systems with energy densities above 20 MeV/fm 3 .
So far we have shown that the IHRG gives reasonable results for µ B = 0 along the T -axis as well as for low T and finite nuclear density. However, the IHRG is defined in the whole (T, µ B ) plane and we may explore the phase diagram e.g. for constant thermodynamical potential, which for a grand-canonical ensemble is the (negative) pressure P . Alternatively, for T → 0 the thermodynamic potential is the energy density E. We show in Fig. 12 the lines of constant pressure P = 63 MeV/fm 3 and of constant energy density E = 0.4 GeV/fm 3 , which are roughly the values of the lQCD equation of state for T c ≈ 155 MeV and vanishing chemical potential. Accordingly, the lines may be interpreted as a QCD phase boundary. Both conditions give similar results but the lines of constant pressure reach further out into the (T, µ B ) plane. When addressing the phase boundary in the context of heavy-ion collisions, an important constraint is strangeness neutrality. In Fig.  12 (a) the lines of constant pressure or energy density were obtained without any constraint on the strange sector, such that we have a finite strangeness N S > 0, while we show in Fig. 12 (b) the lines for a strange neutral medium. However, neither strangeness neutrality nor the two different parametrizations of the IHRG, NLDD1 and NLDD2, have a strong impact on the phase boundaries, which agree close to the axis of the (T, µ B ) plane and cover almost the same area in the phase diagram.
We note in passing that a large uncertainty in the phase boundary stems from the transition at low temperatures (T ≈ 0) and large baryon chemical potentials. As mentioned earlier the nuclear equation of state is known only as a function of the density ρ B , but not of the chemical potential µ B . As mentioned above the strength of the repulsive interaction between the nucleons has a large influence on µ B since the thermodynamics are defined by the effective baryon chemical potential µ * that is shifted by the vector selfenergy (17). Thus models with different interactions can reproduce the same nuclear equation of state as a function of density ρ B , but for different chemical potentials. The grey band at T ≈ 0 in Fig. 12 indicates the uncertainty for small temperatures. The lowest chemical potential is the estimate for a gas of non-interacting nucleons, the largest for the parameter set NL3 from Ref. [70], which has a very strong repulsive interaction but is consistent with the ground-state properties of nuclear matter.
To quantize the effects of strangeness neutrality we show in Fig. 13 the lines of constant energy density calculated for NLDD1 with (solid green line) and without (dashed red line) strangeness neutrality. The lines agree close to the axis of the (T, µ B )-plane because one has the same amount of strange and antistrange particles at µ B = 0 and no strange particles at all at T = 0. These regions of the phase diagram are always strange neutral and therefore have to agree independent of the strange sector while for T = 0 = µ B we find that the phase boundaries for a strange neutral medium are shifted further into the (T, µ B )-plane because strangeness neutrality reduces the amount of particles in the strange sector compared to Lines of constant energy density from the IHRG in the (T, µB)-plane for the parameter set NLDD1 with (solid green) and without (dashed red) strangeness neutrality. The grey band (at T ≈ 0) highlights the ambiguity of the chemical potential at zero temperature for relativistic mean-field theories with different repulsive interactions. The squares are freeze-out points taken from Ref. [79] for 0-5% (red), 30-40% (green) and 60-80% (blue) centrality. The triangles are freeze-out points taken from Ref. [80].
an unconstrained system. This lowers the energy density and the pressure and one needs larger temperatures and baryon chemical potentials to reach the same energy densities and pressures as in an unconstrained medium. We also compare our predictions for the phase boundary to recent results for freeze-out points from the beamenergy-scan (BES) program at RHIC [79] and from the HADES collaboration [80]. The freeze-out data from BES are all at small chemical potentials and located around the predicted phase boundary, the freeze-out data from HADES are at larger baryon chemical potentials and lower than the IHRG boundary. This is expected since in heavy-ion collisions below about 2 AGeV no 'droplets' of QGP can be created and the system entirely stays in the hadronic phase [42].

V. SUMMARY
In this work we have discussed QCD thermodynamics in the regime of hadronic degrees of freedom. At finite temperature T and vanishing chemical potential µ B the interactions are dominated by resonant scatterings and the thermodynamics can well be described by a HRG model where one includes the resonances as non-interacting particles in the partition sum. The approach can reproduce the equation of state, as provided by lQCD, but requires a large amount of interaction strength, i.e. hadronic resonances. The HRG describes only the attractive interactions between the hadrons while short-range repulsive interactions, that are important at large densities, have to be introduced by means of excluded-volume models. These interactions, however, are only thermodynamically motivated and not based on field theory. Hence they are in general not covariant and thus can not be used in transport approaches for non-equilibrium configurations.
A covariant formulation for repulsive (short-range) interactions is provided by relativistic mean-field theories (RMFTs). These approaches describe the interactions by meson exchange and characterize hadronic matter at low temperatures. Consequently, the model is a well-known and suitable approach for nuclear matter at finite density and successful in describing the ground-state properties of the nuclear equation of state. In particular, we have discussed a relativistic mean-field approach for symmetric nuclear matter with density-dependent couplings. This extension allows for a parametrization of realistic nucleon-nucleon interactions -e.g. from Dirac-Brueckner calculations -and introduces a non-trivial density dependence in the scalar and vector interactions.
By combining both hadronic models, the HRG and RMFT, we have constructed an interacting HRG, denoted by IHRG, that reproduces the lQCD equation of state at µ B ≈ 0, T > 0 and the nuclear equation of state at T ≈ 0, µ B > 0. The s-channel interactions are described -as in the HRG -by the inclusion of resonances as non-interacting particles; the t channel, i.e. mesonexchange reactions, are treated by means of RMFT. The repulsive interactions in the IHRG are therefore fundamentally different from those incorporated in excludedvolume models, which are proportional to the total particle number or the pressure. In the IHRG the repulsive interaction is controlled by the net-baryon density and vanishes for µ B = 0. The attractive interactions are realized by different mechanisms and it is therefore not necessary to include a huge amount of hadronic resonances in the IHRG since one can focus on well known hadronic states and resonances; the missing interaction strength is provided by the σ-meson exchange.
We have presented two parametrizations of the IHRG: In the first, denoted by NLDD1, we have only considered the meson-exchange reactions for nucleons. In the second, denoted by NLDD2, we have included additionally the ∆-resonances, which are important for the dynamics of heavy-ion collisions in the pion-nucleon channel. We have defined the interactions of the ∆'s in such a way that the ratio of the effective masses is given by the ratio of the vacuum masses but the effective chemical potentials for nucleons and ∆'s are the same. Using this approximation one can easily extend the IHRG to many interacting baryons, however, the mesons are kept non-interacting.
Since the IHRG can describe the nuclear equation of state at T = 0 and finite density, it should give a better description of the hadronic medium at low temperatures and finite density than the standard HRG. This region of the phase diagram is of particular interest for heavyion collisions at low beam energies that investigate the phase boundary of QCD at large baryon chemical potentials and search for a possible critical point in the phase diagram. We have given an estimate for the QCD phase boundary between the QGP and hadronic configurations based on the assumption that the phase transition occurs at a constant energy density or pressure. We have found that both conditions lead to almost the same transition region for both sets, NLDD1 and NLDD2. The predicted boundary is consistent with the freeze-out analysis from the BES that probed the hot QCD medium at moderate chemical potentials. We, furthermore, note that the phase boundary has a large uncertainty at low temperatures since the nuclear equation of state is only known as a function of the density (canonical ensemble) and not the chemical potential (grand-canonical ensemble).