Physics Overview ---------------- The present code provides a hadron resonance gas (HRG) model of the QCD equation of state (EoS) as a function of temperature *T* and the three chemical potentials: baryon :math:`\mu_B`, strangeness :math:`\mu_S`, and electric charge :math:`\mu_Q`. Because of the hadronic nature of the model, it is expected to be most reliable at zero to low *T* and low–moderate :math:`\mu_B`. Consistent with this expectation, the model reproduces lattice QCD thermodynamics up to :math:`T \lesssim 150\text{-}160\ \mathrm{MeV}` and is commonly used to describe the chemical freeze-out stage of heavy-ion collisions. This code was written following the method described in Ref. [1]_. Ideal HRG ~~~~~~~~~ In the ideal HRG (IHRG), the pressure is given by a sum over hadronic species, .. math:: P_{id}(T,\{\mu\}) = \sum_i p^{id}_i(T,\mu_i), with the chemical potential of species :math:`i` defined as .. _mu-i: .. math:: \mu_i = B_i\mu_B + Q_i\mu_Q + S_i\mu_S The contribution from a single particle species reads .. math:: p^{id}_i(T,\mu_i) = \frac{g_i}{6\pi^2} \int_0^\infty \frac{k^4\,dk}{\sqrt{k^2 + m_i^2}} \left[ \exp\left( \frac{\sqrt{k^2 + m_i^2} - \mu_i}{T} \right) + \eta_i \right]^{-1}. where :math:`g_i` is the spin–isospin degeneracy, :math:`m_i` is the particle mass, and :math:`\eta_i = +1` for fermions, :math:`\eta_i = -1` for bosons, and :math:`\eta_i = 0` in the Boltzmann approximation. Quantum van der Waals HRG ~~~~~~~~~~~~~~~~~~~~~~~~~ Beyond the ideal gas, repulsive and attractive mean-field interactions are included within the multi-component Quantum van der Waals HRG (QvdW-HRG) approach [2]_, [3]_. Short-range repulsion is encoded by an excluded-volume matrix :math:`\tilde b_{ij}`, while attractive interactions are described by coefficients :math:`a_{ij}`, which characterize interactions between species :math:`i` and :math:`j`. The grand-canonical pressure is given by: .. _eq:vdw-p: .. math:: P(T,\{\mu\}) = \sum_i p^{*}_i - \sum_{i,j} a_{ij}\, n_i n_j, where :math:`p^{*}_i \equiv p^{id}_i(T,\mu_i^{*})` denotes the ideal quantum-gas(Fermi/Bose) pressure of species :math:`i`, evaluated at the shifted chemical potentials :math:`\mu_i^{*}`. These effective chemical potentials satisfy the transcendental system: .. _eq:vdw-mu: .. math:: \mu_i^{*} = \mu_i - \sum_j \tilde b_{ij}\, p^{*}_j + \sum_j (a_{ij} + a_{ji})\, n_j. The system is solved numerically using a multidimensional Broyden method. In regions where multiple solutions exist (e.g. near a first-order phase transition), several initial conditions are employed and the physical solution is selected using the Gibbs criterion of maximal pressure. The number densities are obtained from the linear system .. : eq-vdw-n .. math:: \sum_j \left( \delta_{ij} + \tilde b_{ji}\, n^{*}_i \right) n_j = n^{*}_i, \qquad n^{*}_i \equiv n^{id}_i(T,\mu_i^{*}). A commonly used prescription relates :math:`\tilde b_{ij}` to hard-sphere radii :math:`r_i` through the geometric pair excluded volume, .. math:: \tilde b_{ij} = \frac{2\, b_{ii}\, b_{ij}}{b_{ii} + b_{jj}}, \qquad b_{ij} = \frac{2\pi}{3}(r_i + r_j)^3, which reproduces the correct low-density virial limit for multi-component hard spheres [4]_. The *diagonal excluded-volume* special case assigns each species an eigenvolume :math:`v_i = \tfrac{16\pi}{3} r_i^3` and sets :math:`\tilde b_{ij} \equiv v_i` [5]_. Attractive interactions are specified by :math:`a_{ij}`. In practice, only (anti)baryon–(anti)baryon interactions are retained (:math:`a_{BB'}, a_{\bar B \bar B'} \neq 0`), while baryon–antibaryon, meson–baryon, and meson–meson attractions are neglected to avoid double counting of resonant channels already included explicitly. The QvdW-HRG model reduces to the IHRG limit when both :math:`a_{ij}` and :math:`\tilde b_{ij}` tend to zero. The QvdW-HRG reproduces key features of nuclear matter, including the liquid–gas phase transition and the associated critical point. A minimal and widely used implementation assumes that both repulsion and attraction act only between baryons (and, by symmetry, between antibaryons), with universal interaction parameters :math:`a_{BB'} = a_{\bar B \bar B'} = a` and :math:`\tilde b_{BB'} = \tilde b_{\bar B \bar B'} = b` for all (anti)baryon pairs. This two-parameter realization is sufficient to reproduce empirical properties of symmetric nuclear matter at :math:`T = 0`, such as the binding energy per nucleon and the saturation density, with typical values :math:`a = 0.329\,\mathrm{GeV\,fm}^3` and :math:`b = 3.42\,\mathrm{fm}^3` [6]_, [7]_. The model predicts a nuclear matter critical point at :math:`T_c \approx 19.7\,\mathrm{MeV}` and :math:`n_c \approx 0.07\,\mathrm{fm}^{-3}`, in reasonable agreement with constraints from nuclear multifragmentation data [8]_. At vanishing baryon chemical potential, the resulting equation of state also provides a fair description of lattice QCD thermodynamics, including fluctuations of conserved charges [2]_. Thermodynamic Derivatives and Susceptibilities ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ In addition to standard thermodynamic observables such as the pressure, energy density, and entropy density, the module also outputs the Hessian matrix of the pressure in the grand-canonical ensemble, .. math:: H = \begin{pmatrix} \dfrac{\partial^2 P}{\partial T^2} & \dfrac{\partial^2 P}{\partial T\,\partial \mu_Y} \\[6pt] \dfrac{\partial^2 P}{\partial \mu_X\,\partial T} & \dfrac{\partial^2 P}{\partial \mu_X\,\partial \mu_Y} \end{pmatrix} \equiv \begin{pmatrix} \partial s/\partial T & \alpha_Y \\[6pt] \alpha_X & \chi^{XY} \end{pmatrix}. Here :math:`\alpha_X = (\partial n_X / \partial T)_{\mu}` are thermal expansion coefficients, evaluated at fixed chemical potentials, with :math:`X, Y = B, S, Q` denoting the conserved baryon number, strangeness, and electric charge. The lower-right block of the Hessian corresponds to second-order net-charge susceptibilities of the conserved charges, namely :math:`\chi_2^B`, :math:`\chi_2^S`, and :math:`\chi_2^Q`, as well as the off-diagonal correlators :math:`\chi_{11}^{BQ}`, :math:`\chi_{11}^{BS}`, and :math:`\chi_{11}^{QS}`. The net-charge susceptibilities are defined as: .. math:: \chi^{BSQ}_{ijk}(T) = \frac{\partial^{i+j+k} (P/T^4)} {\partial \hat{\mu}_B^{\,i} \partial \hat{\mu}_S^{\,j} \partial \hat{\mu}_Q^{\,k}}. where the dimensionless chemical potentials are given by :math:`\hat{\mu}_X = \mu_X / T`. For the IHRG, the derivatives entering the Hessian matrix are obtained from explicit sums over particle species, since the pressure factorizes as :math:`P(T,\{\mu\}) = \sum_i p_i^{id}(T,\mu_i),` with the single-particle chemical potentials :math:`\mu_i` defined previously. In the interacting QvdW-HRG model, evaluation of the Hessian requires differentiation of the implicit interacting system for the pressure and effective chemical potentials with respect to :math:`T` and the conserved charge chemical potentials :math:`\mu_{B,S,Q}`, as described in the corresponding section. The resulting linear response equations are solved numerically to obtain the thermodynamic derivatives. Additional features ~~~~~~~~~~~~~~~~~~~ The MUSES HRG module includes the possibility to set a finite (anti-)meson excluded-volume parameter :math:`b_M`, in order to suppress hadronic degrees of freedom at high temperatures, where the QGP is expected. This option is particularly useful for a smooth merging of the hadronic EoS with other EoSs in the crossover region. Assigning small but non-zero values, typically :math:`b_{M} \sim 0.1-0.3~\mathrm{fm}^3` (corresponding to :math:`r_M\!\sim\!0.15-0.25~\mathrm{fm}`) reduces the total pressure and entropy density at high :math:`T` without spoiling the existing agreement of the HRG model with lattice data at lower :math:`T`. The HRG module also incorporates an effective mass model for pions [9]_, tuned to reproduce chiral perturbation theory for pions at :math:`T = 0`, which allows for calculations in the pion condensation region up to values of :math:`\mu_Q` currently restricted to :math:`2 \times m_\pi`. In this model, the pion mass is modified by a mean-field interaction term that depends on the pion density, effectively shifting the onset of Bose–Einstein condensation to higher values of :math:`\mu_Q` and allowing for a more realistic description of the system in this regime. References ~~~~~~~~~~ .. [1] V. Vovchenko and H. Stoecker, *Thermal-FIST: A package for heavy-ion collisions and hadronic equation of state*, Comput. Phys. Commun. **244**, 295–310 (2019), doi:10.1016/j.cpc.2019.06.024, arXiv:1901.05249 [nucl-th]. .. [2] V. Vovchenko, M. I. Gorenstein and H. Stoecker, *van der Waals Interactions in Hadron Resonance Gas: From Nuclear Matter to Lattice QCD*, Phys. Rev. Lett. **118**, no.18, 182301 (2017), doi:10.1103/PhysRevLett.118.182301, arXiv:1609.03975 [hep-ph]. .. [3] V. Vovchenko, A. Motornenko, P. Alba, M. I. Gorenstein, L. M. Satarov and H. Stoecker, *Multicomponent van der Waals equation of state: Applications in nuclear and hadronic physics*, Phys. Rev. C **96**, no.4, 045202 (2017), doi:10.1103/PhysRevC.96.045202, arXiv:1707.09215 [nucl-th]. .. [4] M. I. Gorenstein, A. P. Kostyuk and Y. D. Krivenko, *Van der Waals excluded volume model of multicomponent hadron gas*, J. Phys. G **25**, L75–L83 (1999), doi:10.1088/0954-3899/25/9/102, arXiv:nucl-th/9906068 [nucl-th]. .. [5] V. Vovchenko and H. Stoecker, *Examination of the sensitivity of the thermal fits to heavy-ion hadron yield data to the modeling of the eigenvolume interactions*, Phys. Rev. C **95**, no.4, 044904 (2017), doi:10.1103/PhysRevC.95.044904, arXiv:1606.06218 [hep-ph]. .. [6] V. Vovchenko, D. V. Anchishkin and M. I. Gorenstein, *Van der Waals Equation of State with Fermi Statistics for Nuclear Matter*, Phys. Rev. C **91**, no.6, 064314 (2015), doi:10.1103/PhysRevC.91.064314, arXiv:1504.01363 [nucl-th]. .. [7] V. Vovchenko, D. V. Anchishkin, M. I. Gorenstein and R. V. Poberezhnyuk, *Scaled variance, skewness, and kurtosis near the critical point of nuclear matter*, Phys. Rev. C **92**, no.5, 054901 (2015), doi:10.1103/PhysRevC.92.054901, arXiv:1506.05763 [nucl-th]. .. [8] J. B. Elliott, P. T. Lake, L. G. Moretto and L. Phair, *Determination of the coexistence curve, critical temperature, density, and pressure of bulk nuclear matter from fragment emission data*, Phys. Rev. C **87**, no.5, 054622 (2013), doi:10.1103/PhysRevC.87.054622. .. [9] V. Vovchenko, B. B. Brandt, F. Cuteri, G. Endrődi, F. Hajkarim and J. Schaffner-Bielich, *Pion Condensation in the Early Universe at Nonvanishing Lepton Flavor Asymmetry and Its Gravitational Wave Signatures*, Phys. Rev. Lett. **126**, no.1, 012701 (2021), doi:10.1103/PhysRevLett.126.012701, arXiv:2009.02309 [hep-ph].