11institutetext:  Institute of Theoretical Physics, University of Wroclaw, Max
Born place 9, 50-204 Wroclaw, Poland 22institutetext: Helmholtz-Zentrum
Dresden-Rossendorf (HZDR), Bautzner Landstrasse 400, 01328 Dresden, Germany
33institutetext: Center for Advanced Systems Understanding (CASUS), Untermarkt
20, 02826 Görlitz, Germany 44institutetext: Incubator of Scientific
Excellence—Centre for Simulations of Superdense Fluids, University of Wroc?aw,
50-204, Wroclaw, Poland 55institutetext: Institute of Physics, University of
Rostock, Albert-Einstein Str. 23-24, 18059 Rostock, Germany

#  Thermodynamics of quark matter with multiquark clusters in an effective
Beth-Uhlenbeck type approach††thanks: This work has been initiated by Peter
Schuck in discussions during the Workshop on ”Light Clusters in Nuclei and
Nuclear Matter: Nuclear Structure and Decay, Heavy Ion Collisions and
Astrophysics”, at the ECT* Trento, September 2 - 6, 2019.

D. Blaschke 112233 M. Cierniak 11 O. Ivanytskyi 44 G. Röpke 1155

(Received: / Revised version: date)

###### Abstract

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 use simple ansätze 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\leq n\leq 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
$\mathcal{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.

##  1 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 Alford:2007xm ; Anglani:2013gfu ;
Pisarski:2019cvo ; Bauswein:2022vtq  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
Riley:2019yda ; Miller:2019cac  and PSR J0740+6620 Riley:2021pdl ;
Miller:2021qha  (currently the heaviest known neutron star with $2.08\pm
0.07~{}M_{\odot}$, as independent Shapiro–delay based mass measurements
NANOGrav:2019jur ) 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_{\odot}$ LIGOScientific:2017vwq ; LIGOScientific:2018cki ;
LIGOScientific:2018hze . 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 Bauswein:2018bma . In
addition, the hypothesis of a phase transition induced shock revival in
core–collapse supernova Fischer:2017zcr  might hint at an additional
observational avenue of studying the structure of the hot and dense sector of
the QCD phase diagram (cf. Bauswein:2022vtq  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 ($\mu_{B}=0$) sector of the QCD
phase diagram is well explored within this approach, see, for instance, Ref.
Borsanyi:2010bp ; Borsanyi:2012cr ; Bazavov:2014xya ; Bazavov:2018mes .
However, despite recent progress in calculating finite $\mu_{B}$ observables
on the lattice Borsanyi:2021sxv ; Bollweg:2022fqq ; Bollweg:2022rps , 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 $\mu_{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
Nambu:1961tp ; Nambu:1961fr . It is extensively used in studies of dynamical
chiral symmetry breaking and restoration (cf. Buballa:2003qv ). 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 Fukushima:2003fw ;
Ratti:2005jh ; Roessner:2006xn ; Roessner:2007gha ) 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 $\mathcal{O}(\alpha_{s})$ in the
expansion w.r.t. the strong fine structure constant $\alpha_{s}=g^{2}/(4\pi)$
stemming from one-gluon exchange which, following Refs. Blaschke:2020lhm  (for
details, see also Blaschke:2016hzu ; Blaschke:2016fdh ) applies only to
processes with high momentum transfer $p>\Lambda_{\rm 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 $\mu_{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
$\mu_{B}$, where the resulting hadrons are composite quark clusters. In the
spirit of Blaschke:2020lhm , the composite nature of the low density hadronic
degrees of freedom is accounted for by employing a generalized Beth-Uhlenbeck
formula Blaschke:2016hzu ; Blaschke:2016fdh ; Hufner:1994ma ; Zhuang:1994dw ;
Wergieluk:2012gd ; Yamazaki:2012ux ; Blaschke:2013zaa ; Yamazaki:2013yua ;
Torres-Rincon:2016ahl ; Torres-Rincon:2017zbr ; Blaschke:2016sqn , 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 Hagedorn:1965st . 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}$ Borsanyi:2011sw ;
HotQCD:2012fhj ; Bellwied:2015lba ; Bellwied:2013cta  and HIC experiments (cf.
Andronic:2017pug  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 non-interacting 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 Ropke:1986qs ;
Blaschke:2020qrs . 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 ROPKE1982536  and known from solid state and
plasma physics Kraeft:1986 . 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) Horowitz:1985tx ; Ropke:1986qs . In Ref.
Yukalov:1997jk , 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.
Albright:2014gva . The SFM model for the QGP has been considerably improved in
the subsequent work of Khvorostukhin et al. Khvorostukin:2006aw , 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 Kaltenborn:2017hus . 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 Ivanytskyi:2022oxv .

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 $\delta_{i}$ of a hadron species $i$ in the
generalized Beth-Uhlenbeck formula. The ansatz will be similar to the one used
in Blaschke:2016hzu ; Blaschke:2016fdh , in agreement with the Levinson
theorem. The resulting contribution to the thermodynamic observables will be
referred to as the Mott–hadron–resonance–gas (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 $\Phi$-derivable approach which is based on the Luttinger-Ward-
Baym functional Luttinger:1960ua ; Baym:1961zz ; Baym:1962sx  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 Liu:2017qah . 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 $\Phi$-derivable approach to deal with
multi-quark clusters.

As an extension of the previous work Blaschke:2020lhm , 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 continu-um–extrapolated LQCD
calculations of the baryon density and the entropy density as a function of
the temperature for finite values of $\mu_{B}/T$ Borsanyi:2021sxv ;
Bollweg:2022fqq . 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 self–consistently, we make the
assumption of a sudden drop of the quark mass at the $\mu_{B}$-dependent
pseudocritical temperature $T_{c}(\mu_{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_{\rm
Mott}(\mu_{B})$ where all hadronic bound states dissociate into their quark
constituents. This allows us to focus on the cluster contribution to the
Polyakov–loop 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 Polyakov–loop 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
$\mu_{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.

##  2 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 $\mathcal{U}(\phi,\bar{\phi})$ and a
$\mathcal{O}(\alpha_{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 Blaschke:2020lhm , 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 $\Phi$-derivable approach.
Within a Green’s function approach Baym:1961zz ; Baym:1962sx  the correlation
functional $\Phi$ is introduced where the contributions of interaction are
represented by closed-loop diagrams. When for the $\Phi$ 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 Ropke:2012qv . 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.

We introduce an ansatz for the thermodynamic potential which separates the QGP
sector of quark and gluon quasiparticles from that of the MHRG comprised of
hadrons which are understood as quark bound states (multiquark clusters) that
can undergo a Mott dissociation,

$\Omega(T,\mu,\phi,\bar{\phi})=\Omega_{\rm
QGP}(T,\mu,\phi,\bar{\phi})+\Omega_{\rm MHRG}(T,\mu,\phi,\bar{\phi}).$ (1)

In the following subsections we explain these two contributions more in
detail.

###  2.1 QGP thermodynamics of quarks and gluons

The thermodynamics of the deconfined QGP can be separated into a perturbative
and a nonperturbative contribution

$\Omega_{\rm QGP}(T,\mu,\phi,\bar{\phi})=\Omega_{\rm
NP}(T,\mu,\phi,\bar{\phi})+\Omega_{\rm pert}(T,\mu,\phi,\bar{\phi}),$ (2)

where $\Omega_{\rm NP}$ describes the nonperturbative low-energy QCD quark and
gluon degrees of freedom via the mean–field thermodynamic potential

$\Omega_{\rm
NP}(T,\mu,\phi,\bar{\phi})=\Omega_{Q}(T,\mu,\phi,\bar{\phi})+\mathcal{U}(\phi,\bar{\phi})$
(3)

and perturbative corrections are absorbed in $\Omega_{\rm
pert}(T,\mu,\phi,\bar{\phi})$.

We take into account the minimal coupling of quarks to a homogeneous gluon
background mean field in the Polyakov gauge
$A_{4}=\lambda_{3}A_{4}^{3}+\sqrt{3}\lambda_{8}A_{4}^{8}$ while quark
interactions are accounted for by a suitably chosen relativistic energy
density functional $\mathcal{L}_{\rm int}$

$\displaystyle\mathcal{L}_{0}$ $\displaystyle=$
$\displaystyle\bar{q}\left[\gamma_{0}(i\omega_{n}+\mu-
iA_{4})-\vec{\gamma}\cdot\vec{p}-m_{0}\right]q$ (4)
$\displaystyle-\mathcal{U}(\phi,\bar{\phi})-\mathcal{L}_{\rm int}~{},$

with $\mathcal{U}(\phi,\bar{\phi})$ denoting the Polyakov–loop (gluon)
potential for which we will not make use of the traditional form found in
Ratti:2005jh , but the one given in Lo:2013hla  as explained below. The traced
Polyakov loop is defined as

$\phi={\rm Tr}_{c}\exp(i\beta A_{4})/N_{c},$ (5)

where $N_{c}=3$ is the number of colors in the $SU(N_{c})$ gauge theory and
$\beta=1/T$ is the inverse temperature. Abbreviating this expression in
accordance with the standard notation found in the literature, we have
$A_{4}^{3}=\varphi_{3}$ and $A_{4}^{8}=\varphi_{8}$, so that $A_{4}={\rm
diag}(\varphi_{3}+\varphi_{8},-\varphi_{3}+\varphi_{8},-2\varphi_{8})$ and
finally

$\phi=\frac{1}{N_{c}}\left({\rm e}^{-2i\beta\varphi_{8}}+{\rm
e}^{i\beta(\varphi_{8}-\varphi_{3})}+{\rm
e}^{i\beta(\varphi_{8}+\varphi_{3})}\right)~{},~{}\bar{\phi}=\phi^{*}~{}.$ (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\omega_{n}\to
i\omega_{n}-iA_{4}$. We define the distribution function of a quark of given
color $c$ as the Matsubara sum over the quark propagator

$F_{c}(E_{p}-\mu+iA_{4})\equiv\frac{1}{\beta}\sum_{n}\frac{1}{(i\omega_{n}-E_{p}+\mu)\mathbbm{1}_{c}-i(A_{4})_{cc}}~{}.$
(7)

The Polyakov-loop distribution function is then obtained after performing the
color trace

$\displaystyle f^{(1),+}_{\phi}$ $\displaystyle=$ $\displaystyle{\rm
Tr}_{c}F_{c}(E_{p}-\mu+iA_{4})$ (8) $\displaystyle=$
$\displaystyle\sum_{c=1}^{3}f[(E_{p}-\mu)+i(A_{4})_{cc}],$ (9)

where $f[(E_{p}-\mu)]$ is the ordinary Fermi function which can be written as

$f(E_{p}-\mu)\equiv-\frac{1}{\beta}\frac{\partial z^{+}}{\partial E_{p}},$
(10)

with $z^{+}\equiv\ln[1+y^{+}_{1}]$ and $y^{+}_{1}={\rm
e}^{-\beta(E_{p}-\mu)}$. Now we can write

$\displaystyle f^{(1),+}_{\phi}$ $\displaystyle=$
$\displaystyle\frac{1}{\beta}\sum_{c=1}^{3}\frac{\partial\ln\left[1+y^{+}_{1}{\rm
e}^{-i\beta(A_{4})_{cc}}\right]}{\partial E_{p}}$ (11) $\displaystyle=$
$\displaystyle\frac{1}{\beta}\frac{\partial}{\partial
E_{p}}\ln\left[1+3(\bar{\phi}+\phi y^{+}_{1})y^{+}_{1}+{y^{+}_{1}}^{3}\right]$
(12)

and obtain the Polyakov-loop generalized single-quark distribution function

$\displaystyle f^{(1),+}_{\phi}$ $\displaystyle=$
$\displaystyle\frac{(\bar{\phi}+2\phi
y^{+}_{1})y^{+}_{1}+{y^{+}_{1}}^{3}}{1+3(\bar{\phi}+\phi
y^{+}_{1})y^{+}_{1}+{y^{+}_{1}}^{3}}~{}.$ (13)

We note that the distribution function for an antiquark
$\left[f^{(1),-}_{\phi}\right]^{*}$ is obtained from Eq. (13) by the
replacements: $\mu\to-\mu$ and $\phi\to\bar{\phi}=\phi^{*}$. The
nonperturbative part of the QGP thermodynamics (3) consists of a quark
quasiparticle contribution and the Polyakov–loop potential $\mathcal{U}$ for
which we adopt the form defined in Lo:2013hla ,

$\displaystyle\mathcal{U}(T,\phi,\bar{\phi})$ $\displaystyle=$ $\displaystyle
T^{4}\left[-\frac{1}{2}a(T)\bar{\phi}\phi+b(T)\log M_{H}\right.$
$\displaystyle\left.+\frac{1}{2}c(T)\left(\phi^{3}+\bar{\phi}^{3}\right)+d(T)\left(\bar{\phi}\phi\right)^{2}\right]$

with $M_{H}$ being the SU(3) Haar measure

$M_{H}=1-6\bar{\phi}\phi+4\left(\phi^{3}+\bar{\phi}^{3}\right)-3\left(\bar{\phi}\phi\right)^{2},$
(15)

and

$a(T)=\frac{a_{1}+a_{2}\left(T_{0}/T\right)+a_{3}\left(T_{0}/T\right)^{2}}{1+a_{4}\left(T_{0}/T\right)+a_{5}\left(T_{0}/T\right)^{2}}+a_{6}\frac{\mu^{2}T_{0}^{2}}{T^{4}},$
(16)
$b(T)=b_{1}\left(T_{0}/T\right)^{b_{4}}\left[1-e^{b_{2}\left(T_{0}/T\right)^{b_{3}}}\right],$
(17)
$c(T)=\frac{c_{1}+c_{2}\left(T_{0}/T\right)+c_{3}\left(T_{0}/T\right)^{2}}{1+c_{4}\left(T_{0}/T\right)+c_{5}\left(T_{0}/T\right)^{2}},$
(18)
$d(T)=\frac{d_{1}+d_{2}\left(T_{0}/T\right)+d_{3}\left(T_{0}/T\right)^{2}}{1+d_{4}\left(T_{0}/T\right)+d_{5}\left(T_{0}/T\right)^{2}}.$
(19)

The parameters $a_{1}\dots a_{5}$, $b_{1}\dots b_{4}$, $c_{1}\dots c_{5}$,
$d_{1}\dots d_{5}$ are taken from the pure SU(3) lattice gauge fit performed
in Lo:2013hla . The parameter $a_{6}=1.28$ is introduced here to improve the
description of the baryon density in comparison with LQCD results
Borsanyi:2021sxv . 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
Schaefer:2007pw , 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 Borsanyi:2021sxv  discussed below.

The mean field thermodynamics in the quark sector of the PNJL model is
described by

$\displaystyle\Omega_{Q}(T,\mu,\phi,\bar{\phi})$ $\displaystyle=$
$\displaystyle\sum_{f=u,d,s}\frac{\Delta^{2}_{f}(T,\mu)}{8G_{s}}-2N_{f}N_{c}\int_{0}^{\Lambda}\frac{dp\,p^{2}}{2\pi^{2}}E_{p}$
$\displaystyle-$ $\displaystyle
N_{f}N_{c}\int\frac{dp\,p^{2}}{3\pi^{2}}\frac{p}{E_{p}}\left\\{f^{(1),+}_{\phi}+\left[f^{(1),-}_{\phi}\right]^{*}\right\\}.$

Here, $G_{s}$ is the scalar coupling constant, and
$\Delta_{f}(T,\mu)=M_{f}(T,\mu)-m_{f}$ is the dynamically generated quark mass
gap for the flavor $f$, see Tab. 1.

$B_{f}$ | quark |  $M_{f}$ [MeV] |  $m_{f}$ [MeV] | $d_{f}=2N_{c}N_{f}$
---|---|---|---|---
1/3 | q=u,d | $627$ | 5.6 | $12$
1/3 | s | $770$ | 124 | $6$
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$.

The same structure for the mean field thermodynamic potential (2.1) is also
obtained in the confining relativistic density functional (RDF) approach, see
Eqs. (16) and (17) of Ivanytskyi:2022oxv , 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
Ivanytskyi:2021dgq  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 $\Omega_{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 $\Lambda$ is chosen
such, that

$\frac{N_{f}N_{c}}{\pi^{2}}\int^{\Lambda(T,\mu)}_{0}dp\,p^{2}\,E_{p}-\sum_{f=u,d,s}\frac{\Delta^{2}_{f}(T,\mu)}{8G_{s}}=0,$
(21)

simplifying the quasiparticle expression for the quark thermodynamical
potential to

$\Omega_{Q,{\rm
reg}}=-N_{f}N_{c}\int\frac{dp\,p^{2}}{3\pi^{2}}\frac{p}{E_{p}}\left\\{f^{(1),+}_{\phi}+\left[f^{(1),-}_{\phi}\right]^{*}\right\\}.$
(22)

In addition, we introduce a perturbative virial correction contribution in two
loop order, $\Omega_{\rm pert}(T,\mu,\phi,\bar{\phi})$ with a standard
expression (cf. Kapusta:1989tk , see also Turko:2011gw ) of the form

$\displaystyle\Omega_{\rm pert}(T,\mu,\phi,\bar{\phi})=$
$\displaystyle\frac{8}{\pi}\frac{\alpha_{s}N_{f}}{\beta^{4}}$
$\displaystyle\left\\{\frac{1}{6}\left(I^{+}+I^{-}\right)\right.$
$\displaystyle\left.+\frac{1}{4\pi^{2}}\left(I^{+}+I^{-}\right)^{2}\right\\}$

where

$I^{\pm}(T,\mu,\phi,\bar{\phi})=\int^{\infty}_{\Lambda_{\rm
pert}/T}dx\frac{x}{3}f^{(1),\pm}_{\phi}$ (24)

are modified integrals introduced in Blaschke:2020lhm  and extended to finite
chemical potentials with the generalized Polyakov–loop generalized Fermi
distribution functions. The infrared cutoff, $\Lambda_{\rm 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 Blaschke:2005jg ; Shirkov:2002td ;
Peshier:2003ah

$\displaystyle\alpha_{s}(T,\mu)$ $\displaystyle=$
$\displaystyle\frac{12\pi}{11N_{c}-2N_{f}}\left[\frac{1}{\ln\left(x^{2}\right)}-\frac{1}{x^{2}-1}\right],$
(25) $\displaystyle x^{2}$ $\displaystyle=$
$\displaystyle\left(\frac{T}{T_{*}}\right)^{2}+\left(\frac{\mu}{\mu_{*}}\right)^{2},$
(26)

with $N_{c}=3$, $N_{f}=3$, $T_{*}=93.75$ MeV. We choose $\mu_{*}=3\pi T_{*}$
in accordance with the argument of the switch function in Albright:2014gva .

The traced Polyakov loop variables $\phi$ and $\bar{\phi}$ 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 Blaschke:2020lhm  for $\mu=0$. As an
extension of that work, we have introduced finite–$\mu$ 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.

###  2.2 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 Albright:2014gva  both approaches
are matched together, using a switching function for a smooth transition.
There, the crossover region at $\mu=0$ occurs around $T=170$ MeV and goes to
$T=0$ around $\mu=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. Albright:2014gva
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 quark-gluon-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 Ropke:2012qv  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 relevance 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 many-
particle approach and has also been applied in nuclear physics to describe
cluster formation and resonances in nuclear matter Ropke:2014fia ;
Ropke:2020peo .

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 color-neutral
clusters are selected from the Particle Data Group ParticleDataGroup:2022pth
summary tables for mesons in Tab. 3 and for baryons in Tab. 4, according to
the choice made in Ref. Albright:2014gva . 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 Buccella:2020mxi  assuming a
compact three-diquark structure, see also Shahrbaf:2022upc  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 Stancu:2005jv ; Maiani:2023nwj , but those containing only three
light quark flavors are not yet confirmed by experiments.

$B_{i}$ | colored cluster | $d_{i}$ | $M_{i}$ | $M_{{\rm th},i}^{<}$ | $M_{{\rm th},i}^{>}$
---|---|---|---|---|---
|  |  | [MeV] | [MeV] | [MeV]
2/3 |  $D_{1}$ ($uu$/$ud$/$dd$) | $3$ | $783$ | $1254$ | $11.2$
2/3 |  $D_{2}$ ($us$/$ds$) | $6$ | $926$ | $1397$ | $129.6$
4/3 |  $4q_{1}$ ($D_{1}D_{1}$) | $3$ | $1095$ | $2508$ | $22.4$
4/3 |  $4q_{2}$ ($D_{1}D_{2}$) | $6$ | $1238$ | $2651$ | $140.8$
5/3 |  $5q_{1}$ ($ND_{1}$) | $12$ | $1251$ | $3135$ | $28.0$
5/3 |  $5q_{2}$ ($ND_{2}$) | $24$ | $1394$ | $3278$ | $146.4$
5/3 |  $5q_{3}$ ($\Lambda D_{2}$) | $48$ | $1537$ | $3421$ | $264.8$
Table 2: Masses of colored multiquark clusters according to the simple hadron mass formula (52) together with the corresponding multiquark thresholds (53). $B_{i}$ is the baryon number and $d_{i}$ the degeneracy of the state $i$.  PDG | $d_{i}$ | $M_{\rm PDG}$ | $M_{i}$ | $M_{{\rm th},i}^{<}$ | $M_{{\rm th},i}^{>}$
---|---|---|---|---|---
mesons |  | [MeV] | [MeV] | [MeV] | [MeV]
$\pi^{+}$/$\pi^{0}$ | $3$ | $140$ | 140 | $1254$ | $11.2$
$K^{+}$/$K^{0}$ | $4$ | $494$ | $494$ | $1397$ | $129.6$
$\eta$ | $1$ | $548$ | $878$ | $1349$ | $90.1$
$\rho^{+}$/$\rho^{0}$ | $9$ | $775$ | 783 | $1254$ | $11.2$
$\omega$ | $9$ | $783$ | 783 | $1254$ | $11.2$
$K^{*+}$/$K^{*0}$ | $12$ | $895$ | $806^{*)}$ | $2651$ | $140.8$
$\eta^{\prime}$ | $1$ | $960$ | $878$ | $1349$ | $90.1$
$a_{0}$ | $3$ | $980$ | $1095^{*)}$ | $2508$ | $22.4$
$f_{0}$ | $1$ | $980$ | $1095^{*)}$ | $2508$ | $22.4$
$\phi$ | $3$ | $1020$ | 1069 | $1540$ | $248$
$h_{1}$ | $3$ | $1170$ | 1069 | $1540$ | $248$
$a_{1}$/$b_{1}$ | $18$ | $1230$ | $1095^{*)}$ | $2508$ | $22.4$
$K_{1}(1270)$ | $12$ | $1272$ | $1238^{*)}$ | $2651$ | $140.8$
$f_{2}$ | $5$ | $1275$ | $1095^{*)}$ | $2508$ | $22.4$
$f_{1}$ | $3$ | $1282$ | $1095^{*)}$ | $2508$ | $22.4$
$\eta(1295)$ | $1$ | $1294$ | $1190^{*)}$ | $2603$ | $101.4$
$\pi(1300)$ | $3$ | $1300$ | $1095^{*)}$ | $2508$ | $22.4$
$a_{2}$ | $15$ | $1318$ | $1095^{*)}$ | $2508$ | $22.4$
$f_{0}(1370)$ | $1$ | $1350$ | $1095^{*)}$ | $2508$ | $22.4$
$\pi_{1}(1400)$ | $9$ | $1354$ | $1095^{*)}$ | $2508$ | $22.4$
$K_{1}(1400)$ | $12$ | $1403$ | $1238^{*)}$ | $2651$ | $140.8$
$\eta(1405)$ | $1$ | $1409$ | $1190^{*)}$ | $2603$ | $101.4$
$K^{*}(1410)$ | $12$ | $1414$ | $1238^{*)}$ | $2651$ | $140.8$
$\omega(1420)$ | $3$ | $1425$ | $1095^{*)}$ | $2508$ | $22.4$
$K_{0}^{*}(1430)$ | $4$ | $1425$ | $1238^{*)}$ | $2651$ | $140.8$
$K_{2}^{*\pm}(1430)$ | $10$ | $1426$ | $1238^{*)}$ | $2651$ | $140.8$
$f_{1}(1420)$ | $3$ | $1426$ | $1095^{*)}$ | $2508$ | $22.4$
$K_{2}^{*0}(1430)$ | $10$ | $1432$ | $1238^{*)}$ | $2651$ | $140.8$
$\rho(1450)$ | $9$ | $1465$ | $1095^{*)}$ | $2508$ | $22.4$
$a_{0}(1450)$ | $3$ | $1474$ | $1095^{*)}$ | $2508$ | $22.4$
$\eta(1475)$ | $1$ | $1476$ | $1190^{*)}$ | $2603$ | $101.4$
$f_{0}(1500)$ | $1$ | $1505$ | $1095^{*)}$ | $2508$ | $22.4$
$f_{2}^{\prime}(1525)$ | $5$ | $1525$ | $1095^{*)}$ | $2508$ | $22.4$
$\eta_{2}(1645)$ | $5$ | $1617$ | $1190^{*)}$ | $2603$ | $101.4$
$\pi_{1}(1600)$ | $9$ | $1662$ | $1095^{*)}$ | $2508$ | $22.4$
$\omega_{3}(1670)$ | $7$ | $1667$ | $1095^{*)}$ | $2508$ | $22.4$
$\omega(1650)$ | $3$ | $1670$ | $1095^{*)}$ | $2508$ | $22.4$
$\pi_{2}(1670)$ | $15$ | $1672$ | $1095^{*)}$ | $2508$ | $22.4$
$\phi(1680)$ | $3$ | $1680$ | $1381^{*)}$ | $2794$ | $259.2$
$\rho_{3}(1690)$ | $21$ | $1689$ | $1095^{*)}$ | $2508$ | $22.4$
$K^{*}(1680)$ | $12$ | $1717$ | $1238^{*)}$ | $2651$ | $140.8$
$\rho(1700)$ | $9$ | $1720$ | $1095^{*)}$ | $2508$ | $22.4$
$f_{0}(1710)$ | $1$ | $1720$ | $1095^{*)}$ | $2508$ | $22.4$
$K_{2}(1770)$ | $20$ | $1773$ | $1238^{*)}$ | $2651$ | $140.8$
$K_{3}^{*}(1780)$ | $28$ | $1776$ | $1238^{*)}$ | $2651$ | $140.8$
$\pi(1800)$ | $3$ | $1812$ | $1095^{*)}$ | $2508$ | $22.4$
$K_{2}(1820)$ | $20$ | $1816$ | $1238^{*)}$ | $2651$ | $140.8$
$\phi_{3}(1850)$ | $7$ | $1854$ | $1381^{*)}$ | $2794$ | $259.2$
$\pi_{2}(1880)$ | $15$ | $1895$ | $1095^{*)}$ | $2508$ | $22.4$
$f_{2}(1950)$ | $5$ | $1944$ | $1095^{*)}$ | $2508$ | $22.4$
$a_{4}(2040)$ | $27$ | $1996$ | $1095^{*)}$ | $2508$ | $22.4$
$f_{2}(2010)$ | $5$ | $2011$ | $1095^{*)}$ | $2508$ | $22.4$
$f_{4}(2050)$ | $9$ | $2018$ | $1095^{*)}$ | $2508$ | $22.4$
$K_{4}^{*}(2045)$ | $36$ | $2045$ | $1238^{*)}$ | $2651$ | $140.8$
$\phi(2170)$ | $3$ | $2175$ | $1381^{*)}$ | $2794$ | $259.2$
$f_{2}(2300)$ | $5$ | $2297$ | $1095^{*)}$ | $2508$ | $22.4$
$f_{2}(2340)$ | $5$ | $2339$ | $1095^{*)}$ | $2508$ | $22.4$
Table 3: Meson masses and degeneracy factors $d_{i}$ according to the PDG ParticleDataGroup:2022pth , compared to $M_{i}$ from (58) for pseudoscalars or (52) and the continuum thresholds (53). States with ∗) are understood as tetraquark states.  PDG | $d_{i}$ | $M_{\rm PDG}$ | $M_{i}$ | $M_{{\rm th},i}^{<}$ | $M_{{\rm th},i}^{>}$
---|---|---|---|---|---
baryons |  | [MeV] | [MeV] | [MeV] | [MeV]
n/p | $4$ | $939$ | $939$ | $1881$ | $16.8$
$\Lambda$ | $2$ | $1116$ | $1082$ | $2024$ | $135.2$
$\Sigma$ | $6$ | $1193$ | $1082$ | $2024$ | $135.2$
$\Delta$ | $16$ | $1232$ | $1251^{**)}$ | $3135$ | $28$
$\Xi^{0}$ | $2$ | $1315$ | $1225$ | $2167$ | $253.6$
$\Xi^{-}$ | $2$ | $1322$ | $1225$ | $2167$ | $253.6$
$\Sigma(1385)$ | $6$ | $1385$ | $1394^{**)}$ | $3278$ | $146.4$
$\Lambda(1405)$ | $2$ | $1405$ | $1394^{**)}$ | $3278$ | $146.4$
$N(1440)$ | $4$ | $1440$ | $1251^{**)}$ | $3135$ | $28$
$\Lambda(1520)$ | $4$ | $1520$ | $1394^{**)}$ | $3278$ | $146.4$
$N(1520)$ | $8$ | $1520$ | $1251^{**)}$ | $3135$ | $28$
$\Xi^{0}(1530)$ | $4$ | $1532$ | $1537^{**)}$ | $3421$ | $264.8$
$\Xi^{-}(1530)$ | $4$ | $1535$ | $1537^{**)}$ | $3421$ | $264.8$
$N(1535)$ | $4$ | $1535$ | $1251^{**)}$ | $3135$ | $28$
$\Delta(1620)$ | $8$ | $1630$ | $1251^{**)}$ | $3135$ | $28$
$N(1650)$ | $4$ | $1655$ | $1251^{**)}$ | $3135$ | $28$
$\Sigma(1660)$ | $6$ | $1660$ | $1394^{**)}$ | $3278$ | $146.4$
$\Lambda(1670)$ | $2$ | $1670$ | $1394^{**)}$ | $3278$ | $146.4$
$\Sigma(1670)$ | $12$ | $1670$ | $1394^{**)}$ | $3278$ | $146.4$
$\Omega^{-}$ | $4$ | $1673$ | $1368$ | $2310$ | $372$
$N(1675)$ | $12$ | $1675$ | $1251^{**)}$ | $3135$ | $28$
$N(1680)$ | $12$ | $1685$ | $1251^{**)}$ | $3135$ | $28$
$\Lambda(1690)$ | $4$ | $1690$ | $1394^{**)}$ | $3278$ | $146.4$
$\Xi(1690)$ | $4$ | $1690$ | $1537^{**)}$ | $3421$ | $264.8$
$N(1700)$ | $8$ | $1700$ | $1251^{**)}$ | $3135$ | $28$
$\Delta(1700)$ | $16$ | $1700$ | $1251^{**)}$ | $3135$ | $28$
$N(1710)$ | $4$ | $1710$ | $1251^{**)}$ | $3135$ | $28$
$N(1720)$ | $8$ | $1720$ | $1251^{**)}$ | $3135$ | $28$
$\Sigma(1750)$ | $6$ | $1750$ | $1394^{**)}$ | $3278$ | $146.4$
$\Sigma(1775)$ | $18$ | $1775$ | $1394^{**)}$ | $3278$ | $146.4$
$\Lambda(1805)$ | $4$ | $1805$ | $1394^{**)}$ | $3278$ | $146.4$
$\Xi(1820)$ | $8$ | $1823$ | $1537^{**)}$ | $3421$ | $264.8$
$\Lambda(1825)$ | $12$ | $1825$ | $1394^{**)}$ | $3278$ | $146.4$
$N(1875)$ | $8$ | $1875$ | $1251^{**)}$ | $3135$ | $28$
$d$ | $3$ | $1875$ | $1251^{***)}$ | $3762$ | $33.6$
$S(uuddss)$ | $1$ | $1885$ | $1693^{****)}$ | $4048$ | $270.4$
$\Delta(1905)$ | $32$ | $1885$ | $1251^{**)}$ | $3135$ | $28$
$\Lambda(1890)$ | $4$ | $1890$ | $1394^{**)}$ | $3278$ | $146.4$
$N(1900)$ | $8$ | $1900$ | $1251^{**)}$ | $3135$ | $28$
$\Sigma(1915)$ | $18$ | $1915$ | $1394^{**)}$ | $3278$ | $146.4$
$\Delta(1925)$ | $48$ | $1925$ | $1251^{**)}$ | $3135$ | $28$
$\Sigma(1940)$ | $12$ | $1940$ | $1394^{**)}$ | $3278$ | $146.4$
$\Delta(1930)$ | $24$ | $1950$ | $1251^{**)}$ | $3135$ | $28$
$\Xi(1950)$ | $4$ | $1950$ | $1537^{**)}$ | $3421$ | $264.8$
$\Xi(2030)$ | $12$ | $2025$ | $1537^{**)}$ | $3421$ | $264.8$
$\Sigma(2030)$ | $24$ | $2030$ | $1394^{**)}$ | $3278$ | $146.4$
$\Lambda(2105)$ | $14$ | $2105$ | $1394^{**)}$ | $3278$ | $146.4$
$N(2195)$ | $36$ | $2220$ | $1251^{**)}$ | $3135$ | $28$
$\Sigma(2250)$ | $6$ | $2250$ | $1394^{**)}$ | $3278$ | $146.4$
$\Omega^{-}$(2250) | $2$ | $2252$ | $1680^{**)}$ | $3564$ | $383.2$
$N(2250)$ | $20$ | $2275$ | $1251^{**)}$ | $3135$ | $28$
$\Lambda(2350)$ | $10$ | $2350$ | $1394^{**)}$ | $3278$ | $146.4$
$\Delta(2420)$ | $48$ | $2420$ | $1251^{**)}$ | $3135$ | $28$
$N(2600)$ | $24$ | $2600$ | $1251^{**)}$ | $3135$ | $28$
Table 4: Color-singlet hadrons with baryon number $B_{i}=1$, masses $M_{\rm
PDG}$ and degeneracy factor $d_{i}$ according to the PDG
ParticleDataGroup:2022pth . Also shown are the masses $M_{i}$ according to our
simple hadron mass formula (52) and the masses of the continuum thresholds
$M_{{\rm 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 Shahrbaf:2022upc .

The MHRG part takes the form of a cluster decomposition of the thermodynamic
potential for quark matter

$\displaystyle\Omega_{\rm MHRG}(T,\mu,\phi,\bar{\phi})$ $\displaystyle=$
$\displaystyle\sum_{n=2}^{N}\Omega_{n}(T,\mu)+\Phi\left[\left\\{S_{n}\right\\}\right],$
(27) $\displaystyle\Omega_{n}(T,\mu)$ $\displaystyle=$ $\displaystyle
c_{n}\left[{\rm Tr}\ln S_{n}^{-1}+{\rm Tr}(\Pi_{n}S_{n})\right],$ (28)

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
Vanderheyden:1998ph ; Blaizot:2000fc . The functional
$\Phi\left[\left\\{S_{n}\right\\}\right]$ 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$.

The full quark cluster propagators $S_{n}$ fulfill a Dyson–Schwinger equation

$S_{n}^{-1}={S_{n}^{(0)}}^{-1}-\Pi_{n}~{},~{}n=2,\dots,N,$ (29)

where $S_{n}^{(0)}$ is the free $n-$quark cluster propagator and the cluster
selfenergy $\Pi_{n}$ should be obtained from

$\Pi_{n}=\frac{\partial\Phi}{\partial S_{n}}~{}.$ (30)

The free $n-$quark cluster propagator for $n=2,\dots,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. Bastian:2018wfl . 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.

Figure 1:  Feynman diagrams for the $\Phi$ 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.  Figure 2:  Continuation of Fig. 1. Starting from the top, we show the Feynman diagrams for the $\Phi$ 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.  symbol | cluster | quarks | $n$ | $a$ | bipartition | Matsubara frequency | color
---|---|---|---|---|---|---|---
$i$ |  | (antiquarks) |  |  |  | $\omega_{n}^{(a)}$ |
q | quark | 1 (0) | 1 | 1 | – | $((2n+1)\pi T-i\mu)\mathbbm{1}_{c}-A_{4}$ | triplet
D | diquark | 2 (0) | 2 | 2 | (q+q) | $(2n\pi T-i2\mu)\mathbbm{1}_{c}+A_{4}$ | antitriplet
M | meson | 1 (1) | 2 | 0 | (q+$\mathrm{\bar{q}}$) | $(2n\pi T)\mathbbm{1}_{c}$ | singlet
N | nucleon | 3 (0) | 3 | 3 | (q+D) | $((2n+1)\pi T-i3\mu)\mathbbm{1}_{c}$ | singlet
F | 4-quark | 4 (0) | 4 | 4 | (q+N), (D+D) | $(2n\pi T-i4\mu)\mathbbm{1}_{c}-A_{4}$ | triplet
T | tetraquark | 2 (2) | 4 | 0 | (M+M), (D+$\mathrm{\bar{D}}$) | $(2n\pi T)\mathbbm{1}_{c}$ | singlet
Q | 5-quark | 5 (0) | 5 | 5 | (q+F), (D+N) | $((2n+1)\pi T-i5\mu)\mathbbm{1}_{c}+A_{4}$ | antitriplet
P | pentaquark | 4 (1) | 5 | 3 | (M+N), ($\mathrm{\bar{q}}$+F) | $((2n+1)\pi T-i3\mu)\mathbbm{1}_{c}$ | singlet
H | hexaquark | 6 (0) | 6 | 6 | (q+Q), (D+F), (N+N) | $(2n\pi T-i6\mu)\mathbbm{1}_{c}$ | singlet
Table 5: 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. 1 and Fig. 2

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_{n}^{(0)}$ and the choice of the 2CI set of diagrams of the
$\Phi-$functional, built with the set of full cluster propagators
$\left\\{S_{n}\right\\}$. Following the spirit of Refs. Vanderheyden:1998ph ;
Blaizot:2000fc , we consider here the quark density and entropy density as
first derivatives of the thermodynamic potential with respect to the quark
chemical potential $\mu$ and the temperature $T$, respectively. From the
cluster decomposition (27) of the MHRG thermodynamic potential follows the
cluster decomposition of the quark density

$\displaystyle n_{\rm MHRG}(T,\mu)$ $\displaystyle=$
$\displaystyle-\frac{\partial\Omega_{\rm MHRG}(T,\mu)}{\partial\mu}$
$\displaystyle=$ $\displaystyle\sum_{i}a_{i}\,n_{i}(T,\mu),$

where $a_{i}$ is the net quark number in the $n_{i}-$particle state with the
partial density defined as

$\displaystyle n_{i}(T,\mu)$ $\displaystyle=$
$\displaystyle-\frac{\partial}{\partial\mu}\left\\{d_{i}\int\frac{d^{3}q}{(2\pi)^{3}}\int\frac{d\omega}{2\pi}\left[\ln\left(-S_{n_{i}}^{-1}\right)+\right.\right.$
$\displaystyle\left.\left.+\operatorname{Tr}\left(\Pi_{n_{i}}~{}S_{n_{i}}\right)\right]+\sum_{\stackrel{{\scriptstyle
i_{1},i_{2}}}{{i_{1}+i_{2}={n_{i}}}}}\Phi[S_{i_{1}},S_{i_{2}},S_{n_{i}}]\right\\}$

and $d_{i}$ denotes the degeneracy factor for $n_{i}-$particle states. Using
the spectral decomposition rule for an analytic complex function $F(\omega)$

$\displaystyle
F(i\omega_{n})=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\frac{2\,\mathrm{Im}F(\omega)}{\omega-i\omega_{n}}$
(33)

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)

$\displaystyle\frac{1}{3}\sum_{c=1}^{3}\frac{1}{\beta}\sum_{n}\frac{1}{i\omega_{n}^{(a)}-\omega\mathbbm{1}_{c}}=f^{(a),+}_{\phi}(\omega)$
(34)

with the relation $\partial f^{(a)}(\omega)/\partial\mu=-\partial
f^{(a)}(\omega)/\partial\omega$, we get for Eq. (LABEL:eq:n_a) now

$\displaystyle n_{i}(T,\mu)$ $\displaystyle=$ $\displaystyle-
d_{i}c_{a_{i}}\int\frac{d^{3}q}{(2\pi)^{3}}\int\frac{d\omega}{\pi}f^{(a_{i}),+}_{\phi}(\omega)\frac{\partial}{\partial\omega}\left[{\rm
Im}\ln\left(-S_{n_{i}}^{-1}\right)\right.$ (35) $\displaystyle+\left.{\rm
Im}\left(\Pi_{n_{i}}~{}S_{n_{i}}\right)\right]+\sum_{\stackrel{{\scriptstyle
i_{1},i_{2}}}{{i_{1}+i_{2}=n_{i}}}}\frac{\partial\Phi[S_{i_{1}},S_{i_{2}},S_{n_{i}}]}{\partial\mu}~{},$

where a partial integration over $\omega$ has been performed. For two-loop
diagrams of the sunset type holds a cancellation Vanderheyden:1998ph ;
Blaizot:2000fc , which we generalize here for cluster states

$\displaystyle
d_{i}c_{a_{i}}\int\frac{d^{3}q}{(2\pi)^{3}}\int\frac{d\omega}{\pi}$
$\displaystyle
f^{(a_{i}),+}_{\phi}(\omega)\frac{\partial}{\partial\omega}\left({\rm
Re}\Pi_{n_{i}}~{}{\rm Im}S_{n_{i}}\right)$ $\displaystyle-$
$\displaystyle\sum_{\stackrel{{\scriptstyle
i_{1},i_{2}}}{{i+j=n_{i}}}}\frac{\partial\Phi[S_{i_{1}},S_{i_{2}},S_{n_{i}}]}{\partial\mu}=0~{}.$

Using generalized optical theorems Bastian:2018mmc  and the polar
representation of the $n$–quark cluster Green’s function
$S_{n}=|S_{n}|\exp(i\delta_{n})$, where $\delta_{n}$ is the cluster phase
shift, one obtains

$\displaystyle\frac{\partial}{\partial\omega}$ $\displaystyle\left[{\rm
Im}\ln\left(-S_{n}^{-1}\right)+{\rm Im}\Pi_{n}~{}{\rm Re}S_{n}\right]$ (37)
$\displaystyle=$ $\displaystyle 2\,{\rm Im}\left[S_{n}~{}{\rm
Im}\Pi_{n}~{}\frac{\partial}{\partial\omega}S_{n}^{*}~{}{\rm
Im}\Pi_{n}\right]$ $\displaystyle=$
$\displaystyle-2\sin^{2}\delta_{n}\frac{\partial\delta_{n}}{\partial\omega}.$

The density is obtained in the form of a generalized Beth-Uhlenbeck EoS

$\displaystyle n_{\rm MHRG}(T,\mu)$ $\displaystyle=$
$\displaystyle\sum_{i}a_{i}\,d_{i}\,c_{a_{i}}\int\frac{d^{3}q}{(2\pi)^{3}}\int_{0}^{\infty}\frac{d\omega}{\pi}\left\\{f^{(a_{i}),+}_{\phi}\right.$
$\displaystyle\left.-\left[f^{(a_{i}),-}_{\phi}\right]^{*}\right\\}2\sin^{2}\delta_{n_{i}}(\omega,q)\frac{\partial\delta_{n_{i}}(\omega,q)}{\partial\omega}~{},$

where the properties of the distribution function $f^{(a),+}_{\phi}$ and the
phase shift with respect to reflection $\omega\to-\omega$ 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

$\displaystyle f^{(a),\pm}_{\phi}~{}\stackrel{{\scriptstyle\text{(a
even)}}}{{=}}$
$\displaystyle\frac{({\phi}-2\bar{\phi}y_{a}^{\pm})y^{\pm}_{a}+{y_{a}^{\pm}}^{3}}{1-3({\phi}-\bar{\phi}y_{a}^{\pm})y_{a}^{\pm}-{y_{a}^{\pm}}^{3}}~{},$
(39) $\displaystyle f^{(a),\pm}_{\phi}~{}\stackrel{{\scriptstyle\text{(a
odd)}}}{{=}}$ $\displaystyle\frac{(\bar{\phi}+2\phi
y_{a}^{\pm})y^{\pm}_{a}+{y_{a}^{\pm}}^{3}}{1+3(\bar{\phi}+\phi
y_{a}^{\pm})y_{a}^{\pm}+{y_{a}^{\pm}}^{3}}~{},$ (40)

where $y^{\pm}_{a}=e^{-\left(\omega\mp a\mu\right)/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

$\displaystyle s_{\rm MHRG}(T,\mu)$ $\displaystyle=$
$\displaystyle-\frac{\partial\Omega}{\partial T}=\sum_{i}s_{i}(T,\mu)$
$\displaystyle=$
$\displaystyle\sum_{i}d_{i}\,c_{a_{i}}\int\frac{d^{3}q}{(2\pi)^{3}}\int\frac{d\omega}{\pi}\left\\{\sigma^{(a_{i}),+}_{\phi}\right.$
$\displaystyle\left.+\left[\sigma^{(a_{i}),-}_{\phi}\right]^{*}\right\\}2\sin^{2}\delta_{n_{i}}(\omega,q)\frac{\partial\delta_{n_{i}}(\omega,q)}{\partial\omega}~{},$

where $\sigma^{(a)}=f^{(a)}_{\phi}\ln
f^{(a)}_{\phi}(-)^{a}[1(-)^{a}f^{(a)}_{\phi}]\ln[1(-)^{a}f^{(a)}_{\phi}]$ and
$f^{(a)}_{\phi}$ is the cluster distribution function for a net quark number
$a$ modified by the traced Polyakov loop.

The Eqs. (2.2) and (LABEL:eq:s) 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 $\Phi$\- functional Vanderheyden:1998ph ;
Blaizot:2000fc , where it has been shown that the correlation contribution
vanishes, see Eq. (LABEL:cancel). The formula for the pressure as
thermodynamical potential can be obtained from Eq. (2.2) by integration over
the quark chemical potential $\mu$. Analogously, it can be obtained from Eq.
(LABEL:eq:s) by integration over $T$

$p(T,\mu)=\int_{0}^{T}dT^{\prime}s(T^{\prime},\mu)~{}.$ (42)

It is instructive to consider the two limits $\phi=\bar{\phi}=1$
(deconfinement and color singlet states) and $\phi=\bar{\phi}=0$
(confinement). In the situation where $\phi=\bar{\phi}=1$, the distribution
functions reduce to

$\displaystyle f^{(a=0,2,4,\dots),\pm}_{\phi=1}$ $\displaystyle=$
$\displaystyle\frac{y_{a}^{\pm}}{1-y_{a}^{\pm}}~{},$ (43) $\displaystyle
f^{(a=1,3,5,\dots),\pm}_{\phi=1}$ $\displaystyle=$
$\displaystyle\frac{y_{a}^{\pm}}{1+y_{a}^{\pm}}~{},$ (44)

corresponding to the ordinary Bose (Fermi) distribution functions for even
(odd) numbers $a$ of quarks in the cluster. In the case of confinement when
$\phi=\bar{\phi}=0$, the distribution functions take the form

$\displaystyle f^{(a=0,2,4,\dots),\pm}_{\phi=0}$ $\displaystyle=$
$\displaystyle\frac{{y_{a}^{\pm}}^{3}}{1-{y_{a}^{\pm}}^{3}}~{},$ (45)
$\displaystyle f^{(a=1,3,5,\dots),\pm}_{\phi=0}$ $\displaystyle=$
$\displaystyle\frac{{y_{a}^{\pm}}^{3}}{1+{y_{a}^{\pm}}^{3}}~{}.$ (46)

From Eq. (45) one can see how the Polyakov-loop suppression of colored states
acts. Since ${y_{a}^{\pm}}^{3}(T)=y_{a}^{\pm}(T/3)$, in the situation of
confinement when $\phi=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 (2.2) in another form where we split the extended phase
shift $\delta_{i}(\omega,q)$ into its bound state part with a step-up at
$\omega=E_{i}(q)$ and the part related to the elastic scattering phase shift
$\delta_{i}^{\rm sc}(\omega,q)$ depending on the energy $\omega$ of relative
motion of the subclusters and the total momentum $q$ of the cluster. The net
baryon density (2.2) is the difference of contributions from particles and
antiparticles,

$\displaystyle n_{\rm MHRG}(T,\mu)$ $\displaystyle=$ $\displaystyle n^{+}_{\rm
MHRG}(T,\mu)-n^{-}_{\rm MHRG}(T,\mu)~{}.$ (47)

Applying integration by parts to the Beth-Uhlenbeck equation and considering
only the color singlet hadrons for which
$f^{(a),+}_{\phi}=f^{(a),+}_{\phi=1}$, the particle contribution can be
brought to the form

$\displaystyle n^{+}_{\rm MHRG}(T,\mu)$ $\displaystyle=$
$\displaystyle\sum_{i}a_{i}\,d_{i}\,c_{a_{i}}\int\frac{d^{3}q}{(2\pi)^{3}}\left\\{f^{(a),+}_{\phi=1}(E_{i}(q))\right.$
(48) $\displaystyle\left.-f^{(a),+}_{\phi=1}(E_{{\rm
thr},i})+\int\frac{d\omega}{\pi\,T}f^{(a),+}_{\phi=1}\left[\delta^{\rm
sc}_{i}(\omega,q)-\right.\right.$
$\displaystyle\left.\left.\frac{1}{2}\sin(2\delta^{\rm
sc}_{i}(\omega,q))\right]\right\\}~{},$

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 Kraeft:1986 . 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 (LABEL:eq:s), 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 $\pi$ for
every bound state at the bound state energy which is below the continuum edge
of scattering states (so-called ”step-up” phase shift).

##  3 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}(\mu)$, any medium modification of quark masses shall be strongly
suppressed.

Improved microscopic calculations within a rank-2 separable Dyson-Schwinger
equation (DSE) model Blaschke:1999ab  later extended within the Polyakov-loop
DSE model Horvatic:2010md ; Benic:2012ec , and recent results within the
confining density functional approach Ivanytskyi:2022oxv  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.

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_{\rm Mott}(\mu)$, 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}$,

$\displaystyle M_{f}(T,\mu)$ $\displaystyle=$ $\displaystyle
M_{f}\Theta(T_{\rm Mott}(\mu)-T)$ (49) $\displaystyle+m_{f}\Theta(T-T_{\rm
Mott}(\mu)).$

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.

Figure 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 $\mu/T=0$. Figure 4: Temperature dependence of
the Polyakov loop absolute value calculated for several values of $\mu_{B}/T$
indicated in the legend.

We extend the formula (49) to finite $\mu_{B}$ by using the chemical potential
dependence the pseudo-critical temperature $T_{c}(\mu)$ which resulted from
the fitting the lattice QCD results obtained by Taylor expansion techniques
for small chemical potentials Bazavov:2018mes  also for estimating the Mott
temperature,

$T_{\rm Mott}(\mu)=T_{\rm Mott}(0)\left[1-\kappa\left(\frac{\mu}{T_{\rm
Mott}(0)}\right)^{2}\right].$ (50)

We use the values $\kappa=0.012$ and $T_{\rm Mott}(0)=156.5$ MeV from the fit
of the $\mu$\- dependence of the pseudocritical temperature $T_{c}(\mu)=T_{\rm
Mott}(\mu)$ found in Bazavov:2018mes , restricting ourselves to the
$\mathcal{O}(\mu_{B}^{2})$ part of the expansion and neglecting the
$\mathcal{O}(\mu_{B}^{4})$ term because of the smallness of the associated
$\kappa_{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),

$\frac{\partial\Omega}{\partial\phi}=0~{}.$ (51)

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 $\phi$ and $\bar{\phi}$ 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 $\mu_{B}/T$. These results include the
backreaction from colored clusters which turn out to be negligibly small. For
temperatures below $T_{\rm Mott}(\mu)$, 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\sim 200$
MeV, this order parameter for deconfinement stays below $|\phi|\sim 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\sim 350$ MeV, the
order parameter asymptotically approaches $|\phi|=1.0$, the case of
deconfinement.

###  3.2 Hadron mass formulas

####  3.2.1 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_{\rm Mott}(\mu)$ the hadron masses
$M_{i}$ are medium independent and follow the approximate relation
Jankowski:2012ms

$M_{i}=(N_{q}+N_{\bar{q}})M_{q}+(N_{s}+N_{\bar{s}})M_{s}-(n-1)B~{},$ (52)

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_{\bar{q}}$ ($N_{\bar{s}}$) light
(strange) antiquarks. Here, $n=N_{q}+N_{\bar{q}}+N_{s}+N_{\bar{s}}$ is the
total number of quarks and antiquarks and
$a=N_{q}+N_{s}-N_{\bar{q}}-N_{\bar{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$

$\displaystyle M_{\rm thr,i}(T,\mu)$ $\displaystyle=$ $\displaystyle M_{\rm
thr,i}^{<}\Theta(T_{\rm Mott}(\mu)-T)$ (53) $\displaystyle+M_{\rm
thr,i}^{>}\Theta(T-T_{\rm Mott}(\mu))~{},$ $\displaystyle M_{\rm thr,i}^{<}$
$\displaystyle=$
$\displaystyle(N_{q}+N_{\bar{q}})M_{q}+(N_{s}+N_{\bar{s}})M_{s},$ (54)
$\displaystyle M_{\rm thr,i}^{>}$ $\displaystyle=$
$\displaystyle(N_{q}+N_{\bar{q}})m_{q}+(N_{s}+N_{\bar{s}})m_{s}.$ (55)

This behavior of the threshold masses entails that all PDG hadrons undergo the
Mott transition from bound states below $T_{\rm Mott}(\mu)$ to unbound
correlations in the continuum above that temperature, see Tabs. 3 and 4.

For temperatures above $T_{\rm Mott}$, a linear temperature dependence of the
mass and the decay width of the resonance is adopted

$\displaystyle M_{i}(T)$ $\displaystyle=$ $\displaystyle
M_{i,0}+\alpha_{1}\,\Gamma_{i}(T),$ (56) $\displaystyle\Gamma_{i}(T)$
$\displaystyle=$ $\displaystyle\alpha_{2}\,(T-T_{\rm Mott})\,\Theta(T-T_{\rm
Mott})~{},$ (57)

where $M_{i,0}$ stands for the vacuum mass of the resonance $i$ that is taken
from the particle data tables ParticleDataGroup:2022pth  and remains unchanged
for all temperatures below $T_{\rm Mott}$. Since the Mott temperature is
identified with the pseudocritical temperature $T_{c}$ obtained from lattice
QCD simulations Bazavov:2018mes , the ansatz (56), (57) leaves us with only
two free parameters, $\alpha_{1}=11.4$ and $\alpha_{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.2$ MeV) above that for the corresponding ground
state for $T<T_{\rm Mott}$ ($T>T_{\rm Mott}$). While in our calculations the
color singlet hadron masses are taken from the PDG ParticleDataGroup:2022pth ,
the masses of the colored clusters and the continuum thresholds require an
underlying quark model, see tables 2 \- 4.

####  3.2.2 Pseudoscalar meson masses

The pseudoscalar mesons play a special role as the pseudo-Goldstone bosons of
the broken chiral symmetry in low-energy QCD. Therefore, their masses follow a
Gell-Mann–Oakes–Renner relation

$M_{P}^{2}=-\frac{(N_{q}m_{q}+N_{s}m_{s})(N_{q}\langle\bar{q}q\rangle+N_{s}\langle\bar{s}s\rangle)}{(N_{q}-N_{s})f_{\pi}^{2}+2N_{s}f_{K}^{2}},$
(58)

where for the $P=\pi\,(K)$, the pion (kaon) mass of $140$ ($494$) MeV is
obtained with $N_{q}=1(1)$, $N_{s}=0(1)$, $m_{q}=5.6$ MeV, $m_{s}=124$ MeV,
$\langle\bar{q}q\rangle=-(242~{}{\rm MeV})^{3}$,
$\langle\bar{s}s\rangle=-(245~{}{\rm MeV})^{3}$, $f_{\pi}=92.4$ MeV and
$f_{K}=113$ MeV.

####  3.2.3 $\eta$ mesons

For the ground state $\eta$ mesons, we postulate a mass formula composed as a
weighted sum of light and strange quark masses with the average mesonic
binding energy

$M_{\eta}=M_{\eta^{\prime}}=\frac{2}{3}M_{s}+\frac{4}{3}M_{q}-B,$ (59)

with the corresponding threshold masses

$\displaystyle M_{\rm thr,\eta}(T,\mu)$ $\displaystyle=$
$\displaystyle\left(\frac{2}{3}M_{s}+\frac{4}{3}M_{q}\right)\Theta(T_{\rm
Mott}(\mu)-T)$ (60) $\displaystyle+$
$\displaystyle\left(\frac{2}{3}m_{s}+\frac{4}{3}m_{q}\right)\Theta(T-T_{\rm
Mott}(\mu)).$

The higher lying $\eta$ 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.

Figure 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.

###  3.3 Phase shift model

In this work, we want to introduce a generic model ansatz for the phase shift
$\delta_{i}(M,T,\mu)$ that shares all essential features with results that
have been obtained, e.g., within NJL model calculations for the case of the
pion Dubinin:2013yga ; Blaschke:2013zaa  and within the PNJL model for the
pion and the diquark Blaschke:2014zsa  in quark matter.

For the extended phase shifts, which include the bound state contribution as
”step-up” by $\pi$, 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. 1.

The phase shifts jump up from zero to $\pi$ 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 ParticleDataGroup:2022pth , 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. 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 $\pi$ at $E_{{\rm thr},i}$ to the step-up of the
phase shift by $\pi$ at the bound state energy $E_{i}$. In this approximation,
labelled ”SUSD” in the middle panel of Fig. 5, the $\omega$-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_{\rm 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. (LABEL:eq:s) is shown by the red dashed
lines in Fig. 6.

  3. 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_{\rm Mott}$, see the green solid lines in Fig. 6.

Figure 6: Scaled entropy density of hadrons $s_{\rm MHRG}/T^{3}$ as a function
of temperature $T$ calculated for $\mu_{B}/T=0$ (upper panel) and
$\mu_{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 Borsanyi:2021sxv , while the thick solid line represents calculations
within the HRG, see Albright:2014gva .

The simplest model that is in agreement with the Levinson theorem is a SUSD
ansatz Blaschke:2017pvh ; Bastian:2018mmc . It results in the Mott
dissociation of hadrons at high temperatures and chemical potentials

$\displaystyle\delta_{i}(M,T,\mu)$ $\displaystyle=$
$\displaystyle\pi\Theta(M_{{\rm thr},i}-M_{i})$ (61)
$\displaystyle\left\\{\Theta(M-M_{i})-\Theta(M-M_{{\rm thr},i})\right\\}~{}.$

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_{{\rm 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
$\mu=0$ Blaschke:2020lhm . 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\lesssim T[{\rm MeV}]\lesssim 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 state-dependent 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 Wergieluk:2012gd  and extend it
from the original $\pi$ and $\sigma$ channels to include contributions from
each cluster with the appropriate Mott dissociation effect at high
temperatures. This modified phase shift model reads

$\displaystyle\delta_{i}(M)=\delta^{*}_{i}(M)\left[\Theta(M^{2}_{{\rm
thr},i}-M^{2})\right.$ $\displaystyle\left.+\frac{1}{\pi}{\rm
acos}\left(\frac{2M-2M_{{\rm
thr},i}-N_{i}\Lambda_{i}}{N_{i}\Lambda_{i}}\right)\right.$
$\displaystyle\left.\times\Theta(M^{2}-M^{2}_{{\rm thr},i})\Theta((M_{{\rm
thr},i}+\Lambda_{i}N_{i})^{2}-M^{2})\right]$ (62)

with

$\delta^{*}_{i}(M)=\begin{cases}\pi\,\Theta(M^{2}-M_{i}^{2})~{},&T<T_{\rm
Mott}\\\
\pi\frac{\mathrm{atan}\left(\frac{M^{2}-M_{i}^{2}}{\Gamma_{i}M_{i}}\right)-\mathrm{atan}\left(\frac{M_{{\rm
thr},i}^{2}-M_{i}^{2}}{\Gamma_{i}M_{i}}\right)}{\frac{\pi}{2}-\mathrm{atan}\left(\frac{M_{{\rm
thr},i}^{2}-M_{i}^{2}}{\Gamma_{i}M_{i}}\right)}\\\ \times\Theta(M^{2}-M_{{\rm
thr},i}^{2})~{}.&T>T_{\rm Mott}\end{cases}$ (63)

For temperatures above $T_{\rm 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
$\omega=M$) for different temperatures around the Mott transition at $T_{\rm
Mott}=T_{c}=156.5$ MeV.

##  4 Results for QCD thermodynamics

Figure 7: Scaled entropy density $s/T^{3}$ as a function of temperature $T$
calculated at $\mu_{B}/T=0$ (upper panel) and $\mu_{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
Borsanyi:2021sxv  in good agreement with the result of the present model
(thick solid line).  Figure 8: Scaled entropy density $s/T^{3}$ as a function
of temperature $T$ calculated at different values of $\mu_{B}/T$ indicated in
the legend. The shaded regions of the corresponding colors represent the
lattice QCD results Borsanyi:2021sxv . The triangles of corresponding colors
are the values obtained within $\mathcal{O}(\alpha_{s})$ perturbative QCD, see
Tab. 6, and are attained for large $T$. Note the broken temperature axis.
Figure 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.

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 numerical 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 $\mu_{B}/T=3.0$ is presented and presented to the
results of lattice QCD simulations Borsanyi:2021sxv , 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_{\rm
Mott}=T_{c}=156.5$ MeV. For $T<T_{\rm 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. Albright:2014gva , and
added the deuteron as a hexaquark state as well as the hypothetical sexaquark
$S(uuddss)$ which is discussed, e.g., in Shahrbaf:2022upc  and references
therein. At $T=T_{\rm 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\sim 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,
4–quark 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
$\mathcal{O}(\alpha_{s})$ perturbative QCD limits are attained.

Figure 10:  Scaled pressure $p/T^{4}$ as a function of temperature $T$ at
vanishing baryochemical potential $\mu_{B}/T=0$ obtained from integrating the
entropy (42) (blue solid line), in comparison with two sets of lattice QCD
data from Borsanyi et al. (2014) Borsanyi:2013bia  and Borsanyi et al. (2021)
Borsanyi:2021sxv , shown as shaded regions.

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.,
Andronic:2017pug  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\sim 200$ MeV, their entropy fraction in the
hadronic phase below $T_{\rm 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_{\rm Mott}$. The $S(1885)$ contribution is a unique prediction of the
present model. We note that in Ref. Blaschke:2021tul , 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-to-deuteron
ratio is of the order one.

Figure 11: Baryon density as a function of temperature for $\mu_{B}/T=3.0$
(blue solid line) compared to lattice QCD simulations (light blue band)
Borsanyi:2021sxv . 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).

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)
Borsanyi:2013bia  and Borsanyi et al. (2021) Borsanyi:2021sxv  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 $\mu_{B}/T=3.0$ and compare the
resulting total density with the lattice QCD simulation Borsanyi:2021sxv
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_{\rm Mott}$. Fig. 12 demonstrates that the systematics of the chemical
potential dependence is correctly accounted for. In particular, the baryon
density quickly approaches the $\mathcal{O}(\alpha_{s})$ perturbative QCD
values while the lattice QCD data slightly overshoot. At temperatures below
$T_{\rm 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 Bazavov:2014xya . Adding baryonic states from a quark model
description Isgur:1978wd  has been demonstrated to solve this problem
Kaczmarek:2022oiu .

Figure 12: Baryon density as a function of temperature for different values of
$\mu_{B}/T$ compared to lattice QCD simulations Borsanyi:2021sxv . The
triangles of corresponding colors are the values obtained within
$\mathcal{O}(\alpha_{s})$ perturbative QCD, see Tab. 6, which are attained at
large temperature.

##  5 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 low-density limit. This property was taken into
account by using the Polyakov-Nambu-Jona-Lasinio model as an effective
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. Albright:2014gva , 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 $\Phi$. 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 $\Phi$-derivable approach.
When a restriction to closed two-loop diagrams in cluster Green’s functions is
applied to the $\Phi$-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 $\delta_{i}$ is a subject of further investigation. In
addition to the sunset diagrams, other diagrams for the correlation functional
$\Phi$ 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.

### Acknowledgement

D.B. and M.C. acknowledge support from the Polish National Science Centre
(NCN) under grant number
2019/33/B/ST9/03059. The work of O.I. was supported by the program Excellence
Initiative–Research University of the University of Wroc?aw of the Ministry of
Education and Science. The work of G.R. received support via a joint stipend
from the Alexander von Humboldt Foundation and the Foundation for Polish
Science.

##  Appendix 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
single-quark propagators folded with the antisymmetric Gell-Mann matrices
$\lambda_{2}$, $\lambda_{5}$ or $\lambda_{7}$. For a better systematics, we
will replace these matrices by the completely equivalent totally antisymmetric
symbol $\varepsilon_{abc}$ in color space. We introduce the bosonic Matsubara
frequency $\Omega_{m}=\omega_{n}+\omega_{n}^{\prime}$ as the sum of the
fermionic Matsubara frequencies of the quark propagators that constitute the
diquark and obtain the free quark pair propagator.

$\displaystyle S_{2}^{(0)}$ $\displaystyle=$
$\displaystyle\sum_{ab}\sum_{n}\varepsilon_{abc}\frac{1}{(i\omega_{n}-E_{p}+\mu)\mathbbm{1}_{c}-i(A_{4})_{aa}}\cdot$
(64)
$\displaystyle\frac{1}{(i\omega_{n^{\prime}}-E_{p^{\prime}}+\mu)\mathbbm{1}_{c}-i(A_{4})_{bb}}$
$\displaystyle=$
$\displaystyle\sum_{ab}\sum_{n}\frac{\varepsilon_{abc}}{(i\Omega_{m}-E_{p}-E_{p^{\prime}}+2\mu)\mathbbm{1}_{c}-i(A_{4})_{aa}-i(A_{4})_{bb}}$
$\displaystyle\times\left[\frac{1}{(i\omega_{n}-E_{p}+\mu)\mathbbm{1}_{c}-i(A_{4})_{aa}}\right.$
$\displaystyle\left.+\frac{1}{(i(\Omega_{m}-\omega_{n})-E_{p^{\prime}}+\mu)\mathbbm{1}_{c}-i(A_{4})_{bb}}\right]$
$\displaystyle=$
$\displaystyle\frac{-\left.f^{(1),+}_{\phi}\right|_{E_{p}}+\left.\left[f^{(1),-}_{\phi}\right]^{*}\right|_{E_{p^{\prime}}}}{(i\Omega_{m}-E_{p}-E_{p^{\prime}}+2\mu)\mathbbm{1}_{c}+i(A_{4})_{cc}}~{},$
(65)

where $\mathbbm{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}=\varepsilon_{abc}\left[(A_{4})_{aa}+(A_{4})_{bb}\right]$.

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
$\omega=E_{p}+E_{p^{\prime}}$ for the two-particle energy and neglect the
compositeness effect by suppressing the Pauli blocking factor
$-\left.f^{(1),+}_{\phi}\right|_{E_{p}}+\left.\left[f^{(1),-}_{\phi}\right]^{*}\right|_{E_{p^{\prime}}}=1-\left.f^{(1),+}_{\phi}\right|_{E_{p}}-\left.f^{(1),+}_{\phi}\right|_{E_{p^{\prime}}}\longrightarrow
1$. By summation over Matsubara frequencies and color index we then obtain

$\displaystyle f_{\phi}^{(2),+}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\sum_{c=1}^{3}\frac{1}{\beta}\sum_{n}\frac{1}{i\omega_{n}^{(2)}-\omega\mathbbm{1}_{c}}$
(66) $\displaystyle=$
$\displaystyle\frac{1}{3}\sum_{c=1}^{3}\frac{1}{\beta}\frac{\partial\ln\left[1-y_{2}^{+}{\rm
e}^{i\beta(A_{4})_{cc}}\right]}{\partial\omega}$ (68) $\displaystyle=$
$\displaystyle\frac{1}{3\beta}\frac{\partial}{\partial\omega}\ln\left\\{[1-y_{2}^{+}{\rm
e}^{-2i\beta\varphi_{8}}][1-y_{2}^{+}{\rm
e}^{i\beta(\varphi_{8}-\varphi_{3})}]\right.$
$\displaystyle\left.[1-y_{2}^{+}{\rm
e}^{i\beta(\varphi_{8}+\varphi_{3})}]\right\\}$ $\displaystyle=$
$\displaystyle\frac{1}{3\beta}\frac{\partial}{\partial\omega}\ln\left[1-3({\phi}-\bar{\phi}y_{2}^{+})y_{2}^{+}-{y_{2}^{+}}^{3}\right]~{},$
(69)

where $y_{2}^{+}(\omega)=\exp[-\beta(\omega-2\mu)]$ and
$\omega_{n}^{(2)}=(2n\pi T-i2\mu)\mathbbm{1}_{c}+A_{4}$, see Tab. 5. Carrying
out the $\omega-$derivative, we obtain the Polyakov-loop generalized diquark
distribution function

$\displaystyle f_{\phi}^{(2),+}$ $\displaystyle=$
$\displaystyle\frac{({\phi}-2\bar{\phi}y_{2}^{+})y_{2}^{+}+{y_{2}^{+}}^{3}}{1-3({\phi}-\bar{\phi}y_{2}^{+})y_{2}^{+}-{y_{2}^{+}}^{3}}~{}.$
(70)

###  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

$\displaystyle S_{3}^{(0)}$ $\displaystyle=$
$\displaystyle\sum_{ab}\sum_{n}\delta_{ab}\frac{1}{i\omega^{(1)}_{n}-E_{p}\mathbbm{1}_{c}}\cdot\frac{1}{i\omega^{(2)}_{n^{\prime}}-E_{p^{\prime}}\mathbbm{1}_{c}}$
$\displaystyle=$
$\displaystyle\sum_{a}\sum_{n}\frac{1}{(i\omega^{(3)}_{m}-(E_{p}+E_{p^{\prime}})\mathbbm{1}_{c}}\left[\frac{1}{i\omega^{(1)}_{n}-E_{p}\mathbbm{1}_{c}}\right.$
$\displaystyle\left.+\frac{1}{(i(\omega^{(3)}_{m}-\omega^{(1)}_{n})-E_{p^{\prime}}\mathbbm{1}_{c}}\right]$
$\displaystyle=$
$\displaystyle\frac{-\left.f^{(1),+}_{\phi}\right|_{E_{p}}+\left.\left[f^{(2),-}_{\phi}\right]^{*}\right|_{E_{p^{\prime}}}}{i\omega^{(3)}_{m}-(E_{p}+E_{p^{\prime}})\mathbbm{1}_{c}}$
$\displaystyle=$
$\displaystyle\frac{1-\left.f^{(1),+}_{\phi}\right|_{E_{p}}+\left.f^{(2),+}_{\phi}\right|_{E_{p^{\prime}}}}{i\omega^{(3)}_{m}-(E_{p}+E_{p^{\prime}})\mathbbm{1}_{c}}~{}.$
(73)

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 $\omega^{(3)}_{m}=(2m+1)\pi T-i3\mu$ is standard
and gives

$\displaystyle f_{\phi}^{(3)}(\omega)$ $\displaystyle=$ $\displaystyle
f^{(3)}(\omega)=\frac{1}{\beta}\sum_{n}\frac{1}{i\omega_{n}^{(3)}-\omega\mathbbm{1}_{c}}=\frac{y_{3}^{+}}{1+y_{3}^{+}}~{},$
(74)

where $y_{3}^{+}(\omega)=\exp[-\beta(\omega-3\mu)]$ 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 multi-quark distribution functions that are
obtained for multi-quark 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

$\displaystyle f^{(a),\pm}_{\phi}$
$\displaystyle\stackrel{{\scriptstyle\text{(a odd)}}}{{=}}$
$\displaystyle\frac{(\bar{\phi}+2\phi
y_{a}^{\pm})y^{\pm}_{a}+{y_{a}^{\pm}}^{3}}{1+3(\bar{\phi}+\phi
y_{a}^{\pm})y_{a}^{\pm}+{y_{a}^{\pm}}^{3}}~{},$ (75)

which represents the Polyakov–loop–generalized Fermi function, with notable
limiting cases being the single quark distribution function
$f^{(1),\pm}_{\phi}$ and the color singlet Fermi function
$f^{(a),\pm}_{\phi=1}$.

For an even number of net valence quarks, still in the color triplet state, we
obtain

$\displaystyle f^{(a),\pm}_{\phi}$
$\displaystyle\stackrel{{\scriptstyle\text{(a even)}}}{{=}}$
$\displaystyle\frac{(\bar{\phi}-2{\phi}y_{a}^{\pm})y^{\pm}_{a}+{y_{a}^{\pm}}^{3}}{1-3(\bar{\phi}-{\phi}y_{a}^{\pm})y_{a}^{\pm}-{y_{a}^{\pm}}^{3}}~{},$
(76) $\displaystyle~{}y_{a}^{\pm}$ $\displaystyle=$
$\displaystyle\exp[-\beta(\omega\mp a\mu)]~{},$ (77)

which represents the Polyakov–loop–generalized Bose function, with notable
limiting cases being the diquark distribution function
$\left[f^{(2),\pm}_{\phi}\right]^{*}$ (the complex conjugation represents a
diquark in the color–antitriplet state) and the color singlet Bose function
$f^{(a),\pm}_{\phi=1}$.

##  Appendix B High-temperature $\mathcal{O}(g^{2})$ perturbative QCD limit

An important benchmark for the development of an effective model description
that reproduces lattice QCD thermodynamics is the $\mathcal{O}(g^{2})$
perturbative QCD limit at high temperatures. For the readers convenience, we
quote here expressions from Appendix B of Ref. Bollweg:2022fqq , starting with
the pressure at this order, which reads for the three-flavor case ($N_{f}=3$)

$\displaystyle p$ $\displaystyle=$ $\displaystyle p_{\rm id}-g^{2}\,p_{2}~{},$
(78) $\displaystyle p_{\rm id}$ $\displaystyle=$
$\displaystyle\frac{8\pi^{2}}{45}T^{4}\left(1+\frac{63}{32}\right)+\frac{3}{2}\mu^{2}T^{2}+\frac{3}{4\pi^{2}}\mu^{4}~{},$
(79) $\displaystyle p_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{6}T^{4}\left(1+\frac{5}{4}\right)+\frac{3}{4\pi^{2}}\mu^{2}T^{2}+\frac{3}{8\pi^{4}}\mu^{4}~{}.$
(80)

Here, we consider the case of vanishing charge and strangeness chemical
potentials $\mu_{S}=\mu_{Q}=0$, for which
$\mu_{u}=\mu_{d}=\mu_{s}=\mu=\mu_{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,

$\displaystyle n$ $\displaystyle=$ $\displaystyle\frac{\partial
p}{\partial\mu}\bigg{|}_{T={\rm const}}$ (81) $\displaystyle=$ $\displaystyle
3\mu T^{2}+\frac{3}{\pi^{2}}\mu^{3}-g^{2}\left(\frac{3}{2\pi^{2}}\mu
T^{2}+\frac{3}{2\pi^{4}}\mu^{3}\right),$ $\displaystyle s$ $\displaystyle=$
$\displaystyle\frac{\partial p}{\partial T}\bigg{|}_{T={\rm const}}$ (82)
$\displaystyle=$
$\displaystyle\frac{32\pi^{2}}{45}T^{3}\left(1+\frac{63}{32}\right)+3\mu^{2}T$
$\displaystyle-g^{2}\left[\frac{2}{3}T^{3}\left(1+\frac{5}{4}\right)+\frac{3}{2\pi^{2}}\mu^{2}T\right]~{}.$

In table 6, we give the limiting values for the dimensionless baryon density
$n_{B}/T^{3}$ and entropy $s/T^{3}$ to $\mathcal{O}(g^{2})$ for four values of
$\hat{\mu}_{B}=\mu_{B}/T=0$, 1.0, 2.0 and 3.0, for which there are results
from lattice QCD simulations Borsanyi:2021sxv  available.

| $n_{B}/T^{3}=$ | $s/T^{3}=\frac{19\pi^{2}}{9}\left(1-\frac{27g^{2}}{38\pi^{2}}\right)$
---|---|---
$\hat{\mu}_{B}$ | $\left(\frac{\hat{\mu}_{B}}{3}+\frac{\hat{\mu}_{B}^{3}}{27\pi^{2}}\right)\left(1-\frac{g^{2}}{2\pi^{2}}\right)$ | $+\frac{\hat{\mu}_{B}^{2}}{3}\left(1-\frac{g^{2}}{2\pi^{2}}\right)$
0.0 | 0.0 | 17.3072
1.0 | 0.2969 | 17.6008
2.0 | 0.6137 | 18.4816
3.0 | 0.9701 | 19.9497
Table 6: limiting values for the dimensionless baryon density $n_{B}/T^{3}$
and entropy $s/T^{3}$ to $\mathcal{O}(g^{2})$ in perturbative QCD for
$\alpha_{s}=g^{2}/(4\pi)=0.1872$.

These values are approached asymptotically by our model and by the lattice QCD
simulations, as can be seen in Figs. 8 and 12, where they are indicated by
small right arrows of the corresponding colors.

## References

  * (1) M.G. Alford, A. Schmitt, K. Rajagopal, T. Schäfer, Rev. Mod. Phys. 80, 1455 (2008), 0709.4635
  * (2) R. Anglani, R. Casalbuoni, M. Ciminale, N. Ippolito, R. Gatto, M. Mannarelli, M. Ruggieri, Rev. Mod. Phys. 86, 509 (2014), 1302.4264
  * (3) R.D. Pisarski, V.V. Skokov, A.M. Tsvelik, Universe 5(2), 48 (2019), 2202.01036
  * (4) A. Bauswein, D. Blaschke, T. Fischer, _Effects of a strong phase transition on supernova explosions, compact stars and their mergers_ (World Scientific, Singapore, 2022), pp. 283–320, 2203.17188
  * (5) T.E. Riley et al., Astrophys. J. Lett. 887(1), L21 (2019), 1912.05702
  * (6) M.C. Miller et al., Astrophys. J. Lett. 887(1), L24 (2019), 1912.05705
  * (7) T.E. Riley et al., Astrophys. J. Lett. 918(2), L27 (2021), 2105.06980
  * (8) M.C. Miller et al., Astrophys. J. Lett. 918(2), L28 (2021), 2105.06979
  * (9) H.T. Cromartie et al. (NANOGrav), Nature Astron. 4(1), 72 (2019), 1904.06759
  * (10) B.P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 119(16), 161101 (2017), 1710.05832
  * (11) B.P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 121(16), 161101 (2018), 1805.11581
  * (12) B.P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X 9(1), 011001 (2019), 1805.11579
  * (13) A. Bauswein, N.U.F. Bastian, D.B. Blaschke, K. Chatziioannou, J.A. Clark, T. Fischer, M. Oertel, Phys. Rev. Lett. 122(6), 061102 (2019), 1809.01116
  * (14) T. Fischer, N.U. Bastian, D. Blaschke, M. Cierniak, M. Hempel, T. Klähn, G. Martínez-Pinedo, W.G. Newton, G. Röpke, S. Typel, Publ. Astron. Soc. Austral. 34, 67 (2017), 1711.07411
  * (15) S. Borsanyi, Z. Fodor, C. Hoelbling, S.D. Katz, S. Krieg, C. Ratti, K.K. Szabo (Wuppertal-Budapest), JHEP 09, 073 (2010), 1005.3508
  * (16) S. Borsanyi, G. Endrodi, Z. Fodor, S.D. Katz, S. Krieg, C. Ratti, K.K. Szabo, JHEP 08, 053 (2012), 1204.6710
  * (17) A. Bazavov et al., Phys. Rev. Lett. 113(7), 072001 (2014), 1404.6511
  * (18) A. Bazavov et al. (HotQCD), Phys. Lett. B 795, 15 (2019), 1812.08235
  * (19) S. Borsányi, Z. Fodor, J.N. Guenther, R. Kara, S.D. Katz, P. Parotto, A. Pásztor, C. Ratti, K.K. Szabó, Phys. Rev. Lett. 126(23), 232001 (2021), 2102.06660
  * (20) D. Bollweg, D.A. Clarke, J. Goswami, O. Kaczmarek, F. Karsch, S. Mukherjee, P. Petreczky, C. Schmidt, S. Sharma (HotQCD), Phys. Rev. D 108(1), 014510 (2023), 2212.09043
  * (21) D. Bollweg, J. Goswami, O. Kaczmarek, F. Karsch, S. Mukherjee, P. Petreczky, C. Schmidt, P. Scior (HotQCD), Phys. Rev. D 105(7), 074511 (2022), 2202.09184
  * (22) Y. Nambu, G. Jona-Lasinio, Phys. Rev. 122, 345 (1961)
  * (23) Y. Nambu, G. Jona-Lasinio, Phys. Rev. 124, 246 (1961)
  * (24) M. Buballa, Phys. Rept. 407, 205 (2005), hep-ph/0402234
  * (25) K. Fukushima, Phys. Lett. B 591, 277 (2004), hep-ph/0310121
  * (26) C. Ratti, M.A. Thaler, W. Weise, Phys. Rev. D 73, 014019 (2006), hep-ph/0506234
  * (27) S. Roessner, C. Ratti, W. Weise, Phys. Rev. D 75, 034007 (2007), hep-ph/0609281
  * (28) S. Roessner, T. Hell, C. Ratti, W. Weise, Nucl. Phys. A 814, 118 (2008), 0712.3152
  * (29) D. Blaschke, K.A. Devyatyarov, O. Kaczmarek, Symmetry 13(3), 514 (2021), 2012.12894
  * (30) D. Blaschke, A. Dubinin, L. Turko (2016), 1611.09845
  * (31) D. Blaschke, A. Dubinin, L. Turko, Acta Phys. Polon. Supp. 10, 473 (2017), 1612.09556
  * (32) J. Hüfner, S.P. Klevansky, P. Zhuang, H. Voss, Annals Phys. 234, 225 (1994)
  * (33) P. Zhuang, J. Hüfner, S.P. Klevansky, Nucl. Phys. A 576, 525 (1994)
  * (34) A. Wergieluk, D. Blaschke, Y.L. Kalinovsky, A. Friesen, Phys. Part. Nucl. Lett. 10, 660 (2013), 1212.5245
  * (35) K. Yamazaki, T. Matsui, Nucl. Phys. A 913, 19 (2013), 1212.6165
  * (36) D. Blaschke, M. Buballa, A. Dubinin, G. Roepke, D. Zablocki, Annals Phys. 348, 228 (2014), 1305.3907
  * (37) K. Yamazaki, T. Matsui, Nucl. Phys. A 922, 237 (2014), 1310.4960
  * (38) J.M. Torres-Rincon, J. Aichelin (2016), 1601.01706
  * (39) J.M. Torres-Rincon, J. Aichelin, Phys. Rev. C 96(4), 045205 (2017), 1704.07858
  * (40) D. Blaschke, A. Dubinin, A. Radzhabov, A. Wergieluk, Phys. Rev. D 96(9), 094008 (2017), 1608.05383
  * (41) R. Hagedorn, Nuovo Cim. Suppl. 3, 147 (1965)
  * (42) S. Borsanyi, Z. Fodor, S.D. Katz, S. Krieg, C. Ratti, K. Szabo, JHEP 01, 138 (2012), 1112.4416
  * (43) A. Bazavov et al. (HotQCD), Phys. Rev. D 86, 034509 (2012), 1203.0784
  * (44) R. Bellwied, S. Borsanyi, Z. Fodor, S.D. Katz, A. Pasztor, C. Ratti, K.K. Szabo, Phys. Rev. D 92(11), 114505 (2015), 1507.04627
  * (45) R. Bellwied, S. Borsanyi, Z. Fodor, S.D. Katz, C. Ratti, Phys. Rev. Lett. 111, 202302 (2013), 1305.6297
  * (46) A. Andronic, P. Braun-Munzinger, K. Redlich, J. Stachel, Nature 561(7723), 321 (2018), 1710.09425
  * (47) G. Röpke, D. Blaschke, H. Schulz, Phys. Rev. D 34, 3499 (1986)
  * (48) D. Blaschke, H. Grigorian, G. Röpke, Particles 3(2), 477 (2020), 2005.10218
  * (49) G. Röpke, L. Münchow, H. Schulz, Nuclear Physics A 379(3), 536 (1982), ISSN 0375-9474
  * (50) W.D. Kraeft, D. Kremp, W. Ebeling, G. Röpke, _Quantum Statistics of Charged Particle Systems_ (Akademie Verlag, Berlin, 1986), ISBN 978-1461292739
  * (51) C.J. Horowitz, E.J. Moniz, J.W. Negele, Phys. Rev. D 31, 1689 (1985)
  * (52) V.I. Yukalov, E.P. Yukalova, Phys. Part. Nucl. 28, 37 (1997), hep-ph/9709338
  * (53) M. Albright, J. Kapusta, C. Young, Phys. Rev. C 90(2), 024915 (2014), 1404.7540
  * (54) A.S. Khvorostukhin, V.V. Skokov, V.D. Toneev, K. Redlich, Eur. Phys. J. C 48, 531 (2006), nucl-th/0605069
  * (55) M.A.R. Kaltenborn, N.U.F. Bastian, D.B. Blaschke, Phys. Rev. D 96(5), 056024 (2017), 1701.04400
  * (56) O. Ivanytskyi, D. Blaschke, Phys. Rev. D 105(11), 114042 (2022), 2204.03611
  * (57) J.M. Luttinger, J.C. Ward, Phys. Rev. 118, 1417 (1960)
  * (58) G. Baym, L.P. Kadanoff, Phys. Rev. 124, 287 (1961)
  * (59) G. Baym, Phys. Rev. 127, 1391 (1962)
  * (60) S.Y.F. Liu, R. Rapp, Phys. Rev. C 97(3), 034918 (2018), 1711.03282
  * (61) G. Ropke, N.U. Bastian, D. Blaschke, T. Klahn, S. Typel, H.H. Wolter, Nucl. Phys. A 897, 70 (2013), 1209.0212
  * (62) P.M. Lo, B. Friman, O. Kaczmarek, K. Redlich, C. Sasaki, Phys. Rev. D 88, 074502 (2013), 1307.5958
  * (63) B.J. Schaefer, J.M. Pawlowski, J. Wambach, Phys. Rev. D 76, 074023 (2007), 0704.3234
  * (64) O. Ivanytskyi, D. Blaschke, K. Maslov, EPJ Web Conf. 258, 07008 (2022), 2112.09223
  * (65) J.I. Kapusta, _Finite Temperature Field Theory_ , Cambridge Monographs on Mathematical Physics (Cambridge University Press, Cambridge, 1989), ISBN 978-0-521-35155-3
  * (66) L. Turko, D. Blaschke, D. Prorok, J. Berdermann, Acta Phys. Polon. Supp. 5, 485 (2012), 1112.6408
  * (67) D. Blaschke, O. Kaczmarek, E. Laermann, V. Yudichev, Eur. Phys. J. C 43, 81 (2005), hep-ph/0505053
  * (68) D.V. Shirkov, Theor. Math. Phys. 136, 893 (2003), hep-th/0210013
  * (69) A. Peshier, B. Kämpfer, G. Soff, _The QCD equation of state and quark star properties_ , in _NATO Advanced Research Workshop on Superdense QCD Matter and Compact Stars_ (2003), pp. 135–146, hep-ph/0312080
  * (70) G. Röpke, Phys. Rev. C 92(5), 054001 (2015), 1411.4593
  * (71) G. Röpke, Phys. Rev. C 101(6), 064310 (2020), 2004.09773
  * (72) R.L. Workman et al. (Particle Data Group), PTEP 2022, 083C01 (2022)
  * (73) F. Buccella, PoS CORFU2019, 024 (2020)
  * (74) M. Shahrbaf, D. Blaschke, S. Typel, G.R. Farrar, D.E. Alvarez-Castillo, Phys. Rev. D 105(10), 103005 (2022), 2202.00652
  * (75) F. Stancu, AIP Conf. Proc. 775(1), 32 (2005), hep-ph/0504284
  * (76) L. Maiani, A.D. Polosa, V. Riquer, Eur. Phys. J. C 83(5), 378 (2023), 2303.04056
  * (77) B. Vanderheyden, G. Baym, J. Statist. Phys. 93, 843 (1998), hep-ph/9803300
  * (78) J.P. Blaizot, E. Iancu, A. Rebhan, Phys. Rev. D 63, 065003 (2001), hep-ph/0005003
  * (79) N.U.F. Bastian, D. Blaschke, T. Fischer, G. Röpke, Universe 4(6), 67 (2018), 1804.10178
  * (80) N.U.F. Bastian, D.B. Blaschke, Eur. Phys. J. A 57(1), 35 (2021), 1812.11766
  * (81) D. Blaschke, P.C. Tandy, _Mesonic correlations and quark deconfinement_ , in _ECT* International Workshop on Understanding Deconfinement in QCD_ (1999), pp. 218–230, nucl-th/9905067
  * (82) D. Horvatic, D. Blaschke, D. Klabucar, O. Kaczmarek, Phys. Rev. D 84, 016005 (2011), 1012.2113
  * (83) S. Benic, D. Blaschke, M. Buballa, Phys. Rev. D 86, 074002 (2012), 1206.6582
  * (84) J. Jankowski, D. Blaschke, M. Spalinski, Phys. Rev. D 87(10), 105018 (2013), 1212.5521
  * (85) A. Dubinin, D. Blaschke, Y.L. Kalinovsky, Acta Phys. Polon. Supp. 7(1), 215 (2014), 1312.0559
  * (86) D. Blaschke, A. Dubinin, M. Buballa, Phys. Rev. D 91(12), 125040 (2015), 1412.1040
  * (87) D. Blaschke, A. Dubinin, D. Ebert, A.V. Friesen, Phys. Part. Nucl. Lett. 15(3), 230 (2018), 1712.09322
  * (88) S. Borsanyi, Z. Fodor, C. Hoelbling, S.D. Katz, S. Krieg, K.K. Szabo, Phys. Lett. B 730, 99 (2014), 1309.5258
  * (89) D. Blaschke et al., Int. J. Mod. Phys. A 36(25), 2141005 (2021), 2111.03770
  * (90) N. Isgur, G. Karl, Phys. Rev. D 19, 2653 (1979), [Erratum: Phys.Rev.D 23, 817 (1981)]
  * (91) O. Kaczmarek, Lect. Notes Phys. 999, 281 (2022), 2206.15194