\section{Introduction}
\label{sec:intro}
Significant theoretical and experimental effort is dedicated to mapping out the QCD phase diagram in the temperature $T$ and baryon chemical potential $\mu_B$ plane, and to search for the elusive critical point \cite{Aoki:2006we,Aoki:2006br,Borsanyi:2010bp,Bazavov:2011nk,Borsanyi:2013bia,Bazavov:2014pvz,Fodor:2018wul,Busza:2018rrf,Stephanov:2011pb,Stephanov:1998dy,Ratti:2006gh,Fukushima:2010bq,Critelli:2017oub,Luo:2017faz,Parotto:2018pwx,Grefa:2021qvt,Karthein:2021nxe,Adamczyk:2017iwn,Adamczewski-Musch:2020slf,Abelev:2013vea,Fu:2021oaw} (for recent reviews see e.g. \cite{Ding:2015ona,Ratti:2018ksb,Bzdak:2019pkr,Ratti:2021ubw}). Ultrarelativistic heavy-ion collisions are generating the deconfined phase of strongly interacting matter in the laboratory. By systematically decreasing the collision energy of the incoming nuclei, the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory is scanning the phase diagram in the so-called Second Beam Energy Scan (BESII), soon to be followed by even lower collision energies at NICA and at the GSI-FAIR accelerator.

First-principles lattice QCD simulations are available for several thermodynamic quantities, such as the equation of state at zero \cite{Borsanyi:2010cj,Borsanyi:2013bia,Bazavov:2014pvz} and small chemical potential \cite{Guenther:2017hnx,Gunther:2017sxn,Bazavov:2017dus,Borsanyi:2021sxv,Mondal:2021jxk}, QCD transition line \cite{Bellwied:2015rza,Bazavov:2018mes,Borsanyi:2020fev,Borsanyi:2021sxv}, as well as diagonal \cite{Borsanyi:2011sw,Borsanyi:2013hza,Bellwied:2013cta,Borsanyi:2014ewa,Bellwied:2015lba,Noronha-Hostler:2016rpd,Bazavov:2017tot,Borsanyi:2018grb,Bazavov:2020bjn} and off-diagonal \cite{Borsanyi:2018grb} fluctuations of conserved charges. However, they are currently limited to the low baryonic chemical potential regime, due to the fermionic sign problem. Effective models that can reproduce lattice QCD results in certain regimes of temperature and chemical potential are therefore very useful to extend the coverage of the phase diagram beyond the reach of lattice QCD.
The models are also necessary to make a connection to the common heavy-ion observables, such as the measurements of various distributions of identified hadrons.

In the low-temperature regime ($T \lesssim 160$), the bulk thermodynamics of QCD is generally well described by a multi-component gas of free hadrons and resonances \cite{Bazavov:2012jq,Borsanyi:2010bp,Vovchenko:2016rkn,Alba:2017mqu} the so-called Hadron Resonance Gas (HRG) model. This indicates that hadron interactions in this regime may be dominated by the formation of known resonances. The HRG model has been widely used to study the confined phase of QCD matter below the transition line \cite{Karsch:2003vd,Karsch:2003zq,Tawfik:2004sw,Huovinen:2009yb,Ratti:2010kj,Alba:2015iva,Huovinen:2017ogf,Bellwied:2019pxh,Vovchenko:2014pka}. Its remarkable agreement with the equation of state from first-principles lattice calculations has led to its popularity, especially in the study of the chemical freeze-out in HICs \cite{Becattini:2012xb,Alba:2014eba,Vovchenko:2015idt,Andronic:2017pug,Bellwied:2018tkc,Alba:2020jir,Bluhm:2018aei,Flor:2020fdw}. 
However, with the availability of more differential observables like susceptibilities of conserved charges, discrepancies between the predictions of the fundamental theory and the HRG model have been observed \cite{Bazavov:2014xya,Borsanyi:2018grb,Bellwied:2021nrt}, specifically at temperatures $T \sim 150-160$~MeV that characterize the chemical freeze-out in heavy-ion collisions at the highest energies. In particular, lattice QCD results for the partial pressures have shown a need for more resonances in the strange sector than those which are already experimentally well known \cite{Bazavov:2014xya,Alba:2017mqu}. On the other hand, some susceptibility ratios including $\chi_4^B/\chi_2^B$ and $\chi_4^S/\chi_2^S$ suggest the need for repulsive interactions \cite{Bellwied:2015lba,Bazavov:2018mes,Alba:2017mqu}.

Several extensions of the HRG model have thus been proposed to improve the agreement with lattice QCD.
One possibility is extending the hadronic spectrum -- the model input -- to include more states not yet observed \cite{Bazavov:2014xya,Alba:2017mqu,Alba:2020jir}.
Other extensions incorporate additional, non-resonant interactions between hadrons such as excluded volume~\cite{Yen:1997rv,Andronic:2012ut,NoronhaHostler:2012ug,Bhattacharyya:2013oya,Vovchenko:2014pka,Albright:2015uua,Satarov:2016peb,Vovchenko:2017xad,Alba:2017bbr,Vovchenko:2018eod,Motornenko:2020yme}, van der Waals~\cite{Vovchenko:2016rkn,Vovchenko:2017zpj,Samanta:2017yhh,Sarkar:2018mbk,Vovchenko:2020lju}, mean field~\cite{Huovinen:2017ogf,Steinert:2018zni}, or are based on scattering phase shifts~\cite{Venugopalan:1992hy,Friman:2015zua,Vovchenko:2017drx,Lo:2017lym,Dash:2018can,Fernandez-Ramirez:2018vzu,Dash:2018mep}. 
While one of the advantages of the standard HRG model is certainly the lack of free parameters, apart from the uncertainties in the hadronic spectrum, introducing additional interactions unavoidably leads to new free parameters that need to be constrained through comparison with lattice results. We propose a combination of these different corrections to the standard HRG model and investigate its agreement with several lattice results on the susceptibilities. 

In this manuscript, we consider two HRG model extensions: the excluded volume interaction in the baryon sector, and the inclusion of additional particles in the hadronic list, beyond those experimentally observed. We emphasize the complementary nature of these two extensions and identify combinations of susceptibilities of conserved charges that allow to constrain them separately. The resulting HRG model considerably improves the description of lattice QCD results, in particular those for fluctuations of conserved charges.
The model can thus be useful for the analysis of freeze-out in heavy-ion collisions, in particular event-by-event fluctuations.


\section{The HRG Model and its modifications}\label{sec:HRG_ext}
\subsection{Ideal HRG}
The partial pressure for particle species $i$ in the ideal HRG model can be written as:
\begin{equation} \label{eq:press_partitionfxn}
\centering
\begin{split}
    &P_i (T, \mu_B, \mu_Q, \mu_S) = \frac{d_i T}{2\pi^2} \int_0^{\infty} (-1)^{B_i+1} k^2 dk \, \times
    \\ &\times \ln {\left[1 + (-1)^{B_i+1} \lambda_i(T, \mu_i) \exp\left(-\frac{\sqrt{k^2 + m_i^2}}{T}\right)\right]},
\end{split}
\end{equation}
where $d_i$ is the spin degeneracy factor, $B_i$ is the baryon number of species $i$, $k$ is the momentum, $\lambda_i (T,\mu_i) = \exp[(B_i \mu_B + Q_i \mu_Q + S_i \mu_S)/T]$ is the particle fugacity, and $m_i$ is the mass of species $i$. 
The partial pressure can be presented as a series containing the modified Bessel function of the second kind by expanding the logarithm and integrating term-by-term:
\begin{equation} \label{eq:press_Bessfxn}
\centering
\begin{split}
    &P_i (T, \mu_B, \mu_Q, \mu_S) = \\
    &\frac{d_i T^2}{2\pi^2} \sum_{N=1}^{\infty} \frac{[(-1)^{B_i+1}]^{N+1} \lambda_i(T, \mu_i)}{N^2} m_i^2 K_2\left(N\frac{m_i}{T}\right).
\end{split}
\end{equation}
Taking the first term in the expansion corresponds to the Boltzmann approximation.
This approximation is sufficient for our purposes and will be used throughout.
The pressure attains the following convenient form:
\begin{equation} \label{eq:press_phi}
\centering
\begin{split}
    P_i (T, \mu_B, \mu_Q, \mu_S) &= d_i \tilde{\phi}(T,m_i) \lambda_i(T,\mu_i), \\
    \tilde{\phi}(T,m_i) &=  \frac{m_i^2 T^2}{2\pi^2} \, K_2(m_i/T).
\end{split}
\end{equation}
The full pressure in the ideal HRG model corresponds to the sum of partial pressures of all hadronic components.
It is convenient to group the contributions of the various hadrons in accordance with their quantum numbers.
For instance, introducing
\begin{equation} \label{eq:press_sectors}
\centering
    \tilde{\phi}(T) = \sum_{j \in \rm{sectors}} d_j \tilde{\phi}(T,m_j)
\end{equation}
allows us to identify the different sectors of the total HRG pressure broken down by various quantum numbers that are listed in Table \ref{tab:hadron_table}.
Note that all species apart from the $i = 0$ sector~(neutral particles) 
have a corresponding antiparticle sector. 
\begin{table}[h]
\centering
\begin{tabular}{c|c c c|c}
 $i$ & $B$ & $Q$ & $S$ & base hadron \\ \hline
 0 &  0 & 0  & 0 & $\pi^0$ \\ 
 1 &  0 & 1 & 0 & $\pi^+$ \\ 
 2 &  0 & 1 & 1 & K$^+$ \\ 
 3 &  0 & 0 & 1 & K$^0$ \\ 
 4 &  1 & 0 & 0 & n  \\ 
 5 &  1 & 1 & 0 & p  \\
 6 &  1 & 2 & 0 & $\Delta^{++}$  \\
 7 &  1 & -1 & 0 & $\Delta^{-}$ \\
 8 &  1 & 0 & -1 & $\Lambda$ \\
 9 &  1 & 1 & -1 & $\Sigma^+$ \\
 10 & 1 & -1 & -1 & $\Sigma^-$ \\
 11 & 1 & 0 & -2  & $\Xi^0$ \\
 12 & 1 &  -1 &  -2 & $\Xi^-$ \\
 13 & 1 & -1 & -3  & $\Omega^-$
\end{tabular}
\caption{
\label{tab:hadron_table}
The list of hadronic quantum number sets contributing to the pressure of the ideal HRG model.
The last column identifies the lowest mass hadron representing the set of quantum numbers.
}
\end{table}
Therefore, the pressure from  Eq.~(\ref{eq:press_phi}) takes the form of a truncated relativistic expansion in fugacities:
\begin{equation} \label{eq:Pfug}
    \centering
    P (T, \mu_B, \mu_Q, \mu_S) = \tilde{\phi_0}(T) + \sum_{i \neq 0} 2 \, \tilde{\phi_i}(T) \cosh \left(\mu_i/T \right),
\end{equation}
where $\mu_i=B_i \mu_B + Q_i \mu_Q + S_i \mu_S$ is the chemical potential of the corresponding $i$th sector.
Each term in Eq.~\eqref{eq:Pfug} corresponds to the partial pressure associated with the particular set of hadronic quantum numbers.


\subsection{Extended Spectrum}
As a correction to the standard HRG model, our first extension is the incorporation of hadronic states beyond those which are well-known experimentally. When considering the partial pressures from Lattice QCD, it has been shown that the hadronic spectrum in QCD goes beyond what exists in the ordinary version of the hadronic list from the Particle Data Group (PDG) with only the most well-known states \cite{Zyla:2020zbs}. 
The PDG ranks particles by how well established they are, with a rating based on a number of stars (*). 
The **** states are those which are very well established, like e.g. nucleons or $\Delta(1232)$ resonances. 
On the opposite side, the * states are the least established ones, for example $\Delta$(1750) and other high-mass resonances. 
A previous investigation of the agreement between the partial pressures from the lattice and those obtained in the HRG model with different hadronic lists across the various quantum number sectors, found that the PDG2016+ particle list reproduced the largest number of quantities calculated on the lattice without exceeding them~\cite{Alba:2017mqu}. 
This hadronic list incorporates all states from the 2016 version of the PDG~\cite{Patrignani:2016xqp}, starting from the well-known ones, all the way down to those listed as “seen”, with a confidence rating of *. 
Another possibility is to incorporate the states predicted by the Quark Model~\cite{Capstick:1986bm,Godfrey:1985xj,Ebert:2009ub}, which includes an even larger number of hadrons than those contained in the PDG2016+.
As we can see in Eq. (\ref{eq:press_phi}), any additional states seek to increase the pressure of the system, i.e. a larger number of states will lead to a larger overall pressure.

It was recently pointed out in Ref. \cite{Bollweg:2021vqf} that some states calculated within the Quark Model from Ref. \cite{Capstick:1986bm} are overlapping with states later measured by the PDG \cite{Zyla:2020zbs}. 
In light of this, we provide an update to the QM list first published in Ref. \cite{Alba:2020jir} that removes all the duplicate states that remained as an artifact of those early Quark Model calculations \cite{listrepo}.
We note that the calculations here denoted as QM are using this updated list and that we observed a minimal difference between the old and new list when calculating the fluctuations of interest in this study.

In this study, we aim at identifying the most suitable description of the lattice data, with both the excluded volume interactions and additional hadronic states incorporated.
In order to do so, we revisit the comparison with lattice data and investigate several different hadronic lists: 
\begin{itemize}
    \item PDG2016 -- ordinary hadronic list with only the well-known states *** $-$ **** from the 2016 Particle Data Booklet;
    \item PDG2016+ -- the list containing both the established~(***$-$****) and unconfirmed~(*$-$**) states;
    \item Quark Model~(QM) -- the list which incorporates all states predicted by the Quark Model.
\end{itemize}
The latter two lists were introduced and described in detail in \cite{Alba:2017mqu, Alba:2020jir}. 
We also checked that the most recent compilation of the established states from the Particle Data Group -- the PDG2020 list -- yields negligible differences compared to the PDG2016 list, thus we retain the latter list in the analysis for consistency with Refs.~\cite{Alba:2017mqu,Alba:2020jir}.
\subsection{Excluded Volume}

The next extension to the HRG model is the excluded volume model. This corresponds to including repulsive interactions between hadrons. Many versions of the EV-HRG model have been considered in the literature. Here we follow the approach introduced in Refs.~\cite{Vovchenko:2016rkn,Vovchenko:2017xad} where EV interactions are included only for baryon-baryon and antibaryon-antibaryon pairs. 
This corresponds to a minimalistic EV extension that does not affect meson-meson and meson-baryon interactions, which are presumed to be dominated by resonance formation and thus already included in the HRG model.
The pressure is partitioned into contributions of non-interacting mesons and interacting baryons and antibaryons:
\begin{equation} \label{eq:pevhrg}
\centering
p = p_M^{\rm id} + p_B^{\rm ev} + p_{\bar{B}}^{\rm ev},
\end{equation}
where
\begin{equation} \label{eq:pMev}
p_{M}^{\rm id}  = \tilde{\phi}_0(T) ~ + \sum_{i \neq 0, \, i \in M} 2 \, \tilde{\phi}_i(T) \, \cosh(\mu_i/T),
\end{equation}
\begin{equation} \label{eq:pBev}
p_{B(\bar{B})}^{\rm ev}  =  \sum_{i \in B} \tilde{\phi}_i(T) \, \exp(\pm \mu_i/T) \,
\exp\left( \frac{- b \, p_{B(\bar{B})}^{\rm ev}}{T} \right)~.
\end{equation}
Here, $i \in M$ corresponds to mesons~($B_i = 0$), $i \in B$ corresponds to baryons~($B_i = 1$), $b$ is the baryon excluded volume parameter, and $\tilde{\phi}(T)$ is given in Eq. (\ref{eq:press_phi}).
Equation~\eqref{eq:pBev} can be solved in terms of the Lambert W function~\cite{NoronhaHostler:2012ug,Taradiy:2019taz}:
\begin{equation}\label{eq:pBevW}
p_{B(\bar{B})}^{\rm ev}  = \frac{T}{b} \, W[\varkappa_{B(\bar{B})}(T,\mu_B,\mu_Q,\mu_S)],~
\end{equation}
where
\begin{equation}
\varkappa_{B(\bar{B})}(T,\mu_B,\mu_Q,\mu_S) = b \, \sum_{i \in B} \tilde{\phi}_i(T) \, \exp(\pm \mu_i/T).
\end{equation}
The explicit form Eq.~\eqref{eq:pBevW} for the pressure in the EV-HRG model in terms of the Lambert W function allows us to forgo solving the transcendental equation for the pressure.

One should note that different formulations of the excluded volume HRG model exist in the literature. 
In the formulation that we use the EV interactions are introduced only for baryon-baryon and antibaryon-antibaryon pairs with a common EV parameter $b$.
This is consistent with the model used in Refs.~\cite{Vovchenko:2016rkn,Vovchenko:2017xad}.
In other EV models~\cite{Andronic:2012ut,Alba:2017bbr,Motornenko:2020yme} the repulsive interactions are introduced for \emph{all} hadron pairs, and each hadronic species may be characterized by its own value of the EV parameter.
As such, these types of EV models contain many more free parameters, while ours is a minimalistic approach.
The EV effects influence the thermodynamics differently in the latter class of EV models, and the conclusions obtained in this work within the former class of EV models do not necessarily translate.



\section{Susceptibilities}
\label{sec:susc}

The two HRG model extensions introduced above -- the excluded volume corrections and the inclusion of extra states -- are complementary to one another.
This can be seen in the following way.
Adding extra resonances can be interpreted as adding attractive interactions among hadrons that lead to the formation of these resonances.
Given the fact that all the extra states have baryon number equal to either 0 or $\pm 1$, this may correspond to meson-meson and meson-baryon interactions, but not to baryon-baryon interactions. 
On the other hand, the EV corrections considered here correspond to repulsive baryon-baryon interactions but not to any meson-meson or meson-baryon interactions. 
Therefore, the two extensions describe different physics, and thus, they can and should be considered simultaneously.

Both extensions affect the equation of state.
For instance, the inclusion of extra states increases the pressure at a given temperature and chemical potential, while the EV interactions lead to its suppression.
It can be challenging, then, to constrain the two effects separately.
In order to achieve those constraints, we study differential observables that have recently been obtained on the lattice, namely the susceptibilities of conserved charges:
\begin{equation} \label{eq:chis}
\chi^{BQS}_{lmn} = \frac{\partial^{l+m+n}(p/T^4)}{\partial(\mu_{B}/T)^{l} \, \partial(\mu_{Q}/T)^{m} \, \partial(\mu_{S}/T)^{n}}
\end{equation}

The susceptibilities in Eq.~\eqref{eq:chis} can be calculated explicitly in the EV-HRG model with extra states by utilizing Eq.~\eqref{eq:pBevW} and the known properties of the Lambert W function.
One can then construct specific combinations of susceptibilities that are mainly sensitive to either the extra states or baryon excluded volume, but not to both.

\subsection{Extra states from baryon correlators}

Firstly,  we consider the following ratios of second order susceptibilities: $$ \chi_{11}^{BQ} / \chi_2^B, \quad \chi_{11}^{BS} / \chi_2^B.$$
By calculating them explicitly one obtains
\begin{align} \label{eq:chi11BQ}
\frac{\chi_{11}^{BQ}}{\chi_2^B} & = \frac{\sum_{j \in \rm{sectors}} \, B_j \, Q_j \, \tilde{\phi}_j (T) }{\sum_{j \in \rm{sectors}} \, B_j^2 \, \tilde{\phi}_j (T)}, \\
 \label{eq:chi11BS}
\frac{\chi_{11}^{BS}}{\chi_2^B} & = \frac{\sum_{j \in \rm{sectors}} \, B_j \, S_j \, \tilde{\phi}_j (T) }{\sum_{j \in \rm{sectors}} \, B_j^2 \, \tilde{\phi}_j (T)}.
\end{align}
The excluded volume parameter $b$ cancels out in this combination of susceptibilities.\footnote{A similar cancellation has been observed for the $\chi_{11}^{BQ}/\chi_2^B$ ratio in Ref.~\cite{Vovchenko:2017uko} in the framework of the van der Waals HRG model.} 
On the other hand, the ratios are sensitive to the particle list encoded in the ``partial pressures'' $\tilde{\phi}_j$.
Therefore, these ratios can be used to constrain the hadronic spectrum. 
In particular, it follows from Eqs.~\eqref{eq:chi11BQ} and~\eqref{eq:chi11BS} that $ \chi_{11}^{BQ} / \chi_2^B$ and $ \chi_{11}^{BS} / \chi_2^B$ probe the fractions of charged baryons and hyperons, respectively, compared to all baryons.


\begin{figure*}
    \centering
    \includegraphics[width=0.48\linewidth]{chi11BSchi2B.pdf}
    \includegraphics[width=0.48\linewidth]{chi11BQchi2B.pdf}
    \caption{Temperature dependence of second order susceptibility ratios $\chi_{11}^{BS}/\chi_2^B$ and $\chi_{11}^{BQ}/\chi_2^B$. Continuum extrapolated lattice results from Refs. \cite{Bellwied:2019pxh, Bellwied:2015lba} are shown by the black points with error bars, while the  calculations within the EV-HRG model are the curves for different hadronic lists. This combination of susceptibilities leads to the cancellation of the excluded volume parameter, $b$.}
    \label{fig:chi-BS-2B_no_EV_dep}
\end{figure*}

The ratios $\chi_{11}^{BQ}/\chi_2^B$ and $\chi_{11}^{BS}/\chi_2^B$ are shown in Fig.~\ref{fig:chi-BS-2B_no_EV_dep}.
The inclusion of extra states from PDG2016+ and QM leads to the enhancement of $-\chi_{11}^{BS}/\chi_2^B$ and the suppression of $\chi_{11}^{BQ}/\chi_2^B$.
This is driven by the fact that the extra states are mainly hyperons, thus their addition increases the fraction of strange baryons~(probed by $-\chi_{11}^{BS}/\chi_2^B$) and decreases the fraction of non-strange baryons~(probed by $\chi_{11}^{BQ}/\chi_2^B$).
The comparison with continuum extrapolated lattice data \cite{Bellwied:2019pxh, Bellwied:2015lba} suggests the need for additional hyperon states from PDG2016+/QM, as previously discussed in Ref. \cite{Alba:2017mqu}. 
The best agreement with the lattice data is obtained for the QM list.

None of the considered particle lists allow us to describe the lattice data within errors at $T \gtrsim 155$~MeV.
There are several possibilities which might explain these deviations. 
A possible explanation, which would not be captured by our model or any modifications on it, is that this temperature corresponds to the onset of deconfinement, at which new degrees of freedom (quarks) start to be liberated. Otherwise, an improvement in the agreement between lattice results and the model could be obtained through one of the following considerations.
If there are even more strange baryons than predicted by the QM, this could improve the agreement with the lattice data.
For instance, the presence of broad, high-mass Hagedorn states~\cite{Hagedorn:1965st} may have a considerable effect on the susceptibilities as one approaches the Hagedorn temperature $T_H \sim 160-180$~MeV~\cite{Lo:2015cca}. 
If the Hagedorn states contain more strange baryons than non-strange baryons, this may improve the agreement with the lattice data in Fig.~\ref{fig:chi-BS-2B_no_EV_dep}.
However, it might be challenging to preserve at the same time the agreement with the individual susceptibilities rather than in the ratios alone.

Other explanations would go beyond the physics of the model employed in the present paper.
For instance, we have modeled all resonances as free particles with zero width.
On the other hand, many non-strange baryon resonances like $\Delta$'s and $N^*$'s are broad, thus a proper treatment of their spectral functions should be important.
Modeling of broad resonances is challenging, analyses in the literature based on either pion-nucleon scattering phase shifts within the S-matrix approach~\cite{Lo:2017lym} or energy-dependent Breit-Wigner widths~\cite{Vovchenko:2018fmh} indicate that partial pressures of such resonances might be overestimated in the standard HRG model.
This implies that a more involved treatment of broad resonances may lead to a suppressed $\chi_{11}^{BQ}/\chi_2^B$ ratio~(and hence an enhanced $-\chi_{11}^{BS}/\chi_2^B$) and recover the agreement with the lattice data.

Finally, the comparison with the lattice data may be affected if there is a flavor hierarchy in baryon excluded volumes. 
While the excluded volume effects cancel out in $\chi_{11}^{BQ}/\chi_2^B$ and $\chi_{11}^{BS}/\chi_2^B$ ratios when a common EV parameter $b$ is used for all baryons, this would no longer be the case if excluded volumes differ between strange and non-strange baryons.
A smaller EV for strange baryons would lead to a smaller suppression of $\chi_{11}^{BS}$ relative to $\chi_{11}^{BQ}$, thus leading to an improved agreement with the lattice data.
It is possible that a combination of the three effects discussed here is at play, and it would be interesting to study these in more detail in the future.
We indeed find indications for the flavor-dependent excluded volumes in the behavior of fourth-order susceptibilities discussed in the following subsection.

\subsection{Fourth-order cumulants and excluded volume}

In addition to the extra states, we also want to place limits on the excluded volume parameter, $b$. 
This can be done by considering ratios of fourth-to-second order susceptibilities.
The following three ratios are all equal in the EV-HRG model under consideration and sensitive to the EV parameter $b$:
\begin{align} 
\frac{\chi^B_4}{\chi^B_2} = \frac{\chi^{BS}_{31}}{\chi^{BS}_{11}} =  \frac{\chi^{BQ}_{31}}{\chi^{BQ}_{11}} 
& = \frac{1 - 8 \, W(\varkappa_B) + 6 [W(\varkappa_B)]^2}{[1 + W(\varkappa_B)]^4} \nonumber \\
& = 1 - 12 \varkappa_B + O(\varkappa_B^2).
\end{align}    
In the ideal HRG model, i.e. without the EV interactions, $\varkappa_B = 0$ so that $\frac{\chi^B_4}{\chi^B_2} = \frac{\chi^{BS}_{31}}{\chi^{BS}_{11}} = \frac{\chi^{BQ}_{31}}{\chi^{BQ}_{11}} = 1$ regardless of the inclusion of any additional hadronic states. 
The suppression of these ratios relative to unity, on the other hand, is directly sensitive to the EV interactions and can be used to constrain the EV parameter $b$.
Furthermore, the fact that all three ratios are predicted to be equal within the model allows us to probe the limits of validity of the model, which would be signaled by the point where the equality among these three ratios no longer holds in the lattice data.
\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{4th-to-2nd-all.pdf}
    \caption{Temperature dependence of fourth-to-second order susceptibility ratios $\chi_4^B/\chi_2^B$, $\chi_{31}^{BS}/\chi_{11}^{BS}$, and $\chi_{31}^{BQ}/\chi_{11}^{BQ}$, predicted to be equal in the EV-HRG model. The lattice data at finite lattice spacing $N_{\tau}=12$ from Ref. \cite{Borsanyi:2018grb} are shown as grayscale symbols with error bars, while the calculations within the EV-HRG model are shown as bands for a range of excluded volume parameter, $b$, and for different hadronic lists. The ideal HRG result is given by the horizontal line at unity.
    }
    \label{fig:chi-4-2_equiv_EV}
\end{figure}

Figure~\ref{fig:chi-4-2_equiv_EV} depicts the results of the calculation of the ratios $\chi_4^B/\chi_2^B$, $\chi_{31}^{BS}/\chi_{11}^{BS}$ and $\chi_{31}^{BQ}/\chi_{11}^{BQ}$ within the EV-HRG model for range $b = 0.4-1$~fm$^3$ of the EV parameter values. 
The three ratios all coincide with one another, as expected, and exhibit
minimal dependence on the hadronic list utilized. 
Therefore, these particular quantities are indeed sensitive mainly to the excluded volume repulsive interactions rather than to the hadronic spectrum used in the HRG model. 
In Fig. \ref{fig:chi-4-2_equiv_EV} we compare the calculations in the EV-HRG model with various particle lists to lattice data at $N_{\tau} = 12$ from the Wuppertal-Budapest collaboration~\cite{Borsanyi:2018grb}, in the temperature range $T = 135-170$~MeV. 
Since not all of the available lattice data are continuum extrapolated, we choose this larger lattice spacing for the comparisons and avoid constructing ratios that would be a mixture of results at finite $N_{\tau}$ and in the continuum limit. 
The qualitative behavior of the three ratios is very similar in the whole temperature range considered.
Quantitatively, we see that at temperatures below 150 MeV the three ratios sit on top of each other. 
This is expected, although the lattice error bars are relatively sizable at those temperatures.

Statistically significant differences between the three susceptibility ratios in the lattice data emerge at $T \gtrsim$~160~MeV.
On the one hand, this may be a reflection of the transition to the deconfined phase where a hadronic model would be expected to break down.
In the Stefan-Boltzmann limit of massless quarks one has $\chi_4^B/\chi_2^B = \chi_{31}^{BS}/\chi_{11}^{BS} \to 2 / (3\pi^2)$ while $\chi_{31}^{BQ}/\chi_{11}^{BQ} \to 0$.
Thus, the notably smaller values of $\chi_{31}^{BQ}/\chi_{11}^{BQ}$ at $T \gtrsim 160$~MeV compared to the other two ratios might be related to the smaller Stefan-Boltzmann limit for this quantity.
On the other hand, these differences may also reflect a flavor hierarchy in baryon excluded volumes.
The ratios are predicted to be equal in the EV-HRG model if all baryons are assigned a common EV parameter $b$.
However, if for example strange baryons have a different~(smaller) excluded volume, one would expect $\chi_{31}^{BS}/\chi_{11}^{BS}$ to exhibit smaller deviations from the baseline of unity than the other ratios.
We see in Fig.~\ref{fig:chi-4-2_equiv_EV} that the separation between the strangeness and electric charge susceptibilities is such that they lie on the higher and lower ends of the $b = 0.4-1$~fm$^3$ band, respectively.
This is an indication that the strange baryons may indeed have a smaller volume than the non-strange ones.
From a phenomenological point of view, the smaller hyperon volume could reflect the fact that hyperon-hyperon interactions are mediated by the exchange of heavier mesons like $\phi$ compared to nucleon-nucleon interactions which correspond to the exchange of lighter, non-strange mesons like $\sigma$ and $\omega$.
The implementation of smaller excluded volumes for strange particles is possibly not unique. Such extensions of the EV-HRG model have been considered in Refs.~\cite{Alba:2016hwx,Vovchenko:2017zpj,Alba:2017bbr,Motornenko:2020yme} but are not considered in the present study. Of course it can be done in the future in order to model the subtle differences between $\chi_4^B/\chi_2^B$, $\chi_{31}^{BS}/\chi_{11}^{BS}$, and $\chi_{31}^{BQ}/\chi_{11}^{BQ}$, for example at the chemical freeze-out in heavy-ion collisions.

\subsection{Effect of the extra states and excluded volume on other susceptibilities}


Next, we investigate the following two combinations of susceptibilities that are sensitive to the extra strange states: the kurtosis of net-strangeness fluctuations $\chi_4^S / \chi_2^S$ and the correlator $\chi_{11}^{us}$ between net numbers of up and strange quarks.
The effect of extra states on these quantities was investigated in Ref.~\cite{Alba:2017mqu} without excluded volume effects. 
It was shown there that extra states improve the description of $\chi_{11}^{us}$ but in the case of the QM list spoil the agreement with the lattice data for $\chi_4^S/\chi_2^S$. 
Here we investigate how these quantities are affected by the presence of baryon excluded volume in addition to extra states.

Figure~\ref{fig:chi-4S-2S} depicts the temperature dependence of $\chi_4^S/\chi_2^S$.
This quantity does not involve any $\mu_B$ derivatives, thus it does not probe the repulsive baryonic interactions as directly as e.g. $\chi_4^B/\chi_2^B$ that we considered before.
Due to that fact, we calculate over a broader range of $b = 0-1$~fm$^3$ for the EV parameter value, given by the bands in Fig.~\ref{fig:chi-4S-2S}.
Since there are no $\mu_B$ derivatives, both hyperons as well as strange mesons like kaons contribute to this quantity.
Moreover, the values of $\chi_4^S/\chi_2^S$ above unity are due to multi-strange hyperons.
Figure~\ref{fig:chi-4S-2S} shows that the excluded volume suppresses $\chi_4^S/\chi_2^S$. This effect is less pronounced at smaller temperatures, but becomes sizable at $T \gtrsim 150$~MeV.
We see that the lattice data are underestimated when using the standard PDG2016 list, and the excluded volume does not improve the agreement.
For the PDG2016+ list the agreement with the lattice data is obtained at $T \lesssim 160$~MeV for the smallest considered values of the excluded volume~(the upper part of the blue band).
This is consistent with the observation made earlier that strange baryons prefer a smaller excluded volume.
On the other hand, when the QM list is considered, which contains even more extra strange baryons, the best agreement with the lattice data is obtained for the higher value of the EV parameter, $b \simeq 1$~fm$^3$~(the lower part of the red band).
We see similar effects for  $\chi^{us}_{11}$, shown in Fig.~\ref{fig:chi-us}, with a slight over-prediction in the case of the QM list.
The PDG2016 list describes the lattice data at $T \lesssim 150$~MeV, but breaks down at higher temperatures, with no benefit from introducing the excluded volume.
The PDG2016+ and QM lists allow us to extend the agreement with the lattice data for $\chi^{us}_{11}$ to $T = 160-165$~MeV when baryon excluded volumes of up to $b \simeq 0.4$~fm$^3$ and $b \simeq 1$~fm$^3$ are used, respectively.


\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{chi4S2S-all.pdf}
    \caption{Temperature dependence of fourth-to-second order strangeness susceptibility ratio $\chi_4^S/\chi_2^S$. Continuum extrapolated lattice results for this quantity are given by black points with error bars \cite{Alba:2017mqu} and EV-HRG calculations are shown for the full range of the parameter $b$ for different hadronic lists.}
    \label{fig:chi-4S-2S}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{chi11us_cont_ext_b-all.pdf}
    \caption{Up-strange quark susceptibility, $\chi^{us}_{11}$, as a function of the temperature. The continuum extrapolated lattice results are given as black points with error bars \cite{Alba:2017mqu}, while the EV-HRG calculations are shown for a range of the parameter $b$ between 0.4 and 1 $\text{fm}^3$ for different hadronic lists.}
    \label{fig:chi-us}
\end{figure}

Based on the analysis of $\chi_4^S/\chi_2^S$ and $\chi^{us}_{11}$, we observe a correlation between the number of extra strange states and baryon excluded volume: the larger the number of extra states, the larger the excluded volume must be in order to recover agreement with the lattice data for these two quantities.
Thus the attractive interactions via the inclusion of extra states can be balanced by an additional repulsion in the baryon sector in these observables.
One way to break this degeneracy is to consider quantities that probe only a single one of the effects, like the $\chi_{11}^{BQ}/\chi_2^B$ and $\chi_{11}^{BS}/\chi_2^B$ ratios that we considered earlier.
However, those ratios may be sensitive to additional physics like modeling of broad resonances, as discussed above.
Instead, we study here  a couple of additional susceptibilities that probe the strangeness content of baryons.
First we look at $\chi^{BS}_{22}$, which is a doubly strange quantity and thus more sensitive to multi-strange hyperons.
This quantity, shown in the left panel of Fig.~\ref{fig:chi-22BS}, paints a picture consistent with $\chi_4^S/\chi_2^S$ and $\chi^{us}_{11}$: extra states are clearly needed to describe the lattice data, but PDG2016+ prefers a smaller excluded volume while the QM list prefers a larger one.
On the other hand, the mixed $BQS$ susceptibility containing only one strangeness derivative $\chi^{BQS}_{211}$, shown in the right panel of Fig.~\ref{fig:chi-22BS}, exhibits only mild dependence on the excluded volume
but large sensitivity to the number of extra states. 
The minimal dependence on the excluded volume parameter $b$ of this fourth order susceptibility can be understood when one considers that the derivatives with respect to $Q$ and $S$ are only single derivatives. 
Some baryons, e.g. $\Sigma^+$, carry these conserved charges with opposite signs~(see Table~\ref{tab:hadron_table}).
For similar reasons, the contributions of the various baryon-baryon interactions to $\chi^{BQS}_{211}$ can be either positive or negative and as such, each term making up this quantity, as evaluated from Eqs. \eqref{eq:pBevW} and \eqref{eq:chis}, can contribute with different signs. 
This leads to a smaller $b$-dependence of $\chi^{BQS}_{211}$ than that of $\chi^{BS}_{22}$~(the left panel of Fig.~\ref{fig:chi-22BS}), in which case all contributions carry the same sign due to the fact that both the baryon number $B$ and the strangeness $S$ in this observable are squared.
The lattice results tend to lie in between the predictions based on the PDG2016+ and QM lists.
This indicates that the number of extra states might be underestimated in the PDG2016+ list but overestimated in the QM list.
Thus, the difference between the two lists may be taken as a systematic uncertainty in the particle list.


\begin{figure*}
    \centering
    \includegraphics[width=0.48\linewidth]{chi22BS-all.pdf}
    \includegraphics[width=0.48\linewidth]{chi211BQS-all.pdf}
    \caption{Temperature dependence of mixed fourth order susceptibilities $\chi^{BS}_{22}$ and $\chi^{BQS}_{211}$. The results at finite lattice spacing $N_{\tau}=12$ from Ref. \cite{Borsanyi:2018grb} are shown as black points with error bars. Left: Fourth-order baryon-strangeness susceptibility calculated with several different hadronic lists in the EV-HRG model for the full range of the EV parameter $b$. Right: Mixed fourth order $BQS$ susceptibility calculated with several different hadronic lists in the EV-HRG model for a range of the EV parameter $b$ between 0.4 and 1~fm$^3$.}
    \label{fig:chi-22BS}
\end{figure*}


We note that the lattice data on
the susceptibilities $\chi^{BS}_{22}$ and $\chi^{BQS}_{211}$, as well as the ratios $\chi_{31}^{BS}/\chi_{11}^{BS}$ and $\chi_{31}^{BQ}/\chi_{11}^{BQ}$, are not continuum extrapolated but only available up to a lattice spacing of $N_{\tau} = 12$. While continuum extrapolated lattice results for $\chi_4^B/\chi_2^B$ have been shown e.g. in Ref. \cite{Borsanyi:2013hza}, in this manuscript we use $N_{\tau}=12$ also for this quantity, for a consistent comparison with the other ratios.
These data could be subject to slight alterations in the continuum limit, in particular at the lower temperatures.
Additionally, there are many more susceptibilities available from the lattice than those explored in this work. 
However, in this manuscript we focused mostly on baryon and strangeness observables in order to probe the extra strange states and baryon-baryon interactions. 
This allows us to avoid the most severe lattice systematics, for instance taste violation, which mainly affect calculations sensitive to pion degrees of freedom, like e.g. electric charge susceptibilities. 


\section{Conclusions}\label{sec:concl}

We have investigated two common extensions of the Hadron Resonance Gas model 
that implement additional attractive and repulsive interactions among hadrons.
The attractive interactions correspond to adding extra
states exceeding those measured with high confidence by the Particle Data Group, leading to additive corrections to the overall pressure in the HRG model. 
This has been studied with the use of the PDG2016+ and QM particle lists. 
On the other hand, we also apply excluded-volume corrections in the baryon sector, 
which model the presence of repulsive core in (anti)baryon-(anti)baryon interactions. 
We demonstrate that these two extensions are complementary and find support in the available first-principles lattice QCD data.

To constrain the two extensions simultaneously we constructed specific combinations of conserved charge susceptibilities that probe the two effects separately.
We show that the second order ratios, $\chi_{11}^{BS}/\chi_2^B$ and $\chi_{11}^{BQ}/\chi_2^B$, probe only the hadron spectrum but not the excluded volume. 
The inclusion of additional states improves the agreement with the lattice data for these two quantities, with the best description obtained using the QM list.
Even for the QM list, however, deviations from the lattice data emerge at $T \sim 150-155$~MeV.
We argued that these deviations may necessitate a more involved modeling of broad resonances as well the possibility of smaller excluded volumes for strange baryons compared to non-strange ones.
We further studied the constraints on the hadronic spectrum by analyzing the various strangeness susceptibilities, including $\chi_4^S/\chi_2^S$, $\chi^{us}_{11}$, $\chi^{BS}_{22}$ and $\chi^{BQS}_{211}$.
The analysis of the lattice data for these quantities indicates that both the PDG2016+ and QM lists are preferable over the standard list, where the former contains most but not all the extra strange states while the QM list contains too many.
Therefore, the difference between the results using these two lists can be taken as a systematic uncertainty due to the hadron spectrum.

The fourth-to-second order ratios, $\chi_4^B/\chi_2^B$, $\chi_{31}^{BS}/\chi_{11}^{BS}$, and $\chi_{31}^{BQ}/\chi_{11}^{BQ}$ are shown to be mainly sensitive to the excluded volume corrections and, thus, suitable to constrain these corrections.
In the absence of the excluded volume corrections, these three ratios are equal to unity irrespective of the hadronic spectrum.
The excluded volume effects, on the other hand, suppress the ratios and make them behave similarly to the lattice data at $T \sim 155-165$~MeV.
The EV-HRG model that we use assigns a constant EV parameter $b$ for all baryons and predicts $\chi_4^B/\chi_2^B$, $\chi_{31}^{BS}/\chi_{11}^{BS}$, and $\chi_{31}^{BQ}/\chi_{11}^{BQ}$ to be all equal to one another at a given temperature.
The lattice data, on the other hand, reveal statistically significant differences between the three ratios at $T \gtrsim 155$~MeV, which follow a hierarchy $\chi_{31}^{BS}/\chi_{11}^{BS} > \chi_4^B/\chi_2^B > \chi_{31}^{BQ}/\chi_{11}^{BQ}$.
We argued that this hierarchy indicates a flavor dependence in the baryon excluded volumes, namely that strange baryons have generally smaller excluded volumes than non-strange baryons.
Therefore, while using a constant EV parameter in a range $b = 0.4-1$~fm$^3$ may be good enough to capture the general suppression of $\chi_4^B/\chi_2^B$, $\chi_{31}^{BS}/\chi_{11}^{BS}$, and $\chi_{31}^{BQ}/\chi_{11}^{BQ}$ ratios at $T \sim 155-165$~MeV, a more involved model is necessary to describe the subtle differences between the three.

In summary, an extended HRG model incorporating extra states via either the PDG2016+ or QM list as well as baryon excluded volume with parameter $b = 0.4-1$~fm$^3$ significantly improves the description of many lattice QCD susceptibilities at temperatures up to $T \simeq 160-165$~MeV over the standard HRG model. 
We note that the results from these two lists provide an estimate of the theoretical uncertainty on the hadronic list.
This is particularly relevant for the chemical freeze-out conditions realized in heavy-ion collisions and can be used for improved modeling of event-by-event fluctuations measured in the corresponding experiments at the LHC, RHIC, and SPS.
For instance, the EV-HRG model studied here can be directly used in the generalized Cooper-Frye particlization routine developed in Ref.~\cite{Vovchenko:2020kwg}.
We summarize the performance of each list as follows. We find that the QM list performs better for the ratio $\chi_{11}^{BS}/\chi_2^B$, while both lists describe the continuum extrapolated lattice data well for $\chi_{11}^{BQ}/\chi_2^B$ as shown in Fig. \ref{fig:chi-BS-2B_no_EV_dep}.
For ratios of fourth-to-second order susceptibilities as shown in Fig. \ref{fig:chi-4-2_equiv_EV}, we find, as we expect, a small variation in the results within the EV-HRG model for different lists.
From strangeness-sensitive susceptibilities $\chi_4^S/\chi_2^S$, $\chi_{11}^{us}$ and $\chi_{22}^{BS}$, we determined that the PDG2016+ list agrees best with the continuum extrapolated lattice data for a small excluded volume $b \simeq 0.2-0.4$~fm$^3$, while the QM list finds agreement when $b \simeq 1$~fm$^3$.
For the mixed fourth order susceptibility  $\chi_{211}^{BQS}$, the optimal list for treatment of the chemical freeze-out is PDG2016+.
Given that there is still some tension between the PDG2016+ list and the QM list, our results could indicate that more work is still needed to determine the number of states in the hadronic spectrum in heavy-ion collisions.
Further improvements of the model can be achieved by considering differences in excluded volumes of strange and non-strange baryons, as well as a more involved modeling of broad resonances.

\section*{Acknowledgements}
The authors would like to thank Paolo Parotto for discussions and help preparing the updated version of the QM list used in this publication.
This material is based upon work supported by the National Science Foundation under grant no. PHY1654219, by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under contract number DE-AC02-05CH11231 and within the framework of the Beam Energy Scan Theory (BEST) Collaboration. 
V.V. acknowledges the support via the
Feodor Lynen program of the Alexander von Humboldt
foundation. 


