Baryonic susceptibilities, quark-diquark models and quark-hadron duality at finite temperature

Fluctuations of conserved charges such as baryon number, electric charge and strangeness may provide a test for completeness of states in lattice QCD for three light flavours. We elaborate on the idea that the corresponding susceptibilities can be saturated with excited baryonic states with an underlying quark-diquark structure with a linearly confining interaction. Using Polyakov-loop correlators we show that in the static limit, the quark-diquark potential coincides with the quark-antiquark potential in marked agreement with recent lattice studies. We thus study in a quark-diquark model the baryonic fluctuations of electric charge, baryon number and strangeness: $\chi_{BQ}$, $\chi_{BB}$ and $\chi_{BS}$; by considering a realization of the Hadron Resonance Gas model in the light flavor sector of QCD. These results have been obtained by using the baryon spectrum computed within a relativistic quark-diquark model, leading to an overall good agreement with the spectrum obtained with other quark models and with lattice data for the fluctuations.

Quantum Chromo-Dynamics (QCD) is the fundamental non-abelian gauge theory of strong interactions in terms of 2N c N f quarks and antiquarks and 2(N 2 c − 1) gluons with N c = 3 the number of colours and N f = 6 the number of flavour species u, d, s, c, b,t. Quark and gluon confinement requires all physical states to be colour singlet but, are the hadronic states a complete set of eigenstates of QCD spanning the Hilbert space H QCD ? This question is related to the validity and meaning of quark-hadron duality.
In the case of N f = 3 flavours, which will be assumed throughout the paper, the stable and low lying bound states, such as the baryon octet and the pseudoscalar nonet are unambiguously part of the discrete spectrum. 1 The remaining states belong to the continuum spectrum and are experimentally spotted in strong hadronic reactions, interpreted as unstable resonances and characterized by a mass and a width. They have been reported over the years by the PDG booklet [1] as single states rated with *,**,***,**** (for baryons) depending on the increasing confidence on their existence. However, are the PDG states complete and if so in what sense would they be complete?
For such resonance states, the verification of completeness for QCD in the continuum is subtle since within a Hamiltonian perspective they are not proper (normalizable) eigenstates of the QCD Hamiltonian and besides they correspond to unconventional representations of the Poincaré group [2]. Lattice QCD presents the clear advantage that in a finite box all states are discretized and hence become countable; below a certain maximal mass the number of states is finite as long as the volume V remains finite (typically V ∼ (2 − 3 fm) 3 ). Resonances are extracted from those particular energy levels which become insensitive to the box size. The extent to which these states play a key role in the completeness issue is uncertain since, although there is a larger concentration of states around the resonance, in the bulk, the mass separation is ∆M ∼ V −1/3 .
In recent years, the thermodynamic approach to strong interactions pioneered by Hagedorn [3,4] where the vacuum is represented by a non-interacting Hadron Resonance Gas (HRG), has emerged as a practical and viable path to establish completeness of hadronic states on a quantitative level in the hadronic phase. With all the provisos regarding the nature of resonance and bound states, the most impressive and vivid verification has been the recent study of the trace anomaly, (ε − 3P)/T 4 with ε energy density and P the pressure. It was computed directly in lattice QCD by the Wuppertal-Budapest (WB) and the HotQCD collaborations [5,6] and the HRG model using the most recent compilation of the PDG states taking just their masses, i.e. assuming zero widths [1] and for temperatures below T ∼ 170 MeV. Remarkably, a similar degree of success is achieved within uncertainties [7] (see Figs. 1 and 2) using the hadronic spectrum obtained in the Relativized Quark Model of Capstick, Godfrey and Isgur [8,9] which predated the lattice QCD calculations by about 30 years. Finite width effects on the PDG (PDG-Γ) naturally provide a shift towards lower masses [10,11] as a consequence of the mass spectrum spread weighted by the exponentially decreasing Boltzmann factor. This remarkable agreement becomes significantly spoiled when susceptibilities involving conserved charges such as the baryon number, the electric charge and the strangeness are considered [12,13]. In particular, the differences between PDG, PDG-Γ and RQM become more visible and an excess of baryonic states as compared to the Lattice-QCD results is observed [14] illustrating the so-called missing resonance problem.
In fact, since the early days of the quark model the extreme abundance of predicted and experimentally missing baryonic resonances has been a major cause of concern both at the theoretical as well as at the experimental level [15] (see also [16] and [17] for reviews and references therein). Possible ways out of the difficulties have traditionally been attributed either to a weak coupling of the predicted states to the particular production process (photo-production, πN scattering etc.) or to a dynamical reduction of degrees of freedom due to diquark clustering.
The suspicion that the baryonic spectrum can be understood in terms of quark-diquark degrees of freedom is rather old (see e.g. Ref. [18] for a review, but also Ref. [19] for evidence against it). This includes diquark clustering studies [20], non-relativistic [21] and relativistic [22,23] analyses where scalar and axial-vector diquarks have a mass of about 600 MeV and 800 MeV respectively (the diquark mass difference seems quite model independent and about 200 MeV). As expected, in diquark models many states predicted by the quark model do not appear [9]. More specifically, while not strictly forbidden, the relativistic diquark model [22,23] does not predict any missing states below 2 GeV, whereas Isgur and Capstick have predicted 5 unobserved states [9]. Lat-tice QCD has also provided insights, as some evidence on diquarks correlations in the nucleon [24] and the dominance of the scalar diquark channel [25] have been reported. Moreover, the diquark-approximation has been found to work well in Dyson-Schwinger and Faddeev equations studies [26]. The radial Regge behavior found in the relativistic quark-diquark picture from a numerical study of the spectrum of the relativistic two body problem [27] has been confirmed from a direct PDG analysis when the widths of the resonances is implemented in the analysis [28,29].
While the missing resonance problem has attracted a lot of interest both experimentally as well as theoretically, much of the discussion is focused on the individual one-to-one mapping of resonance states which have a mass spectrum and which are produced with different backgrounds. In contrast, the thermodynamic approach offers the possibility to perform a more global analysis where many of the fine details will hopefully be washed out due to the presence of the heat bath. In the present paper we profit from the new perspective provided by Lattice QCD at finite temperature, based on separation of quantum numbers with the study of susceptibilities of conserved charges, where a combination of degeneracy and level density is involved. We try to answer the question whether or not quark-diquark baryonic states saturate the baryonic susceptibilities below the deconfinement crossover as compared to the available lattice QCD calculations [12,13]. Aspects of quark-hadron duality at finite temperature have been discussed in pedagogical way in Ref. [7]. Actually, in a recent work it has been discussed how quark-hadron duality in deep inelastic scattering for baryons suggests an asymptotic quark-diquark spectrum with a linearly rising potential for scaling to hold in the structure functions [28,29]. Motivated by this we want to establish if a similar pattern holds also at finite temperature.
The paper is organized as follows. In Section II we review the relevant aspects of the HRG model from the point of view of the Equation of State and the trace anomaly as compared to lattice QCD. In Section III we analyze the implications for fluctuations of conserved charges at finite temperature in the vacuum. then analyze in Section IV the quark-diquark potential obtained as a correlation function involving Polyakov loops. This allows to define our model and compute the spectrum by diagonalization in Section V where the susceptibilities are analyzed in terms of the free parameters of the theory. Finally in Section VI we come to the conclusions. In the appendices we provide details on the semiclassical determination of the spectrum and also prove a theorem on the sign of susceptibilities which is verified by lattice calculations.

II. COMPLETENESS AND THERMODYNAMIC EQUIVALENCE
A. QCD spectrum and thermodynamics The cumulative number of states may be used as a characterization of the QCD spectrum [31,32]. It is defined as the  [8,9]. Right panel: Trace anomaly as a function of temperature in lattice QCD [5,6] vs HRG using PDG (dashed) and RQM (solid) spectra. We also plot just the contribution of states with M < 0.6 GeV (dotted) and M < 0.8 GeV (dotted-dashed).
number of bound states below some mass M, i.e.
where g i is the degeneracy factor, M i is the mass of the i-th hadron and Θ(x) is the step function, so that the density of states is given by ρ(M) = dN(M)/dM. A way to provide a practical meaning of the previous equation is to use finite box periodic (antiperiodic) boundary conditions for gluon/quark fields. In such a case all states contribute on equal footing. In such setting, stable bound states correspond to eigenvalues which do not depend on the volume of the box for sufficiently large boxes, whereas unstable resonance states correspond to eigenvalues which are volume independent within a given volume interval [33]. Colour neutral eigenstates of the QCD Hamiltonian and in fact excited states have been determined on the lattice in this way [34]. Despite its importance, the cumulative number of states, Eq. (1) has never been evaluated directly in QCD. Within a thermodynamic setup, the issue of completeness acquires a precise meaning where the global aspects of the spectrum rather than individual features are highlighted. The partition function of QCD is given by the standard relation and is the fundamental quantity to study the thermodynamic properties of the theory. Written in terms of the eigenvalues of the QCD Hamiltonian, i.e. H QCD ψ n = E n ψ n , Eq. (2) illustrates the relation between the thermodynamics of the confined phase and the spectrum of QCD. Although Z QCD has been determined independently as an Euclidean path integral on the lattice [12,13], to our knowledge the energy levels contribution to Eq. (2) has not been tested explicitly against the partition function. 2 Actually, the primary quantity is the trace anomaly [12,13] ∆ ≡ ε − 3P 2 As is well known reflection positivity of the Euclidean action on the lattice guarantees the existence of a transfer matrix and hence of a Hamiltonian [35]; the issue is to list the pertinent set of eigenvalues. whence the equation of state (EoS) can be obtained by integration the second equality with suitable boundary conditions and physical quantities such as the entropy density s = ∂ T P or the sound velocity c 2 s = ∂ P/∂ ε ≡ ∂ T P/∂ T ε may be determined.
From this thermodynamic point of view it is fair to say that the completeness of the QCD-spectrum remains to be checked in practice.

B. Hadron Resonance Gas model: the PDG and RQM spectra
The completeness of the listed PDG states [36] is equally a subtle issue. On the one hand they are mapped into the qq and qqq quark model states. On the other hand, most reported states by the PDG are not stable particles but resonances which are produced as intermediate steps in a variety of scattering process.
The HRG model was originally proposed by Hagedorn [4]. In spirit, this is valid under the assumption that physical quantities in the confined phase of QCD admit a representation in terms of hadronic states, which are considered as stable, noninteracting and point-like particles. Within this approach the EoS of QCD is described in terms of a gas of non-interacting hadrons [4,37], and the grand-canonical partition function turns out to be with ζ i = ±1 for bosons and fermions respectively. We consider several conserved charges labeled by the index a, with q i a the charge of the i-th hadron for symmetry a, and µ a the chemical potential associated to this symmetry. The obvious consequence is that a good understanding of the spectrum of QCD turns out to be crucial for a precise determination of the thermodynamic properties of this theory.
In QCD, the quantized energy levels are the masses of lowlying and stable colour singlet states, which are commonly identified with mesons [qq] and baryons [qqq] in the quark model. As already mentioned, unstable hadronic resonances are not proper eigenstates which in the HRG model are regarded as bound states. We reiterate that while there is currently a large scale on-going effort to determine hadronic resonances from lattice-QCD by means of the Lüscher formula [34] (see e.g. [33] for a review) these states are not numerous enough to verify the HRG model.
So far, the states listed by PDG echo the standard quark model classification for mesons [qq] and baryons [qqq]. Then, it would be pertinent to consider also the Relativized Quark Model (RQM) for mesons [8] and baryons [9]. 3 We show in Fig. 1 the hadron spectrum with the PDG compilation (left) and the RQM spectrum (right). The comparison clearly shows that there are further states in the RQM spectrum above some scale M > M min that may or may not be confirmed in the future as mesons or hadrons. In addition there could be exotic, glueballs or hybrids states, predicted by other hadronic models.
The contribution to N(M) obtained just adding the qq (mesons), qqq (baryons) andqqq (antibaryons) components as 4 can be evaluated either from the PDG booklet [1] or alternatively from the RQM of Capstick, Godfrey and Isgur [8,9]. The resemblance of both schemes below 1.7 GeV, see Fig. 2 (left), is noteworthy, particularly if we take into account the 30 years elapsed between the RQM and the current PDG. Finite width effects on the PDG (PDG-Γ) have been estimated [10,11] and naturally provide a shift towards lower masses as a consequence of the mass spectrum spread. The flattening of the curves indicates lack of reported states on the PDG side, as well as a the higher computational cut-off M high = 2.3 GeV in the numerical calculation on the RQM side. This resemblance is also realized at the level of the trace anomaly which in the HRG model is given by and compared with lattice QCD by the Wuppertal-Budapest (WB) and the HotQCD collaborations [5,6] in Fig. 2 (right), suggesting the thermodynamic equivalence of QCD, PDG and RQM below the crossover.

C. Missing states
Given the resemblance of the PDG to the RQM we may speculate on the nature of mesonic and baryonic states separately. While the PDG is a consented compilation of numerous analyses, the RQM corresponds by construction to a solution of the quantum mechanical problem for both qqmesons and qqq-baryons. For colour-singlet states the nparton Hamiltonian takes the schematic form Neglecting spin dependent contributions, not essential in the following argument, and assuming Casimir scaling, the twobody interactions take the form Asymptotic estimates may be undertaken to sidestep the actual numerical evaluation of the eigenvalues. Specifically a semiclassical expansion describes the high mass spectrum for systems where interactions are dominated by linearly rising potentials with a string tension σ , in the range M √ σ . Such an approach relies on a derivative expansion of the cumulative number of states for a given Hamiltonian [38]. 5 At leading order the number of states below a certain mass M takes the form where g n takes into account the degeneracy. For the sake of the argument, let us neglect the Coulomb term in (8), thus v(r) = σ r, as well as the current quark masses. In this case, a dimensional argument, p → M p, r → Mr/σ , gives Using these techniques, one can predict that the large mass expansion of these contributions is The separate contributions of mesons and baryons are presented in Fig. 3 on a log-log scale where we see that again PDG and RQM largely agree and present an approximate linear behaviour on this scale, indicating as expected a power behaviour. However, while in the meson case the N [qq] ∼ M 6 seems to conform with the asymptotic estimate, in the baryon case much lower powers, M 6 − M 8 , than the expected N [qqq] ∼ M 12 are identified. We note that M 6 suggests a similar twobody as in the case of mesons (with a linearly rising potential). We take this feature as a hint that the qqq excited spectrum effectively conforms to a two body system of particles interacting with a linearly growing potential, as we will analyze below in detail in Section V B.

III. FLUCTUATIONS OF CONSERVED CHARGES IN A THERMAL MEDIUM
While mesonic and baryonic contributions can be explicitly distinguished the thermodynamic separation of mesons and baryons in the EoS cannot be done at the QCD level. In order to achieve directly such a separation for baryons in particular we analyze fluctuations containing at least one baryonic charge as well as charge and strangeness, the three of them being conserved charges in strong interactions.
Conserved charges [Q a , H] = 0 play a fundamental role in the thermodynamics of QCD. In the (uds) flavor sector of QCD the conserved charges are the electric charge Q, the baryon number B, and the strangeness S. While their thermal expectation values in the hot vacuum are vanishing, i.e. in the absence of chemical potentials Q a T = 0, where Q a ∈ {Q, B, S}, they present statistical fluctuations characterized by susceptibilities [13,14,39,40] The susceptibilities can be computed from the grandcanonical partition function by differentiation with respect to the chemical potentials, i.e.
where Ω = −T log Z is the thermodynamical potential. 6 One can also work in the quark-flavor basis, Q a ∈ {u, d, s}, where u, d and s is the number of up, down and strange quarks. In this basis B = QCD at high temperature behaves as an ideal gas of quarks and gluons. In this limit the susceptibilities approach to hence for N c = 3 and flavors u, d, s Within the HRG approach, the charges are carried by various species of hadrons, where K 2 (z) refers to the Bessel function of the second kind. 7 This formula will be used to compute the baryonic susceptibilities, namely χ BB , χ BQ and χ BS . Eq. (16) predicts the asymptotic behavior where M 0 is the mass of the lowest-lying state in the spectrum with quantum numbers a and b. Therefore where M p = 938 MeV is the proton mass, M Λ 0 the mass of the Λ 0 baryon, etc. These relations go beyond the HRG model, since for each sector the lightest hadron must certainly saturate the QCD partition function at low enough temperature. This observation makes it appealing to plot the lattice data for the fluctuations in a logarithmic scale. These plots are shown in Fig. 4. Of the six susceptibilities only four are independent, on account of isospin symmetry which requires χ uu = χ dd and χ us = χ ds . On the other hand, the following inequality holds in QCD for degenerate u, d flavors at any temperature, even in the deconfined phase. This is proven in App. B. Because strange and light-quark flavors are not exactly degenerate, the relation χ us (T ) ≤ 0 does not follow as a QCD theorem but it is nevertheless supported by the lattice results. In addition, the following relation holds for any pair of charges, since (∆Q a ± ∆Q b ) 2 is non-negative. In particular it follows The meson contribution to χ ud is always negative while the baryon contribution is always positive. The negative global result indicates that mesons dominate over baryons in χ ud . Going to low temperatures, this provides yet another proof that in QCD with two degenerate light flavors the lightest meson is lighter than the lightest baryon. We display in Fig. 5 some of these relations for the lattice data of the susceptibilities, and check that they are fulfilled.

IV. QUARK-DIQUARK POTENTIAL FROM POLYAKOV LOOP CORRELATORS
In view of the missing resonance problem alluded to previously there is nowadays some discussion about the most probable spatial configuration of quarks inside baryons, and more specifically the structure of excited states. In the quark model such as the RQM [8,9] baryons are qqq states where the interaction is given by a combination of ∆-like pairs of qq interactions and a genuinely Y -like qqq interaction. Remarkably the Y-stringlike behavior of a static baryon energy at finite temperature has been observed in Ref. [41]. An interesting possibility would be that the quarks are distributed according to an isosceles triangle. This is the idea behind an easily tractable class of models, the so-called relativistic quark-diquark (qD) models with non-relativistic [21] and relativistic [22,23] variants where a linearly rising potential is shown to work phenomenologically. Besides this success there was no further reason to invoke such a behaviour. In the next paragraphs we elaborate on this and provide a theoretical argument in favor of the presence of a linear potential.
It should be noted that by definition, a potential can only be evaluated unambiguously (up to an additive constant) in the heavy particles or static limit, since then the position operator is well defined. The calculation of qD static interactions for heavy sources has been addressed on the lattice in several works [42][43][44], where it has actually been found that up to numerical uncertainties the potential is linearly rising and the corresponding string tension is numerically identical with the one of the qq system (see specifically Ref. [44]). In this section we show analytically that under very specific assumptions this must actually be true. To show how this comes about we will rely heavily on our previous work on free energies [45] where a more detailed discussion can be found. Here we just summarize the main issues relevant to the problem at hand.
An operational way of placing static sources in a gauge theory such as QCD, is by introducing in the Euclidean formulation a local gauge rotation Ω realized by the Polyakov loop, We display as dots the lattice data from Refs. [13] (blue) and [12] (red). We also display the HRG model results including the spectrum of the RQM [8,9]   i.e., the gauge covariant operator defined as where P indicates path ordering and A 0 and x are Euclidean. The colour source may be decomposed into irreducible representations, say µ so that one can build the gauge invariant combinations corresponding to the character in that representation, χ µ [Ω(x)]. 8 For instance, the trace in the (anti)fundamental representation corresponds to the character of the SU(3) colour gauge group, 8 χ µ (g) stands here for the character of the element g in the representation µ, not to be confused with the susceptibilities χ ab introduced earlier.
Thus, the qq free energy is given by a thermal expectation value where due to translational invariance the free energy depends only on the separation r = |x x x 1 − x x x 2 |. In this and the following expressions, an ambiguity of an additive constant must be allowed in the free energies coming from the renormalization of the Polyakov loop operator. The potential is obtained as the zero temperature limit of the free energy V qq (r) = F qq (r, 0) .
Likewise the qqq free energy is given by In the limit x x x 3 → x x x 2 we get, for the quark-diquark free energy This limit is singular since at very small distances the interaction is dominated by one gluon exchange, ∼ 1/r, and a self-energy must be added. The renormalization of the new composite operator tr(Ω(x x x 2 ) 2 ) yields a new ambiguity in the form of an additive constant from Using the Clebsch-Gordan series yields the equivalent character relations and therefore The qD potential is obtained in the zero temperature limit. In this limit we expect the energy 6 ⊗ 3 configuration to be larger than that of the qq one, hence

A. The model Hamiltonian
In these models, the baryons are assumed to be composed of a constituent quark, q, and a constituent diquark, D ≡ (qq) [46], and the Hamiltonian writes 9 H qD = p p p 2 + m 2 q + p p p 2 + m 2 D +V qD (r) .
In Sec. IV we have argued that the quark-diquark potential is the same as the quark-antiquark potential, up to an additive constant. For the latter we assume V qq (r) = − τ r + σ r + c, hence with σ = (0.42 GeV) 2 . In addition, we adopt the value τ = π/12, as expected from the Lüscher term in the potential [47]. In this model we can distinguish between two kinds of diquarks: scalar, D, and axial-vector, D AV . When considering the quark content notation of a diquark, we will use [q 1 q 2 ] to denote scalar diquarks, and {q 1 q 2 } for axial-vector diquarks. 9 In this work we are concerned with the overall features of the quark-diquark spectrum, relevant to the thermodynamics of the system, hence we will consider a simplified version of the model, neglecting possible fine interaction terms, like contact terms in the s-channel and spin-dependent interactions, cf. Ref. [46]. Some studies on QCD indicate that the mass difference between these diquarks is [48] ∆m D := m D AV − m D 0.21 GeV .
In what follows we adopt this value for ∆m D . We will assume that the mass parameters of the model are controlled by a constituent quark mass, m cons , and the current quark mass for the strange quark,m s , in the following way: The breaking of flavor SU(3) for diquarks will be modeled as where m D,ns is the scalar mass of non strange diquarks, and n s = 0, 1, 2, is the number of s-quarks in the diquark. The further choice m D,ns = 2m cons (37) is rather natural, but we will not always enforce it. With these assumptions, the only free parameters of the model are m cons ,m s and µ, plus m D,ns if (37) is not enforced. The goal is to reproduce the lattice results for the baryonic fluctuations χ BB , χ BQ and χ BS for temperatures below the crossover. For the baryonic states B = 1 and S = −n s . The electric charges and the degeneracies of the states are summarized in Table I. The spectrum is computed as explained below. With these ingredients the susceptibilities can be evaluated in the model using Eq. (16).

B. Baryon spectrum
In order to obtain the baryon spectrum with the quarkdiquark model we have to diagonalize the Hamiltonian of Eq. (32). Since the problem does not admit a closed analytic solution we will obtain the spectrum numerically after truncation to a model space and diagonalization of the corresponding finite dimensional matrix. This is a variational procedure. Because of the form of the Hamiltonian H = f 1 (p) + f 2 (r) a convenient basis is that of the isotropic harmonic oscillator, with normalized wave functions of the form L l+ 1 2 n−1 (x) are the generalized Laguerre polynomials and the positive parameter b is related to the oscillator mass and frequency. b can be optimized as a variational parameter. The reduced wave functions u nl (r) are normalized to unity The corresponding wave functions in momentum space, R nl (p), have the same form as those in Eq. (38) up to a phase and b → 1/b, namelŷ Then the matrix elements nl|H qD |n l are obtained from with n = 1, . . . , n max , and l = 0, . . . , l max . The value of the parameter b is fixed to minimize the averaged value of the energy levels in the spectrum for each of the multiplets in Table I. Typical values of this parameter are in the range 0.55 fm < ∼ b < ∼ 0.65 fm. We display in Fig. 6 the dependence of the mass of the baryonic state, Λ 0 , as a function of n max , taking l max = n max − 1. We can see that the dependence on n max is very weak already for n max = 4. We have cross-checked many of the results presented subsequently by including larger values of n max and l max , and we find that they do not change appreciably.
After considering a convenient choice of the parameters of the model, for instance and following the procedure mentioned above, we get the spectrum of baryons that is shown in Fig. 7. In this figure we compare this result to the RQM spectrum [9] on a loglog plot where the ∼ M 6 growth of the quark-diquark baryonic spectrum can be clearly identified. It is quite remarkable that below M < 2400 MeV the quark-diquark spectrum is in good agreement with the RQM spectrum. While the authors of Ref. [9] do not compute baryon masses heavier than this, with the present quark-diquark model we have obtained further states up to M ≈ 3400 MeV. Within the HRG picture, these states will contribute to the EoS of QCD as well as to   [9], and the quarkdiquark model used in this paper. We have used the parameters in Eq. (42). We also draw the ∼ M 6 line for illustration. other thermal observables like the fluctuations. The choice of parameters of Eq. (42) is motivated by a comparison with the lattice results for the thermal fluctuations, as we will see in the next section.
Let us mention that, unless otherwise stated, in the following we will use the empirical value of the mass for the nucleon, M n = 938 MeV, and apply the quark-diquark model only for the other baryons. The effect of using the empirical mass for all the baryons in the 1/2 + octet is also discussed. The justification for this is that it is expected that the quarkdiquark picture will be reliable only for excited states. As we will show in the next subsection, this is confirmed from the analysis of the lattice data for the baryonic fluctuations.

C. Baryonic fluctuations
From the spectrum of the quark-diquark model, we can obtain the baryonic fluctuations by using the HRG approach given by Eq. (16). Our goal is to reproduce the lattice results for these quantities, at least for the lowest temperature values. A typical fit of the model prediction with lattice data is shown in Fig. 8. While ideally one would like to have as lowest temperatures as possible so as to determine in a clean way the low-lying states with the proper quantum numbers, in practice we find that at the lowest available temperatures, the contribution of excited states becomes individually small but collectively important. This is a typical problem in intermediate temperature analyses, and this is the reason why all possible constraints on the model, such as the identity of quarkdiquark and quark-antiquark potentials discussed above, are particularly welcome.
In order to perform the best fit to the data, we have chosen to minimize the function Here the T j are the temperatures used in the lattice calculations. The lowest temperature of the data is T 1 = 125 MeV for Ref. [12] and T 1 = 150 MeV for Ref. [13], while j max is the number of data points used in the fit. ∆χ lat ab are the uncertainties of the lattice results.
Since a hadronic model is not expected to reproduce the QCD crossover we fit the data corresponding to the lower temperatures in lattice measurements. Statistical considerations [49] indicate that those data points should be included for which where ν is the number of degrees of freedom. This condition fixes our choice of j max . In our case, ν = n obs j max − n param , with n obs = 3 since three baryonic fluctuations are being fitted, and n param is the number of free parameters of the model. The fits turn out to be very sensitive to the parameter µ in Eq. (33), hence it is convenient to always minimizeχ 2 with respect to this parameter. Whenever we provide or plot a value for the functionχ 2 /ν, it should be understood that this function is already minimized with respect to the parameter µ. Typical values of this parameter are −0.7 GeV < ∼ µ < ∼ 0 GeV. As a first step in the analysis, we show in Fig. 9 a plot of χ 2 /ν as a function of the current quark mass, while fixing the constituent mass to m cons = 0.3 GeV and m D,ns = 2m cons . When including in the fits only the lowest temperature point of the data, i.e. j max = 1, one gets an exceedingly good fit withχ 2 /ν = 1.0 × 10 −4 form s = 0.130 GeV. This is rather surprising, as it means that three central values of the data can be fitted almost exactly with just two parameters:m s and µ. The fit deteriorates but it is still acceptable when the four lowest temperatures are used, j max = 4. We find that the model and the lattice data are compatible with Eq. (45) for temperatures T < ∼ 165 MeV, either if we analyze the lattice data of Ref. [12] or [13]. The best fit in this case withχ 2 /ν = 0.68 is form s = 0.099 GeV and µ = −0.459 GeV. The left panel clearly indicates that the current quark mass takes a value compatible with the PDG, i.e. 80 MeV < ∼m s < ∼ 120 MeV. In addition, the value of the constituent quark mass cannot be determined with precision, but at least we can ensure that it is in the regime 100 MeV < ∼ m cons < ∼ 400 MeV. In the right panel we have fixedm s = 0.10 GeV, and used the scalar mass of non strange diquarks, m D,ns , as a free parameter. The figure shows that the most probable scalar diquark mass is of the order of m D,ns 0.4 − 0.5 GeV and the mass of the constituent quarks m cons 0.3 GeV. However, these values depend of the choice of the current quark mass, so that when increasing the value ofm s , the best fits happen for lower values of m cons . For instance form s = 0.12 GeV the best fit is obtained when m cons < ∼ 0.2 GeV, while form s = 0.90 GeV one obtains values m cons 0.5 − 0.6 GeV.
As mentioned above, for our fits we have taken the value of the PDG for the nucleon mass, M n = 938 MeV, and used the quark-diquark model values for the remaining baryonic states.
One could further adopt the empirical values for the masses of the 1/2 + baryonic octet, that is, M Λ 0 = 1.116 GeV, M Σ = 1.193 GeV, and M Ξ = 1.318 GeV, in view of the fact that they are stable under strong interactions. Such states correspond to 16 of the 18 states of the type [qq]q in Table I [46]. Since the mass of [nn]s is slightly lighter than [ns]n we assign the Λ 0 to this multiplet; the opposite assignation produces very similar results. The effect of using the empirical masses for the whole octet is that the best fits presented in Figs. 9 and 10 (left) are shifted to larger values ofm s , namely ∆m s 24 MeV. The quality of the fits is similar.
The value of the nucleon mass as predicted by our version of the quark-diquark model is M n 1.16 GeV. If this value is used instead of the empirical one, the fits to the susceptibilities worsen. We show in Fig. 11 a plot ofχ 2 /ν in the plane (m s , m cons ) using the nucleon mass as given by the model. One can see that in this case the best fits would correspond to non physical values of the current strange-quark mass, 150 MeV < ∼m s < ∼ 220 MeV. In our model, as in other quark models such as that of Isgur and collaborators [9], we have assumed that the string tension for the light quarks coincides with the one obtained from heavy quarks, and hence for the qD system we have taken σ = (0.42 GeV) 2 . There are models where this value is reduced by a factor of 2 when discussing baryon spectroscopy, as good spectra are obtained for (u, s, d) with σ qD = 2.15 fm −2 [22], for (u, d) σ qD = 1.57 fm −2 [50] and for (c, b) with σ qD = 4.5 fm −2 [27]. We note in passing that this factor of two re-scaling is also needed in the slope of radial Regge trajectories [51]. Motivated by these observations we have repeated our analysis taking σ = (0.42 GeV) 2 /2, with minor modifications but slightly worse fit quality.

VI. CONCLUSIONS
The missing resonance problem, i.e., the apparently overcounting of excited baryonic states by the quark model compared to the experimentally found resonances has been a long standing puzzle which has motivated a wealth of theoretical analysis and experimental work, mainly grounded in the individual identification of resonance states in the production process. This viewpoint demands a good knowledge of the scattering amplitude and its analytical properties in the complex energy plane. Indeed, resonances are uniquely characterized as process independent complex energy poles in unphysical sheets but the extrapolation from the physical axis, where measurements are actually made, to the complex plane is subjected to potentially large uncertainties due to the role played by the background.
The thermodynamic approach to the missing states problem for baryons has several advantages over the more conventional individual states analysis, since it addresses the completeness of states problem from the point of view of quarkhadron duality. It is rather insensitive to resonance energy pro-files as fine details of the the level density are washed by the Boltzmann factor. Lattice QCD has produced thermodynamic quantities, such as the trace anomaly, where with the currently rather small uncertainties one is not able to tell the difference between the current PDG spectrum and a quark model spectrum such as the RQM which was inferred already a few decades ago. Differences in the comparison become more visible when susceptibilities involving baryon number, electric charge and strangeness are considered as different subsets of states are selected.
For zero density, conserved charges such as B, Q, S have zero expectation values but fluctuate statistically in a hot vacuum. Those fluctuations are characterized by susceptibilities that have been determined numerically in QCD on the lattice by the HotQCD and WB collaborations to discriminate among hadronic models with sufficient accuracy, and thus can be used as a benchmark comparison.
We have argued that the asymptotic phase space for qqq systems is much larger than the one actually determined in RQM, which resembles a two body system with a linearly growing potential. This strongly suggests a dominance of quark-diquark dynamics for excited baryons. Therefore we have considered a quark-diquark model with a linearly confining interaction.
By analyzing the free energy of heavy quark sources characterized by Polyakov loops, we have been able to disclose under what conditions the poorly known quark-diquark string tension should coincide with the much familiar quarkantiquark string tension, in agreement with previous and recent lattice results.
Using this a priori fixed quark-diquark potential we have determined the remaining model parameters from conserved charges susceptibilities. The results are reasonable and fall in the bulk of previous intensive studies where a detailed description of the spectrum was pursued. where and and the sum is over mesons or baryons for ζ = ± respectively. For the baryonic susceptibility χ BB , which is the case we are going to consider in the following, ρ BB (M) is equal to the density of states ρ(M), as B = ±1 for (anti)baryons. The density of states can by computed in a derivative expansion [38], an approximation that is closely related to a semiclassical expansion in the high mass spectrum. In this approach it is best to start from the cumulative number whereĤ is the Hamiltonian and the trace is taken in the CM system subspace. The density is then obtained from ρ(M) = dN(M)/dM. In the quark-diquark model, the space of states is divided in sectors, λ , shown in Table I. Each sector contains a tower of multiplets all with degeneracy g λ displayed in the Table. This gives where the sum is over sectors and N (λ ) (M) sums over the tower of states including just one state in each multiplet. For the sake of clarity, in what follows we will drop the label λ . It is understood that the aggregated expressions are obtained by combining the results of the various sectors as in Eq. (A5). Within the semiclassical expansion, the cumulative number can be computed as where H is the (classical) Hamiltonian of the two-body system in the sector λ . The zeroth order term has been made

(B1)
S g (U) is the Euclidean gluonic action and (for simplicity here we use a notation of QCD in the continuous formulation), D a (U) = γ µ (∂ µ + iA µ ) + m a , where the gluon field A µ (x) is a hermitian matrix and m a is the mass of the flavor a. The Dirac matrices are hermitian.
As is well-known [52], the identity γ 5 D a (U)γ 5 = D a (U) † (B3) implies that the eigenvalues of the Dirac operator are either real or come in conjugated pairs, hence det D a (U) is real. Because u and d are degenerated, the weight det D u (U) det D d (U) is positive. The weight det D s (U) will be assumed to be positive too as required to be able to apply importance sampling Monte Carlo with dynamical quarks. This allows to define the real action so that Z = DU e −S g (U)−S q (U) .
The operator counting the flavor a is Q a = d 3 xψ a (x)γ 0 ψ a (x) a = u, d, s.
Because this quantity is conserved, we can use equivalently Q a = T d 4 xψ a (x)γ 0 ψ a (x) where T = 1/β is the temperature. Its expectation value can be obtained using Wick's theorem Of course this expectation value vanishes due to charge conjugation. Namely, using D a (U c ) = CD a (U) T C −1 with U c ≡ U * (or A c µ = −A T µ ), and Cγ T µ C −1 = −γ µ , it follows that the measure including the actions are even under U → U c , whereas Tr(γ 0 D −1 a (U)) is odd. For the ud correlation, again applying Wick contractions, This quantity is negative definite because Tr(γ 0 D −1 a (U)) is purely imaginary, as follows from (γ 0 D a (U)) † = −(γ 5 γ 0 )γ 0 D a (U)(γ 5 γ 0 ) −1 . (B10) One observation is that the fact that Tr(γ 0 D −1 a (U)) is purely imaginary provides another proof of Q a = 0 since this quantity is real, because Q a is hermitian in the real time formulation.
Another observation is that χ ud ≤ 0 will always hold in a Monte Carlo calculation, since Tr(γ 0 D −1 a (U)) 2 ≤ 0 for every configuration, even if the assumption det D s (U) > 0 were violated. On the other hand, for Q 2 u there are two Wick contractions Eq. (B10) again implies that the second term is real and the inequality χ uu ≥ |χ ud | implies that this second term is not only positive but at least twice as large (in average) as minus the first one. Since there is no reason to expect that Tr(γ 0 D −1 u (U)) 2 − Tr(γ 0 D −1 u (U)γ 0 D −1 u (U)) is definite positive for each gauge configuration, the condition χ uu ≥ |χ ud | could fail to hold for non positive definite det D s (U), thereby providing a test on this.