Thermal-FIST: A package for heavy-ion collisions and hadronic equation of state

TL;DR

Thermal-FIST is a versatile C++ toolkit for analyzing hadron resonance gas models, enabling detailed studies of particle production, fluctuations, and interactions in heavy-ion collisions with an interactive GUI.

Contribution

It introduces a comprehensive, user-friendly software package that combines advanced statistical analysis, event generation, and a graphical interface for heavy-ion collision research.

Findings

Supports various statistical ensembles including grand canonical and canonical.

Includes features for fluctuations, correlations, and hadronic interactions.

Provides a fast event generator and a graphical user interface.

Abstract

Thermal-FIST (Thermal, Fast and Interactive Statistical Toolkit) is a C++ package designed for a convenient general-purpose physics analysis within the family of hadron resonance gas (HRG) models. This mainly includes the statistical analysis of particle production in heavy-ion collisions and the phenomenology of hadronic equation of state. Notable features include fluctuations and correlations of conserved charges, effects of probabilistic decay, chemical non-equilibrium, and inclusion of van der Waals hadronic interactions. Calculations are possible within the grand canonical ensemble, the canonical ensemble, as well as in mixed-canonical ensembles combining the canonical treatment of certain conserved charges with the grand-canonical treatment of other conserved charges. The package contains a fast thermal event generator, which generates particle yields in accordance with the HRG chemistry, and particle momenta based on the Blast Wave model. A distinct feature of this package is the presence of the graphical user interface frontend -- QtThermalFIST -- which is designed for fast and convenient general-purpose HRG model applications.

Figures

Tables

Table 1. Table 1: List of HRG model parameters contained in the ThermalModelParameters structure.

Parameter	Unit	Notes
Temperature ($T$)	GeV	–
Baryon chemical potential (${\mu }_B$)	GeV	Not in CE
Electric chemical potential (${\mu }_Q$)	GeV	Not in CE
Strangeness chemical potential (${\mu }_S$)	GeV	Not in CE/SCE
Charm chemical potential (${\mu }_C$)	GeV	Not in CE/SCE/CCE
${\gamma }_q$	–	–
${\gamma }_S$	–	–
${\gamma }_C$	–	–
Volume ($V$)	fm3	–
Canonical volume ($V_c$)	fm3	In CE, SCE, CCE
Baryon charge ($B$)	–	CE only
Electric charge ($Q$)	–	CE only
Strangeness ($S$)	–	CE, SCE
Charm ($C$)	–	CE, CCE

Table 2. Table 2: List of hadron properties contained in an input file.

Column	Value type	Description
PDG ID	long long	PDG code
Name	string	Particle name
Stability flag	bool	Stability w.r.t decays
Mass [GeV]	double	Particle mass
Degeneracy	double	Internal degeneracy factor
Statistics	int	Fermi (+1), Bose (-1), Boltzmann (0)
$B$	int	Baryon number
$Q$	int	Electric charge
$S$	int	Strangeness
$C$	int	Charm
$n_{\left|S\right|}$	double	Absolute strange quark content
$n_{\left|C\right|}$	double	Absolute charm quark content
Width [GeV]	double	Resonance width
Threshold [GeV]	double	Threshold mass

Table 3. Table 3: List of the ThermalModelBase methods to set HRG parameters.

Method	Parameter	Notes
SetTemperature	$T$	–
SetBaryonChemicalPotential	${\mu }_B$	–
SetElectricChemicalPotential	${\mu }_Q$	–
SetStrangenessChemicalPotential	${\mu }_S$	–
SetCharmChemicalPotential	${\mu }_Q$	–
SetGammaq	${\gamma }_q$	–
SetGammaS	${\gamma }_S$	–
SetGammaC	${\gamma }_C$	–
SetVolume	$V$	–
SetVolumeRadius	$R$	Sets $V=\left(4/3\right)\pi R^3$
SetCanonicalVolume	$V_c$	CE, SCE, CCE
SetCanonicalVolumeRadius	$R_c$	Sets $V_c=\left(4/3\right)\pi R_c^3$
SetBaryonCharge	$B$	For CE only
SetElectricCharge	$Q$	For CE only
SetStrangeness	$S$	CE, SCE
SetCharm	$C$	CE, CCE

Table 4. Table 4: List of classes inheriting from the ThermalModelBase class and the corresponding HRG model.

Class	HRG model variant	Ensemble
ThermalModelIdeal	Id-HRG	GCE
ThermalModelDiagonalEV	Diagonal EV-HRG	GCE
ThermalModelCrosstermsEV	Non-diagonal EV-HRG	GCE
ThermalModelVDW	QvdW-HRG	GCE
ThermalModelCanonical	Id-HRG	CE/GCE
ThermalModelCanonicalStrangeness	Id-HRG	SCE
ThermalModelCanonicalCharm	Id-HRG	CCE

Table 5. Table 5: List of fit parameters in a thermal model fit, their default values and fit flags

Parameter	Default
	value	fit flag
$T$	0.150 GeV	true
${\mu }_B$	0	true
${\mu }_Q$	0	false
${\mu }_S$	0	false
${\mu }_C$	0	false
${\gamma }_q$	1	false
${\gamma }_S$	1	false
${\gamma }_C$	1	false
$R$	8 fm	true
$R_c$	8 fm	true

Equations

$$p(T,\mu )=\sum _i{p}_i^{\mathrm{id}}(T,{\mu }_i),$$

$$p_i^{\mathrm{id}}(T,{\mu }_i)=\frac{d_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},$$

$$n_i^{\mathrm{id}}(T,{\mu }_i)=\frac{\partial p_i^{\mathrm{id}}(T,{\mu }_i)}{\partial {\mu }_i}=\frac{d_i}{2{\pi }^2}{\int }_0^{\infty }{k}^{2}dk{\left[\exp \left(\frac{\sqrt{k^2+m_i^2}-{\mu }_i}{T}\right)+{\eta }_i\right]}^{-1},$$

$${\epsilon }_i^{\mathrm{id}}(T,{\mu }_i)=\frac{d_i}{2{\pi }^2}{\int }_0^{\infty }{k}^{2}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}.$$

$${\mu }_i=B_i{\mu }_B+S_i{\mu }_S+Q_i{\mu }_Q+C_i{\mu }_C$$

$$e^{{\mu }_i/T}\to e^{{\mu }_i/T}{\gamma }_q^{|q_i|}{\gamma }_S^{|S_i|}{\gamma }_C^{|C_i|}\text{[or }{\mu }_i\to {\mu }_i+T\log ({\gamma }_q^{|q_i|}{\gamma }_S^{|S_i|}{\gamma }_C^{|C_i|})\text{]},$$

$$p(T,\mu )=\sum _i{p}^{\mathrm{id}}(T,{\mu }_i^{\ast }),{\mu }_i^{\ast }={\mu }_i-v_{i}p,$$

$$n_i(T,\mu )$$

$$s(T,\mu )$$

$$\epsilon (T,\mu )$$

$$p(T,\mu )=\sum _i{p}_i(T,\mu ),$$

$$p_i(T,\mu )=p_i^{\mathrm{id}}(T,{\mu }_i^{\ast }),{\mu }_i^{\ast }={\mu }_i-\sum _j{\tilde{b}}_{ij}{p}_j,i=1,\dots ,f.$$

$$\sum _j[{\delta }_{ij}+{\tilde{b}}_{ji}{n}_i^{\mathrm{id}}(T,{\mu }_i^{\ast })]n_j=n_i^{\mathrm{id}}(T,{\mu }_i^{\ast }),i=1\dots f.$$

$$s(T,\mu )$$

$$\epsilon (T,\mu )$$

$${\tilde{b}}_{ij}=\frac{2b_{ii}{b}_{ij}}{b_{ii}+b_{jj}}\text{with}{b}_{ij}=\frac{2\pi }{3}(r_i+r_j{)}^3.$$

$$p(T,\mu )=\sum _i{p}_i^{\mathrm{id}}(T,{\mu }_i^{\ast })-\sum _{i,j}{a}_{ij}{n}_i{n}_j.$$

$$\sum _j[{\delta }_{ij}+{\tilde{b}}_{ji}{n}_i^{\mathrm{id}}(T,{\mu }_i^{\ast })]n_j=n_i^{\mathrm{id}}(T,{\mu }_i^{\ast }),i=1\dots f,$$

$${\mu }_i^{\ast }+\sum _j{\tilde{b}}_{ij}{p}_j^{\ast }-\sum _j(a_{ij}+a_{ji})n_j={\mu }_i,i=1,\dots ,f.$$

$$s(T,\mu )$$

$$\epsilon (T,\mu )$$

$$\int dk\to {\int }_{m_i^{\min }}^{m_i^{\max }}dm{\rho }_i(m)\int dk,$$

$${\rho }_i(m)=A_i\frac{2m{m}_i{\Gamma }_i(m)}{(m^2-m_i^2{)}^2+m_i^2[{\Gamma }_i(m){]}^2},$$

$${\rho }_i(m)=A_i\frac{1}{(m-m_i{)}^2+[{\Gamma }_i(m){]}^2/4}.$$

$$⟨N_i⟩=⟨N_i^{\ast }⟩+\sum _R⟨n_i{⟩}_R⟨N_R^{\ast }⟩,$$

$${\chi }_{lmnk}^{BSQC}=\frac{{\partial }^{l+m+n+k}p/T^4}{\partial ({\mu }_B/T{)}^l\partial ({\mu }_S/T{)}^m\partial ({\mu }_Q/T{)}^n\partial ({\mu }_C/T{)}^k}.$$

$${\chi }_2^i$$

$${\chi }_3^i$$

$${\chi }_4^i$$

$$\frac{1}{V{T}^3}⟨\Delta N_i^{\ast }\Delta N_j^{\ast }⟩=\frac{{\partial }^2(p/T^4)}{\partial ({\mu }_i/T)\partial ({\mu }_j/T)},$$

$$⟨(\Delta N_i{)}^2⟩$$

$$+2\sum _R⟨\Delta N_i^{\ast }\Delta N_R^{\ast }⟩⟨n_i{⟩}_R+\sum _{R,R^{\prime }}⟨\Delta N_R^{\ast }\Delta N_{R^{\prime }}^{\ast }⟩⟨n_i{⟩}_R⟨n_i{⟩}_{R^{\prime }}.$$

$$⟨(\Delta N_i{)}^3⟩$$

$$+\sum _R⟨N_R^{\ast }⟩⟨(\Delta n_i{)}^3{⟩}_R,$$

Full text

††journal: Computer Physics Communications

Program summary

Program Title: Thermal-FIST, version 1.2

Licensing provisions: GPLv30

Programming language: C++

Computer: any with a C++ compiler and, optionally, the Qt5 framework for the GUI frontend

Operating system: cross-platform, tested on Linux Ubuntu 16.04, 18.04; Mac OS X Yosemite, Microsoft Windows 10; Android 8, 9

External routines: Eigen template library for the linear algebra routines [1], MINUIT2 package from CERN ROOT [2], Mersenne Twister random number generator [3], Qt5 framework [4] (for the GUI only), QCustomPlot Qt widget [5] (for the GUI only)

Nature of problem:

The HRG model and its various modifications constitute a common framework used for modeling of the hadronic equation of state and particle production in heavy-ion collisions. Even the simplest versions of the HRG model require careful considerations of the many details, including the resonance decay feed-down, implementation of charge conservation constraints relevant for heavy-ion collisions, chemical non-equilibrium effects. A notable extra effort is required in order to treat the fluctuations and correlations of various charges charges, which presently are being extensively studied in the heavy-ion collision experiments and lattice QCD calculations. The inclusion of hadronic interactions, modeled by an excluded-volume (EV) or a van der Waals (vdW) type framework, additionally requires a numerical solution to a system of many transcendental equations.

Solution method:

The Thermal-FIST package contains a class-based library which calculates relevant HRG observables for a specified setup. The setup includes a particle list, usually to be supplied with an external file, an HRG model specification (statistical ensemble, van der Waals interaction parameters, etc.), a set of thermal parameters, and conservation laws constraints. Whenever necessary, the systems of transcendental equations are solved numerically with the Broyden’s method. The package includes a fitter for extracting thermal parameters from hadron yield data through the $\chi^{2}$ minimization. The HRG model based Monte Carlo event generator is a complementary feature to analytic calculations. A general-purpose thermal analysis is made maximally convenient with QtThermalFIST – a GUI frontend based on the Qt framework where all typical calculations, such as the properties of the equation of state or the thermal fits, can be straightforwardly performed.

*Additional comments:

*If the EV/vdW interactions are present, exact analytic calculations are presently only possible within the grand canonical ensemble. Approximate calculations are possible for the strangeness-canonical ensemble on the condition that strange particles form a small subsystem relative to the total system. Effects of probabilistic decays on fluctuation observables are generally included only up to the moments of the 2nd order. The only exception is the ideal HRG model in the grand canonical ensemble, where these effects are included up to the moments of the 4th order. On the other hand, the Monte Carlo event generator is not constrained by the above restrictions.

*Running time:

*Depending on the specific task, calculation times may vary from milliseconds (an evaluation of thermodynamic functions within the ideal HRG model in the GCE) to several minutes (a thermal fit with sophisticated features such as energy-dependent widths, van der Waals interactions, or exact charge conservation). The thermal event generator takes few milliseconds to generate an event with about a thousand particles in the GCE. Event generation in the CE is considerably slower and depends on the rejection sampling rate realized for the given setup.

1 Introduction

The abundant hadron production in heavy-ion collision reactions has long been treated in the framework of the thermal-statistical models [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. Such a description assumes emission of particles from a thermally and chemically equilibrated source created in these reactions. Fitting the observed yields of stable hadrons allows to determine the thermal parameters, corresponding to the so-called chemical freeze-out stage of the collision. In most cases, the ideal hadron resonance gas (Id-HRG) model has been used, and a surprisingly good description of many experimental hadron yield data from heavy-ion collisions have been achieved within this simple approach for a broad range of collision energies, ranging from moderate energies at the SchwerIonen-Synchrotron (SIS) to the highest energy of the Large Hadron Collider (LHC) (see, e.g., Refs. [11, 12, 13] for an overview).

HRG-type models also play an important role in the phenomenology of the QCD equation of state. The Id-HRG model is a popular choice for describing the low temperature, confined phase of QCD. At temperatures between $T\sim 100-150$ MeV and at zero chemical potential, the Id-HRG model indeed appears to reproduce many lattice QCD observables [14, 15, 16, 17].

Several implementations of the Id-HRG model exist on the market, including SHARE [18, 19, 20], THERMUS [21], and THERMINATOR [22, 23].

In many applications, significant deviations from the ideal gas picture can be expected. Extensions of the ideal gas picture have been discussed mostly within the excluded volume (EV) HRG model [24, 25, 26], where the effects of repulsive hadronic interactions at short distances are introduced (see, e.g., [27] for recent developments). Another extension is the quantum van der Waals (QvdW) model [28, 29, 30], which allows to include both the repulsive and attractive interactions between hadrons. Recently, repulsive interactions have received renewed interest in the context of lattice QCD data on fluctuations of conserved charges. It was indicated that large deviations of several fluctuation observables from the Id-HRG baseline are captured by HRG models with repulsive baryon-baryon interactions [29, 31, 32, 33]. The Thermal-FIST package presented here allows to calculate thermodynamic features and fluctuation observables within a HRG with arbitrary attractive and repulsive QvdW parameters characterizing interactions between each type of hadron species.

Possible effects of chemical non-equilibrium for light, strange, and/or charm quarks in heavy-ion collisions have been studied in the framework of the Id-HRG model within the SHARE package [11, 34]. The Thermal-FIST package allows to study the chemical non-equilibrium effects simultaneously with the effects of EV/vdW interactions or exact charge conservation within the canonical ensemble.

Recent thermal model applications also include the description of the measurements of multiplicity fluctuations [35, 36, 37, 38]. Treatment of probabilistic decays [39] and hadronic interactions [36] is important in such applications, and it is included into the Thermal-FIST package.

The applications of the Thermal-FIST package which can be found in the literature include chemical freeze-out analysis in proton-proton [40, 41] and nucleus-nucleus collisions [42, 43, 44], influence of the EV interactions on thermal fits [45, 46, 47, 27] and equation of state [48, 49, 32, 33], Monte Carlo analysis of EV effects in the canonical ensemble [50], effects of the vdW interactions between baryons and nuclear liquid-gas transition on various observables [29, 51, 30, 52, 53], and thermal production of light nuclei [54, 55].

2 Hadron resonance gas

2.1 Ideal HRG

In the simplest setup, the thermodynamics of hadronic phase is described by a multi-component, ideal gas of point-like hadrons – the Id-HRG model. In the grand canonical ensemble (GCE) formulation of the Id-HRG model there are no correlations between different hadronic species. Thus, the pressure is given by

[TABLE]

where the sum goes over all hadron species included in the model, $p^{\rm id}_{i}(T,\mu_{i})$ is the pressure of the ideal Fermi or Bose gas at the corresponding temperature and chemical potential for species $i$:

[TABLE]

where $d_{i}$ and $m_{i}$ are, respectively, the spin degeneracy factor and mass of hadron species $i$, and where $\eta_{i}$ equals +1 for fermions, -1 for bosons, and 0 for the Boltzmann approximation.

Other thermodynamic quantities are given by expressions similar to (1): a sum over the corresponding ideal gas quantities for all hadron species. The particle density of hadron species $i$ is $n_{i}^{\rm id}(T,\mu_{i})$, i.e. it is simply given by the ideal gas relation for species $i$:

[TABLE]

while the energy density is $\varepsilon(T,\mu)=\sum_{i}\varepsilon_{i}^{\rm id}(T,\mu_{i})$ with

[TABLE]

Within the GCE formulation, all conserved charges, such as baryonic number $B$, electric charge $Q$, strangeness $S$, and charm $C$, are conserved on average. For the chemical equilibrium case these average values are regulated by the corresponding independent chemical potentials: $\mu_{B}$, $\mu_{Q}$, $\mu_{S}$, and $\mu_{C}$111Other charges, such a bottomness, can be considered additionally, in the same manner.. The chemical potential of the $i$th hadron species is thus determined as

[TABLE]

with $B_{i}=0,\,\pm 1$, $S_{i}=0,\,\pm 1,\,\pm 2,\,\pm 3$, $Q_{i}=0,\,\pm 1,\,\pm 2$, and $C_{i}=0,\,\pm 1,\,\pm 2$, being the corresponding conserved charges of the $i$th hadron: baryon number, strangeness, electric charge, and charm.

2.2 Chemical non-equilibrium and fugacity factors

The assumption of full chemical equilibrium can be relaxed in the heavy-ion collision applications of the HRG model. This is usually done by introducing additional fugacity factors which regulate the absolute abundances of quarks of different flavor [56, 11], initially introduced as phenomenological parameters [57, 58]. In this case, the Boltzmann factors in Eqs. (2)-(4) are modified to

[TABLE]

where $|q_{i}|$, $|S_{i}|$, and $|C_{i}|$ correspond, respectively, to the absolute light, strange, and charm quark content of hadron $i$, and where $\gamma_{q}$, $\gamma_{S}$, and $\gamma_{C}$ are parameters which regulate deviations from chemical equilibrium in the light, strange, and charm quark sectors, respectively. $\gamma_{i}=1$ corresponds to the chemical equilibrium scenario for the corresponding quark flavor sector.

2.3 Excluded-volume corrections

The repulsive interactions between hadrons can be modeled by an excluded volume correction of the van der Waals type, whereby the volume available for hadrons to be created and move in is reduced by the sum of all their eigenvolumes. Such a correction was first studied in the hadronic equation of state in Refs. [59, 60, 61, 62, 63]. A thermodynamically consistent procedure for a single-component gas was first formulated in Ref. [24].

It should be noted that an excluded-volume correction is only an effective approach to treat repulsive interactions between hadrons. A common (and simplest) assumption is the constant eigenvolume parameter for all hadronic species [64]. A more realistic approach is to allow the possibility for different hadrons to have different eigenvolumes [25, 26]. Such an approach, however, is incomplete. For instance, it is not possible to take into account the expected differences between baryon-baryon and baryon-antibaryon interactions [65] within such a model. Therefore, each pair of hadron species can be characterized by its own “excluded volume” parameter [66] in the most general case.

All of the above options are considered and implemented in the Thermal-FIST package.

2.3.1 Diagonal EV model

In the Diagonal EV-HRG (DEV-HRG) model [25, 26], each hadron is assigned an excluded-volume parameter $v_{i}$. It is common to characterize the excluded-volume parameter $v_{i}$ with the effective hard-core radius $r_{i}$, by using the classical relation $v_{i}=(16\pi/3)\,r_{i}^{3}$. The excluded volume correction leads to the transcendental equation for the system pressure,

[TABLE]

which is solved numerically. The numerical solution is obtained in Thermal-FIST using the Broyden’s method [67].

Other thermodynamic functions are then obtained from the standard thermodynamic relations:

[TABLE]

2.3.2 Non-diagonal EV model

The repulsive interactions in the Non-Diagonal EV-HRG (NDEV-HRG) model [66, 46, 27] are characterized by the matrix $\tilde{b}_{ij}$ of the excluded volume type parameters, which characterize the repulsive interactions for each pair of particle species. The total pressure is partitioned into the sum of “partial” pressures,

[TABLE]

which are determined by the following system of transcendental equations:

[TABLE]

The solution to Eq. (12) at given $T$ and $\mu$ is obtained numerically in Thermal-FIST, using the multi-dimensional Broyden’s method.

The particle number densities $n_{i}\equiv(\partial p/\partial\mu_{i})_{T}$ are found as the solution to the system of linear equations

[TABLE]

The entropy and energy densities are given by

[TABLE]

Parameters $\tilde{b}_{ij}$ can be chosen arbitrarily for each pair of particle species. They need not to be symmetric. One possibility is the classical picture of a multi-component gas of hard spheres222It should be noted that the classical picture is not necessarily valid at the nuclear scale, see [33]., where each species is characterized by the hard-core radius $r_{i}$ and where $\tilde{b}_{ij}$ are given by

[TABLE]

The NDEV-HRG model reduces to the DEV-HRG model in the partial case $\widetilde{b}_{ij}\equiv v_{i}$.

2.4 Quantum van der Waals model

The presence of both, the short-range repulsive and the intermediate/long range attractive interactions between hadrons can be treated in the framework of the Quantum van der Waals (QvdW) equation [28], extended to multiple components present in a HRG [29, 30]. This QvdW-HRG model is defined by the following pressure function

[TABLE]

The particle number densities, $n_{i}$, satisfy the system of linear equations

[TABLE]

while the shifted chemical potentials, $\mu_{i}^{*}$, satisfy the system of transcendental equations

[TABLE]

The pressure, $p(T,\mu)$, at a given temperature $T$ and chemical potentials $\mu$ is determined by first solving numerically the system of equations (19) for $\mu_{i}^{*}$333Note that particle number densities $n_{j}$ in Eq. (19) are given as a function of $\{\mu_{i}^{*}\}$ via Eq. (18). and then plugging in the result into Eq. (17).

The entropy and energy densities are given as

[TABLE]

The parameters $\tilde{b}_{ij}$ correspond to the repulsive vdW interactions, and they have the same physical meaning as in the NDEV-HRG model. The parameters $a_{ij}$ correspond to the attractive vdW interactions444When introducing terms corresponding to the attractive vdW interactions one should take care not to double count those interactions that lead to the formation of resonances. between hadron species $i$ and $j$, modeled in the mean-field approximation. The QvdW-HRG model reduces to the NDEV-HRG model for the case $a_{ij}\equiv 0$.

The QvdW-HRG model permits the inclusion into the HRG model of the basic features of nuclear matter, in particular, the nuclear liquid-gas phase transition and the associated criticality. The critical point of nuclear matter was shown to be important for the fluctuation observables in heavy-ion collisions [68, 29, 69, 53].

2.5 Finite resonance widths

The finite widths of the resonances may play an important role for some applications of the HRG model. These can be taken into account in a simplified way, by adding an additional integration into Eqs. (2)-(4) over the resonance masses [70, 21]

[TABLE]

where $\rho_{i}(m)$ is the mass distribution for resonance $i$, taken either in the relativistic Breit-Wigner form,

[TABLE]

or the nonrelativistic Breit-Wigner form,

[TABLE]

Here $A_{i}$ is determined from the normalization condition $\int_{m_{i}^{\rm min}}^{m_{i}^{\rm max}}dm\,\rho_{i}(m)=1$.

There are different possibilities for the choice of $m_{i}^{\rm min}$ and $m_{i}^{\rm max}$, as well as for the mass-dependence of the width $\Gamma_{i}(m)$. Two options are implemented in Thermal-FIST: (1) integration in the truncated interval, defined by $m_{i}^{\rm min}=\max(m_{i}-2\Gamma_{i},m_{i}^{\rm thr})$ and $m_{i}^{\rm max}=m_{i}+2\Gamma_{i}$, with energy independent widths [70, 21], and (2) integration in the full interval, defined by $m_{i}^{\rm min}=m_{i}^{\rm thr}$ and $m_{i}^{\rm max}=\infty$, with energy dependent widths. Here $m_{i}^{\rm thr}$ is the threshold mass. More details about the different options and their implementation can be found in Ref. [43].

2.6 Feeddown from resonance decays

Applications of the HRG model to particle production must take into account the feeding from decays of unstable particles (resonances) to the final hadron yields measured in experiments. Feeding from resonance decays gives dominant contribution in many cases, as much as 70% of all final pions may come from resonance decays [21]. This is illustrated in Fig. 1, which shows different feeddown contribution to the total yields of pions and protons, evaluated in the Id-HRG model at $T=155$ MeV, $V=4000$ fm3, and zero chemical potentials.

The final mean multiplicity $\langle N_{i}\rangle$ of $i$th particle species is calculated in the HRG model as a sum of the primordial mean multiplicity $\langle N^{*}_{i}\rangle\equiv n_{i}\,V$ and resonance decay contributions as follows

[TABLE]

where $\langle n_{i}\rangle_{R}$ is the average number of particles of type $i$ resulting from decay of resonance $R$. $\langle n_{i}\rangle_{R}$ includes contribution from both, the direct decays of resonance $R$ resulting in the production of hadron $i$, as well as the contributions resulting from the chain of decays via lower-mass resonances. Note that for energy-dependent resonance widths the relation (25) is modified to reflect energy dependent branching ratios, see Ref. [43] for details.

The exact decay feeddown chain for each hadron yield should be matched to the relevant experimental conditions. By default, the hadrons which decay strongly or electromagnetically are marked unstable in Thermal-FIST. It is also possible, for each hadron species, to evaluate final yields which include contributions from strong decays only, from strong and electromagnetic decays, from strong, electromagnetic and weak decays, or in accordance with manually set particle stability flags that define the feeddown.

2.7 Correlations, fluctuations and susceptibilities

2.7.1 Conserved charges

Fluctuations and correlations of conserved charges are characterized by the corresponding susceptibilities

[TABLE]

These susceptibilities are useful, as they carry information about finer details of the equation of state, they can be obtained in the lattice QCD simulations, and they can also be measured in heavy-ion collisions, although in many cases proxy observables have to be used.

In the most general case, the derivatives defining a particular $\chi_{lmnk}^{BSQC}$ can be evaluated numerically by calculating the pressure function at different values of chemical potentials. The diagonal susceptibilities of an arbitrary conserved charge (which can be, e.g., some linear combination of $B$, $S$, $Q$, and $C$), can be calculated within Thermal-FIST analytically up to the 4th order. This analytic procedure is described in [30] for the QvdW-HRG model, all other models included are obtained as partial cases of the QvdW-HRG model.

2.7.2 Hadron yield fluctuations and probabilistic decays

Fluctuations and correlations of hadron numbers can also be considered. A hadron yield is not a conserved quantity, and, therefore, it cannot be calculated from first-principles using lattice QCD. On the other hand, in contrast to the conserved charges, hadron yields are quantities which are directly measured in heavy-ion collision experiments, and they are readily available for calculation in a HRG approach.

The cumulants $\chi_{k}^{i}$ of the primordial hadron number distributions, $N_{i}^{*}$, can be calculated by taking the derivatives with respect to the chemical potential $\mu_{i}$ of hadron specie $i$. In particular

[TABLE]

where $\chi_{k}^{i}$ are expressed here through the central moments of the primordial hadron number distribution. $\chi_{k}^{i}$ are evaluated in Thermal-FIST in the same way as the fluctuations of conserved charges.

The two-particle correlator $\langle\Delta N_{i}^{*}\Delta N_{j}^{*}\rangle$ for the primordial hadron yields $N_{i}^{*}$ and $N_{j}^{*}$ is given by the mixed derivative

[TABLE]

while the higher-order correlations are obtained in a similar fashion.

The correlations and fluctuations in final hadron yields are affected by the probabilistic decays of the resonances. The variance of the final hadron yield $N_{i}$ is calculated as follows [71, 72, 73, 36]

[TABLE]

Here $\langle(\Delta n_{i})^{2}\rangle_{R}$ is the variance of the number of hadrons of type $i$ which result from the chain of probabilistic decays of the resonance $R$.

The higher-order fluctuations of final hadron yields are considered in Thermal-FIST for the Id-HRG model only. The third and fourth cumulants in this model read [39]

[TABLE]

Here $\langle(\Delta n_{i})^{3}\rangle_{R}$ and $\langle(\Delta n_{i})^{4}\rangle_{R,c}=\langle(\Delta n_{i})^{4}\rangle_{R}-3\,\langle(\Delta n_{i})^{2}\rangle_{R}^{2}$ are, respectively, the third and fourth cumulants of the distribution of the number of hadrons of type $i$ which result from the chain of probabilistic decays of the resonance $R$.

2.8 Canonical ensemble

The exact conservation of conserved charges can be enforced within the canonical thermodynamics. The effects of exact charge conservation influence hadron yields in small systems, typically when the number of hadrons carrying a particular conserved charge is of the order of unity or smaller. The exact charge conservation also affects strongly the fluctuations of hadron numbers, even in the thermodynamic limit.

Analytic calculations in the canonical ensemble are restricted in Thermal-FIST to the Id-HRG model.

2.8.1 Full canonical ensemble

The canonical ensemble treatment implies exact conservation of the corresponding conserved charges. The partition function of the Id-HRG in the canonical ensemble for fixed baryon number $B$, electric charge $Q$, strangeness $S$, and charm $C$, reads [70, 74]

[TABLE]

Here

[TABLE]

with $V_{c}$ being the canonical correlation volume, i.e. the volume over which the exact conservation is enforced. The mean primordial multiplicity for species $j$ is given by

[TABLE]

Similar expressions exist for other thermodynamic functions, as well as for fluctuations and correlations of particle numbers [72]. Note that all the above relations include the effects of quantum statistics, given by the summation over $n$. This sum is truncated at a sufficiently large $n_{\rm max}$ value in calculations. The calculations proceed by evaluating the partition functions, given by Eq. (2.8.1), through a numerical integration. For the case when both, quantum statistical effects are neglected for baryons and when there are no particles in the list with $|B_{j}|>1$, the integration over $\phi_{B}$ is performed analytically, using the method of Ref. [75]. This allows to significantly speed up the calculations.

The inclusion of quantum statistics and fluctuations is a new element compared to the functionality of the presently available other open source packages. Most packages do not include the canonical ensemble at all. One exception is THERMUS [21], where calculations are, however, restricted to the case of Boltzmann statistics and where it is not possible to consider particles which have a baryon number $|B_{j}|>1$.

It is also possible to selectively treat certain conserved charges grand-canonically, while preserving the canonical treatment for other conserved charges. To remove the canonical treatment of a certain conserved charge, for instance the baryon charge $B$, one removes the integration over $\phi_{B}$ in Eq. (2.8.1), sets $\phi_{B}$ to zero in the integrand, and adds the baryon number fugacity factor into the expression (36) for the one-particle partition function. The selective canonical treatment of conserved charges can clarify the role of canonical effects for different conserved charges in a given setup.

The canonical ensemble implementation works reliably and fast when the system volume or conserved charges are sufficiently small. This covers the typical cases where the effects of exact charge conservation on mean hadron multiplicities are significant.

2.8.2 Strangeness/charm canonical ensemble

The abundances of hadrons carrying strangeness produced in intermediate-energy heavy-ion collisions are notably smaller than those of light flavored hadrons. The abundances of charm hadrons are even smaller. Therefore, in many cases it is sufficient to implement the canonical treatment of strangeness or charm only, while preserving the grand canonical treatment of baryon number and electric charge. The canonical treatment of strangeness is achieved in the so-called strangeness-canonical ensemble (SCE). In principle, the full canonical ensemble implementation as described above provides such a functionality. In some cases, however, the system volume can be large, and the above implementation can be inefficient. Therefore, here we consider an alternative approach to calculate the partition functions in Eq. (2.8.1). Assuming the Boltzmann statistics for strange particles, the SCE partition function is expressed as the sum over the Bessel functions [76]

[TABLE]

Here $\mathcal{Z_{\rm NS}}$ is the grand-canonical partition function of the subsystem consisting of non-strange particles, and

[TABLE]

with $S_{i}=\sum_{j\in(S_{j}=i)}\langle N_{j}^{\rm gce}\rangle$ being the cumulative mean multiplicity of all particle carrying strangeness $i$ calculated in the grand canonical ensemble.

The canonical treatment of charm here assumes that there are no multi-charmed particles in the particle list and that the system is net charm free, which is sufficient for most applications. The yields of charmed hadrons in the charm-canonical ensemble (CCE) are then calculated as follows

[TABLE]

2.9 Thermal fits

Perhaps the most common application of the HRG model is fitting the hadron yield data from relativistic heavy-ion collisions – the thermal fits. Such an approach assumes thermal and (partial) chemical equilibrium between all stable hadrons and all resonances at the so-called “chemical freeze-out” stage of a heavy-ion reaction. The HRG model fits are performed by minimizing the value

[TABLE]

where $N_{i}^{\rm exp}$ and $N_{i}^{\rm HRG}$ are the experimental and calculated in the HRG hadron multiplicities, respectively; $N_{\rm dof}$ is the number of degrees of freedom, that is the number of the data points minus the number of fitting parameters; and $\sigma_{i}^{2}=(\sigma_{i}^{syst})^{2}+(\sigma_{i}^{stat})^{2}$ is the sum of the squares of the statistical and systematic experimental errors. Note that $N_{i}^{\rm HRG}$ is the total hadron yield, including the resonance feeddown, calculated in accordance with Eq. (25) and using the appropriate feeddown flags. $N_{i}$ in Eq. (40) can also represent a ratio of two yields. Thermal-FIST employs the MINUIT2 package [77] for the $\chi^{2}$ minimization procedure.

In the simplest setup, corresponding to the full chemical equilibrium in the grand canonical ensemble, there are only three fit parameters: the temperature $T$, the baryonic chemical potential $\mu_{B}$, and the system volume parameter $V$. The electric charge and strangeness chemical potentials $\mu_{Q}$ and $\mu_{S}$ are not fitted. Instead, at each fixed temperature $T$ and baryochemical potential $\mu_{B}$, the $\mu_{Q}$ and $\mu_{S}$ are determined in a unique way in order to satisfy two conservation laws given by the “initial” conditions: the electric-to-baryon charge ratio of $Q/B=0.4$, and the vanishing net strangeness $S=0$. These two conditions are relevant if pre-freezeout radiation is neglected. Otherwise, $\mu_{S}$ and/or $\mu_{Q}$ can also be considered as additional fit parameters. For completeness, Thermal-FIST also allows to constrain the baryochemical potential $\mu_{B}$ to a fixed entropy-per-baryon ratio, $S/B$.

The chemical potentials disappear in the canonical ensemble formulation. In this case, the total baryon, electric charge, strangeness, and charm numbers appear instead, these are fixed by the “initial” conditions.

Some modifications, such as chemical under- or over-saturation of the light, strange, and/or charm quarks can also introduce additional parameters, $\gamma_{q}$, $\gamma_{S}$, and $\gamma_{C}$ (see, e.g., Refs. [11, 78]).

3 Monte Carlo event generator

Thermal-FIST implements Thermal Event Generator (TEG) – a Monte Carlo generator of hadronic microstates that correspond to a particular formulation of the HRG model in grand canonical or canonical ensemble. The effects of radial flow are included in the framework of the blast wave model. Monte Carlo implementation of the probabilistic decays of the primordial resonances is also provided.

The TEG is useful for heavy-ion collision applications, as well as for calculation of observables which are otherwise problematic using the analytic methods. These include e.g., effects of radial flow and momentum cuts. Another example is the higher-order fluctuations and correlations of final hadron yields, which may include the effects of the probabilistic decays of resonances and of the residual EV/QvdW interactions between the primordial hadrons. Monte Carlo formulation also allows to study the simultaneous effects of the EV interactions and the global charge conservation, which presently cannot be done using the analytic methods.

The TEG assumes that hadrons stem from a thermally equilibrated source, i.e. they are emitted from a particular “freeze-out” hypersurface. Each point at the hypersurface is characterized by the same values of all thermal parameters of a HRG. The generation of each event consists of three steps:

The multiplicities of all primordial hadrons are generated from the probability distribution which corresponds to the partition function of a particular HRG model. 2. 2.

Momenta of all primordial hadrons are generated, independently for each hadron. The momentum distribution is given by the spherical or longitudinal Blast Wave model, with a kinetic freeze-out temperature parameter $T_{\rm kin}$ that can be different from the chemical freeze-out temperature $T_{\rm ch}$. 3. 3.

Probabilistic decay chain of all resonances in the system is simulated until only the stable hadrons remain (optional).

It should be noted that the TEG is a rather simplistic event generator and has certain limitations. First, the present implementation is restricted to a HRG in the Boltzmann approximation, thus, the quantum statistical effects are omitted. Second, it is also assumed that the momenta and coordinates of the thermal hadrons at freeze-out are uncorrelated. The validity of such approximation can be questioned if two-particle correlations at freeze-out, e.g. due to the EV interactions, are non-negligible. Third, the modeling of decays assumes isotropic two-body and three-body decay kinematics, but provides only an approximate treatment of the isotropic many-body (four or more) decay kinematics. Therefore, an analysis of the observables which may be sensitive to the kinematics of many-body decays should be done with care.

3.1 Multiplicity sampling

The implementation of the multiplicity sampling is done according to the procedure described in Ref. [50]. Here only the basic details are described.

3.1.1 Poisson distribution

In the simplest case one has the Id-HRG model in the grand canonical ensemble. In this case the multiplicity distribution for each hadron species is given by the Poisson distribution, i.e. all multiplicities are described by the multi-Poisson distribution:

[TABLE]

where $\langle N_{i}\rangle\equiv\phi_{i}(T)\,e^{\mu_{i}/T}\,V$ is the mean number of hadron species $i$ in the GCE.

In the sampling procedure, first $\langle N_{i}\rangle=n_{i}\,V$ are calculated analytically, and then $\{N_{i}\}$ are generated for each event independently for each hadron species $i$ from the Poisson distribution.

3.1.2 Exact charge conservation

In the canonical ensemble all globally conserved charges are conserved exactly. In the Id-HRG model this condition modifies the multi-Poisson distribution as follows:

[TABLE]

where the index $k$ runs through all conserved charges in the system.

The implementation of exact charge conservation proceeds by combining the sampling from the multi-Poisson distribution with a rejection sampling. To speed up the process, we use the multi-step sampling procedure described in Ref. [79].

3.1.3 EV/QvdW interactions

A presence of the EV interactions leads to an appearance of a reduced volume in the partition function. In addition, the number of particles in the system is restricted from above, such that the total cumulative eigenvolume of all hadrons does not exceed the system volume. In the QvdW-HRG model, the presence of vdW interactions leads to the following modification of the microstate probability in the GCE (see Ref. [50] for details):

[TABLE]

The Monte Carlo procedure in this case uses the importance sampling technique. Namely, the $\{N_{i}\}$ multiplicities for each event are still generated from the multi-Poisson distribution, but each generated event is assigned a weight $w=\Pi_{\rm vdW}/\Pi$. Therefore, any observable is calculated as a weighted average. The theta function in Eq. (3.1.3) is taken into account via a rejection sampling as well: all $\{N_{i}\}$ configurations for which $\Theta(\{N_{i}\};V)$ evaluates to zero are rejected. The average multiplicities $\langle N_{i}\rangle$ used for the sampling from the auxiliary multi-Poisson distribution (41) are evaluated analytically using the GCE formulation of the QvdW-HRG model.

The multiplicity sampling for the QvdW-HRG model in the CE contains an additional rejection sampling step, as described in the previous subsection. The DEV-HRG and NDEV-HRG models follow from the QvdW-HRG model as partial cases, as elaborated in Sec. 2.4.

3.2 Momentum sampling

The thermal momenta of all the primordial hadrons are generated independently for each hadron. A possibility of collective motion at freeze-out is included in the framework of the blast-wave model. Two options are available: the spherically symmetric blast wave scenario [80], appropriate for intermediate collision energies, and the cylindrically symmetric blast wave scenario [81], which is more appropriate for high collision energies.

3.2.1 Spherically symmetric blast wave model

In the spherically symmetric blast wave model the momentum distribution is given by the Siemens-Rasmussen formula [80]:

[TABLE]

where $\gamma=(1-\beta^{2})^{-1/2}$, $\alpha=\gamma\,\beta\,p/T_{\rm kin}$, $\beta\in[0,1)$ is the radial flow velocity parameter, and $\tilde{N}$ is the normalization constant. $T_{\rm kin}$ is the kinetic freeze-out temperature, which can be taken different from the chemical freeze-out temperature that determines the multiplicity distribution. The momentum vector $\mathbf{p}$ is considered in the spherical coordinate basis, i.e.

[TABLE]

The spherical coordinate angles $\varphi_{p}$ and $\theta_{p}$ are generated assuming the isotropic distribution, i.e. $\varphi_{p}$ and $\cos\theta_{p}$ are uniformly and independently distributed in the $[0,2\pi)$ and $[-1,1]$ intervals, respectively. The distribution function for $p$ follows from Eq. (45):

[TABLE]

For practical purposes, it is convenient to perform the following variable change

[TABLE]

The new variable $\xi$ takes values in the finite interval $0\leq\xi<1$, which is convenient for a numerical implementation. The distribution function of $\xi$ reads

[TABLE]

The function $dN/d\xi$ goes to zero at the edges of the interval $0\leq\xi<1$, and a has a maximum in-between. Therefore, the value of $\xi$ is generated using the rejection sampling technique. The value of the maximum of $dN/d\xi$ is determined using the ternary search. The $p_{0}=1$ GeV value is used in the present implementation. After the value of $\xi$ is generated, the absolute value of momentum is calculated via Eq. (47).

The above procedure is repeated for all particles in the event.

3.2.2 Cylindrically symmetric blast wave model

The momentum distribution in the cylindrically symmetric blast wave model is given by [81]:

[TABLE]

where $m_{T}=\sqrt{p_{T}^{2}+m^{2}}$ is the transverse mass, $y=\displaystyle\frac{1}{2}\log\frac{\omega_{p}-p_{z}}{\omega_{p}+p_{z}}$ is the longitudinal rapidity, $\rho=\tanh^{-1}\beta_{r}$, and $\beta_{r}=\beta_{s}\tilde{r}^{n}$ is the transverse flow velocity profile, $\eta_{\rm max}$ is the longitudinal rapidity cutoff. $\beta_{s}$ is the transverse flow velocity at the surface. The mean transverse flow velocity is $\langle\beta_{T}\rangle=\frac{2}{2+n}\,\beta_{s}$. The particle three-momentum is parameterized as $\mathbf{p}=(p_{T}\,\cos\varphi_{p},p_{T}\sin\varphi_{p},m_{T}\sinh y)$.

Due to the azimuthal symmetry, the azimuthal angle $\varphi_{p}$ is distributed uniformly in the interval $[0,2\pi)$, and it is independent of $p_{T}$ and $y$ .

The distribution function for the transverse momentum $p_{T}$ reduces to [81]:

[TABLE]

The $p_{T}$ value is generated using the change of variable $p_{T}=-p_{0}\,\log\xi$, and the rejection sampling technique for the generation of $\xi$, using the same method as described in the previous subsection for the Siemens-Rasmussen formula.

The rapidity distribution at a fixed value of $p_{T}$ follows from Eq. (3.2.2):

[TABLE]

The value of $y$ is generated from Eq. (3.2.2) using the rejection sampling technique.

4 Thermal-FIST structure and implementation

Thermal-FIST is implemented as a C++ library of classes and functions. The applications of the package proceed by linking the ThermalFIST library within a C++ program (macro), which performs the needed task.

The library documentation is available online [82]. Here only the basic structure is described.

4.1 Base classes

4.1.1 IdealGasFunctions module

Calculation of thermodynamic properties of an ideal gas in the GCE forms the basis of all HRG calculations in Thermal-FIST. These are given as functions of the temperature $T$, the chemical potential $\mu$, the particle mass $m$, the internal degeneracy factor $d$, and the statistics $\eta$. Here $\eta=+1$ corresponds to the Fermi-Dirac statistics, $\eta=-1$ corresponds to the Bose-Einstein statistics, and $\eta=0$ corresponds to the Maxwell-Boltzmann statistics.

These functions are implemented in the IdealGasFunctions module, and include pressure, particle number density, energy density, entropy density, scalar density, and the leading four particle number susceptibilities $\chi_{k}$.

The calculations within the Maxwell-Boltzmann statistics ($\eta=0$) use known analytical expressions involving the modified Bessel functions. For example, the pressure reads

[TABLE]

Similar expressions exist for all other ideal gas functions.

There are two options for the numerical calculations of the ideal quantum gas functions ($\eta=\pm 1$): (i) using the cluster expansion technique, (ii) using the numerical integration with Gauss-Laguerre quadratures (the default method).

Cluster expansion

The first method uses the cluster expansion of an ideal quantum gas. The calculations proceed by calculating the truncated cluster expansion, e.g., the pressure is calculated as

[TABLE]

similar expressions are used for other thermodynamic functions. The value of $k_{\rm max}$ can be regulated, it should be sufficiently large for accurate calculations. Setting $k_{\rm max}=1$ one recovers the Maxwell-Boltzmann statistics. Note that cluster expansion formally diverges if $\mu>m$. Therefore, this method should not be used for calculating the Fermi-Dirac functions for $\mu>m$, and extra care taken for the case $\mu\lesssim m$.

Numerical integration

The numerical integration method proceeds by evaluating the integrals numerically, with the help of quadratures. First, the integrals are written in the dimensionless form,

[TABLE]

where $\tilde{k}=k/T$, $\tilde{m}=m/T$, and $\tilde{\mu}=\mu/T$.

For $\mu<m$, the calculations proceed by applying the 32-point Gauss-Laguerre quadrature to the dimensionless integral in Eq. (4.1.1).

For $\mu>m$, the Bose-Einstein integrals are divergent. The Fermi-Dirac integrals, however, are convergent and can be computed. The numerical accuracy of the above-described method, however, is usually unsatisfactory in such a case. Therefore, for calculating the Fermi-Dirac integrals at $\mu>m$, Eq. (4.1.1) is rewritten as follows [83, 84]:

[TABLE]

Here $\tilde{p}_{F}=\sqrt{\tilde{\mu}^{2}-\tilde{m}^{2}}$ and

[TABLE]

The first integral in Eq. (4.1.1) is evaluated using the 32-point Gauss-Legendre quadrature, while the second integral in Eq. (4.1.1) is evaluated using the 32-point Gauss-Laguerre quadrature.

Expressions similar to Eq. (4.1.1) are obtained for other quantities using the standard thermodynamic relations.

4.1.2 ThermalModelParameters

The HRG model parameters are kept in the ThermalModelParameters structure. Each implementation of the HRG model has this structure, which contains the values of thermal parameters used in calculations. List of all parameters is presented in Table 1.

4.1.3 ThermalParticle class

The ThermalParticle class contains all information about the properties of a particular hadron species in the hadron list. This includes the Particle Data Group (PDG) code, the mass, the degeneracy, the type of quantum statistics, the quantum numbers, absolute quark contents, and, if applicable, the decay channels and the width.

The class also provides various methods. The calculation of the ideal gas functions is performed through the Density method. These calculations take into account a possible additional integration due to a finite width. Other methods provide setting of various options, including whether the integrals are evaluated using the cluster expansion or the numerical integration, the number of terms in the cluster expansion (if this method is used), the shape and the prescription used to treat the finite resonance widths. More details can be found in the online documentation [82], in particular in the annotated ThermalParticle.h source file.

4.1.4 ThermalParticleSystem class

The ThermalParticleSystem class contains the list of particles to be used in calculations, represented by a vector of ThermalParticle objects, one per each particle species.

The particle list is usually provided from an external file. The LoadList method of the class loads the particle list from an external file. The external file which contains the particle list has a format of a table, where each row correspond to a single particle specie. Each column in the table corresponds to a particular hadron property. List of all columns is provided in Table 2.

4.1.5 ThermalModel classes

Implementation of the different HRG models is contained in classes starting with the prefix ThermalModel.

ThermalModelBase

ThermalModelBase is the base class for all ThermalModel classes. This is an abstract class containing fields and methods common for an arbitrary variant of HRG model. The exact implementation of some of these methods can be different in different versions of HRG, and is given in the corresponding derived classes. The ThermalModelBase class instance contains a pointer to a ThermalParticleSystem object containing the particle list, and a ThermalModelParameters structure with all the thermal model parameters. Both are provided during the creation of the ThermalModelBase instance. Some of the important methods of the class include:

Methods to set the various thermal model parameters. The list is given in Table 3.

• 2.

SetUseWidth method to specify the treatment of finite resonance widths.

• 3.

ConstrainMuB(bool), ConstrainMuQ(bool), ConstrainMuS(bool),

and ConstrainMuC(bool) methods to define whether the values of the corresponding chemical potentials should be constrained by, respectively, a fixed entropy per baryon ratio $S/B$, electric-to-baryon charge ratio $Q/B$, zero strangeness and zero charm.

• 4.

ConstrainChemicalPotentials() procedure constrains, if required, the chemical potentials $\mu_{B}$, $\mu_{Q}$, $\mu_{S}$, and/or $\mu_{C}$.

• 5.

CalculatePrimordialDensities() procedure calculates the primordial hadronic densities for current values of the thermal parameters. This is the main method were most calculations take place.

• 6.

CalculateFeeddown() procedure calculates various feeddown contributions to final yields (densities) of hadrons after the primordial densities were computed.

• 7.

double GetDensity(int PDGID, Feeddown::Type feeddown) returns the calculated density for a particle with the given PDG code and feeddown. The feeddown parameter can be 0 (primordial density), 1 (with feeddown according to the stability flags), 2 (with strong/electromagnetic/weak decay feeddown), 3 (with strong/electromagnetic decay feeddown), and 4 (with strong decay feeddown). This method should only be called after the primordial and final densities were calculated.

Classes which implement the different variants of the HRG model inherit from the ThermalModelBase class and are listed in Table 4. More details about the usage of the ThermalModel classes can be found in the corresponding source files and in the annotated sample routines in the src/examples folder.

4.2 ThermalModelFit class

The ThermalModelFit class implements routines related to the thermal fitting of hadron yields and/or yield ratios. The class instance contains a pointer to the ThermalModelBase class instance, provided on object initialization, as well as the set of parameters for fitting, and the data to fit.

The list of all possible fit parameters is shown in Table 5. If the parameters $\mu_{Q}$, $\mu_{S}$, or $\mu_{C}$ are not fitted, then they are constrained by the corresponding $Q/B$ ratio and strangeness/charm neutrality conditions. Obviously, (some of) these parameters are not used in the fit also if CE/SCE/CCE ensemble is employed where these parameters do not appear at all.

The data points to fit can be either read from an external file with the loadExpDataFromFile(string filename) method, or added manually one-by-one through the AddQuantity method. See the annotated

src/examples/BagModelFit/BagModelFit.cpp macro for an example of using the ThermalModelFit class.

4.3 Event generator

The thermal event generators with momentum distribution given by the spherical or longitudinally symmetric Blast-Wave model are implemented in SphericalBlastWaveEventGenerator and LongitudinalBlastWaveEventGenerator classes, respectively. The underlying thermal model configuration, as well as the momentum spectrum parameters described in Sec. 3, are provided through the constructor. A new event can be generated with the GetEvent() method and written into file with the GetEvent().writeToFile() method.

4.4 Graphical user interface (GUI)

The package contains QtThermalFIST module – a GUI frontend for the Thermal-FIST library, where some general purpose HRG model tasks can be performed in a convenient way. The GUI is written using the cross-platform framework Qt5. QtThermalFIST consists of five tabs.

In the Thermal model tab a comprehensive HRG model calculation can be performed at given values of thermal parameters. It is possible to analyze the primordial and total yields of all hadron species, particle number fluctuations, and the equation of state properties. Calculations can be performed in the grand canonical or canonical ensemble, within the Id-HRG, DEV-HRG, NDEV-HRG, or QvdW-HRG model, using a number of common eigenvolume parametrizations. Inclusion of quantum statistics is optional. There is also a possibility to constrain electric, strangeness, and/or charm chemical potentials from conservation laws, and to specify different prescriptions for the treatment of finite resonance widths.

Thermal fitting of the hadron production data is performed in the Thermal fits tab. It is possible to specify which parameters are fitted (and in which range) and which are fixed. The experimental data can be loaded from an external file. Some samples are provided in the input/data folder. Possibility to input the data directly inside the GUI is provided as well. It is also possible to view the fit results in a form of thermal fit plots. Another useful feature is calculation of the $\chi^{2}$ profiles.

The Equation of state tab offers a possibility to study the temperature dependence at a fixed $\mu_{B}$ of some common equation of state observables, conserved charges susceptibilities, and particle number densties (primordial or with feeddown). It is also possible to consider ratios of any pair of these observables. At finite $\mu_{B}$ the chemical potentials $\mu_{Q}$ and $\mu_{S}$ can be fixed from a fixed $Q/B$ ratio and zero net strangeness, or set to zero otherwise. For a number of observables at $\mu_{B}=0$ the published lattice QCD data of the Wuppertal-Budapest and/or HotQCD collaborations is plotted along with the calculation results for convenience.

The Event generator tab can be used for generating events with the TEG. The GUI program shows histograms with momentum distributions of the different generated hadrons. It is also possible here to write the generated events to a specified file.

The Particle list editor tab provides a user-friendly editor of the particle list. The edited list can be saved to a file and/or used on the fly in HRG model calculations in the other tabs.

In general, most of the options in QtThermalFIST are self-explanatory, a more detailed description of the usage can be found in the Quick Start Guide which comes with the package.

5 Installation

The current version of the source code of Thermal-FIST can be obtained from https://github.com/vlvovch/Thermal-FIST

The package is platform-independent. The preferred method to configure the package is to use cmake. An example for a Linux system:

⬇

Clone the repository from GitHub

git clone https://github.com/vlvovch/Thermal-FIST.git

cd Thermal-FIST

Create a build directory, configure the project with cmake

and build with make

mkdir build

cd build

cmake ../

make

Run the GUI frontend

./bin/QtThermalFIST

Run the test calculations from the paper

./bin/examples/cpc1HRGTDep

./bin/examples/cpc2chi2

./bin/examples/cpc3chi2NEQ

./bin/examples/cpc4mcHRG

The above commands will build the Thermal-FIST package in the build directory within the root folder of the package, run the GUI frontend program, QtThermalFIST, if it was built, and various test calculations which are presented below. More specifically, the library will be located in build/lib directory, the QtThermalFIST program in build/bin, and the sample macros for using the library in build/bin/examples. Note that in order to build the QtThermalFIST GUI one may need to install first the freely available Qt5 framework [85].

Automated tasks can be performed by writing C++ macros which link to the Thermal-FIST library and perform the necessary calculations. Sample macros in src/examples can be used as templates.

6 Test results

6.1 Thermodynamic properties at zero chemical potential

Here we consider the temperature dependence of thermodynamic functions, namely the pressure and the energy density, calculated at $\mu=0$ within different variants of the HRG model. The Thermal-FIST results are compared to the calculations performed using the THERMUS-2.3 package. In order to ensure a consistent comparison, the Thermal-FIST calculations here employ the hadron list from THERMUS-2.3. Quantum statistics and the energy independent Breit-Wigner scheme for resonance widths are used in both codes.

Three different variants of the HRG model are considered. The first variant is the standard Id-HRG model. The second variant is the EV-HRG model with a constant radius parameter $r=0.3$ fm assigned to all hadron species. Equation of state within this particular EV-HRG model was earlier considered in Refs. [65, 49]. The third variant is the QvdW-HRG model, which includes vdW interactions between baryon-baryon and antibaryon-antibaryon pairs, with vdW parameters $a=329$ MeV fm3 and $b=3.42$ fm3, common for all pairs. This model was formulated in Ref. [29].

The temperature dependence of the scaled pressure $p/T^{4}$ and the scaled energy density $\varepsilon/T^{4}$ calculated for the three different HRG models described above are depicted in Fig. 2 by lines (FIST) and symbols (THERMUS-2.3). One can see that the EV/vdW effects kick in at higher temperatures, where hadronic densities become large. The Thermal-FIST and THERMUS-2.3 results for Id-HRG and EV-HRG models are quantitatively consistent with each other. It is unfortunately not possible to present such a comparison for the QvdW-HRG model since THERMUS-2.3 does not contain an implementation of this model.

The annotated macro, which performs the above calculations, can be found in the following location: src/examples/cpc/cpc1-HRG-TDep.cpp

6.2 Thermal fits to heavy-ion hadron yield data

6.2.1 Excluded volume/van der Waals effects

In this subsection we consider the thermal fits to hadron yield data. It is illustrated how different EV/vdW parametrizations influence the fits. For testing purposes, we consider here those EV/vdW configurations which have a strong influence on thermal fits. In particular, we redo here some of the calculations which were presented in Ref. [45]. We take the hadron yield data measured by the ALICE collaboration in 0-5% most central Pb-Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV. The actual data used in the fits can be found in Table 1 of Ref. [45]. We assume $\mu=0$, and consider the temperature dependence of $\chi^{2}$ of the fit. To calculate this temperature dependence, we minimize $\chi^{2}$ at each temperature by fitting the system volume parameter.

Four different variants of the HRG model are considered here.

The ideal HRG model. 2. 2.

The DEV-HRG model with mass-proportional eigenvolumes, $v_{i}=v_{p}\,(m_{i}/m_{p})$, where $v_{p}=(16/3)\pi r_{p}^{3}$ with $r_{p}=0.5$ fm is the assumed proton eigenvolume, and $m_{p}=0.938$ GeV/$c^{2}$ is the proton mass. This is the bag model HRG considered in Ref. [45]. 3. 3.

The DEV-HRG model with point-like mesons, and finite-sized (anti)baryons with $v_{B}=(16/3)\pi r_{B}^{3}$, $r_{B}=0.3$ fm. This HRG model was considered in Refs. [45, 65]. 4. 4.

The QvdW-HRG model which includes the QvdW interactions for baryon-baryon and antibaryon-antibaryon pairs, with QvdW parameters $a=329$ MeV fm3 and $b=3.42$ fm3, common for all pairs. This model was formulated in Ref. [29].

As in the previous example, the calculations are performed using both Thermal-FIST and THERMUS-2.3. Therefore, Thermal-FIST calculations use the hadron list from THERMUS-2.3. Quantum statistics and the energy independent Breit-Wigner modeling of finite resonance widths are used in both codes. In order to correctly calculate the feeddown contribution to the yield of the unstable $\phi$ meson, we apply the fix [86] to THERMUS-2.3.

The temperature dependence of $\chi^{2}$ of the fit to the ALICE data, performed within the four models described above, is depicted in Fig. 3 by lines (Thermal-FIST) and symbols (THERMUS-2.3). The Id-HRG model shows a single $\chi^{2}$ minimum at $T\sim 155$ MeV, the result reported in the literature by many groups [34, 87, 88, 89]. Calculations within the non-ideal HRG models show peculiar two-minimum structures in Fig. 3, the mathematical origin of these structures was discussed in Refs. [45, 27]. Note that physical interpretation of the high-temperature 2nd minima should be done with care, lattice QCD calculations suggest that the crossover transition to quarks and gluons may already be completed at lower temperatures [90, 91]. The calculations are presented here merely in the context of code testing and comparison. The calculations within the Id-HRG and the two EV-HRG models are quantitatively consistent between Thermal-FIST and THERMUS-2.3. There are no THERMUS-2.3 calculations for the QvdW-HRG model available due to a lack of its implementation in THERMUS-2.3.

The annotated macro, which performs the above calculations, can be found in the following location: src/examples/cpc/cpc2-chi2-vs-T.cpp

6.2.2 Chemical non-equilibrium fits

Previous considerations were restricted to the chemical equilibrium HRG model, i.e. for $\gamma_{q}=\gamma_{S}=1$. The chemical non-equilibrium scenario, $\gamma_{q}\neq 1$, $\gamma_{S}\neq 1$, for hadron production in heavy-ion collisions was advocated in Refs. [11, 34, 92]. Fits within the chemical non-equilibrium scenario lead to significantly smaller $\chi^{2}/dof$ values. The resulting deviations from chemical equilibrium results are significant, e.g. the extracted temperature is about 15-20 MeV lower, the extracted $\gamma_{q}$ values are in the 1.6-1.7 range, very close to the Bose-Einstein condensation singularity for pions. Therefore, the improvement in fit quality comes at the cost of abandoning the chemical equilibrium scenario for heavy-ion collisions.

Here we do not discuss which physical scenario should be preferred. Instead, we verify whether the Thermal-FIST package can reproduce the previously published results obtained in the chemical non-equilibrium scenario. To our knowledge, the chemical non-equilibrium calculations were previously restricted to the usage of SHARE package.

We consider the Id-HRG model fits in Thermal-FIST within chemical equilibrium ($\gamma_{q}=\gamma_{S}=1$) and chemical non-equilibrium ($\gamma_{q}\neq 1$, $\gamma_{S}\neq 1$) scenarios. The data fitted include the 4$\pi$ multiplicities reported by the NA49 collaboration [93, 94, 95, 96, 97, 98] for most central Pb-Pb collisions at $\sqrt{s_{NN}}=7.6,\,8.8,\,12.3$, and $17.3$ GeV, and the midrapidity yields measured by the ALICE collaboration in 0-5% most central Pb-Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV. These data sets are consistent with the ones analyzed with SHARE in Refs. [11, 34].

There are some technical differences between the current analysis and and previously published SHARE analyses. First, the electric charge and strangeness chemical potentials $\mu_{Q}$ and $\mu_{S}$ are not fitted here, but are fixed from the conditions of the electric-to-baryon charge ratio $Q/B=0.4$ and strangeness neutrality $S=0$. Second, there are differences in the hadron list employed in Thermal-FIST and SHARE, and there are also differences in the treatment of the finite resonance widths. Therefore, some differences between the two codes are expected, the important question is whether the qualitative features reported in Refs. [11, 34] can be reproduced.

The $T$-$\mu_{B}$ values extracted from the thermal fits to NA49 and ALICE data within the chemical non-equilibrium and chemical equilibrium scenarios are depicted in Fig. 4a by full and open red squares, respectively. The SHARE non-equilibrium results from Refs. [11, 34] are shown by the blue circles. The Thermal-FIST and SHARE results are consistent with each other to a good precision. The $\gamma_{q}$ values extracted from fits within Thermal-FIST are in the range $\gamma_{q}=1.6-1.7$, close to the onset of pion Bose-Einstein condensation, and consistent with the values reported in Refs. [11, 34]. Figure 4b depicts the $\chi^{2}/dof$ values obtained in chemical (non-)equilibrium fits within Thermal-FIST. The reduced $\chi^{2}$ is significantly smaller in the chemical non-equilibrium scenario, as reported in Refs. [11, 34]. The presented calculation shows that Thermal-FIST is able to reproduce previously published systematics of the chemical non-equilibrium scenario.

The annotated macro, which performs the above calculations, can be found in the following location: src/examples/cpc/cpc3-chi2NEQ.cpp

6.3 Analytic calculations vs Monte Carlo

The last section here illustrates applications of the Monte Carlo TEG. The TEG is most useful when an analytic approach is unavailable or problematic. One such example is a simultaneous inclusion of the CE and EV effects. This procedure is described in detail in Ref. [50]. Another possibility is a study of fluctuations and correlations of various hadron yields. Various effects, such as radial flow and kinematic cuts can be naturally included in the Monte Carlo approach whereas the whole procedure of the event-by-event analysis resembles closely the experimental situation.

In this section we consider a comparison of the analytic and Monte Carlo results for various 2nd order fluctuations and correlations of hadron yields. Such a comparison serves as an important cross-check of consistency between analytic and Monte Carlo methods in general, and of the probabilistic decays implementation in Sec. 2.7.2 in particular.

More specifically, we consider the mean-to-variance ratio of the net-kaon distribution $M_{k}/\sigma_{2}^{k}$, the ratio $\sigma_{11}^{Qk}/\sigma_{2}^{k}$ of the correlator between the total electric charge and final net kaon number over the variance of the final net kaon number, the ratio $\sigma_{11}^{pQ}/\sigma_{2}^{p}$ of the correlator between the final proton number and the total electric charge over the variance of the final net proton number, and the ratio $\sigma_{11}^{pk}/\sigma_{2}^{k}$ of the correlator between the final proton number and the final net kaon number over the variance of the final net kaon number, in the framework of the Id-HRG model. The beam energy dependence of $M_{k}/\sigma_{2}^{K}$ measured in Au-Au collisions was recently published by the STAR collaboration [99], while the experimental analysis of the correlation observables is ongoing.

The quantities $\sigma_{2}^{k}$, $\sigma_{2}^{p}$, $\sigma_{11}^{Qk}$, $\sigma_{11}^{pQ}$, and $\sigma_{11}^{pk}$ can be expressed as linear combinations over the correlators (30) of the number of final state hadrons, which are calculated within the Thermal-FIST. The Thermal-FIST calculations here are performed along the chemical freeze-out curve of Ref. [100], the resulting beam energy dependences are shown in Fig. 5 by red lines. The red lines in Fig. 5 correspond to analytic treatment of fluctuations, including the effects of probabilistic decays.

The calculation results for the same observables using the Monte Carlo TEG are shown by the red symbols. These are consistent with the analytic results, thus verifying the accuracy of the probabilistic decay treatment implemented in Thermal-FIST. In contrast, the simplified “average” decay treatment procedure [39] leads to a markedly different result for $M_{k}/\sigma_{2}^{K}$ [dashed black curve in Fig. 5(a)] in comparison to probabilistic decays and to Monte Carlo result. Full probabilistic treatment of resonance decays thus appears to be important for interpretation of the corresponding data within the HRG approach.

The quantities $\sigma_{11}^{Qk}/\sigma_{2}^{K}$, $\sigma_{11}^{pQ}/\sigma_{2}^{p}$, and $\sigma_{11}^{pk}/\sigma_{2}^{k}$ are considered as possible proxies for the corresponding ratios of susceptiblities of conserved charges, $\chi_{11}^{QS}/\chi_{2}^{S}$, $\chi_{11}^{BQ}/\chi_{2}^{B}$, and $\chi_{11}^{BS}/\chi_{2}^{S}$. The Thermal-FIST calculations for the latter are shown in Fig. 5 by the black lines (analytical) and symbols (Monte Carlo), and these differ substantially from the measurable correlator ratios. These significant differences should be considered when interpreting the experimental data in the context of the equation of state of QCD matter.

The annotated macro, which performs the above calculations, can be found in the following location: src/examples/cpc/cpc4-mcHRG.cpp

7 Concluding remarks

The Thermal-FIST package provides the ability to perform an analysis of both, the hadronic part of the QCD equation of state and of the statistical description of hadron production in heavy-ion collisions.

The powerful graphical user interface in the package will be useful for intepreting the future hadron yield data which will come in the ongoing and future heavy-ion experiments.

Generalization and development of different variants of the HRG model is also presently an active research topic. In that regard, the flexible modular structure of Thermal-FIST is particulary suitable for continuous improvement of the package and implementation of various new features to test different new ideas.

Acknowledgments

We would like to thank P. Alba, D. Anchishkin, V. Begun, E. Bratkovskaya, B. Dönigus, M. Gorenstein, I. Kisel, M. Lorenz, D. Oliinychenko, and L. Satarov for fruitful discussions and physics suggestions. We also acknowledge useful feedback from F. Flor, R. Hensch, A. Motornenko, R. Poberezhnyuk, P. Parotto, M. Puccio, and J. Stumm. H.St. acknowledges the support through the Judah M. Eisenberg Laureatus Chair at Goethe University, and the Walter Greiner Gesellschaft, Frankfurt.