{'fraction_non_alphanumeric': 0.03495622885406078, 'fraction_numerical': 0.013821668613736706, 'mean_word_length': 4.3556512553831155, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quark Hadron Duality at Finite Temperature

At low temperatures we expect that all QCD observables are defined in terms of hadrons. This includes the partition function as well as the Polyakov loop in all representations. We analyze the physics underlying a microscopic derivation of the hadron resonance gas.


Introduction
The promise of a phase transition at finite temperature in QCD from the confined hadronic phase to a deconfined quark-gluon phase has triggered a lot of activity both theoretically as well as experimentally. This original prediction by early lattice calculations [1][2][3] has turned into a smooth crossover after many years of accumulated experience [4]. The determination and understanding of the QCD equation of state (EoS) from the lattice [5] plays a crucial role in the current analysis of ultra-relativistic heavy ion collisions [6].
While in the limit of massless quarks the QCD Lagrangian is scale invariant, implying a vanishing trace of the energy momentum tensor, the symmetry is broken explicitly by the finite quark masses and anomalously by the necessary renormalization which introduces a energy scale and generates a trace anomaly. The most recent up-to-date results for 2+1 flavor lattices have been obtained by the Wuppertal-Budapest(WB) [7] and HotQCD [8] collaborations and after continued discrepancies there is a final consensus that maximal violation of scale invariance occurs at T s ∼ 200 MeV where the trace anomaly reaches its maximum value.
At very high temperatures, the typical momentum scales or µ ∼ 2πT are large and finite quark mass effects can be neglected. Due to asymptotic freedom the strong and running coupling constant becomes small and thus interactions can be neglected. Thus, one effectively has a gas of free and massless 4N f N c fermions (quarks and anti-quarks) and 2(N 2 c − 1) bosons (gluons), and scale invariance is restored. At the same time colour is delocalized corresponding to a deconfined phase. This allows in principle to count the number of 2N f N c quark and 2(N 2 c − 1) gluon elementary species by means of the Stefan-Boltzmann law. Current analyses show that this happens for temperatures much larger than T s . Yet, there is the firm belief that because of confinement, hadron states, composite, extended and most often unstable bound states made of quarks and gluons build a complete basis at very low temperatures, where they effectively behave as point-like, stable, non-interacting and structureless particles. This quark-hadron duality resembles to a large extent the similar duality between different degrees of freedom as the one found in deep inelastic scattering (for a review see [9]). It is remarkable that in the large N c limit some of these requirements are indeed fulfilled (for a review see e.g. [10]), except and most notably the point-like character. Of course, as the temperature is raised we expect many excited states to contribute, but also that finite width and finite size effects play a role. How this hadronization happens in detail has not really been understood so far. Lattice calculations suggest that there is a smooth transition or crossover from the purely hadronic phase to the purely quark-gluon phase [4]. Actually, it is not obvious when this hadronic picture fails in practice or what is the main mechanism behind quark and gluon liberation at the lowest possible temperature. The present contribution tries to address this problem guided by our own experience on the field.

The hadron resonance gas
In the opposite limit of very low temperatures T T s one expects an interacting gas of the lightest colour neutral particles (for N f = 3 pions and kaons) [11]. Because of spontaneous chiral symmetry breaking, the lowlying pseudoscalar particles are the lightest Goldstone bosons made of u, d, s quarks. For the typical low momenta in the heat bath p ∼ 2πT m π , m K Goldstone bosons interact weakly through derivative couplings and hence interactions are strongly suppressed. Thus scale invariance violations of the non-interacting gas are due to the finite pion and kaon masses. Therefore, in the chiral limit of massless pseudoscalars scale invariance is also exact at sufficiently low temperatures. Because of the small signal, current lattice calculations of the trace anomaly are just above the edge of this pion and kaon gas, which is expected to work for T ∼ m π /2π. When the temperature is raised, hadronic interactions among pions start playing a role and two-and more particle states contribute to the thermodynamic properties. The calculation may be organized according to the quantum virial expansion [12] where there are two kinds of contributions. The excluded volume corrections come from repulsive interactions which prevent particles to approach each other below a certain distance. The resonance contributions stem from attractive interactions which generate states living long enough to produce pressure in the system, meaning that the resonance can hit the wall of the container before it decays. A well known example is ππ scattering where one has attractive and resonating states in the isospin I = 0, 1 corresponding to the σ and ρ resonances whereas one has a repulsive core in the I = 2 exotic channel [13,14] providing a measure of the finite pion size. Once a ρ meson is created, it may interact with another pion and produce a resonance, ρπ → A 1 (which is a 3π state) and so on. For baryons, the situation is similar where N π → ∆ is a good example of a resonance contribution. This separation between attractive and repulsive contributions leaves out the residual interaction stemming from the background scattering. The Hadron Resonance Gas (HRG) corresponds to assume that all interactions among the lightest and stable particles can be described by an ideal gas of non-interacting resonances which are effectively pictured as stable, point-like and elementary particles, hence the trace anomaly is given by where the sum is over all hadronic states including spin-isospin and antiparticle degeneracies, η n = ∓1 for mesons and baryons respectively, E n (p) = p 2 + M 2 n is the energy and M n are the hadron masses. This collection of masses is usually called the hadron spectrum and most often identified with the PDG compilation which represents an established consensus among particle physicists [15]. The states classification echoes the non-relativistic quark model, namely mesons are qq states and baryons qqq states. Those falling outside this category are classified as "further states". The hadron spectrum obtained from the PDG is depicted in Fig. 1 as well as its separation into mesonic and baryonic spectra.  Trace anomaly for the HRG using the PDG [15] (dashed) and the RQM [16,17] (full) compared to the continuum extrapolated results of the WB [7] and HotQCD [8] collaborations. We also show the contributions from states with M ≤ 600 MeV (dotted) and with M ≤ 800 MeV (dot-dashed).
The result of using Eq. (1) with the PDG states [15] compared to the continuum extrapolated results of the Wuppertal-Budapest (WB) [7] and HotQCD [8] collaborations is shown in Fig. 2. It is amazing that this exceedingly simple picture works accurately almost below T s (the maximal scaleviolating temperature) at about T 170MeV. Note that the lowest lattice data points at T = 120 MeV are first saturated when states with masses above the vector mesons are included. Higher temperatures start feeling excited hadronic states, which by their nature embody relativistic dynamics. The quantum statistical bosonic and baryonic character of the states accounts by less than 1% correction to the classical Maxwell-Boltzmann distribution (corresponding to take η n → 0) at about T = 200MeV, already well beyond the range where lattice and HRG agree and much smaller than the lattice 10% uncertainties. The transition from hadrons to quarks has been analyzed on the light of strangeness [18] and an observable which vanishes for the HRG has been proposed.
Actually, the HRG model has arbitrated the discrepancies between the different collaborations in the past and the lattice community have had a long struggle until agreement between them and with the HRG model has been declared. While this has made the HRG model a sort of holy grail (see e.g. [19] and references therein), it is good to remind that despite of the phenomenological success it is not a theorem, nor a well defined approximation from the original QCD Lagrangian. The most compelling argument is that it is to date unclear how corrections to this simple approach should be implemented nor what the error estimate for the HRG should be. For a comparison, we also show the result of Eq. (1) using instead the Relativized Quark Model (RQM) [16,17] which essentially combines two basic elements, the static energy among the constituents and a relativistic form of the kinetic energy which does not consider the spin of the particles but does not contain more parameters than QCD itself. The likewise impressive agreement of the RQM trace anomaly with the lattice data not only illustrates our point on the lack of uniqueness of what is actually being checked by these comparisons, but also that the RQM may provide information not listed in the PDG booklet.
For instance, the PDG lists the quantum numbers, decay modes, masses and widths of the resonances building the hadron spectrum, but no information on their size is provided. While for unstable particles this is a problematic issue (see e.g. [20] and references therein), within the heavyion literature the assumption of a constant volume is frequently made (see e.g. [21]). This information can be accessed by means of quark models or lattice calculations. On the other hand, as discussed in Ref. [22], the purely resonance character makes the very definition of the mass ambiguous, and this allows to generate an error band on the PDG prediction of the HRG model which gives a spread for the trace anomaly about half the lattice uncertainties.
3. QCD at finite temperature As a guideline, let us provide some of the main features of QCD at finite temperature emphasizing some relevant aspects for the discussion. Many gaps in this sketchy presentation may be filled by consulting e.g. [5] and references therein. The QCD Lagrangian (in Euclidean spacetime) is given by The QCD Lagrangian is invariant under colour gauge transformations The QCD thermodynamics is obtained from the partition function with the periodic or anti-periodic boundary conditions for gluons and quarks respectively where the Matsubara frequencies are w n = 2nπT for gluons and w n = (2n + 1)πT for quarks. Preservation of the quark antiperiodic boundary conditions implies that only periodic gauge transformations are allowed, namely Within the convenient Polyakov gauge (A 0 (x) stationary and everywhere diagonal [23]) the most general remaining transformation is either stationary and diagonal or of the type where λ = diag(n 1 , · · · , n Nc ) and Tr λ = 0. Large Gauge Invariance implies periodicity in A 0 with period 2πT /g This is the finite temperature version of the Gribov copies, i.e., the fact that there is no complete gauge fixing in a non-abelian gauge theory. A drastic consequence of this periodicity property is that it becomes explicitly Broken in perturbation theory [24,25]. Thus, we may consider this as a signal of the relevance of non-perturbative finite temperature gluons.
In the limit of massless quarks (m f = 0) the QCD Lagrangian is scale invariant, i.e. x −→ λx. This symmetry is broken by quantum corrections due to the necessary regularization. To see this consider the partition function dependence on the coupling constant g after the rescaling of the gluon fieldĀ µ = gA µ (and ignoring renormalization issues) Note that the only dependence on g is the one shown explicitly. Thus, where in the last equation we have assumed a vacuum space time independent configuration, with V the volume of the system. On the other hand the free energy and internal energy are given by the thermodynamic relations and the trace anomaly becomes Generally, a renormalization scale µ has to be introduced to handle both IR and UV divergences. Thus, on purely dimensional grounds one has Physical results should not depend on the renormalization scale, thus using that dP/dµ = 0 we get Introducing the beta function the trace anomaly becomes then where Here we have implemented, in full harmony with standard lattice practice (see, e.g., [26]), a subtraction to renormalize the vacuum contribution [27,28]. The vanishing of the trace anomaly at zero temperature is consistent with quark-hadron duality, see Eq. (1). Also in the massless quark limit, the QCD Lagrangian is invariant under Chiral symmetry is spontaneously broken by the chiral condensate in the vacuum down to SU V (N f ) The Gell-Mann-Oakes-Renner relation, which in the simplest N f = 2 case becomes relates the quark condensate and the quark mass to physically measurable quantities such as the pion mass m π and the pion weak decay constant f π . This relation exemplifies the quark-hadron duality, namely the fact that, in the confined phase of QCD, we expect quark observables to be representable by hadronic observables.
In the opposite limit that all quarks become infinitely massive, m f → ∞, the quark determinant becomes a constant which can be factored out of the path integral and the resulting action corresponds to a pure Yang-Mills theory This purely gluonic action exhibits a larger symmetry: 't Hooft Center Symmetry Z(N c ), namely, invariance under gauge transformations which are periodic modulo a center element of SU(N c ) An example for z = e i2πk/Nc and in the Polyakov gauge is given by the choice with λ = diag(1 − N c , 1, . . . , 1). The Polyakov loop is an order parameter of the center symmetry which is related to the free energy of a colour charge in the medium. In the Polyakov gauge, the Polyakov loop in the fundamental representation reads The vacuum expectation value of the Polyakov loop transforms under the previous gauge transformation as Therefore unbroken symmetry (L T = zL T )) implies L T = 0 and hence F q = ∞. The divergence of the free energy of a colour charge in the fundamental representation is interpreted as a signal for confinement [29]. The renormalization of the Polyakov loop is a subtle issue addressed in the lattice in Ref. [30].
In gluodynamics the center symmetry is spontaneously broken above a critical temperature and L T /N c approaches unity (or any other central element) as the temperature increases. In fact, at high temperatures A 0 T and one may expand [31], Note that while this formula suggests that L T ≤ N c , the renormalization overshoots this value at a perturbative level in a tiny but visible way [31][32][33].
In full QCD the center symmetry is explicitly broken, which results in L T > 0 at all temperatures. In the limit of heavy quarks and low temperature one has L T = O(e −mq/T ) 1 or F q → m q + + . . . , where is the smallest residual binding energy of doubly heavy QQ meson. Likewise for light dynamical quarks L T ∼ e −∆/T where ∆ denotes the ground state mass of a heavy-light Qq meson (the mass of the heavy quark excluded) [34].
One can also define Polyakov loops in higher representations which have been subject of attention only a few times despite its very interesting properties [35][36][37][38].
Scale invariance is also broken in gluodynamics. For instance in perturbation theory to two loops one has [39], where at lowest order Thus, taking µ = 2πT , we expect to have a free gas of gluons and quarks in the high temperature limit. In the simple case of non-interacting particles the partition function is given by where η = −1 for bosons, η = +1 for fermions and g i is the number of species. From the partition function we have the thermodynamic identities For the massless quark and gluon gas the pressure is given by which is the Stefan-Boltzmann law. Because in the high temperature limit the particles are effectively massless, scale invariance is restored and hence the (reduced) trace anomaly vanishes These expectations have been checked in Ref. [40] by a lattice study in a wide temperature window, T = 0.7 . . . 10 3 T c .
Thus, the quark condensate qq and the Polyakov loop in the fundamental representation L T become true order parameters in quite opposite situations. While qq signals spontaneous chiral symmetry breaking in the massless quark limit, L T signals confinement for infinitely heavy quarks. The real situation is somewhat intertwined, and can be summarized as follows.
• Order parameter of chiral symmetry breaking (m q = 0) • Order parameter of deconfinement (m q = ∞) Polyakov loop: Center symmetry Z(N c ) broken • In the real world m q is finite but inflection points nearly coincide (accidental?) The chiral-deconfinement crossover is a unique prediction of lattice QCD. Whether or not this result is accidental, could be answered by computing the (connected) crossed correlator, qq tr c e igA 0 /T − qq tr c e igA 0 /T = ∂L T ∂m q , which corresponds to the quark mass dependence of the Polyakov loop. Finally, the correlation function between Polyakov loops in an arbitrary representation R exhibits Casimir scaling (quenched approximation) [41] Tr R Ω( It can be shown that this correlation function, which for the fundamental representation is related to the singletqq free energy, F 1 (r, T ), can also be written as a ratio of partition functions betweenQQ sources placed at a distance | x 1 − x 2 | and the vacuum [42], hence satisfying a spectral decomposition with integral weights w n and positive energies E n (| where the singlet free energy F 1 (r, T ) has been introduced. One important property is that at large distances (for the unquenched full QCD case), 4. Relativized Quark-Gluon models

Relativity and thermodynamics
One of the most troublesome aspects of hadron binding is that it makes relatively heavy particles from massless ones, hence most of the mass comes from the interaction. A prominent example is the Glueball in gluodynamics, where the lightest 0 ++ state [43] has a mass M 0 ++ / √ σ ∼ 4.5 while gluons are massless. Fully relativistic few body equations are not only hard to handle but encounter many difficulties regarding cluster decomposition properties [44][45][46]. This feature is particularly interesting as it is related to the compositeness nature of relativistic particles which build the hadrons, and strictly speaking we have to deal with relativistic statistical mechanics of interacting particles, a subject which has a long history [47].
Unfortunately, as we have shown, the physics of finite temperature QCD below the phase transition involves the excited hadronic spectrum, and thus relativity becomes an essential ingredient in the game. Relativized quark models (RQM) combine two basic elements, the static energy among the constituents and a relativistic form of the kinetic energy which does not consider the spin of the particles [16,17].

The linear potential
In the Born-Oppenheimer approximation the object to be analyzed is the interaction between heavy sources A and B. In perturbation theory one has one gluon exchange which yields a Coulomb like interaction, where α s = g 2 /(4π) and λ A and λ B are the generators 1 of the SU (3) colour group corresponding to the representation of the source. This form of the colour interaction exhibits Casimir scaling, a property that is violated only at three loops in perturbation theory [48] and appears to hold non-perturbatively on the lattice with an additional linear potential contribution [49]. Thus, to a good approximation, the interaction between heavy sources on the lattice reads, 1 These are more customarily denoted by T , while λ is twice the generator.
Thus, for quark-antiquark or gluon-gluon pairs coupled to a singlet state or a quark-quark pair coupled to the antifundamental representation (diquark) the following relations are obtained In what follows we use σ to denote the string tension σ F . As a consequence of Casimir scaling the ratio between the fundamental QQ ≡ (3×3) 1 , adjoint By making simplifying assumptions, easy relations can be derived from Casimir scaling. For instance, consider the lowest glueball state by analyzing two massless spin-1 particles in the CM system assuming spin independent interactions, and similarly for the ρ meson as composed of two massless quarks. Neglecting the Coulomb term, the respective mass operators appearing in the Salpeter equation read Simple dimensional considerations imply that the eigenmasses are proportional to the square root of the string tension, thus Here, as it is customary, we have matched the scales of gluodynamics and QCD by assuming a common value of σ F in both theories. A rough estimate of the mass follows from using the uncertainty principle for the ground state, namely, by taking pr ∼ 1. For the glueball this yields

The cumulative number
The spectrum of the RQM model of Isgur and collaborators forqq in the case of mesons and qqq for baryons [16,17] is concisely shown in Fig. 3. A detailed comparison to individual states unveils a rather good description of the data. Of course, we do not expect this or any quark model Hamiltonian to describe accurately the individual levels. This, however, poses an interesting problem on how two different spectra including many excited states can quantitatively be compared, beyond eyesight and subjective impression just based on contemplating Figs. 1 and 3. One way suggested by Hagedorn in the early 60's is by means of the cumulative number of states The question is then to decide to what extent N QCD (M ) or N PDG (M ) coincide. For our discussion we will consider a simplified version of that model where hyperfine splittings due to spin and three-body forces are ignored and the Hamiltonian for n constituents (we restrict to qq and qqq systems) is taken to be We will consider explicitly some cases of interest below, but already at this level some important observations can be made on the growth of the cumulative number of states. To this end, let us adopt a semiclassical approximation, which should be reliable when the number of states is large. The number of states in the CM system and at rest, below a certain mass 1.00 0. 50 2.00 0. 30 3.00 1.

M takes the form
where g n takes into account the degeneracy. For the sake of the argument, let us neglect the Coulomb term, thus v(r) = κ r, as well as the current quark masses. In this case, a dimensional argument, p → M p, r → M r/κ, gives It is not hard to show that lifting any of the above approximations only modifies this result by subleading powers of M . Thus, for a finite number of degrees of freedom, the leading contribution to the cumulative number scales as a power of the mass. The qualitative power behaviour can be clearly identified as straight lines in the log-log plot of Fig. 4 for the cumulative number, where we compare the resulting cumulative number both for the PDG and the RQM separating the mesonic and baryonic contributions.
As we can also see, the PDG and RQM spectra look very much alike below 1.7 − 1.8 GeV. We note, however, that the spectrum for the RQM saturates sharply at 2.5 GeV which is about the cut-off mass where Isgur, Godfrey and Capstick stopped to compute states. On the other hand, the PDG states saturates at lower energy values in a softer fashion. Note that the RQM looks like a linear extrapolation (mind the log scale) of the PDG spectrum. The agreement at lower masses is not highly surprising since RQM parameters were tuned to reproduce low lying states on the one side and the listed PDG states fit into the quark model scheme on the other side. On the other hand, there are only σ and the u, d, s quark masses as fitting parameters, so we should regard the agreement as a further big success of the RQM from a global point of view.

Bound states vs resonances
One issue which is problematic from the start is that if N (M ) really counts the number of bound states, then we expect it to be, for N f = 2 flavors, just pions (π + , π − , π 0 ), nucleons (p, n) and anti-nucleons (p,n), so The cumulative number is by its own stair-case nature a piecewise discontinuous function, but as we go to higher states the discontinuities smooth out. This becomes particularly visible in the RQM in the range 2 GeV < M < 2.5 GeV.
In addition, while in the PDG we have an issue regarding completeness of states, in the RQM case this is not the case; besides angular and spin flavor quantum numbers, the radial number can just be checked with the oscillation theorem, in theqq case, or simply by diagonalizing in a complete basis of normalizable states; no state will be left out in the process. In the PDG however, it is unclear if there are missing or redundant states although the consensus of listing states fitting into the quark model makes this argument into a circular one.
Of course, the meaning of M PDG n and M RQM n is different. While in the PDG listing we usually encounter a Breit-Wigner resonance parameterization characterized by a mass and a width, in the RQM we have just the mass of a bound qq or qqq state. We expect that when we couple these bound states obtained from the RQM to the continuum there will be, besides a width, a mass shift M RQM n → M RQM n +∆M n , whence the cumulative number for bound states will generally differ as the one for Breit-Wigner resonances, but the sign of ∆M n is a priori unknown.

Hadron sizes
The solution of the multiparton Hamiltonian Eq. (42) onto colour singlet states yields the corresponding hadron wave functions. To estimate the hadron size we can make use of the virial theorem based on stationarity of eigenstates under unitary coordinate scaling r i → λ r i around the λ = 1 value, Ψ( x 1 , . . . , x n ) → λ 3n/2 Ψ(λ x 1 , . . . , λ x n ). For light quarks we may, for simplicity, take the massless quark limit m q → 0. Thus, for the various terms in the Hamiltonian p i → p i /λ, likewise 1/r ij → 1/r ij /λ and r ij → λ r ij where r ij is the mean distance between particles i and j. Due to the virial theorem we have for a multiparton state In particular, for mesons and baryons respectively. So, in this model the size of a hadron grows linearly with the mass per constituent. For instance, from the relations m ρ = 2σ rq q ρ and M N = 3σ r qq N we have rq q ρ = 0.42 fm and r qq N = 0.52 fm. From M N /m ρ = 3/2 we would get rq q ρ = r qq N so we can identify the constituent quark mass as M q = σ r . Note that the above virial relation includes the Coulomb like contribution −α s /r, since this term scales exactly as the kinetic piece. In the case of hadrons with one heavy quark, the virial theorem yields Mq Q = 2σ r qQ + m Q and M qqQ = 2σ r qQ + σ r qq + m Q .

String breaking
Confinement is often attributed to this ever-linear growing of the energy with the distance. This is true on the lattice only in the quenched approximation, where quark-antiquark creation is suppressed. In full QCD however, the string breaks, a fact that has been observed by lattice calculations at a distance r c = 1.25 fm [50]. This happens when a lightqq pair is created in between the heavy quark and anti-quark sources QQ, thus two colour singletqQ andQq mesonic states can be created. On the other hand, charge conjugation implies that the binding energy of theqQ and theQq is the same and equals the residual energy of a heavy-light meson with total mass Mq Q . Thus the string breaking distance corresponds to Good approximations to these states exist in nature for q = u, d and Q = c, b. From heavy-quark QCD we expect a universal (independent of the heavy quark spin and flavor) spectrum ∆ h,α of hybrid hadron masses (heavyquark mass subtracted). For the lightest, pseudoscalar, hybrid meson, the following sequence should approach a value of ∆, for increasingly heavier quarks and using the PDG values in the MS-scheme, which gives the estimate for r c = 2 fm. Another estimate can be made based on a constituent quark model picture where the total mass of the quark is M q = M 0 + m q , with M 0 the constituent quark mass and m q the current quark mass. Spontaneous breaking of chiral symmetry implies that in the chiral limit (massless current quark masses, m q → 0) the total mass is non vanishing, thus M 0 = 0. Actually, for light u, d mesons current masses can be neglected. Then the lightqq meson has a mass Mq q = 2M q and the mass of the light qqq baryon is M qqq = 3M q . The mass of a heavy-light meson would then be Mq Q = M q + M Q = 2M 0 + m q + m Q hence σr c = 4M 0 + 2m q , which for √ σ = 420−440 MeV and M 0 = 300−350 MeV yields the estimate r c = 1.2 − 1.5 fm for the string breaking distance, a quite reasonable value.

Avoided crossings
The observation of string breaking for aQQ system requires taking into account the mesonicM M channels into which the system may decay after qq pair creation from the vacuum. This coupled channel dynamics spans the Hilbert space H = HQ Q + HQ qqQ and features the avoided crossing phenomenon familiar from molecular physics in the Born-Oppenheimer approximation [51]. For two channels, sayQQ andM M =QqqQ, one can compute the direct correlators yielding the lowest energies These two channels are orthogonal for all r. For simplicity we have disregarded the Coulomb piece −4α s /3r as well as the residual interaction between the two heavy-light mesons M =qQ andM =Qq which is of van der Waals type and corresponds to meson exchange. Note that these curves cross when σr c = 2∆. Thus, the spectrum in the Hilbert space withQQ andM M =QqqQ components reads, where now the statesQQ andQqqQ are piecewise orthogonal. The point r = r c corresponding to degenerate states is singular. A linear combination ofQQ andQqqQ involves a crossed correlator between the channels and representing a variational improvement. The avoided crossing occurs because the finite energy of the non-diagonalQQ →M M interaction lifts the degeneracy, a feature called level repulsion. The adiabatic potential curves E 0 (r) and E * 0 (r) appear as avoided crossings on the lattice [50] with a finite and small energy repulsion and a narrow transition region of about 0.1fm, see Fig. 5, which resemble the simple shape of Eq. (51). Therefore, as long as the size of the system remains small, we may ignore the string breaking effect. Otherwise, one has to consider a coupled channel dynamics withQQ andQqqQ states. Excited states potential curves should follow a similar pattern as Eq. (51) but with suitable modifications. Before mixing one has the crossing among the energy levels up toqq pair creation, where the double excitation character of the adiabatic potential curve is displayed explicitly and a universal string tension is assumed. The crossings must happen at σr qQ . Avoided crossings take place when mixing among different sectors is allowed yielding energy curves sketched in Fig. 5 when using the spectrum from the RQM for c-hadrons [16,17].   [50]. Right panel: Complete double heavy-light spectrum of the RQM for c-hadrons [16,17] and included a tiny offset to make clear the structure. We also draw the envelope σr with σ = (420MeV) 2 .

Limitations in Counting states
As we have said, when the size of the system is large enough the string breaks and the assumption of a linear potential becomes invalid, see Eq. (51). Because of the linear dependence of the size with the mass, Eq. (46), this provides a maximum mass value beyond which the RQM becomes inapplicable. A rough estimate for mesons can be made by taking max Mq q = 2σr c which for r c = 1.25 fm gives max Mq q ∼ 2.2 GeV. For baryons the situation is more complex. The value max M qqq = 3σr c ∼ 3.4 GeV is actually an upper bound for an equilateral triangular configuration, however, the string may break more economically when just two constituents are sufficiently far apart, r ij ∼ r c . This corresponds to an elongated isosceles triangle (quarkdiquark configuration) such that max M qqq = 2σr c + σ r 12 ∼ 2σr c + M q ∼ 2.5 GeV. The RQM departs from the PDG at Mq q ∼ 2 GeV (see Fig. 4). While this poses the problem on the validity of the RQM for masses beyond the PDG saturation, it also suggests that higher mass states break up into weakly bound molecular systems with a small net contribution to the cumulative number. Actually, as argued in [52], counting hadronic states implicitly averages over some scale, and so states such as the deuteron generate fluctuations in a smaller scale. 2 For systems with a heavy quark, we have that for mesons ∆ Qq = M Qq − m Q = 2σ r Qq and baryons ∆ Qqq = M Qqq − m Q = 2σ r Qq + σ r qq string breaking occurs when ∆ Qq , ∆ Qqq ∼ 2σr c .

The total cumulative number and equation of state in the confined phase
The relativized quark model describes all states as bound states ofqq for mesons and qqq for baryons [16,17]. The total cumulative number is then defined as This counts the number of bound states below M which is depicted in Fig. 8 (note the log-scale) where a clear straight line is observed. By quark-hadron duality, in the limit of very low temperatures we expect to have a gas of pions (the lightest hadrons) which due to spontaneous breaking of chiral symmetry interact weakly at low energies through derivative couplings. In the chiral limit the pions would become massless resulting in a small trace anomaly in the temperature regime where heavier hadrons are suppressed. For a gas of hadrons the pressure reads where the sum is over all hadronic states including spin-isospin and antiparticle degeneracies. From here, and using the cumulative number Eq. (53) obtained with the RQM [16,17], it is straightforward to compute the trace anomaly. The comparison, shown in Fig. 2, with the continuum extrapolated results of the WB [7] and HotQCD [8] collaborations is remarkable. As a side remark let us mention that in QCD 80% of the trace anomaly stems from the gluonic part of the operator (right-hand of Eq. (13)) [8].

Polyakov loop correlators
A straightforward consequence of Eq. (52) is that in the confined phase the correlator between Polyakov loops in the fundamental representation at large distances becomes, according to Eq. (32), with w 0 = 1 and E 0 (r) = σr, where ∆ n = ∆ where we have replaced σr → VQ Q (r) = σr−π/(12r). The result is depicted in Fig. 6 for T = 50, 100, 150, 200, 250, 300 MeV using the RQM in the case of c-quarks for F 1 (∞, T ). We see that at small temperatures there is little change in qualitative agreement with lattice calculations [53]. Corrections to Eq. (55) are expected, as it is based on a sharp string breaking transition and the fact that we truncated the spectrum to one lightqq pair creation. The avoided crossing structure shown in Fig. 5 is modified by the finite string breaking transition region, which on the lattice and for the ground state was found to be about ∆r = 0.1fm. Moreover, for ∆q Q = σr c we expect theqQ system to break up intoqQ andqq.

The Polyakov loop
We can likewise consider the spectrum of a system with one heavy quark such as c, b, t quarks. An equivalent cumulative number can also be defined with similar features. Because of the heavy mass it is more convenient to subtract the heavy quark mass, m Q , from the hadron mass, M = ∆ + m Q . Using again the RQM for hadrons with one heavy quark,qQ for mesons and Qqq for baryons [16,17], the total cumulative number is defined as The result is depicted in Fig. 7, where similar patterns as the light quark systems are encountered, namely power behaviour for large ∆ for individual Nq Q (∆) and N Qqq (∆) contributions, and Hagedorn type spectrum for the combined result, with T H ∼ 210MeV [54]. The Polyakov loop is obtained by computing the corresponding partition function of Nq Q (∆) and N Qqq (∆) with the pertinent statistics since from Eq. (33) and Eq. (55) we get The result can be seen in Fig. 7 and compared to lattice data from [55] for the HISQ/tree action and [56] for the continuum extrapolated stout result.  [16,17] and from the PDG [15]. Right: Polyakov loop as a function of temperature (in MeV). Lattice data from [55] for the HISQ/tree action and [56] for the continuum extrapolated stout result. We compare lowestlying charmed hadrons from PDG [15], the RQM spectrum with one b or one c-quark [34].

Hagedorn spectrum
As we have mentioned, the HRG has been successfully applied in many situations below the phase transition, such as the EoS, and quark number susceptibilities, but little has been achieved with regard to understanding the emergence of hadronization in the low temperature regime. So, what is the complete hadron spectrum? Of course, we expect QCD to give the answer to this question, but this requires a knowledge of all multiparticle states of stable particles, and most of them are in the continuum. The PDG tries to answer this question by filling in the expected quark model statesqq for mesons and qqq for baryons. While states falling outside this category are listed as further states, whether or not the list contains redundant states is difficult to say.
In Fig. 8 we show the total cumulative number for the PDG [15] compared to the RQM [16,17]. As can be seen there is an exponential growth for M ≤ 2 GeV for the PDG states [15] and M ≤ 2.5 GeV for the RQM states [16,17]. This remarkable feature of the cumulative number, first noted by Hagedorn [57], has fascinated theoreticians for decades. The exponential growth is of the form N (M ) ∼ e M/T H where T H is the so-called Hagedorn temperature. This pattern was predicted with very few states and an ever increasing exponential spectrum was anticipated as new states entered the hadronic list. Different upgrading analyses have confirmed this exponential growth [58,59] and even two different Hagedorn temperatures have been reported (see however the discussion around Fig. 4). The onset of the Hagedorn spectrum has been questioned more recently [60]. The growth happens at about M ∼ 1.5 GeV and continues with the same slope until M ∼ 2.1 GeV. The pattern is even more clear in the RQM, where the trend stretches up to M ∼ 2.5 GeV. The consequence of a truly Hagedorn temperature is that the partition function, and hence all thermodynamics quantities develop a pole at the Hagedorn temperature, This form is not observed on the lattice although good fits in the range 50 MeV < T < 180 MeV for the trace anomaly yield a value of about T H ∼ 220 MeV. One important consequence of such a favorable comparison between PDG and RQM is that it may provide a clue to the onset of the Hagedorn spectrum. We can estimate the asymptotic behaviour for bothqq and qqq spectra in the RQM as the model Hamiltonian is known. As we have shown we expect a power behaviour for both Nq q (M ) and N qqq (M ), a fact confirmed in Fig. 4 where the separated contributions in a log-log scale exhibit this power like behaviour. The question is on how is it possible that summing two polynomials do we end up with an exponential?
Using the MIT Bag model [61] this question was answered by Kapusta [62] since then the cumulative number can be computed as n N n (M )/n!g n from Eq. (43) by using a constant Bag volume and evaluating the phase-space integral for free particles. We have checked this explicitly by summing the Bag modes [37].

Finite width corrections
Of course, the discrete summation involved in the cumulative number definition requires considering all masses of states to be interpreted as bound states. In reality most of the states listed in the PDG are resonances, i.e. unstable particles which have instead a mass spectrum characterized by a distribution ρ(µ) which is peaked at µ = M n with a certain width Γ n . This idea has been implemented by the half-width rule, where resonance masses are regarded as random variables with an uncertainty of half the width [63]. This interpretation has been fruitful in describing the Regge spectrum of mesons [64,65], hadronic form factors [66] as well as providing an error estimate for the cumulative number itself [22]. This also provides a quantitative way of defining a figure of merit for the quark model taking the width as a genuine uncertainty.
In a quantum mechanical picture, where the resonance decay can be viewed as a tunneling process, the mass shift of a unstable state is negative as the infinite barrier becomes finite making the energy shift negative. As a consequence, the cumulative number would increase for a fixed mass, N resonance (M ) > N bound (M ). This effect goes in the opposite direction of making N RQM (M ) and N PDG (M ) to agree in the upper part of the spectrum, and that the outnumbering ofqq and qqq states in the RQM would be a genuine one.

Excluded volume condition
The reason for the failure of the HRG in describing the trace anomaly at T ∼ 170 MeV may be sought by questioning any of the assumptions involved. One of them, the finite size of hadrons may be tested by computing the excluded volume.
In an ideal world where the hadrons live forever they still have a finite size. This is the case for instance in the large N c limit, where one has Γ/M = O(N −1 c ) but their size is r = O(N 0 c ). Due to the finite hadronic size, when hadrons overlap their underlying composite nature becomes relevant, particularly regarding the Pauli principle as applied to the constituent quarks.
The finite volume corrections have often been addressed along the lines of statistical mechanics for real gasses where hadron volumes are usually assumed to be similar (see e.g. [21]). Taking into account how these corrections originate in the quantum virial expansion as repulsive contributions through negative phase shifts it is unclear what are the actual values one should take for the volume without actually carrying the phase shift analysis. In a bound states picture, which is the spirit of the HRG since finite width effects are disregarded, the size of the bound state is expected to increase with the mass of the state.
One simple estimate of the size of the hadron can be made by using the MIT bag model [61] where the volume V i of the hadron is a natural concept since in that model hadrons with a mass M i have a sharp edge, r i , where and thus the volume grows with the mass.
In the case of the RQM where particles are interacting through a confining σr ij potential, and hadrons are localized but have a diffuse edge, the volume is not a well defined concept. To estimate the value of the volume occupied by the hadron in the RQM, we take as the meson radius rq q /2 and for the baryon the radius of the equilateral triangle r qq / √ 3. This yields the volume estimate where c i = 1/4 for mesons and c i = 1/(3 √ 3) for baryons. With this prescription, baryons occupy a smaller volume for the same mass. These are crude estimates, which assume a sharp edge of the hadron, but they show that the volume does depend on the mass.
Rather than trying to model finite volume corrections, we will analyze the quite natural condition that the excluded volume cannot be negative, or equivalently that the occupied volume is smaller than the total volume. This means In Fig. 9 we illustrate the situation by depicting the occupied volume condition using the MIT Bag model volume, Eq. (60), and our estimate for the RQM hadronic volume, Eq. (61), using both the RQM spectrum and the PDG spectrum. As one can see the condition is violated for temperatures T ∼ 150 − 180 MeV where the HRG departs from the lattice data. This is an interesting fact, since it relates the crossover to a volume effect, which deserves further investigation. For comparison we also show in the right panel the results for a constant hadronic volume taking the radius values r = 0.5, 0.75, 1 fm. The choice r ∼ 0.75 fm resembles the RQM and bag excluded volume condition.

Insights from gluodynamics
Matters get simpler in the case of pure gauge theories where gluons are the only degrees of freedom. In QCD this corresponds to taking the limit m q → ∞ for all quark species. As mentioned, the Polyakov loop in the fundamental representation becomes a true order parameter of the deconfinement phase transition since at low temperatures it vanishes as L T = O(e −mq/T ) in this limit. The Polyakov loop in the adjoint representation does not vanish below T c . Gluon-Glueball duality means in this case that any observable defined as an expectation value of a gauge invariant and hence colour singlet operator can be determined in terms of purely colour singlet states in the confined phase. As we will illustrate, this means in practice that the EoS can be determined in terms of glueballs (bound states of gluons) whereas the Polyakov loop in the adjoint representation can be be determined from gluelumps which are bound states of gluons in the presence of a colour octet source [67][68][69][70].

The glueball gas model
In pure gauge theories the glueball spectrum has been determined in the lattice [71]. (For full unquenched QCD glueball studies see e.g. [72]. A general perspective including phenomenology has been reviewed in Ref. [73].) Gluons are massless particles for which the helicity formalism [74] is the proper one. For spin 1 particles this means that just the ±1 projections are possible. This issue was first discussed by Barnes [75] within the context of gluonium (see [76][77][78] for an upgraded discussion). In the helicity formalism two-gluon states J P C = 1 −+ , 1 ++ , 3 −+ , 5 −+ are forbidden in accordance to lattice results [71,79]. The glueball spectrum of two gluons might be obtained from a full fledged solution of the Bethe-Salpeter equation for two spin-1 particles [80][81][82][83]. In practice, solving these equations requires an ansatz for the kernel which is usually obtained by making reasonable approximations. For our discussion we will make some drastic approximations which actually illustrate one important point regarding gluon-glueball duality at finite temperature. We will assume the following rotational invariant and spin independent Hamiltonian in the CM system In this scheme the total spin of the glueball J is obtained by composing the gluons spin and the relative orbital angular momentum. Due to Bose statistics the total glueball wave function must be, besides a colour singlet state, fully symmetric. Thus, the spin and orbital part must be both either symmetric or antisymmetric. This gives a total degeneracy of (2l + 1)g(g ± 1)/2 for even/odd l, where g is the number of spin states of the gluon. The lowest one is the 0 ++ glueball [43]. Then the cumulative number becomes where ν l = g(g ± 1)/2 for even/odd l. A simple estimate using the uncertainty principle for the ground state is made by taking pr ∼ 1, thus For excited states, an estimate can be obtained by using the Bohr-Sommerfeld quantization condition (WKB spectrum). For s-wave, this gives 2 a 0 dr p(r) = 2π(n + α), where p(r) = (M −σ A r +3α s /r)/2, a is the classical turning point, p(a) = 0, and α depends on the boundary conditions. This produces (neglecting the Coulomb term) More systematically, the eigenvalues of the CM Hamiltonian, Eq. (63), can be computed by diagonalizing it in the harmonic oscillator wave functions basis where L l+ 1 2 n−1 (x) are associated Laguerre polynomials. The advantage of the harmonic oscillator basis is that the matrix elements of the pseudodifferential operator | p| are simply related to those of r. 4 The reduced wave functions u nl (r) are normalized to unity and satisfy the equation Here b has dimensions of length. The single-particle energies are where the oscillator frequency is ω = 1/(M b 2 ). Thus, we can compute the matrix elements H nl,n l up to some maximum values of n and l and diagonalize the truncated Hamiltonian. A list of eigenvalues (for a linear potential) for n ≤ 8 and l ≤ 13 can be looked up in Ref. [84]. The quoted accuracy is right provided we take the dimension as large as N = 100.
A derivative expansion [85] can be used to evaluate the cumulative number of the Hamiltonian in Eq. (63). In the present context this is closely related to a semiclassical expansion and a large mass expansion. A direct application of the results in [85] gives 5 H(p, r)) + . . .
Beware that the Fourier-Bessel transform introduces an additional phase (−1) n to the momentum space wave function P nl (p) ≡ ∞ 0 j l (pr)R nl (r)r 2 dr = (−1) n R nl (r = pb 2 ). 5 For simplicity here we disregard the effect of the different degeneracies in even and odd partial waves. It has no effect at leading order in . By restoring the 's in the Hamiltonian, σ A → σ A / and α s → α s , we can see that the first quantum corrections enter in the last term quoted. As advertised, the derivative expansion corresponds to a power expansion in M , with a general form N (M ) = n a n M n .
We show in Fig. 10 the comparison of the full diagonalization in the Harmonic Oscillator basis to the semiclassical large mass expansion. As we see, the agreement is rather good for a relatively low value of N (M ). It is now straightforward to compute the corresponding contribution to the trace anomaly from the expression 6 where M n are eigenvalues of H with the corresponding spin-orbit-colour degeneracies compatible with the colour singlet and Bose character of the twogluon states. The most recent lattice data from the WB collaboration [40] include 7 points for A(T ) below T c which being a dimensionless quantity depends just on the dimensionless ratio T /T c . In our case we only have the dimension √ σ and thus the dimensionless A 2g (T ) depends on T / √ σ. Thus, we can determine T c / √ σ from a fit of A 2g (T ) to the lattice data [40].
For g = 3 we get χ 2 = 3.2 for the first 6 lattice data but χ 2 /ν = 56/(7−1) when the last point is added. Actually, the value of σ is already determined by the first point, although its error decreases as new points are included in the fit. This provides some robustness to the analysis, and indicates that for T ≤ 0.9T c the trace anomaly is mainly saturated by two-gluon glueballs. The fit gives a quite reasonable result. We can proceed further in the analysis by rewriting the trace anomaly for bosons as and making use of the formula For instance, in the two-gluon sector we get Applying the results in Ref. [85] for the coefficients a i yields .
The interesting feature is that the effect of the full quantized spectrum can be very well described by the semiclassical expansion for the temperatures available on the lattice. Thus, we can use instead the semiclassical expansion to fit the parameters and sidestep the diagonalization.
As a further remark let us note that a direct attempt to fit the polynomial formula fails since there are too few data points and the curve is smooth. The correlations implicit in the 2g-Hamiltonian among the different coefficients are however enough to guide the fit. The last data point is intriguing and the most obvious candidate to fill the gap is by looking at three-gluon glueballs. The three gluon potential on the lattice has been analyzed [86] and it has been found that the triangle is preferred to the starfish configuration. In a semirelativistic framework three gluon glueballs have been addressed in Ref. [87]. Thus we can take the partonic Hamiltonian result which sets a scale separation between 2g-WKB and 3g glueballs, Whether or not 3g saturate the missing contribution, or even multigluonic states are needed remains to be seen. On the other hand, let us note that using the lowest glueballs on the lattice we found [88] that while the cumulative number could, after come coarse graining, be described as exponentially growing with T H = 2.8T c , the effect on the trace anomaly was tiny. The High-Precision Thermodynamics in connection to the Hagedorn Density of States has been discussed in [89]. A satisfactory fit with T H = 1.024(3)T c was obtained by adding a string motivated Hagedorn spectrum to the lowest 0 ++ and 2 ++ glueballs (see also [90] and [40] ). The fate of glueballs above the phase transition has been analyzed in Ref. [91].

The adjoint Polyakov loop and the gluelump spectrum
According to duality the adjoint Polyakov loop can be computed in the confined phase in terms of the gluelump spectrum. In the simplest case, this system corresponds to one massless spin-1 particle and one gluon source (infinitely heavy) which is the CM system. The Salpeter equation for the Hamiltonian operator (the Coulomb term is omitted) reads which by rescaling the coordinate can be brought to the glueball spectrum, Thus, the smallest mass gap in the pure gauge theory is the gluelump and not the glueball! This scaling relation of the spectrum implies that the partition function fulfills Z gluelumps (T ) = (2/g)Z glueballs (T / √ 2). Likewise, for the adjoint Polyakov loop at low temperatures implies,

The powers of deconfinement
One straightforward application of the glueball gas is that because they are heavy, the pressure in the confined phase is tiny, where M G T c . On the other hand, at high temperatures the pressure is due to a free gluon gas, When the trace anomaly for gluodynamics was first evaluated on the lattice [26] across the phase transition, the common wisdom was the standard textbook explanation of the deconfinement based on the MIT bag model, where hadrons are bubbles in the strongly interacting vacuum with an energy density B ∼ M H /V ∼ 1GeV/fm 3 ∼ (0.3GeV) 4 [61]. Thus the total pressure was and continuity of the pressure implies P gluons (T c ) = B for M G T c , and thus yields the trace anomaly However, this behaviour is in complete disagreement with the old [26] and newest [40] lattice data. Ten years ago we looked into the Polyakov loop and found that, contrary to expectations, one has inverse temperature power corrections above the phase transition of the form of Eq. (24) [31], with a remarkable good description of the data in a wide range of temperatures. These power corrections were quite surprising and completely unexpected, 7 since it indicated the breakdown of perturbation theory for temperatures as large as 5T c , but have been verified on the lattice [35,36] and other models [93].
One possibility to explain the trace anomaly data is to assume instead that there is a temperature dependence in the bag constant (a fuzzy bag [94], see also [88,[95][96][97] for a unified setup). Tc T 2 θ(T − T c ) with the lattice data from the WB collaboration [40].
An inspiring consequence from the power corrections of Eq. (24) is to assume which gives The result is presented in Fig. 11 and compared with the most recent lattice data from the WB collaboration [40]. As we see the agreement for T > 1.5 T c of such a simple model is impressive. These features are confirmed for any N c [98,99]. Further analyses involve renormalization group resumations [88], hard thermal loops [100] as well as holographic methods [101][102][103], but a clear physical picture of the undoubtful but mysterious powers is still lacking.

Chiral quark models at finite temperature
Up to now, we have assumed that duality holds without specifying how this might occur from a microscopic point of view. Here we want to illustrate how duality arises from quark models and the relevant provisos to achieve this goal. A full QCD argument for the hadronization of the Polyakov loop was advanced and elaborated in Refs. [34,37]. Chiral quark models have been used (e.g. [104] and references therein) for a long time to describe the chiral phase transition. It was then realized that they lead to a wrong N c counting of thermal corrections, since they are described as one quark loop and hence are O(N c ). In reality, they correspond to zero point energy of meson states which are O(N 0 c ). For instance, in the quark condensate one has qq 0 = O(N c ) but qq T − qq 0 = O(N 0 c ) (wrongly counted by the CQM as O(N c )). This observation was well known although surprisingly nothing was done about it. The requirement of large gauge invariance motivated us the introduction of the Polyakov line as an independent variable [105] which has a local and quantum character. Previous works had already dealt with this coupling [106,107] initiating the Polyakov-Nambu-Jona-Lasinio (PNJL) saga [108][109][110][111][112][113][114][115][116][117] where the Polyakov line has been taken global and classical. This allowed to study the interplay between chiral symmetry restoration and breakdown of the center symmetry although, as we have repeatedly pointed out in our previous works [105,[118][119][120][121], this generates an undesirable ambiguity of group coordinates as well as a non-vanishing value for the Polyakov loop in the adjoint-representation, contradicting lattice simulations. 8 These difficulties may be overcome [105,118] by recognizing the local and quantum nature of the Polyakov loop. Using this interpretation we have shown that a gateway to the hadron resonance gas may be established if the motion of quarks in the field generated by the Polyakov loop is quantized.
We note that on the lattice [124] the HRG has been deduced in strong coupling and for large N c by considering hadrons at rest.

The Quantum and local Polyakov loop
As we have mentioned above, large gauge invariance is an important restriction at finite temperature which breaks down in perturbation theory. In the Polyakov gauge one can automatically implement large gauge invariance by considering the Polyakov loop line Ω( x) as an independent variable, which in the Polyakov gauge becomes a diagonal unitary matrix in colour space. This is equivalent to a minimal coupling scheme in the time derivative of a dynamical quark or gluon field. We refer to Refs. [105,[118][119][120] for further motivation.
Following this prescription, the partition function of the Chiral Quark Model (CQM) at finite temperature can be written as where Ω = e igA 0 /T and DΩ is the invariant SU(N c ) Haar group integration measure, for each SU(N c ) variable Ω( x) at each point x. Here the action is The fermionic contribution depends on the quarks (and anti-quarks), and it is obtained from the corresponding fermion determinant. Assuming massdegenerated quarks for simplicity the effective action reads where c.c. stands for the complex conjugated contribution stemming from the anti-quarks and 2 is the spin degeneracy factor for spin 1/2 particles.
Here E q (p) = p 2 + M 2 q is the energy of a quark with total mass M q = M 0 + m q , and M 0 is the constituent mass. As one can see the diagonal part of the Polyakov loop corresponds to consider chemical potentials for different color species. Large color gauge invariance is implemented by just averaging over group elements. We will check whether this minimal coupling scheme for the Polyakov line complies with known QCD properties [34,37], such as the fact that the expectation value of the Polyakov loop admits a hadronic representation.

Hadronic representation of the Polyakov loop
In order to see this, consider a system with N f dynamical quarks and an extra heavy quark (not antiquark) of an arbitrarily large mass m H at rest located at a fixed point and with fixed spin and colour a = 1, . . . , N c . From Eq. (88) the change in the effective action is yielding the partition function after averaging over color degrees of freedom implied by DΩ. Thus we get To evaluate the r.h.s. we explicitly separate in the corresponding HRG hadrons composed of hadrons made of N f dynamical quarks and one extra heavy quark. The mass of such hadrons can be written as where by definition in ∆ it has been subtracted the mass of the heavy quark. On the other hand, the HRG partition function with N f + 1-flavors with this extra heavy quark H can be separated into hadrons containing it or not. To do this, we reinstate the finite box quantization conditions on the momentum p to make sense of the limit m H → ∞ and V → ∞ (the Compton wavelength of the heavy quark is shorter than the box size) and get Quark-Hadron duality at this level implies Z HRG = Z CQM so that we get, 1 N c tr c Ω = 1 2 α g α e −∆α/T , providing confidence on the assumed minimal coupling of the Polyakov line to quarks.

From chiral quark models to the hadron resonance gas
The previous models can be pictured as multiquark states which are created or annihilated at point x and momentum p with factors Ω( x)e −Ep/T and Ω( x) † e −Ep/T . This corresponds to classical particles with internal quantum numbers and statistic, i.e. to a second quantized but not first quantized formalism.
At low temperatures quark Boltzmann factors are small e −Ep/T 1, and the quark contribution to the action becomes small where we have quantized the Heavy-light ground state system as H(p, x)ψ n = p 2 + m 2 q + σr ψ n = ∆ n ψ n .
In previous works the quantization of the quark motion was not considered, and as a consequence we failed to see the connection to the HRG. In particular, we found that L T ∼ e −M 0 /T so we proposed to determine the constituent quark mass from an analysis of the Polyakov loop on the lattice at low temperatures. 9 While fits along these lines turned out to provide too large a constituent mass we preferred to keep the L T ∼ e −M 0 /T behaviour with a suitable proportionality factor [119]. As we have seen here the Boltzmann factor contains the gap ∆ (the heavy-light meson mass) which according to our discussion above corresponds to twice the constituent quark mass, so that L T ∼ e −2M 0 /T and now the value for M 0 from the lattice is phenomenologically acceptable. A similar argument holds for the correlation function between Polyakov loops. This requires an assumption for a four-point correlator in the pure gluonic theory which at low temperature we assume to be tr c Ω( x) † tr c Ω(0)tr c Ω( x 1 ) † tr c Ω( x 2 ) G = e −σr/T e −σr 12 /T + e −σ| x− x 2 |/T e −σr 1 /T (104) and satisfies cluster decomposition properties. This yields a disconnected piece due to the effect of the quark determinant, tr c Ω( x)tr c Ω † (0) = e −σr/T [1 + Zq q + . . . ] + | tr c Ω | 2 1 + Zq q + . . .
This reproduces the free energy formula Eq. (55) based on a avoided crossing structure of theQQ energy levels, see Fig. 5, and yields the result for F 1 (r, T ) sketched in Fig. 6.

Gluon models with Polyakov loop
We may introduce gluon fields besides those involved by the Polyakov line variable. A particularly interesting case is gluodynamics. At one loop the effective potential has been computed in [126]. More recently, one loop gluon actions with Polyakov loops in the adjoint representation were suggested in Refs. [127][128][129][130], where the background gauge for the classical gluon-field was assumed and a Polyakov gauge for the quantum fluctuations where H 2 is the two-gluon Hamiltonian whose spectrum provides the two gluon glueballs discussed in the previous section, and Tr 2 is the trace in the corresponding two-gluon Hilbert space. Note that the factor 1/2 corresponds to the correct Boltzmann counting which originates in the quantum statistical mechanics to avoid the Gibbs paradox in the high-temperature limit [131]. Following similar steps as before we can obtain the gluelump representation of the Polyakov loop in the adjoint representation. One can analyze what happens for higher Fock states since the glueball gas corresponds to having log Z = − n Tr n log 1 − e −Hn/T , where here Tr = n Tr n represents the trace over the whole multigluon Hilbert space. Pursuing the low temperature expansion to higher orders we have checked that up to three gluons included, the mapping with the glueball mass works but ambiguities arise for 4-gluon states, where the possibility of forming separate and weakly interacting 2-gluon glueballs first arises. We are facing again the concept of colour irreducible clusters inside colour neutral states [37] and the very definition of a hadron.

Conclusions
Quark-Hadron Duality at finite temperature is the statement that at low temperatures hadrons can be considered a complete basis of states. The naive hadron resonance gas, while simple minded, works well enough at sufficiently high temperatures as to deserve dedicated attention on why becomes this picture invalid. This success remains a mystery over the years since Hagedorn first proposed it.
The listed PDG states incorporate currently just the qq or qqq states which fit into the conventional quark model, but what is the nature of states that are needed when approaching the crossover from below? As we have seen, saturating at subcritical temperatures requires many hadronic states, and so the excited spectrum involves relativistic effects even for heavy quarks. This point makes relativistic quark models a potential source for investigation of the hadron spectrum from a global and thermodynamic perspective since the number of needed excited states challenges any lattice QCD calculation. Moreover, the hadronic mapping of Polyakov loops in fundamental and higher SU(N c ) colour group representations allows to deduce multiquark states, gluelumps and hybrid states, containing one or several heavy quark or gluon sources. This goes beyond the models and opens up the possibility of a Polyakov loop spectroscopy including exotics. While the question on what is the complete hadronic particle spectrum remains still open, we envisage the possibility of grasping the yet unknown physics of the phase transition by enquiring this question at the lowest possible temperature where the hadron resonance gas picture fails.
One of us (E.R.A.) warmly thanks Michal Prasza lowicz for the invitation and for providing a stimulating atmosphere. We also thank Wojtek Broniowski for many discussions and Marco Panero for clarifications. This work has been supported by Plan Nacional de Altas Energías