Equation of state of neutron star matter and its warm extension with an interacting hadron resonance gas

We propose an interpolating equation of state that satisfies phenomenologically established boundary conditions in two extreme regimes at high temperature and low baryon density and at low temperature and high baryon density. We confirm that the hadron resonance gas model with the Carnahan-Starling excluded volume effect can reasonably fit the empirical equation of state at high density up to several times the normal nuclear density. We identify the onsets of strange particles and quantify the strangeness contents in dense matter. We finally discuss the finite temperature effects and estimate the thermal index $\Gamma_{\rm th}$ as a function of the baryon density, which should be a crucial input for the core-collapse supernova and the binary neutron star merger simulations.


Introduction
Strongly interacting matter has been investigated well in two extremes of the temperature T and baryochemical potential µ B . The heavy-ion collision experiments realize matter at high-T and low-µ B , which can also be studied by the Monte-Carlo simulation of quantum chromodynamics (QCD) on the lattice [1,2]. Another extreme matter at low-T and high-µ B is found in cores of neutron stars, for which the Monte-Carlo lattice-QCD simulation is of no use due to the notorious sign problem.
Confronting the lattice-QCD thermodynamics, the hadron resonance gas (HRG) model has become a common tool to approximate the equation of state (EoS) below the critical temperature T c ∼ 150 MeV (see, e.g., Ref. [3] for a recent review). The HRG model is a parameter-free description of thermal properties with experimentally observed particle and resonance spectra. In the HRG model, nonperturbative hadronic interactions are assumed to be dominated by poles and branch cuts corresponding to resonances [4]. The ideal HRG (IHRG) model consists of noninteracting hadrons and resonances, and yet, its validity has been endorsed by the successful thermal fit of particle abundances in heavy-ion collision experiments. The IHRG EoS is recently extended to the terrain of finite baryochemical potential in a region of µ B /3 T [5], guided Email addresses: fujimoto@nt.phys.s.u-tokyo.ac.jp (Yuki Fujimoto), fuku@nt.phys.s.u-tokyo.ac.jp (Kenji Fukushima), hidaka@post.kek.jp (Yoshimasa Hidaka), a.hiraguchi@nycu.edu.tw (Atsuki Hiraguchi), iida@kochi-u.ac.jp (Kei Iida) by the conserved charge susceptibilities from lattice-QCD calculations [6][7][8]. This region may include a possibility of Quarkyonic matter [9][10][11][12][13] and the associated triple-point structure in the QCD phase diagram [14]. Furthermore, owing to its parameter-free nature, it can also be applied to a wide variety of problems such as a reference for baryon number fluctuations [15], a shift of the chemical freezeout due to the inverse magnetic catalysis [16], and the rotational effect on the deconfinement temperature [17], to mention a few.
Thermodynamic quantities from the IHRG model, however, blow up above T c . We expect that the validity region of the HRG model could be extended to higher temperature by introducing the interaction effect. The excluded volume (EV) effect is the simplest way to implement the interacting HRG and the formulation was given in Ref. [18]; the short-range repulsive interaction was modeled via a hard-core correction following the thermodynamically consistent way as developed in Ref. [19] (see also Ref. [20]). By supplementing the repulsive interaction with an additional attractive interaction, the interacting HRG model amounts to the van der Waals (VDW) EoS [21,22] (see also Ref. [23] and references therein). Albeit a few known theoretical problems [24], this crude model approach is successful in describing the lattice-QCD results including the EoS and the cumulants (see, e.g., Refs. [18,21]). Indeed the nuclear interaction has an intricate structure, but the nuclear matter properties can be captured by two parameters (see, e.g., Refs. [25,26] for the two-parameter nuclear matter description by the scattering length and the effective range).
Up to here, we have reviewed the interacting HRG model in the context of high-T and low-µ B studies. Then, it is a natural anticipation that the same machinery of the interacting HRG model can also work at low-T for astrophysical applications; interestingly enough, this idea of applying the HRG description for neutron stars can be traced back to the pioneering work by Hagedorn [27] more than half a century ago even before the discovery of pulsars. One might think that this anticipation is immediately falsified by the nuclear liquid-gas phase transition of first order in the low-T and high-µ B region. Although the IHRG model is unable to describe the first-order phase transition, the interacting VDW-HRG model can properly account for the liquid-gas phase transition as studied in Refs. [21,28].
In the present work, we will look into the EoS from the VDW-HRG model by making a quantitative comparison with the Chiral Effective Field Theory (χEFT). We will show that the hard-core EV effect violates the causality, but the VDW-HRG model with refined repulsive interaction leads to more reasonable behavior at low T and high µ B . The modified repulsive interaction is incorporatedà la Carnahan-Starling (CS) [29]. The CS excluded volume has been adopted in the hadron physics [30][31][32], and astrophysics [33,34]. Within the HRG model with the CS-type EV, i.e., the CS-HRG model, we will construct the EoS, p(n B ) (the pressure vs. the baryon number density), for dense and warm matter at T < 60 MeV. Even though the parameters in the CS-HRG model are fixed by the χEFT, the CS-HRG can provide us with useful insights. One example is that we can diagnose the strangeness contents in dense matter hadron by hadron. Another useful application is the finite-T extension of the EoS; in many simulations of binary neutron star mergers the thermal component of the EoS is treated as an ideal gas for convenience [35,36], and there is a recent discussion in the astrophysical context [37] that thermal pions could be important.
At the conceptual level, the interacting HRG model embodies the idea of Quarkyonic matter that claims continuous duality between baryonic and quark matter (see, however, Ref. [38] for the earlier analysis based on the VDW model indicating the large-N c transition between N c = 3 and N c → ∞). The idea is consonant to the crossover EoS construction based on quark-hadron continuity [39][40][41]. There are theoretical attempts to build a Quarkyonic model [42][43][44][45][46][47][48][49], and as argued in a quantum percolation picture [50] quark degrees of freedom emerge from interactions. In this sense, the interacting HRG model could be regarded as a concrete modeling that exhibits the quark-hadron duality from the hadronic side. This theoretical argument also justifies the extended validity of the interacting hadronic model for high-T and/or high-µ B matter in which quark degrees of freedom could be mixed together.
We note that, throughout this work, we use natural units; c = = k B = 1.

Interacting hadron resonance gas model
The HRG encompasses the numerous contributions from experimentally measured states of hadrons and resonances. The total thermodynamic quantities, such as the pressure, can be decomposed into three pieces: where M , B, andB denote contributions from mesons, baryons, and anti-baryons, respectively, and µ = (µ B , µ Q , µ S ) are the chemical potentials conjugate to the net baryon number B, the electric charge Q, and the strangeness S. We here limit ourselves to B = 1 baryons and B = −1 anti-baryons, and we discard composite baryons with B > 1 such as the deuteron, light nuclei, hypernuclei, etc.; the deuteron, for example becomes unbound in neutron-rich matter. Moreover, since matter itself is a gigantic nucleus (if it is self-bound), light cluster contributions may lead to double-counting problems. For our HRG model in this work, we have adopted the particle data group list of particles and resonances (where the resonances are handled in the zero-width approximation; see Ref. [51] for a refined S-matrix treatment of resonances) contained in the THERMUS-V3.0 package [52].

Ideal hadron resonance gas (IHRG)
In the IHRG model, each particle is treated as the Bose/Fermi ideal (i.e., non-interacting) gas. In this model, the mesonic, the baryonic, and the anti-baryonic contributions in Eq. (1) are given by where the pressure function p id i in the grand canonical ensemble is We note that the overall + (−) sign corresponds to the fermion (boson). The energy dispersion is with the spin degeneracy factor, g i = 2s i +1, and the chemical potential, The other thermodynamic quantities such as the number density n i and the energy density ε i can also be derived accordingly as For cold matter at T = 0 we can analytically carry out the fermion integrals (and mesons are irrelevant) to find the pressure, , and the number and the energy densities are n id

Van der Waals hadron resonance gas (VDW-HRG)
In this work, we incorporate the interaction effect based on the VDW construction, which comprises the repulsive interaction by the EV effect and the attractive interaction. The VDW EoS was originally formulated in classical systems in the canonical ensemble with the fixed number of particles, and thus the following two extensions were necessary: the reformulation in the grand canonical ensemble was given in Ref. [28], and the quantum statistics was taken into account in Refs. [53,54]. Hereafter, we will employ the formulation in Ref. [55]. At this point, we shall make a comment on the interactions. We will consider the interaction effect only in the B sector. At zero and low temperatures, theB sector is simply negligible. The reason why we neglect the interaction in the M sector is that mesons are only weakly interacting in the limit of large colors. Indeed, the substantial mesonic EV effect leads to the discrepancy between the EV-HRG and the lattice-QCD data due to too strong suppression of thermodynamic quantities, which was already implied in Figure 4 of Ref. [18].
To introduce the repulsive interaction through the EV effect, we replace the volume V in the partition function by f (η)V with η being the packing fraction 1 . The function f (η) ∈ (0, 1] measures the volume fraction in which particles with hard sphere can move around. The concrete expression of f (η) will be given soon in Eqs. (11), and (12). The packing fraction, η, can be related to the hardsphere radius, R, as η = 4π 3 R 3 n (see Ref. [55] for more description).
The VDW-HRG model in the B sector can be obtained from the pressure defined by where the last term ∝ a represents the attractive interaction effect, and the shifted chemical potentialμ i iŝ 3 being the eigenvolume of baryons. From Eq. (8), the energy density is immediately derived as 1 Strictly speaking, the EV effect has to be introduced through the canonical partition function as formulated originally.
The physical meaning of this pressure expression should be transparent from the density, that is, the µ B -derivative of the pressure: For the conventional VDW model, the standard choice By definition the range of f VDW (η) must be limited within 0 < f VDW (η) ≤ 1, where f VDW = 1 refers to the no EV limit and f VDW = 0 corresponds to the maximal packing. When the maximal packing is reached at η = 1/4, the thermodynamic quantities become singular. Singularity in the thermodynamic quantities at the maximal packing leads to the acausal behavior in the speed of sound c s > 1, where the speed of sound is defined as c 2 s ≡ ∂p/∂ε. In the neutron star environment, using the hard-sphere radius R B = 0.511 fm (which will be determined in Sec. 4.1), we have a rough estimate for the maximal packing density as to n B 2.8 n 0 and the density at which the superluminal speed of sound (c s > 1) is reached as to n B 2.2 n 0 ; n 0 = 0.16 fm −3 is the normal nuclear density. These values are lower than the density in central cores of typical-mass neutron stars, and the above simple choice (11) is obviously inappropriate.

Carnahan-Starling (CS) refinement of the excluded volume term
In the VDW-type EV treatment, the interaction sphere is too rigid and the maximal packing occurs at unphysically low density. It would be sensible to smear the interaction clouds so that the EV effects can mildly set in. For the astrophysical application, therefore, we should improve the function (11) and a promising candidate is the Carnahan-Starling (CS)-type EV [29].
For the CS-HRG model [55], the choice of f (η) is When η 1, we can easily confirm that f CS (η) ≈ 1−4η = f VDW (η) up to the linear order in η. Here again, the range of f CS must be limited within 0 < f CS ≤ 1. The maximal packing occurs at η = 1. The corresponding density, n B = 4η/b = 3/(4πR 3 B ), then reaches n B = 11.2 n 0 with R B = 0.511 fm. Because of a longer tail up to η ∼ 1/2 in Eq. (12) as compared to a sharp drop at η = 1/4 in Eq. (11), the maximal packing density in the CS-type EV is about four times larger than that of the VDW-type EV, which extends the validity range. This is a naïve estimate, and the serious calculation concerning the causality gives a rather smaller value of the limiting density, i.e., n B 3.7 n 0 for R B = 0.511 fm.

Need for the attractive interaction
In the spirit of the HRG model, the attractive interaction leading to the resonance formation is implicitly incorporated through added resonances; however, not all the attractive interactions are taken into account by resonances. The nuclear force has attractive regions, which are dominated by one-pion exchange, as well as repulsive regions by the exchange of heavy mesons and multi-pions, as seen also in the first-principles lattice-QCD calculation [56]. The VDW model is not the direct description of the nuclear force itself but it can emulate such characteristics of the nuclear force; the EV effect captures the repulsive core nature at short range, while the intermediate and long-range parts are reasonably captured by resonances and the attractive interaction term in the VDW model.
One might wonder how that the attractive term ∝ a affects physical observables. Only to avoid singular behavior of thermodynamic quantities, the minimal model with either VDW-or CS-type EV would be enough. We would, however, stress that the attractive term is indispensable for quantitative analysis. We will come back to this point when we discuss the M -R relation using the EoS from the CS-HRG model later.

EoS construction
We are ready to construct the EoS for neutron stars at T = 0 as well as for compact-binary mergers and supernovae at T > 0. Under these circumstances, the system is β-equilibrated via the weak processes such as n → p + +ν , p + → n + ν , and neutral in the electric charge. These conditions of βequilibrium and electric charge neutrality amount to setting the values of the chemical potentials, µ Q and µ S .   Owing to large asymmetry of baryons and anti-baryons at high µ B , we can safely neglect the anti-baryonic contribution to the EoS, and we can also drop the mesonic contribution at T = 0. The meson condensation would cause subtlety; in particular, negatively charged pions would form a condensate when µ > m π − , where µ is the leptochemical potential and m π = 140 MeV is the pion mass. We assume no condensate and this will be justified selfconsistently as shown in Fig. 3. The values of µ at various temperatures in Fig. 3 do not exceed m π .
We impose the following conditions on µ Q and µ S : where n Q is the electric charge density in the interacting HRG model, whose expression is given by For the expression of the lepton density, n , we substitute the electron mass, m e = 0.511 MeV, and the muon mass, m µ = 106 MeV for the mass in the free particle expression (4): We neglected the neutrino contribution by assuming that neutrinos quickly escape in neutron star systems both for T = 0 and T > 0, which is justified for ordinary neutron stars. Now we can solve these conditions to find For the actual procedures, we fix the value of n B , and solve Eqs. (10)-(8), (15), and (16) in terms of three variables, µ B , µ Q (= −µ ), and ∆µ B . Once three variables are obtained, we can compute the thermodynamic quantities correspondingly. We note in passing that µ B and ∆µ B always appear in the combination ofμ B ≡ µ B + ∆µ B , so that we can separately treat Eq. (8) and solve two coupled equations to determineμ B and µ Q . With bearing the application to neutron stars in mind, we calculate the EoS at T = 0 following the procedures outlined above. We plot our main results for the EoS from the CS-HRG model in Fig. 1. The two parameters, R B and a, are fitted to reproduce the N 3 LO χEFT EoS, which is taken from Ref. [58] based on the pure neutron and symmetric nuclear matter calculations [59]. It is evident from the comparison between our results (by the thick curve) and the IHRG results (by the dashed curve) that the interaction is crucial for neutron star matter. We also observe that our fitted EoS stays close to the phenomenological nuclear EoS, for which we choose SLy4 [57]. It should be noted that the VDW-HRG model with f VDW (η) in Eq. (11) can by no means fit the N 3 LO χEFT nor the SLy4 EoS at all. In Fig. 1 there appears a bump around the n B = 1.8 n 0 ; it is related to the onset of hyperons as discussed later in Fig. 5. Here, too rapid stiffening is tamed by softening induced by the appearance of hyperons. The interacting HRG model has originally been used in the finite-T circumstances, so it should be a reasonable setup for generalizing the EoS construction to the finite-T case. In this way, the finite-T EoS has been obtained from the CS-HRG model and we plot the results at T = 0, 20, 30, 40, and 50 MeV in Fig. 2. Overall, the pressure becomes higher with increasing T as expected, but we find an exception; the EoS at T = 20 MeV goes below the T = 0 EoS around n B 1.8 n 0 . We can understand this as follows. At T = 0 new particles cannot appear until the density exceeds the mass threshold. Below the onset density, the EoS becomes stiffer as the density rises up. On the other hand, at finite T , the EoS can accommodate any massive particles according to the thermal weights, which makes a qualitative difference from the T = 0 EoS and explains the change of the ordering around n B 1.8 n 0 .
Finally, we show the leptochemical potential µ in Fig. 3. As we already mentioned before, µ never exceeds the pion mass for any temperature and density. This selfconsistently justifies our assumption of no meson condensate.

Astrophysical applications
We will showcase the numerical results for the EoS and the phenomenological implications at T = 0 and finite T in order. The density and temperature reached in supernovae are as high as n B ∼ 2-3 n 0 and T ∼ 100 MeV while those in neutron star mergers are n B 5 n 0 and T 50 MeV.

Neutron star matter and structure at T = 0
With the EoS from the CS-HRG model, we calculate the mass-radius (M -R) relation by solving the Tolman-Oppenheimer-Volkoff (TOV) equation. As shown in Fig. 4, our EoS incidentally matches with the phenomenologically accepted EoSs such as SLy4, even in the crust region of neutron stars, and so we use our EoS down to the surface of the star when solving the TOV equation; see Ref. [64]. We compare the cases with/without the attractive interaction. If we turn off the attractive interaction at a = 0, the radius at M = 1.4 M is too large and out of the NICER observation as seen in Fig. 4. Usually, the lower density part of the EoS is responsible for the radius of stars, and a larger radius is favored for a stiffer EoS at low density. This tendency can be confirmed in Fig. 1; the IHRG EoS is constructed without the attractive interaction, and the lower density part (n B /n 0 1) from the IHRG model exceeds that from the CS-HRG model. This behavior is consistent with such an interpretation that the EoS without attractive interaction gives too stiff EoS at low density and leads to a too large radius.   As mentioned earlier, for R B = 0.511 fm, the maximum packing density, above which the model breaks down, is n B = 11.2 n 0 for the CS-HRG model. On the M -R relation, the maximum mass of M = 2.56 M is attained at n B = 9.32 n 0 , which certainly lies within the validity range of the CS-HRG model. Regarding the maximum mass, some controversies are unavoidable. Combining the GW170817 event with the accompanying electromagnetic observation, the maximum mass could be constrained to be at most 2.3 M [65][66][67][68]. Meanwhile, a compact object with ∼ 2.6 M has been observed in the GW190814, which may be identified as a massive neutron star. Near the maximum mass region, another subtlety arises from a possible transition to quark matter [69,70]. It is a nontrivial question where the validity bound of our model should be. If we locate the validity bound at n B 3.7 n 0 (see explanations in Sec.   Fig. 2. For reference, the ab initio calculations for pure neutron matter [72] are overlaid in the grey color with the same line style as our results. The canonical value of the adiabatic index for nonrelativistic ideal gas, Γ = 5/3, is shown by a horizontal line.
to M 1.5 M .
The HRG model provides us with a convenient picture to probe the particle abundances. In Fig. 5 we show the fraction Y i = n i /n B of the particle species i. At small density, the neutron, n, is dominant with a small fraction of the proton, p, which slowly increases with increasing density. The onset of the hyperons is observed slightly below 2 n 0 . As is consistent with the conventional scenario (see, e.g., Sec. 5 of Ref. [71]), Σ − is activated first as we increase the density, and then Λ is produced afterwards.

Thermal index
In the applications for astrophysical phenomena such as supernovae and binary neutron star mergers, the thermal corrections to the EoS are often modeled by an ideal gas approximation [35,36,73]. In order to define the thermal part of the EoS, which is parametrized by the thermal index, Γ th , we introduce the rest-mass density of baryons as ρ B = m B n B with the nucleon mass, m B = 939 MeV. We can decompose the energy density ε as ε = (1 + e)ρ B , where e is the specific internal energy. We can add the thermal corrections to the pressure and the energy on top of the T = 0 parts as In the simulations of neutron star mergers, the cold (T = 0) component is used before shock heating associated with the stellar collision sets in. The relation between the thermal pressure and the energy should be supplemented with the additional constraint, that is commonly parametrized by In the phenomenological studies, Γ th is a free parameter. It is customary to choose Γ th around ∼ 1.7. For example, Γ th = 1.8 was adopted in Ref. [74]. If Γ th is too small (like ∼ 1.3), the thermal pressure is not large enough to sustain matter, resulting in a rapid proto-neutron star contraction for supernovae, and in a faster collapse of merger remnants to black holes for binary merger. In this way, smaller values of Γ th may have impact on core-collapse supernova and binary neutron star merger simulations [75]. In contrast, a larger Γ th (like ∼ 2.0) would elongate the life-time of the post-merger dynamics. Thus, for reliable theoretical predictions, it is of utmost importance to constrain Γ th . Moreover, although it is often assumed to be constant, Γ th may depend on the density and temperature [72].
We can make use of our EoS to infer Γ th , which can be represented in terms of the thermodynamic quantities as Here, ε th = ρ B e th is the thermal part of the energy density. In Fig. 6, we show our estimate for the thermal index, Γ th , as a function of the density. We find that Γ th becomes less sensitive to the density as T gets larger; e.g., Γ th is almost constant around ∼ 1.4 at T = 50 MeV. The preceding ab initio calculation based on χEFT [72] is also overlaid on Fig. 6. Our Γ th and the ab initio Γ th [72] differ qualitatively; it may be partially because the slope of the ab initio EoS at T = 0 is gentle compared with the stateof-the-art χEFT EoS [58,59] to which our model is fitted to. We note that Γ th in Fig. 6 is computed under the assumption of neutrinoless β-equilibrium. The values of Γ th shown in this figure may deviate from effective values realized in dynamical astrophysical phenomena, particularly when neutrinos are trapped inside hot neutron stars.

Summary and outlooks
We demonstrated a successful construction of the equation of state based on the van der Waals prescription of the hadron resonance gas model with the Carnahan-Starling refinement of the excluded volume term. The ideal hadron resonance gas description does not work at high density as shown in Fig. 1; at T = 0 both the repulsive and the attractive interactions are crucial. In our proposal a very simple parametrization can fit in with the empirical EoSs. This provides us with an intuitive picture of the quarkyonic regime as well as a paractically useful tool for phenomenological studies. It is then a straightforward extension to include the finite-T corrections. We quantified the thermal index, Γ th , as a function of the density for various temperatures. Our estimated Γ th favors a rather smaller value than the conventionally adopted one.
Our present modeling is the simplest one, and we can consider improvements in many respects. For example, the interaction should be species dependent, and in particular, the strangeness sector could behave differently. In some channels involving the strangeness, the Pauli blocking is relaxed and the short-range repulsive core is absent. As neutron stars contain hyperons as inevitable physical degrees of freedom, the flavor-dependent interactions should be considered along the line of Ref. [22]. It is also an interesting question to think of meson interaction effects. Moreover, since the validity region of our proposed model should cover the high-T and low-µ B regime, we can apply our predicted EoS for the heavy-ion collision experiments at lower collision energies. These await to be investigated as future extensions.