Constraining the hadronic spectrum through QCD thermodynamics on the lattice

TL;DR

This study uses lattice QCD simulations with imaginary strangeness chemical potentials to investigate the hadronic spectrum, revealing the need for additional strange resonances beyond established states to match thermodynamic data.

Contribution

It introduces a novel lattice QCD approach at imaginary chemical potentials to constrain the hadronic spectrum, highlighting missing states in the HRG model.

Findings

Established states are insufficient to match lattice results.

Additional PDG or Quark Model states are needed for agreement.

Interactions beyond HRG may be required in the strange meson sector.

Abstract

Fluctuations of conserved charges allow to study the chemical composition of hadronic matter. A comparison between lattice simulations and the Hadron Resonance Gas (HRG) model suggested the existence of missing strange resonances. To clarify this issue we calculate the partial pressures of mesons and baryons with different strangeness quantum numbers using lattice simulations in the confined phase of QCD. In order to make this calculation feasible, we perform simulations at imaginary strangeness chemical potentials. We systematically study the effect of different hadronic spectra on thermodynamic observables in the HRG model and compare to lattice QCD results. We show that, for each hadronic sector, the well established states are not enough in order to have agreement with the lattice results. Additional states, either listed in the Particle Data Group booklet (PDG) but not well established, or predicted by the Quark Model (QM), are necessary in order to reproduce the lattice data. For mesons, it appears that the PDG and the quark model do not list enough strange mesons, or that, in this sector, interactions beyond those included in the HRG model are needed to reproduce the lattice QCD results.

Equations

$${\chi }_{lmn}^{BQS}={\left(\frac{{\partial }^{l+m+n}P(T,{\mu }_B,{\mu }_Q,{\mu }_S)/T^4}{\partial ({\mu }_B/T{)}^l\partial ({\mu }_Q/T{)}^m\partial ({\mu }_S/T{)}^n}\right)}_{\mu =0}.$$

$$P_{tot}(T,\mu )=\sum _k{P}_k(T,{\mu }_k)=\sum _k(-1{)}^{B_k+1}\frac{d_{k}T}{(2\pi {)}^3}\int d^3\vec{p}$$

$$\ln \left(1+(-1{)}^{B_k+1}\exp \left[-\frac{(\sqrt{{\vec{p}}^2+m_k^2}-{\mu }_k)}{T}\right]\right),$$

$${\left(\frac{{\mu }_S}{{\mu }_B}\right)}_{LO}=-\frac{{\chi }_{11}^{BS}}{{\chi }_2^S}-\frac{{\chi }_{11}^{QS}}{{\chi }_2^S}\frac{{\mu }_Q}{{\mu }_B}.$$

$$P({\hat{\mu }}_B,{\hat{\mu }}_S)$$

$$P_{01}^{BS}={\chi }_2^S-{\chi }_{22}^{BS},$$

$$\mathrm{Im}{\chi }_1^B=-P_{11}^{BS}\sin ({\mu }_I)-P_{12}^{BS}\sin (2{\mu }_I)$$

$$-P_{13}^{BS}\sin (3{\mu }_I),$$

$${\chi }_2^B=P_{10}^{BS}+P_{11}^{BS}\cos ({\mu }_I)+P_{12}^{BS}\cos (2{\mu }_I)$$

$$+P_{13}^{BS}\cos (3{\mu }_I),$$

$$\mathrm{Im}{\chi }_1^S=(P_{01}^{BS}+P_{11}^{BS})\sin ({\mu }_I)+2P_{12}^{BS}\sin (2{\mu }_I)$$

$$+3P_{13}^{BS}\sin (3{\mu }_I),$$

Full text

Introduction

The precision achieved by recent lattice simulations of QCD thermodynamics allows to extract, for the first time, quantitative predictions which provide a new insight into our understanding of strongly interacting matter. Recent examples include the precise determination of the QCD transition temperature Aoki et al. (2006, 2009); Borsanyi et al. (2010a); Bazavov et al. (2012a), the QCD equation of state at zero Borsanyi et al. (2010b, 2014a); Bazavov et al. (2014a) and small chemical potential Borsanyi et al. (2012a); Gunther et al. (2016); Bazavov et al. (2017) and fluctuations of quark flavors and/or conserved charges near the QCD transition Borsanyi et al. (2012b); Bazavov et al. (2012b); Bellwied et al. (2015a). The latter are particularly interesting because they can be related to experimental measurements of particle multiplicity cumulants, thus allowing to extract the freeze-out parameters of heavy-ion collisions from first principles Karsch (2012); Bazavov et al. (2012c); Borsanyi et al. (2013, 2014b); Noronha-Hostler et al. (2016). Furthermore, they can be used to study the chemical composition of strongly interacting matter and identify the degrees of freedom which populate the system in the vicinity of the QCD phase transition Koch et al. (2005); Bazavov et al. (2013); Bellwied et al. (2013).

The vast majority of lattice results for QCD thermodynamics can be described, in the hadronic phase, by a non-interacting gas of hadrons and resonances which includes the measured hadronic spectrum up to a certain mass cut-off. This approach is commonly known as the Hadron Resonance Gas (HRG) model Dashen et al. (1969); Venugopalan and Prakash (1992); Karsch et al. (2003a, b); Tawfik (2005). There is basically no free parameter in such a model, the only uncertainty being the number of states, which is determined by the spectrum listed in the Particle Data Book. It has been proposed recently to use the precise lattice QCD results on specific observables, and their possible discrepancy with the HRG model predictions, to infer the existence of higher mass states Noronha-Hostler et al. (2009); Majumder and Muller (2010); Bazavov et al. (2014b), not yet measured but predicted by Quark Model (QM) calculations Capstick and Isgur (1986); Ebert et al. (2009) and lattice QCD simulations Edwards et al. (2013). This leads to a better agreement between selected lattice QCD observables and the corresponding HRG curves. However, for other observables the agreement with the lattice gets worse, once the QM states are included.

Amongst experimentally measured hadronic resonances within the Particle Data Group (PDG) list, there are different confidence levels on the existence of individual resonances. The most well-established states are denoted by **** stars whereas * states indicate states with the least experimental confirmation. Furthermore, states with the fewest stars often do not have the full decay channel information known nor the branching ratios for different decay channels.

In Fig. 1 we compare, for several particle species, the states listed in the PDG$2016$ (including states with two, three and four stars) Patrignani et al. (2016), in the PDG$2016+$ (including also states with one star) Patrignani et al. (2016) and those predicted by the original Quark Model Capstick and Isgur (1986); Ebert et al. (2009) and a more recent hypercentral version (hQM) Ferraris et al. (1995). The latter contains fewer states than the ones found in Refs. Capstick and Isgur (1986); Ebert et al. (2009), due to inclusion of an interaction term between the quarks in the bound state, and the decay modes are listed for most of the predicted states. No mass cut-off has been imposed. The total number of measured particles and anti-particles, excluding the charm and the bottom sector, increases from the $2016$ to the $2016+$ listing: considering particles and antiparticles and their isospin multiplicity we get 608 states with two, three and four stars and 738 states when we also include the one star states. In the QM description the overall increase is much larger: in total there are 1446 states in the non-relativistic QM Capstick and Isgur (1986); Ebert et al. (2009) and 1237 in the hQM Ferraris et al. (1995); Bijker et al. (2016). The QM predicts such a large number of states because they arise from all possible combinations of different quark-flavor, spin and momentum configurations. However, many of these states have not been observed in experiments so far; besides the basic QM description does not provide any information on the decay properties of such particles. As already mentioned, the hQM reduces the number of states by including an interaction term between quarks in a bound state. A more drastic reduction can be achieved by assuming a diquark structure Ferraris et al. (1995); Santopinto (2005); Santopinto and Ferretti (2015) as part of the baryonic states, although experiments and lattice QCD may disfavor such a configuration Edwards et al. (2012).

In this paper, we perform an analysis of several strangeness-related observables, by comparing the lattice QCD results to those of the HRG model based on different resonance spectra: the PDG 2016 including only the more established states (labeled with two, three and four stars), the PDG 2016 including all listed states (also the ones with one star), and the PDG 2016 with the inclusion of additional Quark Model states. This is done in order to systematically test the results for different particle species, and get differential information on the missing states, based on their strangeness content. The observables which allow the most striking conclusions are the partial pressures, namely the contribution to the total pressure of QCD from the hadrons, grouped according to their baryon number and strangeness content. The main result of this paper is a lattice determination of these partial pressures. This is a difficult task, since the partial pressures involve a cancellation of positive and negative contributions (see the next section), and they span many orders of magnitude, as can be seen in Fig. 2. From this analysis a consistent picture emerges: all observables confirm the need for not yet detected or at least not yet fully established strangeness states. The full PDG2016 list provides a satisfactory description of most observables, while for others the QM states are needed in order to reproduce the lattice QCD results. Besides, it appears that the PDG and the Quark Model do not list enough strange mesons or that, in this sector, interactions beyond those included in the HRG model are needed to reproduce the lattice QCD data Alba (2017); Cabrera et al. (2014).

HRG and the strangeness sectors

The HRG model provides an accurate description of the thermodynamic properties of hadronic matter below $T_{c}$. This is especially true for global observables such as the total pressure and other collective thermodynamic quantities. However, it was recently noticed that more differential observables which are sensitive to the flavor content of the hadrons show a discrepancy between HRG model and lattice results Bazavov et al. (2014b). An example of such discrepancy is shown in Fig. 3 and will be explained below. Such observables involve the evaluation of susceptibilities of conserved charges in the system at vanishing chemical potential:

[TABLE]

Cumulants of net-strangeness fluctuations and correlations with net-baryon number and net-electric charge have been evaluated on the lattice in a system of $(2+1)$ flavours at physical quark masses and in the continuum limit Borsanyi et al. (2014c, 2013); Bellwied et al. (2015a). The same quantities can be obtained within the HRG model. In this approach, the total pressure in the thermodynamic limit for a gas of non-interacting particles in the grand-canonical ensemble is given by:

[TABLE]

where the sum runs over all the hadrons and resonances included in the model. Here the single particle chemical potential is defined with respect to the global conserved charges (baryonic $B$, electric $Q$ and strangeness $S$) as $\mu_{k}=B_{k}\mu_{B}+Q_{k}\mu_{Q}+S_{k}\mu_{S}$. More details on the HRG model used here can be found in Ref. Alba et al. (2014). In order to describe the initial conditions of the system occurring during a heavy-ion collision, we require strangeness neutrality and the proper ratio of protons to baryons given by the colliding nuclei, $n_{Q}=\tfrac{Z}{A}n_{B}\simeq 0.4n_{B}$. These conditions yield $\mu_{S}$ and $\mu_{Q}$ as functions of $\mu_{B}$; their specific dependence on $\mu_{B}$ is affected by the amount of strange particles and charged particles included in the model. To leading order in $\mu_{B}$, the ratio $\mu_{S}/\mu_{B}$ reads Borsanyi et al. (2013); Bazavov et al. (2012c):

[TABLE]

The inclusion of a larger number of heavy hyperons, such as $\Lambda$ and $\Xi$, and the constraint of strangeness neutrality are reflected by a larger value of the strange chemical potential $\mu_{S}$ as a function of temperature and baryo-chemical potential. In Fig. 3 this ratio is shown as a function of the temperature: our new, continuum extrapolated lattice results are compared to the HRG model calculations based on the 2012 version of the PDG and on the Quark Model states (as done in Ref. Bazavov et al. (2014b)). One should expect agreement between HRG model and lattice calculations up to the transition temperature which has been determined independently on the lattice to be $\sim 155$ MeV Aoki et al. (2006, 2009); Borsanyi et al. (2010a); Bazavov et al. (2012a). The HRG model based on the QM particle list yields a better agreement with the lattice data within error bars, while the HRG results based on the PDG2012 spectrum underestimate the data. However, for other observables such as $\chi_{4}^{S}/\chi_{2}^{S}$ and $\chi_{11}^{us}$ (see the two panels of Fig. 4), the agreement between HRG model and lattice results is spoiled when including the QM states. The QM result overestimates both $\chi_{4}^{S}/\chi_{2}^{S}$ and $\chi_{11}^{us}$; $\chi_{4}^{S}/\chi_{2}^{S}$ is proportional to the average strangeness squared in the system: the fact that the QM overestimates it, means that it contains either too many multi-strange states or not enough $|S|=1$ states. Moreover, the contribution to $\chi_{11}^{us}$ is positive for baryons and negative for mesons: this observable provides the additional information that the QM list contains too many (multi-)strange baryons or not enough $|S|=1$ mesons.

In this paper, we try to solve this ambiguity, even though we are aware that it might be difficult to resolve the contribution of high mass particles in our simulations. We separate the pressure of QCD as a function of the temperature into contributions coming from hadrons grouped according to their quantum numbers. This is done by assuming that, in the low temperature region we are interested in, the HRG model in the Boltzmann approximation yields a valid description of QCD thermodynamics. If this is the case, the pressure of the system can be written as Bazavov et al. (2013); Noronha-Hostler et al. (2016):

[TABLE]

where $\hat{\mu}_{i}=\mu_{i}/T$, and the quantum numbers can be understood as absolute values. These partial pressures are the main observables we study. Notice that we do not distinguish the particles according to their electric charge content.

Assuming this ansatz for the pressure, the partial pressures $P^{BS}_{ij}$ can be expressed as linear combinations of the susceptibilities $\chi^{BS}_{ij}$. An example of one such formula is:

[TABLE]

which gives the strange meson contribution to the pressure. This means that in principle one could determine these partial pressures directly from $\mu=0$ simulations, by evaluating linear combinations of the $\chi^{BS}_{ij}$ directly. This can be done on the lattice, by calculating fermion matrix traces, that can be evaluated with the help of random sources Allton et al. (2002); Bellwied et al. (2015a).

This is not the approach we pursue here, since the noise level in the calculation would be too high, certainly for the $S=2$ or $3$ sectors, but as Fig. 6(bottom) shows, probably already for $S=1$. The higher order fluctuations are already quite noisy, because they involve big cancellations between positive and negative contributions Bellwied et al. (2015a). In addition, when we take linear combinations to calculate the partial pressures, we introduce extra cancellations between the susceptibilities. Therefore we propose to use an imaginary strangeness chemical potential and extract the partial pressures from the $\operatorname{Im}\mu_{S}$ depdendence of low order susceptibilities. For earlier works exploiting imaginary chemical potentials, see de Forcrand and Philipsen (2002); D’Elia and Lombardo (2003); Wu et al. (2007); D’Elia et al. (2007); Conradi and D’Elia (2007); de Forcrand and Philipsen (2008); D’Elia and Sanfilippo (2009). A more recent work, that uses imaginary chemical potentials to estimate higher order susceptibilities is D’Elia et al. (2016).

Lattice method

Our lattice simulations use the same 4stout staggered action as  Bellwied et al. (2015a, b); Borsanyi et al. (2016); Gunther et al. (2016). We generate configurations at $\mu_{B}=\mu_{S}=\mu_{Q}=0$ as well as $\mu_{B}=\mu_{Q}=0$ and $\operatorname{Im}\mu_{S}>0$, in the temperature range $135\rm{MeV}\leq T\leq 165\rm{MeV}$. All of our lattices have an aspect ratio of $LT=4$. We run roughly $1000-2000$ configurations at each simulation point, separated by 10 HMC trajectories. For the determination of the strangeness sectors we use the HRG ansatz of equation (4) for the pressure. With the notation $\mu_{S}=i\mu_{I}$ we get by simple differentiation:

[TABLE]

and similar terms for $\chi^{S}_{2}$, $\chi^{S}_{3}$ and $\chi^{S}_{4}$. An advantage of the imaginary chemical potential approach is that we do not have to make any extra assumptions, as compared to the direct evaluation of the linear combinations, see equation (5), as the linear combinations there were already derived from the HRG ansatz. On the other hand, our method reduces the errors considerably, as the lower derivatives already contain the information on the higher strangeness sectors. Simulations at imaginary chemical potential are not hampered by the sign problem, so the evaluation of the lower order susceptibilities at $\operatorname{Im}\mu_{S}>0$ is not any harder than at $\mu=0$.

To obtain the Fourier coefficients $(P^{BS}_{01}+P^{BS}_{11})$, $P^{BS}_{12}$ and $P^{BS}_{13}$ we perform a correlated fit with the previous ansatz for the observables $\chi^{S}_{1}$, $\chi^{S}_{2}$, $\chi^{S}_{3}$ and $\chi^{S}_{4}$ at every temperature. To obtain $P^{BS}_{10}$ we fit $\chi^{B}_{2}$, to obtain $P^{BS}_{11}$ we fit $\chi^{B}_{1}$. We note that $P^{BS}_{11}$ could be deduced from $\chi^{B}_{2}$ as well, but with considerably higher statistical errors. To get $P^{BS}_{01}$ we just take the difference $(P^{BS}_{01}+P^{BS}_{11})-P^{BS}_{11}$. As an illustration that the HRG based ansatz fits our lattice data we include an example of the correlated fit for $\chi^{S}_{1}$, $\chi^{S}_{2}$, $\chi^{S}_{3}$ and $\chi^{S}_{4}$ in Fig. 5

For statistical errors, we use the jackknife method. For the continuum limit we use $N_{t}=10,12$ and $16$ lattices. To estimate the systematic errors we repeat the anaylsis in several different ways: to connect the lattice parameters to physical temperatures we use two different scale settings, based on $w_{0}$ and $f_{\pi}$. More details on the scale setting can be found in Bellwied et al. (2015a). For each choice of the scale settings we use two different spline interpolations for the temperature dependence of the $P^{BS}_{ij}$. Both of these describe the data well. For each of these four choices we do the continuum limit in four different ways, by applying tree level improvement or not, and by using a straight $a+b/N_{t}^{2}$ and a rational function $1/(a+b/N_{t}^{2})$ ansatz for the continuum limit. These 16 different results are then weighted with the AKAIKE information criterion Akaike (1974) using the histogram method Durr et al. (2008). Examples of linear continuum limit extrapolations are included in Fig. 6 (top). For the $B=0$, $S=1$ sector, which includes a large contribution from kaons, the continuum extrapolation could not be carried out using these lattices. For this case we obtain a continuum estimate, based on the assumption that only $N_{t}=10$ is not in the scaling regime, therefore using only the $N_{t}=12$ and $N_{t}=16$ lattices, and the same sources of systematic error as before, but now with uniform weights.

Finally, as a comparison we show in Fig. 6 (bottom) one of the partial pressures determined with both methods for $N_{t}=16$, using the same action. The figure shows that using imaginary chemical potential improved the accuracy drastically already in the $S=1$ sector. In the $S=2,3$ sectors, the direct method would be too noisy to plot, while the imaginary $\mu$ method allows for a quite accurate determination of the strangeness sectors.

Results and their interpretation

We evaluate the contributions to the total QCD pressure from the following sectors: strange mesons, non-strange baryons, and baryons with $|S|=1,~{}2,~{}3$. For each sector, we compare the lattice QCD results to the predictions of the HRG model using the PDG2016, PDG2016+, hQM and QM spectra.

In Figs. 7 and 8 we show our results. Fig. 7 shows the contributrion of strange mesons, while Fig. 8 shows the contribution of non-strange baryons (upper left), $|S|=1$ baryons (upper right), $|S|=2$ baryons (lower left) and $|S|=3$ baryons (lower right).

We observe that, in all cases except the non-strange baryons, the established states from the most updated version of the PDG are not sufficient to describe the lattice data. For the baryons with $|S|=2$, a considerable improvement is achieved when the one star states from PDG2016 are included. The inclusion of the hQM states pushes the agreement with the lattice results to higher temperatures, but one has to keep in mind that the crossover nature of the QCD phase transition implies the presence of quark degrees of freedom in the system above $T\simeq 155$ MeV, which naturally yields a deviation from the HRG model curves. Notice that, in the case of $|S|=1,~{}3$ baryons, it looks like even more states than PDG2016+ with hQM are needed in order to reproduce the lattice results: the agreement improves when the resonances predicted by the QM Capstick and Isgur (1986); Ebert et al. (2009) are added to the spectrum. Fig. 2 shows the relative contribution of the sectors to the total pressure. Notice that three orders of magnitude separate the $|S|$=1 meson contribution from the $|S|=3$ baryon one. The method we used for this analysis, namely simulations at imaginary $\mu_{S}$, was crucial in order to extract a signal for the multi-strange baryons.

As for strange mesons, we point out that the PDG2016 and 2016+ coincide since there is no star ranking for mesons. In this sector, it was not possible to perform a continuum extrapolation for the data, since apparently they are not in the scaling regime. However, there is clear trend in the $N_{t}=10,~{}12,~{}16$ data that makes it very natural to assume that the continuum extrapolated results will lie above the HRG curves. We also include a continuum estimate of this quantity, based on only the $N_{t}=12$ and $16$ lattices, which is clearly above the HRG curves. This might mean that, for strange mesons, the interaction between particles is not well mimicked by the HRG model in the Boltzmann approximation, or that we need even more states than the ones predicted by the QM. This was already suggested in Ref. Lo et al. (2015), based on a different analysis. In general, one should keep in mind that here we use a version of the HRG model in which particles are considered stable (no width is included). Any width effects on the partial pressures can be considered in future work. Analogously, our previous lattice QCD results did not show indications of finite volume effects for the total pressure. These effects have not been checked for the partial pressures presented here.

Our analysis shows that, for most hadronic sectors, the spectrum PDG2016, does not yield a satisfactory description of the lattice results. All sectors clearly indicate the need for more states, in some cases up to those predicted by the original Quark Model. One has to keep in mind that using the QM states in a HRG description will introduce additional difficulties in calculations used in heavy ion phenomenology, as the QM does not give us the decay properties of these new states. The HRG model is successfully used to describe the freeze-out of a heavy-ion collision, by fitting the yields of particles produced in the collision and thus extracting the freeze-out temperature and chemical potential Cleymans et al. (2006); Manninen and Becattini (2008); Andronic et al. (2011), which are known as “thermal fits”. To this purpose, one needs to know the decay modes of the resonances into the ground state particles which are reaching the detector. As of yet, the QM decay channels are unknown so predictions for their decay channels are needed first, before one can use them in thermal fits models.

In conclusion, we re-calculate the two observables which triggered our analysis, namely $(\mu_{S}/\mu_{B})_{LO}$ and $\chi_{4}^{S}/\chi_{2}^{S}$, with the updated hadronic spectra. They are shown in the two panels of Fig. 9. The upper panel shows $(\mu_{S}/\mu_{B})_{LO}$ as a function of the temperature: the lattice results are compared to the HRG model curves based on the PDG2016, PDG2016+ and PDG2016+ with the inclusion of the states predicted by the hQM. The two latter spectra yield a satisfactory description of the data up to $T\simeq 145$ MeV. In the case of $\chi_{4}^{S}/\chi_{2}^{S}$, all three spectra yield a good agreement with the lattice results. Our analysis shows that the original QM overestimates these quantities because it predicts too many $|S|=2$ baryons and not enough $|S|=1$ mesons. In the context of future experimental measurements this study gives guidance to the RHIC, LHC and the future JLab experiments on where to focus their searches for as of yet undetected hadronic resonances.

Acknowledgements

We acknowledge fruitful discussions with Elena Santopinto, Igor Strakovsky, Moskov Amaryan and Mark Manley. This project was funded by the DFG grant SFB/TR55. This material is based upon work supported by the National Science Foundation under grant no. PHY-1513864 and by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, within the framework of the Beam Energy Scan Theory (BEST) Topical Collaboration. The work of R. Bellwied is supported through DOE grant DEFG02-07ER41521. An award of computer time was provided by the INCITE program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. The authors gratefully acknowledge the Gauss Centre for Supercomputing (GCS) for providing computing time for a GCS Large-Scale Project on the GCS share of the supercomputer JUQUEEN juq at Jülich Supercomputing Centre (JSC), and at HazelHen supercomputer at HLRS, Stuttgart. The authors gratefully acknowledge the use of the Maxwell Cluster and the advanced support from the Center of Advanced Computing and Data Systems at the University of Houston.