Thermodynamics of quark matter with multiquark clusters in an effective Beth-Uhlenbeck type approach

. We describe multiquark clusters in quark matter within a Beth-Uhlenbeck approach in a background gluon field coupled to the underlying chiral quark dynamics using the Polyakov gauge. An effective potential for the traced Polyakov loop is used to establish the center symmetry of the SU(3) color which suppresses colored states and its dynamical breaking as an aspect of the confinement/deconfinement transition. Quark confinement is modeled by a large quark mass in vacuum which is motivated by a confining density functional approach. A multiquark cluster containing n quarks and antiquarks is described as a binary composite of smaller subclusters n 1 and n 2 ( n 1 + n 2 = n ). It has a spectrum consisting of a bound state and a scattering state continuum. For the corresponding cluster-cluster phase shifts we discuss simple ans¨atze that capture the Mott dissociation of clusters as a function of temperature and chemical potential. We go beyond the simple ”step-up-step-down” model that ignores continuum correlations and introduce an improved model that includes them in a generic form. In order to explain the model, we restrict ourselves here to the cases where the cluster size is 1 ≤ n ≤ 6. A striking result is the suppression of the abundance of colored multiquark clusters at low temperatures by the coupling to the Polyakov loop and their importance for a quantitative description of lattice QCD thermodynamics at non-vanishing baryochemical potentials. An important ingredient are Polyakov-loop generalized distribution functions of n -quark clusters which are derived here for the first time. Within our approach we calculate thermodynamic properties such as baryon density and entropy. We demonstrate that the limits of a hadron resonance gas at low temperatures and O ( g 2 ) perturbative QCD at high temperatures are correctly reproduced. A comparison with lattice calculations shows that our model is able to give a unified, systematic approach to describe properties of the quark-gluon-hadron system


Introduction
The phase diagram of quantum chromodynamics (QCD) is a subject of intense theoretical and experimental research.Exploratory probes of the hot and increasingly dense nuclear matter are being performed in heavy-ion collision (HIC) experiments at the Relativistic Heavy-Ion Collider (RHIC) at BNL Brookhaven, the Large Hadron Collider (LHC) and the Super Proton Synchrotron (SPS) at CERN in Geneva with plans for many additional accelerator experiments in the near future.In the low-temperature range of the phase diagram, the possibility of a phase transition to a cold and dense quark-gluon plasma (QGP) in neutron stars is widely considered (see [1,2,3,4] for recent reviews).Many ongoing observational efforts aim at detecting signs of a phase transition inside neutron star cores.The current iteration of these multi-messenger observations include the NICER mission, and the LIGO-VIRGO collaboration with their gravitational wave detectors.The efforts of the NICER team have produced vital data on the masses and radii of pulsars PSR J0030+0451 [5,6] and PSR J0740+6620 [7,8] (currently the heaviest known neutron star with 2.08 ± 0.07 M ⊙ , as independent Shapirodelay based mass measurements [9]) while the analyses of the gravitational wave signal from the binary neutron star merger event GW170817 have put constraints on the tidal deformability of neutron stars in their typical mass range of 1.4 M ⊙ [10,11,12].Future measurements of postmerger gravitational waves may contribute to deciding the question whether a strong phase transition to deconfined quark matter may occur in those systems [13].In addition, the hypothesis of a phase transition induced shock revival in core-collapse supernova [14] might hint at an additional observational avenue of studying the structure of the hot and dense sector of the QCD phase diagram (cf.[4] and references therein).
All of these efforts require a good theoretical understanding of the phase diagram of strongly interacting matter and its structure.The only ab-initio solutions of QCD, which could be applied to the study of the phase diagram are lattice-QCD (LQCD) calculations, where Monte-Carlo simulations of the QCD partition function are performed on a discretized Euclidean space-time lattice.The vanishing baryochemical potential (µ B = 0) sector of the QCD phase diagram is well explored within this approach, see, for instance, Ref. [15,16,17,18].However, despite recent progress in calculating finite µ B observables on the lattice [19,20,21], the available density range is still insufficient to provide theoretical data applicable to both HIC and neutron star studies due to the numerical sign problem for finite µ B calculations.
Therefore, in order to explore the QCD phase diagram at arbitrary temperatures and densities one has to use effective models based on the underlying symmetries of QCD.One such model is the Nambu-Jona-Lasinio (NJL) model [22,23].It is extensively used in studies of dynamical chiral symmetry breaking and restoration (cf.[24]).However, this model does not capture the aspect of quark confinement and leads to an unphysical dominance of colored quark degrees of freedom already at relatively low temperatures, below 100 MeV.Therefore, an extension of the NJL model by an additional coupling between the fermionic degrees of freedom and a gluon background field in the Polyakov gauge (or PNJL model [25,26,27,28]) has been introduced which addresses the aspect of quark confinement by a statistical suppression of the contribution of free quarks to thermodynamic observables at low temperatures and chemical potentials.In this work we further extend the PNJL-based QGP model by a perturbative QCD correction of O(α s ) in the expansion w.r.t. the strong fine structure constant α s = g 2 /(4π) stemming from one-gluon exchange which, following Refs.[29] (for details, see also [30,31]) applies only to processes with high momentum transfer p > Λ pert .This extension of the PNJL model leads to an improvement of the comparison with LQCD thermodynamics at high temperatures.
The model is defined to allow for a thermodynamic description of matter at finite temperature T and baryochemical potential µ B under conditions of isospin symmetry and strangeness neutrality.
For a realistic description of strong interacting matter, the effective model should take into account quark hadronization in the domain of low T and µ B , where the resulting hadrons are composite quark clusters.In the spirit of [29], the composite nature of the low density hadronic degrees of freedom is accounted for by employing a generalized Beth-Uhlenbeck formula [30,31,32,33,34,35,36,37,38,39,40], which provides an interface between the low density and low temperature hadronic matter and the high density and high temperature PNJL-based QGP.This approach to quark clustering closely resembles the idea of the hadron-resonance-gas (HRG) model, which is crucial to the interpretation of LQCD calculations at low temperatures.
The HRG model assumes that the strong interactions are saturated by hadron resonance formation [41].Indeed, this assumption produces a description of low energy QCD matter as a free HRG well in agreement with both LQCD studies below the pseudocritical temperature T c [42,43,44,45] and HIC experiments (cf.[46] and reference therein).
Unsurprisingly, however, at increasing temperatures this agreement falters due to the rising importance of the quark substructure of hadrons leading to repulsive interactions, which the unmodified HRG model does not take into account.One fact is that highly excited resonances in the continuum of scattering states cannot be treated like bound states.These weak correlations are seen in the phase shifts and give only a small contribution to the virial expansion of the thermodynamic functions, as known, for instance from the Beth-Uhlenbeck formula discussed below.
Another fact is that the model of ideal mixture of noninteracting clusters fails when density is increasing.Various methods of solving this problem exist, for instance by a hadronic excluded volume.One physical mechanism for the repulsive interactions in a hadron gas is seen in the quark substructure of hadrons which, following the Pauli exclusion principle, leads to a repulsive interaction by quark Pauli blocking among hadrons [47,48].These effects of exchange symmetry are precursory to the delocalization of the hadron wave functions in the hot, dense hadron gas which entails the dissociation of hadrons in the analogue of a Mott effect introduced in nuclear physics [49] and known from solid state and plasma physics [50].In order to provide a description of the QGP on the same microscopic footing with the gas of hadronic resonances as bound states of a confining interquark potential, the principle of saturation of the color charge within a range of nearest neighbors has been applied within the string-flip model (SFM) [51,47].In Ref. [52], an improvement of the SFM has been presented which includes relativistic kinematics for the quarks and discusses a mixed phase between the SFM plasma and a schematic HRG model containing multi-quark clusters.The results were compared with the early lattice QCD simulations of that time.However, a consistent formulation of medium effects on bound state formation is absent in this approach so that this mixed phase construction can only be viewed as an early realization of the switch function concept of Ref. [53].The SFM model for the QGP has been considerably improved in the subsequent work of Khvorostukhin et al. [54], where a confining quasiparticle model of quarks and gluons has been introduced in the form of a relativistic density functional (RDF) approach.The form of this density functional has later been generalized to include the coupling to a nonlinear vector mean field that allowed to discuss heavy hybrid neutron stars and mass twin solutions for them [55].The most recent development of the SFM is a lagrangian RDF formulation that obeys chiral symmetry and contains vector meson and diquark mean fields and is shown to obey modern constraints from multi-messenger observations of neutron stars and their mergers [56].
However, to provide qualitative insight into the transition from an HRG to a QGP, in the present work we abstain from a detailed microscopic model for solving the few-quark problem in a dense medium.For our study, we will incorporate the effect of the repulsive interaction as well as the eventual temperature and density dissociation of hadrons (the Mott effect) directly into a model for the phase shift δ i of a hadron species i in the generalized Beth-Uhlenbeck formula.The ansatz will be similar to the one used in [30,31], in agreement with the Levinson theorem.The resulting contribution to the thermodynamic observables will be referred to as the Mott-hadron-resonancegas (MHRG) contribution.The full MHRG+PNJL model embodies the main consequence of quark chiral symmetry restoration, i.e. the lowering of the continuum thresholds for multiquark scattering states, which entails the dissociation of multiquark bound states (Mott effect).A general approach to strong correlations in many-particle systems is the Φ-derivable approach which is based on the Luttinger-Ward-Baym functional [57,58,59] and has recently been applied to a T -matrix description of the QGP thermodynamics as well as spectral and transport properties by Liu and Rapp [60].However, these authors did not extend the approach to include bound state formation and the limiting case of the hadron resonance gas phase.This is the main goal of the present work, where we generalize the Φ-derivable approach to deal with multi-quark clusters.
As an extension of the previous work [29], we have considered clusters beyond the color-singlet states, and we allow for any combination of quarks and antiquarks up to 6 valence degrees of freedom.In this way, we also take into account all colored multi-quark cluster states which, however, are strongly suppressed in the region below the Mott temperature due to their coupling to the Polyakov loop.
As a benchmark for defining the model, we use continuum-extrapolated LQCD calculations of the baryon density and the entropy density as a function of the temperature for finite values of µ B /T [19,20].These data allow us to test the PNJL model with the MHRG virial correction terms.We perform exploratory calculations accounting for a limited number of resonances from the Particle Data Group (PDG) up to a mass threshold of about 2600 MeV.
Furthermore, instead of deriving the quark mass selfconsistently, we make the assumption of a sudden drop of the quark mass at the µ B -dependent pseudocritical temperature T c (µ B ) from the constituent quark mass value to the current quark mass value.This assumption entails that the pseudocritical temperature can be identified with the Mott temperature T Mott (µ B ) where all hadronic bound states dissociate into their quark constituents.This allows us to focus on the cluster contribution to the Polyakovloop calculation, without taking into account the complicated interplay between the singlet clusters, the scalar mean field and the Polyakov loop.Such an analysis is definitely warranted, especially due to the discrepancy between the pseudocritical temperature of the chiral crossover transition obtained in LQCD simulations and Polyakovloop generalized chiral quark models.But this issue is beyond the scope of the present work.
This work is organized as follows: In Sect. 2 we present the theoretical aspects of our approach, i.e., the QGP thermodynamics based on the PNJL thermodynamic potential, and the treatment of clusters in the context of the MHRG approach.A major issue is the determination of the scattering phase shifts, which follow in principle from the solution of a few-quark in-medium wave equation, but here are postulated in a generic, medium-dependent form in Sect. 3 along with other parameters of the model.Results for the thermodynamic quantities baryon density and entropy as a function of T and µ B are presented and discussed in Sect. 4. In Sect.5, a summary is given and conclusions are drawn.Technical details are collected in Appendices A and B.

Generalized Beth-Uhlenbeck approach to quark clustering in the Polyakov gauge
In order to satisfactory reproduce the high temperature lattice QCD predictions, we shift now towards an effective Polyakov-loop quark-gluon plasma (QGP) model based on the quark quasiparticle description combined with the Polyakov-loop potential U(ϕ, φ) and a O(α s ) perturbative virial correction.In an effort to maintain the description of hadrons as nonperturbative correlations of quarks and gluons, we couple this model to the Mott-hadron resonance gas (MHRG), in analogy to [29], and extended it to finite chemical potentials.
In this work we develop the technique of the cluster virial expansion approach for multi-quark clusters at finite temperatures and chemical potentials on the basis of a cluster generalization of the so-called Φ-derivable approach.Within a Green's function approach [58,59] the correlation functional Φ is introduced where the contributions of interaction are represented by closed-loop diagrams.When for the Φ functional a restriction to closed two-loop diagrams in cluster Green's functions is applied, this approach is equivalent to the generalized Beth-Uhlenbeck approach to clustering in hot, dense Fermi-systems [61].The intricacy of color confinement in low-density quark matter is considered by coupling the quarks and their clusters to the Polyakov-loop background field which serves to suppress the appearance of colored clusters in the region of quark confinement.We demonstrate a satisfactory comparison of our results for the thermodynamics of clustered quark matter with recent lattice QCD simulations at finite temperature and chemical potentials, where they are available.
(1) In the following subsections we explain these two contributions more in detail.

QGP thermodynamics of quarks and gluons
The thermodynamics of the deconfined QGP can be separated into a perturbative and a nonperturbative contribution where Ω NP describes the nonperturbative low-energy QCD quark and gluon degrees of freedom via the mean-field thermodynamic potential and perturbative corrections are absorbed in Ω pert (T, µ, ϕ, φ).We take into account the minimal coupling of quarks to a homogeneous gluon background mean field in the Polyakov gauge A 4 = λ 3 A 3 4 + √ 3λ 8 A 8 4 while quark interactions are accounted for by a suitably chosen relativistic energy density functional L int with U(ϕ, φ) denoting the Polyakov-loop (gluon) potential for which we will not make use of the traditional form found in [26], but the one given in [62] as explained below.
The traced Polyakov loop is defined as where N c = 3 is the number of colors in the SU (N c ) gauge theory and β = 1/T is the inverse temperature.Abbreviating this expression in accordance with the standard notation found in the literature, we have A 3 4 = φ 3 and (6) The minimal coupling of the quark to the homogeneous gluon background field in the Polyakov gauge leads to a shift of the Matsubara frequency iω n → iω n − iA 4 .We define the distribution function of a quark of given color c as the Matsubara sum over the quark propagator The Polyakov-loop distribution function is then obtained after performing the color trace where f [(E p − µ)] is the ordinary Fermi function which can be written as with z + ≡ ln[1 + y + 1 ] and y + 1 = e −β(Ep−µ) .Now we can write and obtain the Polyakov-loop generalized single-quark distribution function We note that the distribution function for an antiquark is obtained from Eq. ( 13) by the replacements: µ → −µ and ϕ → φ = ϕ * .The nonperturbative part of the QGP thermodynamics (3) consists of a quark quasiparticle contribution and the Polyakov-loop potential U for which we adopt the form defined in [62], with M H being the SU(3) Haar measure and The parameters a 1 . . .a 5 , b 1 . . .b 4 , c 1 . . .c 5 , d 1 . . .d 5 are taken from the pure SU(3) lattice gauge fit performed in [62].The parameter a 6 = 1.28 is introduced here to improve the description of the baryon density in comparison with LQCD results [19].The parameter T 0 is related to the critical temperature for deconfinement in the case of pure gluodynamics, where the value T 0 = 270 MeV is taken from pure gauge theory simulations on the lattice.One has to invoke a flavor dependence of this parameter [63], which for our applications to the realistic case of N f = 2 + 1 flavors is chosen to be T 0 = 205 MeV in order to provide the best agreement with the LQCD data [19] discussed below.
The mean field thermodynamics in the quark sector of the PNJL model is described by Here, G s is the scalar coupling constant, and ∆ f (T, µ) = M f (T, µ) − m f is the dynamically generated quark mass gap for the flavor f , see Tab. 1.The same structure for the mean field thermodynamic potential (20) is also obtained in the confining relativistic density functional (RDF) approach, see Eqs. ( 16) and ( 17) of [56], since it can be mapped to an NJL model with medium dependent couplings.Solutions of the gap equation for the medium dependent quark mass in this case are shown, e.g., in [64] and exhibit a significantly larger constituent quark mass than NJL models, realising the effective quark confinement in the RDF approach.In this work, we will use a sudden switch model that entails a corresponding behavior of the thresholds for multi-quark continuum states as discussed in more detail in subsection 3.1 below.
The divergence of the integral within Ω Q is a known deficiency of the (P)NJL model.We regularize the divergent part of this integral by introducing a 3-dimensional momentum cutoff.The cutoff Λ is chosen such, that simplifying the quasiparticle expression for the quark thermodynamical potential to In addition, we introduce a perturbative virial correction contribution in two loop order, Ω pert (T, µ, ϕ, φ) with a standard expression (cf.[65], see also [66]) of the form where are modified integrals introduced in [29] and extended to finite chemical potentials with the generalized Polyakovloop generalized Fermi distribution functions.The infrared cutoff, Λ pert = 222 MeV, represents the momentum range below which nonperturbative physics dominates.This value is adjusted in order to provide agreement with the data on the LQCD thermodynamics (see Section 4).We use a temperature and chemical potential dependent regularized running coupling [67,68,69] with N c = 3, N f = 3, T * = 93.75MeV.We choose µ * = 3πT * in accordance with the argument of the switch function in [53].
The traced Polyakov loop variables ϕ and φ are obtained by solving the corresponding gap equations derived by minimizing the total thermodynamic potential.This includes the PNJL and perturbative parts of the QGP model, which were derived in [29] for µ = 0.As an extension of that work, we have introduced finite-µ descriptions of each of the contributions to the thermodynamic potential.Additionally, in the MHRG part we have permitted the contributions of color-triplet and color-antitriplet clusters, which, while strongly suppressed in the HRG phase where the traced Polyakov-loop is close to zero, have an impact on the Polyakov-loop behaviour and on the thermodynamics in the phase where the approximate chiral symmetry is restored and bound states underwent Mott dissociation.The results of these calculations will be shown in the following subsection.

Clusters in quark matter and the MHRG approach
Below the deconfinement transition, quark matter appears predominantly as color-neutral clusters, the hadrons.The quark substructure of the hadronic system determines the interaction between hadrons, but also intrinsic excitations.The inclusion of cluster formation, including excited states and resonances is an indispensable ingredient to the thermodynamics of quark-gluon matter.
In contrast to plasma physics where we start from a known hamiltonian and use for the evaluation of the partition function well-elaborated methods such as partial summation of ladder diagrams to introduce bound state formation, this first principle approach is presently not in reach in QCD.Data from Lattice gauge theory simulations are not available for nuclear matter in the hadronic phase.Usually, this problem is circumvented by interpolation between two different approaches to describe high energy density matter, the QGP at high temperatures as described above, and the hadron-resonance gas (HRG) below the deconfinement transition.For instance, in [53] both approaches are matched together, using a switching function for a smooth transition.There, the crossover region at µ = 0 occurs around T = 170 MeV and goes to T = 0 around µ = 1.25 GeV.In addition, the HRG is modified by an excluded-volume ansatz to mimic the interaction between hadrons in the hadronic phase.Adapting parameters for this semi-empirical approach, in Ref. [53] good agreement with lattice data is obtained.The introduction of a switching function "by hand", however, does not provide any insights to the microphysical processes that underlie the hadron-to-QGP transition.
We intend to find a unified description for the quarkgluon-hadron matter system including cluster formation and dissociation.The hadrons are obtained as bound states, but their disappearance at high densities is described as melting owing to medium effects.Instead of a phenomenological excluded volume concept, a quark exchange process may be considered as origin of repulsion.Thus, a main point is the account for antisymmetrization of the fermionic quark states what leads to Pauli blocking and the dissolution of bound states, denoted above as Mott effect.In the absence of a fundamental hamiltonian of the system, we will use multiquark correlation functions to encode and model the dynamical properties of the system.
In this work we are interested in the systematic description of correlations in the quark system which are seen from the bound state formation, but also in the continuum correlations.For the two-particle system, continuum correlations are related to the scattering phase shifts as known from the Beth-Uhlenbeck formula.A generalization of this Beth-Uhlenbeck formula has been discussed [61] where the scattering of two clusters was considered.
When considering more than two elementary constituents of a cluster, different channels are possible where this cluster appears.Here, we deal with binary processes for the formation and the decay of the clusters into two subclusters.We consider two-component scattering phase shifts which indicate the occurrence of a resonance, but reflect also the possible existence of a bound state of these constituents.
In general, we have different possibilities for the decomposition of a cluster into subclusters.This is of rel- evance for the decays into various channels with corresponding branching ratios, but not for the thermodynamic properties which are determined by the density of states.
For the thermodynamic properties it is sufficient to consider only a specific binary decomposition.This approach is a generalization of the two-particle problem to a manyparticle approach and has also been applied in nuclear physics to describe cluster formation and resonances in nuclear matter [70,71].Within our approach, we consider all possible clusters formed from light quarks of flavors f = u, d, s, including also colored clusters.This is necessary to describe the transition to the QGP state.The set of clusters considered in the present work is shown in Tabs.2-4.The colorneutral clusters are selected from the Particle Data Group [72] summary tables for mesons in Tab. 3 and for baryons in Tab. 4, according to the choice made in Ref. [53].We added the deuteron as a hexaquark state as well as the hypothetical sexaquark S(uuddss), for which a mass of 1885 MeV was motivated by Buccella [73] assuming a compact three-diquark structure, see also [74] and references therein.The contribution of heavier clusters is discussed below.It should be mentioned that also further clusters such as penta-quarks are discussed recently [75,76], but those containing only three light quark flavors are not yet confirmed by experiments.
The MHRG part takes the form of a cluster decomposition of the thermodynamic potential for quark matter where n denotes the total number of valence quarks and antiquarks in the cluster; c n = 1/2 for bosonic and c n = −1 for fermionic clusters [77,78].The functional Φ [{S n }] contains all two-cluster irreducible (2CI) closed-loop diagrams that can be formed with the complete set of cluster Green's functions S n .We will restrict ourselves to a maximum number of N = 6 quarks in the cluster and to the class of two-loop diagrams of the "sunset" type which are shown in Fig. 1 and Fig. 2 for the case of N = 6.[72].Also shown are the masses Mi according to our simple hadron mass formula (52) and the masses of the continuum thresholds M th,i for the decay into their quark constituents.
States denoted by * * ) are considered as pentaquark states; * * * ) the deuteron is a six-quark state with molecular structure of two almost unbound nucleons; * * * * ) the sexaquark is a hypothetical deeply bound six-quark state [74].
The full quark cluster propagators S n fulfill a Dyson-Schwinger equation where n is the free n−quark cluster propagator and the cluster selfenergy Π n should be obtained from The free n−quark cluster propagator for n = 2, . . ., N in general is obtained from the free quark propagator by a product ansatz and subsequent summation over n − 1 Matsubara frequencies in order to obtain a one-frequency function, see, e.g., Ref. [79].Here we perform a cluster virial expansion which we restrict to the second virial coefficients for interactions of pairs of clusters, including also quark pairs as the simplest example.In that case we obtain the free a−quark cluster propagator by performing one Matsubara summation in the product ansatz for the bipartition of the cluster as illustrated in Appendix A. The generalized Matsubara frequencies which include the coupling to the Polyakov-gauge gluon background field, are given in Tab. 5.
The equations ( 27), ( 28), ( 29) and (30) form a closed set of equations that for its solution requires the knowledge of the free n−quark cluster propagator S (0) n and the choice of the 2CI set of diagrams of the Φ−functional, built with the set of full cluster propagators {S n }.Following the spirit of Refs.[77,78], we consider here the quark density and entropy density as first derivatives of the thermodynamic potential with respect to the quark chemical potential µ and the temperature T , respectively.From the cluster decomposition ( 27) of the MHRG thermodynamic potential follows the cluster decomposition of the quark density where a i is the net quark number in the n i −particle state with the partial density defined as and d i denotes the degeneracy factor for n i −particle states.Using the spectral decomposition rule for an analytic com- plex function F (ω) and the Matsubara summation in the case of many-quark Green's functions coupled to the gluon background field (see Tab. 5 and Appendix A.2) with the relation ∂f (a) (ω)/∂µ = −∂f (a) (ω)/∂ω, we get for Eq.(32) now  where a partial integration over ω has been performed.For two-loop diagrams of the sunset type holds a cancellation [77,78], which we generalize here for cluster states Using generalized optical theorems [80] and the polar representation of the n-quark cluster Green's function , where δ n is the cluster phase shift, one obtains The density is obtained in the form of a generalized Beth-Uhlenbeck EoS where the properties of the distribution function f where y ± a = e −(ω∓aµ)/T and a is the net number of valence quarks present in the cluster.See Appendix A for a detailed derivation.
In an analogous manner follows for the MHRG entropy density where ϕ ] and f (a) ϕ is the cluster distribution function for a net quark number a modified by the traced Polyakov loop.
The Eqs. ( 38) and ( 41) are consistent with each other because they follow from the functional for the thermodynamic potential (27) in the case of two-loop skeleton diagrams for the Φ-functional [77,78], where it has been shown that the correlation contribution vanishes, see Eq. (36).The formula for the pressure as thermodynamical potential can be obtained from Eq. ( 38) by integration over the quark chemical potential µ.Analogously, it can be obtained from Eq. ( 41) by integration over T p(T, µ) It is instructive to consider the two limits ϕ = φ = 1 (deconfinement and color singlet states) and ϕ = φ = 0 (confinement).In the situation where ϕ = φ = 1, the distribution functions reduce to corresponding to the ordinary Bose (Fermi) distribution functions for even (odd) numbers a of quarks in the cluster.In the case of confinement when ϕ = φ = 0, the distribution functions take the form From Eq. ( 45) one can see how the Polyakov-loop suppression of colored states acts.Since y ± a 3 (T ) = y ± a (T /3), in the situation of confinement when ϕ = 0, the generalized distribution functions of odd (even) colored clusters become Fermi (Bose) functions with the temperature T being replaced by T /3, so that in this case of full Polyakov-loop suppression the thermal occupation is the same as for the unsuppressed case at one third of the temperature.
We bring expressions (38) in another form where we split the extended phase shift δ i (ω, q) into its bound state part with a step-up at ω = E i (q) and the part related to the elastic scattering phase shift δ sc i (ω, q) depending on the energy ω of relative motion of the subclusters and the total momentum q of the cluster.The net baryon density (38) is the difference of contributions from particles and antiparticles, Applying integration by parts to the Beth-Uhlenbeck equation and considering only the color singlet hadrons for which f ϕ=1 , the particle contribution can be brought to the form and a similar relation for the antiparticles.Here, the bound state part can be addressed to the first two contribution, similar to the Planck-Larkin-Brillouin partition function for the hydrogen plasma [50].Within the generalized cluster Beth-Uhlenbeck formula, we have for the a-quark correlation the continuum contribution as expressed by the phase shifts, and, in addition, the contribution of bound states, if they exist.For the entropy density (41), a similar decomposition as ( 47) and ( 48) can be made.As well known, the subdivision of contribution from bound states and from the continuum of scattering states is arbitrary.Because of the Levinson theorem we can introduce an extended scattering phase shift which jumps by π for every bound state at the bound state energy which is below the continuum edge of scattering states (so-called "step-up" phase shift).

Exploratory model calculation for multiquark correlations 3.1 Quark masses and continuum thresholds
The temperature and chemical potential dependent quark masses determine the behavior of the continuum thresholds.We want to point out that in non-confining NJL models a smooth temperature dependence of quark masses has been obtained which results from the consideration of a medium of uncorrelated quarks in these models.Because in reality the quarks are strongly correlated (confined) inside hadrons even for temperatures up to the chiral crossover T c (µ), any medium modification of quark masses shall be strongly suppressed.[81] later extended within the Polyakov-loop DSE model [82,83], and recent results within the confining density functional approach [56] motivate a steep medium dependence (switch) from a large constituent quark mass in cool and dilute matter to the current quark mass of the QCD Lagrangian in a hot and dense medium.

Improved microscopic calculations within a rank-2 separable Dyson-Schwinger equation (DSE) model
Therefore, we assume in the present work a "sudden switch" scenario for the quark masses.It consists of a constant and large quark mass as a reflection of the quark confinement within the hadronic phase of the model that is delimited by the Mott dissociation temperature of the hadrons (and multi-quark states) T Mott (µ), and a tail function for the medium-dependent quark mass in the QGP phase, which for simplicity we assume here to coincide with the constant current quark masses m f , Here the subscript f = q, s denotes the light (q = u, d) and strange (f = s) quark flavors with their corresponding current quark masses m f and the constituent quark masses M f = M f (0, 0) in the vacuum.This quark mass model is depicted for light quarks as the red solid line in Fig. 3.We are aware that such a "sudden switch" model is a strong idealization of the reality and could lead to a kinky behaviour of the resulting thermodynamic functions, in particular when they are derivatives of the thermodynamic potential.
Consistent with the quark confinement property we make in our model the idealizing assumption of constant hadron masses.This captures, however, the fact that hadrons are almost pointlike, strongly localized states which "feel" the medium only when it provides conditions very close to the Mott dissociation.In this case, hadron-hadron interactions induce finite mass shifts and collisional broadening.
We extend the formula (49) to finite µ B by using the chemical potential dependence the pseudo-critical temperature T c (µ) which resulted from the fitting the lattice QCD results obtained by Taylor expansion techniques for small chemical potentials [18] also for estimating the Mott temperature, We use the values κ = 0.012 and T Mott (0) = 156.5 MeV from the fit of the µ-dependence of the pseudocritical temperature T c (µ) = T Mott (µ) found in [18], restricting ourselves to the O(µ 2 B ) part of the expansion and neglecting the O(µ 4 B ) term because of the smallness of the associated κ 4 parameter found in that work.
With the medium dependence of the quark mass being defined, one can solve the gap equation for obtaining the value of the traced Polyakov loop which would correspond to an extremum of the thermodynamic potential (1), This can be done first at the quark mean field level, including the perturbative QCD contribution.Then, the effect of colored clusters can be included, which are depending on ϕ and φ via the Polyakov-loop generalized multiquark distribution functions ( 76) and (75).In Fig. 4, we show the temperature dependence of the absolute value of the Polyakov loop for different values of chemical potentials µ B /T .These results include the backreaction from colored clusters which turn out to be negligibly small.For temperatures below T Mott (µ), because of the large quark mass, the values of the traced Polyakov loop stay close to zero, reflecting confinement.Also in the transition region above the Mott temperature, up to T ∼ 200 MeV, this order parameter for deconfinement stays below |ϕ| ∼ 0.5, indicating strong color correlations induced by the gluon background field.It can also be seen from Fig. 4 that for temperatures above T ∼ 350 MeV, the order parameter asymptotically approaches |ϕ| = 1.0, the case of deconfinement.

Multiquark clusters
In the model exploited in the present work, we do not obtain the hadron masses as solutions of the corresponding equations of motion.Instead, we assume that below the Mott dissociation temperature T Mott (µ) the hadron masses M i are medium independent and follow the approximate relation [84] where M q(s) = 627 (770) MeV is the light (strange) quark mass in vacuum and B = 471 MeV is the binding energy per quark bond, i.e. per relative Jacobi coordinate of which we have n − 1 in a multiquark state composed of N q (N s ) light (strange) quarks and N q (N s) light (strange) antiquarks.Here, n = N q + N q + N s + N s is the total number of quarks and antiquarks and a = N q + N s − N q − N s is the net quark number in the multiquark state (hadron).
From the medium-dependent quark mass spectrum follows immediately the temperature and chemical potential dependent threshold mass of the continuum states for a multiquark bound state species i This behavior of the threshold masses entails that all PDG hadrons undergo the Mott transition from bound states below T Mott (µ) to unbound correlations in the continuum above that temperature, see Tabs. 3 and 4.
For temperatures above T Mott , a linear temperature dependence of the mass and the decay width of the resonance is adopted where M i,0 stands for the vacuum mass of the resonance i that is taken from the particle data tables [72] and remains unchanged for all temperatures below T Mott .Since the Mott temperature is identified with the pseudocritical temperature T c obtained from lattice QCD simulations [18], the ansatz ( 56), ( 57) leaves us with only two free parameters, α 1 = 11.4 and α 2 = 1.9, which determine the slopes of masses and widths, respectively.Their values are determined such that the behaviour of the total entropy density and the total baryon density at T c are continuous.This behavior is illustrated in Fig. 3 for the nucleon and the pion.For the higher-lying meson (baryon) states which all decay into their ground state and a meson (mostly pion), we adopt a mass formula that assumes a tetraquark (pentaquark) structure, i.e. adding 2M q − 2B = 312 MeV to the ground state mass.The continuum threshold for these states is then 2M q = 1254 MeV (2m q = 11.2MeV) above that for the corresponding ground state for T < T Mott (T > T Mott ).While in our calculations the color singlet hadron masses are taken from the PDG [72], the masses of the colored clusters and the continuum thresholds require an underlying quark model, see tables 2 -4.

η mesons
For the ground state η mesons, we postulate a mass formula composed as a weighted sum of light and strange quark masses with the average mesonic binding energy with the corresponding threshold masses The higher lying η states are understood as tetraquarks, i.e. adding a light quark-antiquark pair as in the case of the higher lying meson and baryon states.

Phase shift model
In this work, we want to introduce a generic model ansatz for the phase shift δ i (M, T, µ) that shares all essential features with results that have been obtained, e.g., within NJL model calculations for the case of the pion [85,36] and within the PNJL model for the pion and the diquark [86] in quark matter.
For the extended phase shifts, which include the bound state contribution as "step-up" by π, we can use experimental determined ones, including the binding energies, or results from a potential model calculation.
Three simple approximations which are illustrated in the three panels of Fig. 5 shall be considered: 1.The phase shifts jump up from zero to π at the bound state mass and then remain constant.This simplest step-up model is labelled "SU" and shown in the upper panel of Fig. 5.For the full set of hadronic states in the tables 3, 3 extracted from the listing of the PDG [72], this model reproduces the HRG thermodynamics, see the dotted lines in Fig. 6.The SU model violates the Levinson theorem and does not account for the dissociation of hadrons at high temperatures and densities.It overestimates the contribution of clusters.2. The simplest way to account not only for the formation of bound states but also for their dissociation by simultaneously fulfilling the Levinson theorem is to add a step down by π at E thr,i to the step-up of the phase shift by π at the bound state energy E i .In this approximation, labelled "SUSD" in the middle panel of Fig. 5, the ω-integral in the Beth-Uhlenbeck equation (48) for the hadron contribution to the net baryon density can be trivially performed and results in the absence of this contribution for temperatures T > T Mott .By the definition of this SUSD model, resonances in the continuum are absent.In this case the entropy of the MHRG model which follows from Eq. ( 41) is shown by the red dashed lines in Fig. 6. 3. Include the possibility to describe resonances in the continuum using a Breit-Wigner ansatz and postulate a structureless negative background phase shift which assures the Levinson theorem, namely that at asymptotic large energies the total phase shift in each channel goes to zero.The resonance masses and width as parameters of the Breit-Wigner ansatz are temperature dependent, a simple linear dependence is assumed.In this case, labelled "SUC" and shown in the bottom panel of Fig. 5, there is a finite, nonvanishing contribution to the thermodynamics from hadronic resonances in the continuum for temperatures above T Mott , see the green solid lines in Fig. 6.
The simplest model that is in agreement with the Levinson theorem is a SUSD ansatz [87,80].It results in the Mott dissociation of hadrons at high temperatures and chemical potentials Here, the only ingredients are the medium-dependent hadron masses M i and the sum of the masses of their quark constituents which determine the threshold mass M thr,i of the corresponding continuum of unbound states.The main drawback of the simple ansatz of Eq. ( 61) is that it neglects any contribution from the continuum of scattering states.This ansatz has recently been applied in the description of lattice QCD thermodynamics at µ = 0 [29].In the absence of continuum correlations, the requirement of a smooth transition region between the pure HRG and pure QGP phases for the temperature interval 160 ≲ T [MeV] ≲ 200 was fulfilled by adopting a smooth drop of the quark masses, viz. the continuum threshold over that region of temperatures.This resulted in a sequential dissociation of hadrons with the consequence of a staircase-like fine structure of the thermodynamic potential which entailed a fuzzy behavior of its derivatives.Besides this caveat, if the Mott dissociation is related to the chemical freeze-out, this scenario would entail a statedependent smearing of freeze-out parameters for which there are no indications in heavy-ion collision experiments.
For the present work we take the ansatz of [34] and exit from the original π and σ channels to include contributions from each cluster with the appropriate Mott dissociation effect at high temperatures.This modified phase shift model reads For temperatures above T Mott , a linear temperature dependence of the mass and the decay width of the resonance is adopted according to Eqs. ( 56) and (57), respectively.
The encoded mass dependence of the phase shift is shown in the bottom panel of Fig. 5 for the generic case of the pion at rest in the medium (where ω = M ) for different temperatures around the Mott transition at T Mott = T c = 156.5 MeV.Fig. 6.Scaled entropy density of hadrons sMHRG/T 3 as a function of temperature T calculated for µB/T = 0 (upper panel) and µB/T = 3.0 (lower panel) calculated for the SU (blue dotted curves), SUSD (red dashed curves) and SUC (green solid curves) models of the hadronic phase shifts discussed in the text.The shaded regions correspond to the lattice QCD results [19], while the thick solid line represents calculations within the HRG, see [53].

Results for QCD thermodynamics
After having defined the Beth-Uhlenbeck type approach to hadron-quark-gluon matter in Sect. 2 and making plausible assumptions about the key ingredients for its numer- ical evaluation in Sect.3, we are now in the position to present results for the entropy density and the net baryon density.In Fig. 7, we present the scaled entropy s/T 3 as function of the temperature.In the upper panel, we show the results for vanishing baryochemical potential, while in the lower panel the case µ B /T = 3.0 is presented and presented to the results of lattice QCD simulations [19], which are given as a grey band.There is a very good agreement.On this basis, we can now give an interpretation of the QCD thermodynamics.Within the sudden drop model for the quark mass, there is a perfect confinement of quark and gluon degrees of freedom for temperatures below the pseudocritical one, which in our case coincides with the Mott temperature of hadron dissociation T Mott = T c = 156.5 MeV.For T < T Mott , the entropy is saturated by the contributions from the ideal hadron resonance gas.We have implemented here the same set of hadronic states from the PDG as in Ref. [53], and added the deuteron as a hexaquark state as well as the hypothetical sexaquark S(uuddss) which is discussed, e.g., in [74] and references therein.At T = T Mott , the hadron contribution drops by about half of its value as a result of the temperature dependence of the hadronic correlations encoded in their phase shift adopted in this model.The full spectrum of hadronic states contributes to the entropy as continuum correlations modeled by their phase shifts.With increasing temperature, these contributions die out and vanish above T ∼ 200 MeV.We note that it is important for obtaining a continuous behaviour of the entropy at the transition temperature that the contributions in the deconfined phase, stemming from quark quasiparticles, gluons, perturbative virial corrections and colored clusters conspire so as to compensate the sudden drop in the hadronic entropy.One may attribute this to quark-hadron duality.A remarkable finding is that at finite chemical potential the colored clusters (diquarks, 4quark and 5-quark states) play an important role.In Fig. 8, we demonstrate that the chemical potential dependence of the lattice QCD data is correctly reproduced and that for large temperatures the O(α s ) perturbative QCD limits are attained.
In Fig. 9, we demonstrate the very power of the present approach.Namely, that we can identify the composition of the matter on both sides of the transition.The fractional contributions of hadron species reflect the composition of the HRG within the statistical model which gives an excellent understanding of hadron production in heavy-ion collision experiments (see, e.g., [46] and references therein).We want to underline that the red and turquoise solid lines in Fig. 9 standing for the quark and gluon fractions of the entropy, respectively, shows the degree of realization of quark and gluon confinement in this model.While quarks and gluons dominate the high temperatuure region above T ∼ 200 MeV, their entropy fraction in the hadronic phase below T Mott is below the per cent level.This realization of sudden quark-gluon confinement is due to the combination of the Polyakov-loop suppression of the colored state and the large constituent quark mass in the hadronic phase.
We want to comment also on the changes that the entropy fractions undergo at the Mott transition.The contributions of all multiquark clusters get suppressed by about a factor two because of the dramatic change in the generic phase shift function, when they change their character from a bound state to a correlation in the continuum, see Fig. 5.This change can be directly read-off from Fig. 9 for color singlet hadrons (pions, kaons, nucleons, ...) while for colored clusters there is the effect of Polyakov loop suppression which acts on their distribution functions in the hadronic phase, where the traced Polyakov loop is close to zero, see Fig. 4. In the case of diquarks, both effects compensate each other, so that the entropy fraction of diquarks is continuous at the Mott transition, while for the 4-quark and 5-quark states, the Polyakov loop suppression overcompensates the effect of Mott dissociation at T Mott .The S(1885) contribution is a unique prediction of the present model.We note that in Ref. [89], predictions have been made for the relative yields of S(1950)/d and S(1700)/d in heavy-ion collisions using a thermal statistical model with excluded volume corrections.Similar to the entropy fraction shown in Fig. 9, the sexaquark-todeuteron ratio is of the order one.A quantity of central interest for the QCD thermodynamics is the pressure because it plays the role of the thermodynamic potential of the grand canonical ensemble.We use Eq. ( 42) to obtain the pressure by integrating the entropy density over temperature.The result is shown in Fig. 10.The comparison to two sets of lattice QCD data by Borsanyi et al. (2014) [88] and Borsanyi et al. (2021) [19] shows an excellent agreement.
Finally, in Figs.11 and 12 we show the results for the scaled baryon density as a function of the temperature.In Fig. 11 we show the different components that contribute to the baryon density for µ B /T = 3.0 and compare the resulting total density with the lattice QCD simulation [19] shown as a blue band.We find it remarkable that the colored clusters contribute more than the color singlet hadrons in the transition region above T Mott .Fig. 12 demonstrates that the systematics of the chemical potential dependence is correctly accounted for.In particular, the baryon density quickly approaches the O(α s ) perturbative QCD values while the lattice QCD data slightly overshoot.At temperatures below T Mott , where the baryon density is given by the color singlet baryons only, there is a discrepancy with the lattice QCD results which eventually can be related to the known problem of missing strange baryon states in the PDG HRG description [17].Adding baryonic states from a quark model description [90] has been demonstrated to solve this problem [91].

Summary and Conclusions
We have developed a systematic approach to the thermodynamics of the quark-gluon-hadron system that studies quark-quark correlations in the form of bound states and continuum correlations.The challenge is to describe, in contrast to conventional plasma systems, the confinement property that allows only color-neutral clusters in the lowdensity limit.This property was taken into account by using the Polyakov-Nambu-Jona-Lasinio model as an ef- fective approach for the many-quark system, where the gluon degrees of freedom were eliminated.We obtained results for arbitrary temperatures and for finite baryonic chemical potentials.Two limiting cases are reproduced: hadronic matter described as hadron resonance gas at low temperatures, for which experimental data are available, and the quark-gluon plasma at high temperatures, for which results from perturbation theory are known.In contrast to semi-empirical approaches such as Ref. [53], which use a switch function to interpolate between the two limiting cases, we have outlined a microscopic approach that involves the formation of a condensate, chiral symmetry, and deconfinement.In addition to colored single-quark quasiparticles, colored clusters are also included.
A main topic is the treatment of correlations, for which we employed the correlation functional Φ.Within a cluster decomposition, the special class of "sunset" diagrams was singled out and a generalized Beth-Uhlenbeck equation was used.This approach allows not only a systematic inclusion of continuum correlations via the scattering phase shifts, but also the inclusion of in-medium effects describing the dissolution of bound states with increasing energy density.
In this work, we have developed the technique of the cluster-virial expansion approach for multi-quark clusters at finite temperatures and chemical potentials based on a cluster generalization of the Φ-derivable approach.When a restriction to closed two-loop diagrams in cluster Green's functions is applied to the Φ-functional, this approach is equivalent to the generalized Beth-Uhlenbeck approach to clustering in hot and dense Fermi systems.
The complexity of color confinement in low-density quark matter is accounted for by coupling the quarks and their clusters to the Polyakov-loop background field, which serves to suppress the appearance of colored clusters in the quark confinement region.We have compared our results for the thermodynamics of clustered quark matter with recent lattice QCD simulations at finite temperature and small chemical potentials where they are currently available and found satisfactory agreement.
Our approach contains some simplifying approximations, especially with respect to the in-medium phase shifts.A self-consistent solution of the gap equation in a correlated medium would be of interest to obtain consistent results.We expect that such a solution would exhibit a behaviour which justifies the assumption of a sudden switch model for the quark masses that we made in this work.The few-quark wave equation describing the shift in binding energies and the phase shifts δ i is a subject of further investigation.In addition to the sunset diagrams, other diagrams for the correlation functional Φ can be considered to improve this approximation.Baryon densities and entropies have been compared with lattice QCD simulations to verify the quality of our approach.Consistent expressions for the various equations of state for the corresponding thermodynamic variables are needed to simulate matter under extreme conditions, such as those encountered in ultrarelativistc heavy-ion collisions and in dense astrophysical objects.

A Matsubara summations for multiquark Polyakov-loop propagators A.1 Diquark propagator and distribution function
In this work we restrict ourselves to the consideration of the scalar, color antitriplet diquark.The free color antitriplet diquark propagators are obtained after color and Matsubara sum of the product of two color triplet singlequark propagators folded with the antisymmetric Gell-Mann matrices λ 2 , λ 5 or λ 7 .For a better systematics, we will replace these matrices by the completely equivalent totally antisymmetric symbol ε abc in color space.We introduce the bosonic Matsubara frequency Ω m = ω n + ω ′ n as the sum of the fermionic Matsubara frequencies of the quark propagators that constitute the diquark and obtain the free quark pair propagator.
where 1 c is the diagonal matrix in color space and we have employed the result of the summation over color and Matsubara frequency of the single quark propagators given by Eq. ( 13) and used the identity −(A 4 ) cc = ε abc [(A 4 ) aa + (A 4 ) bb ].
The diquark distribution function in the presence of the Polyakov-loop background field can now be obtained analogously to the case of the single quark distribution.For that we introduce the abbreviation ω = E p + E p ′ for the two-particle energy and neglect the compositeness effect by suppressing the Pauli blocking factor − f A.2 Three-quark (nucleon) propagator and distribution function We obtain higher order cluster Green's functions as a result of their bipartition into lower lying cluster states which are already known.In this way, we combine the color antitriplet diquark (anticolor a) with a color triplet quark (color a) in order to obtain a color singlet nucleon, see Tab. 5.In this way, we obtain the free three-quark (nucleon=quark+diquark) propagator as Now it is obvious that the nucleon is a color singlet state and no summation over colored degrees of freedom need to be done.The summation over the fermonic Matsubara frequancy ω where y + 3 (ω) = exp[−β(ω − 3µ)] and we have neglected the phase space occupation factor in the numerator of Eq. (73).

A.3 Polyakov-loop generalized multi-quark distribution functions
On the basis of the derivations presented in Appendices A.1 and A.2, we can now summarize the generalized multiquark distribution functions that are obtained for multiquark states which are dimers of other multiquark states.To this end it is sufficient to know the Matsubara frequencies for the a−quark cluster and its color structure (triplet, antitriplet or singlet) as given in Tab. 5.
For the color triplet state with an odd number of net valence quarks, we obtain B High-temperature O(g 2 ) perturbative QCD limit An important benchmark for the development of an effective model description that reproduces lattice QCD thermodynamics is the O(g 2 ) perturbative QCD limit at high temperatures.For the readers convenience, we quote here expressions from Appendix B of Ref. [20], starting with the pressure at this order, which reads for the three-flavor case (N f = 3) In table 6, we give the limiting values for the dimensionless baryon density n B /T 3 and entropy s/T 3 to O(g 2 ) for four + e iβ(φ8−φ3) + e iβ(φ8+φ3) , φ = ϕ * .

Fig. 1 .
Fig. 1.Feynman diagrams for the Φ functional in the cluster two-loop approximation.A n-quark cluster Green's function is shown as n solid lines bound with one arrow, a meson Green's function is a dashed line with an arrow.Starting from the top, the figure shows the diquark as two-quark state, the meson as quark-antiquark state, the nucleon as quark-diquark state and the 4-quark cluster as nucleon-quark or diquark-diquark state.
q+Q), (D+F), (N+N) (2nπT − i6µ)1c singlet Table5.Matsubara frequencies for an n−particle cluster with a net valence quarks in the Polyakov-gauge quark model, obtained as composite of a pair of lower order clusters, see Fig.1and Fig.2

Fig. 2 .
Fig. 2. Continuation of Fig. 1.Starting from the top, we show the Feynman diagrams for the Φ functional in the cluster two-loop approximation whith the tetraquark as meson-meson or diquark-antidiquark state, the 5-quark cluster as 4-quark-quark state or nucleon-diquark state, the pentaquark as nucleon-meson state or 4-quark-antiquark state and the hexaquark channels: 5-quark-quark, 4-quark-diquark and nucleon-nucleon states.
shift with respect to reflection ω → −ω have been used and the "no sea" approximation has been employed which removes the divergent vacuum contribution.The Polyakov-loop modified distribution functions are defined as f (a),± ϕ (a even)

Fig. 3 .Fig. 4 .
Fig. 3. Mass spectrum of light and strange quarks used in this work, together with mass and decay width of the Breit-Wigner model for the pion and the nucleon as generic examples for hadrons as a function of temperature for vanishing chemical potential µ/T = 0.

Fig. 5 .
Fig. 5. Approximations for the hadronic phase shift for the example of the pion at rest in the medium.Upper panel: step-up (SU), middle panel: step-up-step-down (SUSD), bottom panel: step-up-continuum (SUC) models.

Fig. 7 .
Fig.7.Scaled entropy density s/T 3 as a function of temperature T calculated at µB/T = 0 (upper panel) and µB/T = 3.0 (lower panel).Partial contributions of hadrons, quarks, gluons as well as perturbative correction and total s/T 3 are represented by different curves indicated in figures.The shaded regions correspond to the lattice QCD results[19] in good agreement with the result of the present model (thick solid line).

Fig. 8 .
Fig. 8. Scaled entropy density s/T 3 as a function of temperature T calculated at different values of µB/T indicated in the legend.The shaded regions of the corresponding colors represent the lattice QCD results [19].The triangles of corresponding colors are the values obtained within O(αs) perturbative QCD, see Tab. 6, and are attained for large T .Note the broken temperature axis.

Fig. 9 .
Fig.9.Fractions of entropy density carried by different species as a function of temperature T calculated at vanishing baryochemical potential.Notably, the red and turquoise solid lines stands for the entropy fractions of quarks and gluons, resp., which are dominant in the QGP phase and strongly suppressed below the Mott transition temperature, because of the effective confinement mechanism in the present work.

Fig. 11 .
Fig.11.Baryon density as a function of temperature for µB/T = 3.0 (blue solid line) compared to lattice QCD simulations (light blue band)[19].Shown are also the contributions from partial densities of quarks (orange dashed line), colored clusters (red long dashed line), hadrons (green dotted line) and perturbative QCD (magenta very long dashed line).

Fig. 12 .
Fig.12.Baryon density as a function of temperature for different values of µB/T compared to lattice QCD simulations[19].The triangles of corresponding colors are the values obtained within O(αs) perturbative QCD, see Tab. 6, which are attained at large temperature.
consider the case of vanishing charge and strangeness chemical potentialsµ S = µ Q = 0, for which µ u = µ d = µ s = µ = µ B /3.From the pressure as the thermodynamic potential, one can derive the other equations of state, such as the density and the entropy density,

Table 1 .
Properties of quarks: B f -baryon number, d f -degeneracy, M f -constituent mass and m f -current mass of quark flavor f .

Table 2 .
(53)es of colored multiquark clusters according to the simple hadron mass formula(52)together with the corresponding multiquark thresholds(53).Bi is the baryon number and di the degeneracy of the state i.

Table 4 .
Color -singlet hadrons with baryon number Bi = 1, masses MPDG and degeneracy factor di according to the PDG