††thanks: e-mail<EMAIL_ADDRESS>e-mail:
<EMAIL_ADDRESS>e-mail<EMAIL_ADDRESS>
# Quark cluster expansion model for interpreting finite-T lattice QCD
thermodynamics

David Blaschke  Institute of Theoretical Physics, University of Wroclaw,
50-204 Wroclaw, Poland  Bogoliubov Laboratory of Theoretical Physics, JINR
Dubna, 141980 Dubna, Russia  National Research Nuclear University (Moscow
Engineering Physics Institute), 115409 Moscow, Russia  Kirill A. Devyatyarov
Bogoliubov Laboratory of Theoretical Physics, JINR Dubna, 141980 Dubna, Russia
National Research Nuclear University (Moscow Engineering Physics Institute),
115409 Moscow, Russia  Olaf Kaczmarek  Fakultät für Physik, Universität
Bielefeld, Germany  Key Laboratory of Quark & Lepton Physics (MOE) and
Institute of Particle Physics, Central China Normal University, Wuhan 430079,
China

###### Abstract

In this work, we present a unified approach to the thermodynamics of
hadron–quark–gluon matter at finite temperatures on the basis of a quark
cluster expansion in the form of a generalized Beth–Uhlenbeck approach with a
generic ansatz for the hadronic phase shifts that fulfills the Levinson
theorem. The change in the composition of the system from a hadron resonance
gas to a quark–gluon plasma takes place in the narrow temperature interval of
150–190 MeV, where the Mott dissociation of hadrons is triggered by the
dropping quark mass as a result of the restoration of chiral symmetry. The
deconfinement of quark and gluon degrees of freedom is regulated by the
Polyakov loop variable that signals the breaking of the $Z(3)$ center symmetry
of the color $SU(3)$ group of QCD. We suggest a Polyakov-loop quark–gluon
plasma model with $\mathcal{O}(\alpha_{s})$ virial correction and solve the
stationarity condition of the thermodynamic potential (gap equation) for the
Polyakov loop. The resulting pressure is in excellent agreement with lattice
QCD simulations up to high temperatures.

Polyakov quark–gluon plasma; hadron resonance gas; Beth–Uhlenbeck approach;
lattice QCD thermodynamics

###### pacs:

11.10.Wx,12.38.Gc, 12.38.Mh, 25.75.Nq

††preprint: draft

##  I Introduction

Recently, continuum extrapolated lattice QCD (LQCD) thermodynamics results for
physical quark masses have become available Borsanyi:2010bp ; Bazavov:2011nk ;
Borsanyi:2013bia ; Bazavov:2014pvz . It has now become a major goal to
construct an effective low-energy QCD model that would reproduce these results
in the finite temperature and low chemical potential domain to high accuracy.
Such a description could form a basis for extrapolations to the region of low
temperatures and high baryochemical potentials where the sign problem still
prevents LQCD obtaining benchmark solutions.

QCD Dyson–Schwinger equations Roberts:2000aa ; Fischer:2018sdj  and functional
renormalization group (FRG) methods Tripolt:2013jra ; Gao:2020qsj  are
promising tools to investigate the nonperturbative aspects of the QCD phase
diagram in the vicinity of the chiral and deconfinement transitions from first
principles. However, to date, they do not self-consistently account for the
bound state formation. Therefore, despite a satisfactory description of the
temperature dependence of the renormalized light chiral condensate Gao:2020qsj
, a question arises regarding the relationship to the chiral perturbation
theory limit, where the onset of chiral condensate melting is solely due to
light pseudoscalar meson excitations in the medium Gerber:1988tt . The account
for bound states in the FRG approach to low-energy QCD is still under
development Alkofer:2018guy  and not yet applicable to address the role of
hadronic bound states for chiral condensate melting.

Complementary to these first-principle approaches is the Polyakov–quark–meson
(PQM) model (see, e.g., Herbst:2013ufa  and references therein), which is
suitable to study the interrelation of chiral and deconfinement transitions.
Both transitions are closely correlated, but in comparison with LQCD, the
chiral restoration temperature result is too high. It could be suspected that
the inclusion of further mesonic and baryonic degrees of freedom could improve
the situation. Furthermore, the dynamic quark self-energy effects due to
dressing by the meson cloud should be taken into account. In the PQM model,
mesons are not composites but elementary degrees of freedom.

In order to properly account for the composite nature of mesons and baryons,
as well as their Mott dissociation in a hot and dense medium, a Beth–Uhlenbeck
approach can be employed Hufner:1994ma ; Zhuang:1994dw ; Wergieluk:2012gd ;
Yamazaki:2012ux ; Blaschke:2013zaa ; Yamazaki:2013yua ; Torres-Rincon:2016ahl
; Torres-Rincon:2017zbr ; Dubinin:2016wvt ; Blaschke:2016fdh ;
Blaschke:2016hzu . This approach has been limited to the application of low-
lying mesons only and did not take into account perturbative corrections,
meaning that a quantitative description of LQCD thermodynamics has not yet
been possible.

To overcome these limitations, the generic behavior of the scattering phase
shifts in the hadronic channels has been constructed in the spirit of a
cluster expansion model which reproduces the full hadron resonance gas at low
temperatures and the quark–gluon plasma (QGP) with $\mathcal{O}(\alpha_{s})$
virial corrections at high temperatures Blaschke:2016fdh ; Blaschke:2016hzu .
The model embodied the main consequence of chiral symmetry restoration in the
quark sector: the lowering of the thresholds for the two and three-quark
scattering state continuous spectrum, which triggers the transformation of
hadronic bound states to resonances in the scattering continuum. In the early
version of this model, the gap equation for the Polyakov-loop was incomplete,
and the model of the phase shifts was a rather complicated one.

We suggest here a Polyakov-loop quark–gluon plasma model with
$\mathcal{O}(\alpha_{s})$ virial correction in order to obtain a satisfactory
agreement with lattice QCD simulations up to high temperatures and solve the
complete stationarity condition of the thermodynamic potential (gap equation)
for the Polyakov loop. The phase shift model employed in this work is simpler
than in previous works and universally applicable for all hadronic species. It
is in accordance with the Levinson theorem and results in the vanishing of
hadronic contributions to the thermodynamics at high temperatures.

##  II Cluster Virial Expansion to Quark-Hadron Matter

The main idea behind unifying the description of the quark–gluon plasma (QGP)
and the hadron resonance gas (HRG) phase of low-energy QCD matter is the fact
that hadrons are strong, nonperturbative correlations of quarks and gluons. In
particular, mesons and baryons are bound states (clusters) of quarks and
should therefore emerge in a cluster expansion of interacting quark matter as
new, collective degrees of freedom.

For the total thermodynamic potential of the model, from which all other
equations of state can be derived, we make the following ansatz:

$\displaystyle\Omega_{\rm total}(T;\phi)=\Omega_{\rm QGP}(T;\phi)+\Omega_{\rm
MHRG}(T)~{},$ (1)

where $\Omega_{\rm QGP}(T;\phi)=\Omega_{\rm PNJL}(T;\phi)+\Omega_{\rm
pert}(T;\phi)$ describes the thermodynamic potential of the quark and gluon
degrees of freedom with a perturbative part $\Omega_{\rm pert}(T;\phi)$ and a
nonperturbative mean field part $\Omega_{\rm PNJL}(T;\phi)=\Omega_{\rm
Q}(T;\phi)+\mathcal{U}(T;\phi)$ that can be decomposed into the quark
quasiparticle contribution $\Omega_{\rm Q}(\phi;T)$ and the gluon contribution
that is approximated by a mean field potential $\mathcal{U}(T;\phi)$. Note
that all these contributions to the QGP thermodynamic potential are
intertwined by the traced Polyakov loop $\phi$ as the order parameter for
confinement. The correlations beyond the mean field approximation which
correspond to the hadronic bound states and their scattering state continuum
are described by the Mott–HRG pressure $P_{\rm MHRG}(T)$. This is an HRG
pressure that takes into account the dissociation of hadrons by the Mott
effect, when their masses would exceed the mass of the corresponding continuum
of unbound quark states. A detailed description and numerical evaluation of
these contributions is given below.

###  II.1 Beth–Uhlenbeck Model for HRG with Mott Dissociation

For the MHRG part of the pressure of the model, we have $P_{\rm
MHRG}(T)=-\Omega_{\rm MHRG}(T)$,

$\displaystyle P_{\rm MHRG}(T)=\sum_{i=M,B}{P_{i}(T)}~{},$ (2)

where the sum extends over all mesonic (M) and baryonic (B) states from the
particle data group (PDG), comprising an ideal mixture of hadronic bound and
scattering states in the channel $i$ that are described by a Beth–Uhlenbeck
formula. Then, the partial pressure of the hadron species $i$ is

$P_{i}(T)=\mp
d_{i}\int_{0}^{\infty}\frac{dp~{}p^{2}}{2\pi^{2}}\int_{0}^{\infty}\frac{dM}{\pi}~{}T\ln\left(1\mp{\rm
e}^{-\sqrt{p^{2}+M^{2}}/T}\right)\frac{d\delta_{i}(M;T)}{dM}~{},$ (3)

where $d_{i}$ is the degeneracy factor. For the phase shift of the bound
states of $N_{i}$ quarks in the hadron $i$, we adopt the simple model that is
in accordance with the Levinson theorem:

$\delta_{i}(M;T)=\pi\left[\Theta(M-M_{i})-\Theta(M-M_{{\rm
thr},i}(T))\right]\Theta(M_{{\rm thr},i}(T)-M_{i}).$ (4)

Inserting (4) into (3) results in

$P_{i}(T)=\mp
d_{i}\int_{0}^{\infty}\frac{dp~{}p^{2}}{2\pi^{2}}~{}T\left[\ln\left(1\mp{\rm
e}^{-\sqrt{p^{2}+M_{i}^{2}}/T}\right)-\ln\left(1\mp{\rm
e}^{-\sqrt{p^{2}+M_{{\rm thr},i}(T)^{2}}/T}\right)\right]~{}\Theta(M_{{\rm
thr},i}(T)-M_{i}).$ (5)

The temperature dependent threshold mass of the two (or three) quark continuum
for mesonic (baryonic) bound state channels $i$ is

$M_{{\rm thr},i}(T)=\sqrt{2}\left[(N_{i}-N_{s})m(T)+N_{s}m_{s}(T)\right]~{},$
(6)

where $N_{s}=0,1,...,N_{i}$ is the number of strange quarks in hadron $i$. The
factor $\sqrt{2}$ originates from quark confinement in the following way. In
the confining vacuum, the quarks are not simple plane waves with an
arbitrarily long wavelength, but due to the presence of bag-like boundary
conditions, their wavelength shall not exceed a certain length scale.
Therefore, a minimal quark momentum applies to the quark dispersion relations
$E_{q,{\rm min}}(T)=\sqrt{m_{q}^{2}(T)+p^{2}_{q,{\rm min}}}$, which for the
choice $p_{q,{\rm min}}=m_{q}(T)$ results in $E_{q,{\rm
min}}(T)=\sqrt{2}m_{q}(T)$. For details, see Dubinin:2013yga . The chiral
condensate is defined as

$\langle\bar{\psi}\psi\rangle_{q,T}=-\frac{\partial\Omega(T)}{\partial
m_{q}}~{},~{}~{}q=u,d,s,$ (7)

where $m_{l}$ ($m_{s}$) is the current-quark mass in the light (strange) quark
sector, $l=u,d$. It is an order parameter for the dynamical breaking of the
chiral symmetry that is reflected in the corresponding temperature dependence
of the dynamical quark masses $m_{q}(T)$.

In our present model, we do not treat the dynamical quark mass as an order
parameter that should follow from the solution of an equation of motion (gap
equation) that minimizes the thermodynamic potential, as in the case of the
Polyakov-loop variable $\phi$, but we use the quantity $\Delta_{l,s}(T)$ from
simulations of 2 + 1 flavor lattice QCD as an input. This quantity has been
introduced in Cheng:2007jq  with the definition

$\Delta_{l,s}(T)=\frac{\langle\bar{\psi}\psi\rangle_{l,T}-(m_{l}/m_{s})\langle\bar{\psi}\psi\rangle_{s,T}}{\langle\bar{\psi}\psi\rangle_{l,0}-(m_{l}/m_{s})\langle\bar{\psi}\psi\rangle_{s,0}}~{},$
(8)

and was used later on, e.g., in Borsanyi:2010bp ; Bazavov:2011nk . Further, we
assume the following for the temperature-dependent light quark mass,

$m(T)=m(0)\Delta_{l,s}(T)+m_{l}~{},$ (9)

where $m_{l}=5.5$ MeV is the current-quark mass, and for the strange quark
mass, we adopt

$m_{s}(T)=m(T)+m_{s}-m_{l}=m(0)\Delta_{l,s}(T)+m_{s}~{},$ (10)

with $m_{s}=100$ MeV. The LQCD result for the temperature dependence of the
chiral condensate Borsanyi:2010bp ; Bazavov:2011nk  can be fitted by

$\Delta_{l,s}(T)=\frac{1}{2}\left[1-\tanh\left(\frac{T-T_{c}}{\delta_{T}}\right)\right]~{},$
(11)

where $T_{c}=154$ MeV is the common pseudocritical temperature of the chiral
restoration transition of both LQCD Collaborations and $\delta_{T}=26$ MeV is
its width for the data from Borsanyi:2010bp , while $\delta_{T}=22.7$ MeV for
those from Bazavov:2011nk , see Figure 1. For our present applications
modeling QCD thermodynamics, we use the fit of the chiral condensate (11), but
with the modern value of $T_{c}=156.5\pm 1.5$ MeV Bazavov:2018mes . We have
checked that the results for the total pressure of our model are practically
inert against a changing value of $\delta_{T}$ within the above range of
variation. Inserting (9) and (10) into (6), we get

$\displaystyle M_{{\rm thr},i}(T)$ $\displaystyle=$
$\displaystyle\sqrt{2}\left[N_{i}m(T)+N_{s}(m_{s}(T)-m(T))\right]$ (12)
$\displaystyle=$
$\displaystyle\sqrt{2}\left[m_{s}N_{s}+m_{l}(N_{i}-N_{s})+m(0)N_{i}\Delta_{l,s}(T)\right]~{},$

and using (9) results in

$M_{{\rm
thr},i}(T)=\sqrt{2}\left\\{m_{s}N_{s}+m_{l}(N_{i}-N_{s})+m(0)N_{i}\left[\frac{1}{2}-\frac{1}{2}\tanh\left(\frac{T-T_{c}}{\delta_{T}}\right)\right]\right\\}~{}.$
(13)

Figure 1: Comparison of the fit (11) for the temperature dependence of the
chiral condensate $\Delta_{l,s}(T)$ and the LQCD data for it from the
Wuppertal–Budapest Collaboration Borsanyi:2010bp  and the hotQCD Collaboration
Bazavov:2011nk .

In Fig. 2 we show the pressure (2) as a function of temperature for the hadron
resonance gas (HRG) model with stable hadrons (red line) and for the HRG with
Mott dissociation of hadrons (MHRG) according to the simple phase shift model
(4) employed in the present work. These results are compared to the LQCD data
from the HotQCD Collaboration Bazavov:2014pvz  (green band) and the
Wuppertal–Budapest Collaboration Borsanyi:2013bia  (blue band). We want to
point out that due to the Mott dissociation effect hadrons completely vanish
from the system at $T\approx 190$ MeV while the pressure of the ideal HRG
model is misleadingly still in perfect agreement with LQCD data at this
temperature!

Figure 2: Pressure as a function of temperature for the HRG model with stable
hadrons (red line) and for the MHRG model with Mott dissociation of hadrons
according to the simple phase shift model (4) employed in the present work.
These results are compared to the LQCD data from the HotQCD Collaboration
Bazavov:2014pvz  (green band) and the Wuppertal–Budapest Collaboration
Borsanyi:2013bia  (blue band).

###  II.2 Polyakov-Loop Improved Nambu–Jona-Lasinio (PNJL) Model

The underlying quark and gluon thermodynamics are divided into a perturbative
contribution $\Omega_{\rm pert}(T)$ which is treated as virial correction in
two-loop order following Turko:2011gw  and a nonperturbative part described
within a PNJL model in the form

$P_{\rm PNJL}(T;\phi)=P_{Q}(T;\phi)+\mathcal{U}(T;\phi)~{},$ (14)

where the quark quasiparticle contribution is given by

$P_{Q}(T;\phi)=4N_{c}\sum_{q=u,d,s}\int\frac{dp~{}p^{2}}{2\pi^{2}}\frac{T}{3}\ln\left[1+3\phi(1+Y_{q})Y_{q}+Y_{q}^{3}\right],~{}Y_{q}=e^{-\sqrt{p^{2}+m_{q}^{2}(T)}/T},$
(15)

and the Polyakov-loop potential $\mathcal{U}(T;\phi)$ takes into account the
nonperturbative gluon background in a meanfield approximation using the
polynomial fit of Ratti:2005jh

$\mathcal{U}(T;\phi)=\frac{b_{2}(T)}{2}\phi^{2}+\frac{b_{3}}{3}\phi^{3}-\frac{b_{4}}{4}\phi^{4}~{},$
(16)

where the temperature-dependent coefficient $b_{2}(T)$ is given by

$b_{2}(T)=a_{0}+a_{1}\left(\frac{T_{0}}{T}\right)+a_{2}\left(\frac{T_{0}}{T}\right)^{2}+a_{3}\left(\frac{T_{0}}{T}\right)^{3}~{},$
(17)

and the coefficients are given in Table 1.

Table 1: Set of values for the Polyakov-loop potential $\mathcal{U}(T;\phi)$ Ratti:2005jh .  $a_{0}$ |  $a_{1}$ |  $a_{2}$ |  $a_{3}$ |  $b_{3}$ |  $b_{4}$
---|---|---|---|---|---
6.75 |  $-$1.95 |  2.625 |  $-$7.44 |  0.75 |  7.5

###  II.3 Perturbative Contribution

It is well known that the lattice QCD thermodynamics at high temperatures of
$T\sim 1$ GeV follow a Stefan–Boltzmann like behavior $\propto T^{4}$, but
with a $15-20\%$ reduction of the effective number of degrees of freedom. It
has been observed, e.g., in Turko:2011gw , that this deviation can be
described by the virial correction to the pressure due to the quark–gluon
scattering at $\mathcal{O}(\alpha_{s})$ shown in Figure 3. Here, we modify the
standard expression Kapusta:1989tk  of the form

$\Omega_{\rm
pert}(T;\phi)=-\frac{8}{\pi}\alpha_{s}T^{4}\left[I(T;\phi)+\frac{3}{\pi^{2}}(I(T;\phi))^{2}\right]$
(18)

by introducing the modified integral

$I(T;\phi)=\int_{\Lambda/T}^{\infty}{\rm d}x~{}xf_{\phi}(x),$ (19)

where the generalized Fermi distribution function of the PNJL model for the
case of vanishing quark chemical potential considered here is defined as

$f_{\phi}(x)=[\phi(1+2Y)Y+Y^{3}]/[1+3\phi(1+Y)Y+Y^{3}],~{}Y=\exp(-x)$ (20)

and $\Lambda=m_{l}(T)$ is the momentum range below which nonperturbative
physics dominates and is accounted for by the dynamically generated quark
mass. Here, we use a temperature-dependent, regularized running coupling
Peshier:2003ah ; Blaschke:2005jg ; Shirkov:2002td

$\displaystyle\alpha_{s}=\frac{g^{2}}{4\pi}=\frac{12\pi}{11N_{c}-2N_{f}}\left(\frac{1}{\ln(r^{2}/c^{2})}-\frac{c^{2}}{r^{2}-c^{2}}\right),$
(21)

where $r=3.2T$, $c=350$ MeV and $N_{c}=N_{f}=3$.

Figure 3: Two-loop diagram for the contribution of the one-gluon exchange
interaction to the thermodynamic potential of quark matter.

##  III Stationarity Condition for the Polyakov Loop

The pressure follows from the thermodynamic potential under the condition of
stationarity with regard to variations of the order parameters. Since the
chiral condensate is fixed by the fit (11) to the numerical result from
lattice QCD, the Polyakov loop $\phi$ is the only free-order parameter in the
system to be varied; this condition means

$P_{\rm QGP}(T)=-\min_{\phi}\left\\{\Omega_{\rm QGP}(T;\phi)\right\\}~{}.$
(22)

This is realized by demanding

$\frac{d\Omega_{\rm
QGP}(T;\phi)}{d\phi}=\frac{dU(T;\phi)}{d\phi}+\frac{d\Omega_{Q}(T;\phi)}{d\phi}+\frac{d\Omega_{\rm
pert}(T;\phi)}{d\phi}=0~{},$ (23)

where the separate contributions come from the variations of the Polyakov loop
potential

$\frac{dU(T;\phi)}{d\phi}=b_{2}(T)\phi+b_{3}\phi^{2}-b_{4}\phi^{3}~{},$ (24)

and the quark quasiparticle pressure

$\frac{d\Omega_{Q}(T;\phi)}{d\phi}=4N_{c}\sum_{q=u,d,s}\int\frac{dp~{}p^{2}}{2\pi^{2}}\frac{(1+Y_{q})Y_{q}}{1+3\phi(1+Y_{q})Y_{q}+Y_{q}^{3}}~{},$
(25)

with $Y_{q}=\exp[-\sqrt{p^{2}+m_{q}^{2}(T)}/T]$, and the
$\mathcal{O}(\alpha_{s})$ quark loop contribution

$\frac{d\Omega_{\rm
pert}(T;\phi)}{d\phi}=-\frac{8}{\pi}\alpha_{s}T^{4}\left[\frac{dI(\phi,T)}{d\phi}+\frac{6}{\pi^{2}}I(\phi,T)\frac{dI(\phi,T)}{d\phi}\right]~{},$
(26)

where

$\frac{dI(T;\phi)}{d\phi}=\int_{\Lambda/T}^{\infty}{\rm
d}x~{}x\frac{df_{\phi}(x)}{d\phi},$ (27)

and

$\frac{df_{\phi}(x)}{d\phi}=\frac{Y+2Y^{2}-2Y^{4}-Y^{5}}{\left(1+3\phi(1+Y)Y+Y^{3}\right)^{2}}=\frac{(1+2Y)Y-(2+Y)Y^{4}}{\left(1+3\phi(1+Y)Y+Y^{3}\right)^{2}}~{},~{}Y=\exp(-x)~{}.$
(28)

The equation resulting from the stationarity condition (23) can be dubbed a
“gap equation” for $\phi$ since it has a similar structure to the quark mass
gap equation, known from Nambu–Jona-Lasinio models. In previous work
Blaschke:2016fdh ; Blaschke:2016hzu , the contribution to this gap equation
from the $\mathcal{O}(\alpha_{s})$ quark loop diagram was omitted. Since this
perturbative contribution is calculated with the Polyakov-loop generalized
quark distribution functions $f_{\phi}$, it has to be included to the
generalized gap equation for the traced Polyakov loop. This has been done in
the present paper for the first time.

Moreover, an infrared cutoff is placed at the loop integrals in the
perturbative contribution which is set at the medium-dependent quark mass.
Therefore, the full perturbative contribution in accordance with
Kapusta:1989tk  is restored only at high temperatures, where $m(T)\to m_{l}$
and $\phi\approx 1$, while in the vicinity of the chiral and deconfinement
transition, the effects of both the quark mass and Polyakov loop are taken
into account. The solution of this gap equation gives the temperature
dependence of the traced Polyakov loop $\phi$. This is discussed in the next
section.

##  IV Results

###  IV.1 Polyakov Loop

We performed the numerical solution of the gap Equation (23) for the traced
Polyakov loop as a function of the temperature which enters via the
coefficient $b_{2}(T)$ of the Polyakov-loop potential (16) and the Boltzmann
factors $Y_{q}$ and $Y$ of the distribution functions in the integrals (26)
and (27). The result is shown in Figure 4 along with a comparison to the
lattice QCD data for the renormalized Polyakov loop from the TUMQCD
Collaboration Bazavov:2016uvm  and the Wuppertal–Budapest Collaboration
Borsanyi:2010bp .

Figure 4: The traced Polyakov loop $\phi$ from the solution of the
stationarity condition (23) on the thermodynamical potential as a function of
temperature (magenta solid line) compared with the lattice results for the
renormalized Polyakov loop the TU Munich QCD (TUMQCD) Collaboration
Bazavov:2016uvm  (green band) and the Wuppertal–Budapest Collaboration
Borsanyi:2010bp  (blue symbols).

In Figure 5, we compare the Polyakov-loop susceptibility $d\phi/dT$ with the
chiral susceptibility $d\Delta_{l,s}/dT$ and obtain a strong synchronization
effect of the chiral and Polyakov-loop crossover transitions. This is
demonstrated by the almost coincident vertical lines indicating the peak
positions of these transitions at $T_{\chi}=156.5$ MeV and $T_{\phi}=159.0$
MeV, respectively.

Figure 5: The temperature derivatives of the chiral condensate (chiral
susceptibility $d\Delta_{l,s}/dT$, red solid line) and of the Polyakov loop
(Polyakov-loop susceptibility $d\phi/dT$, ) as functions of temperature. The
vertical lines indicate their almost coincident peak positions at
$T_{\chi}=156.5$ MeV and $T_{\phi}=159.0$ MeV, respectively.

###  IV.2 Pressure

The main result of this work is a unified approach to the pressure of
hadron–quark–gluon matter at finite temperatures that is in excellent
agreement with lattice QCD thermodynamics (see Figure 6). The nontrivial
achievement of the presented approach is that the Mott dissociation of the
hadrons described by the MHRG model pressure conspires with the quark–gluon
pressure described by the Polyakov-loop quark–gluon model with
$\mathcal{O}(\alpha_{s})$ corrections in such a way that the resulting
pressure as a function of temperature yields a smooth crossover behavior. By
virtue of the Polyakov-loop-improved perturbative correction, the agreement
with the lattice QCD thermodynamics extends to the high temperatures of
$T=1960$ MeV reported in Bazavov:2017dsy ; see Figure 7.

Figure 6: The temperature dependence of the total scaled pressure (red solid
line) and it’s constituents: MHRG (coral dotted line), quark (dashed magenta
line), Polyakov-loop potential $\mathcal{U}(T;\phi)$ (dash–dotted green line)
and perturbative QCD contribution (dash-dotted blue line) compared to the
lattice QCD data: HotQCD Collaboration Bazavov:2014pvz  (green band) and
Wuppertal–Budapest Collaboration Borsanyi:2013bia  (blue band).  Figure 7: The
temperature dependence of the total scaled pressure (red solid line) and its
constituents: MHRG (coral dotted line), quark (dashed magenta line), Polyakov-
loop potential $U(\phi,T)$ (dash-dotted green line) and perturbative QCD
contribution (dash-dotted blue line) compared to the lattice QCD data: HotQCD
Collaboration Bazavov:2014pvz  (green band) and Wuppertal–Budapest
Collaboration Borsanyi:2013bia  (blue band), and the high-temperature result
Bazavov:2017dsy  (magenta band).

###  IV.3 Quark Number Susceptibilities

In the present work, we did not consider the generalization of the approach to
finite chemical potentials which would then allow us to evaluate the
(generalized) susceptibilities as derivatives of pressure with respect to the
corresponding chemical potential in appropriate orders. On that basis, ratios
of susceptibilities could be formed as they indicate different aspects of the
QCD transition between the limiting cases of a HRG and a QGP. Here, we would
like to discuss, as an insight into these extensions of the approach, one of
the simplest susceptibility ratios, namely the dimensionless ratio of quark
number density to quark number susceptibility:

$R_{12}(T)=\frac{n_{q}(T)}{\mu_{q}\,\chi_{q}(T)}\bigg{|}_{\mu_{q}=0}~{},$ (29)

where $n_{q}(T)=\partial P(T,\mu_{q})/\partial\mu_{q}|_{\mu_{q}=0}$ and
$\chi_{q}(T)=\partial^{2}P(T,\mu_{q})/\partial\mu_{q}^{2}|_{\mu_{q}=0}$. This
ratio (29) has two well-known limits. At low temperatures, in the hadron
resonance gas phase, it is given by

$R_{12}^{HRG}(T)=\frac{T}{3\mu_{q}}\tanh\left(\frac{3\mu_{q}}{T}\right)~{},$
(30)

while in the QGP phase for massless quarks it approaches

$R_{12}^{QGP}(T)=\frac{1+(1/\pi^{2})(\mu_{q}/T)^{2}}{1+(3/\pi^{2})(\mu_{q}/T)^{2}}~{}.$
(31)

An evaluation of (29) for the present model for the QCD pressure would require
its extension to a finite $\mu_{q}$, which we will perform in a subsequent
work. In the present model, we used our knowledge of the composition as a
function of temperature to define a proxy for (29) by interpolating between
the two known limits (30) and (31) with the partial pressure of the HRG,
$x_{HRG}(T)=P_{\rm MHRG}(T)/P_{\rm tot}(T)$, as

$R_{12}(T)=x_{HRG}(T)R_{12}^{HRG}(T)+[1-x_{HRG}(T)]R_{12}^{QGP}(T)~{}.$ (32)

The result is shown in Figure 8 for two values of $\mu_{q}/T$, for which
lattice QCD results in the two-flavor case PhysRevD.71.054508  are shown for a
comparison. The fact that the present approach reproduces the transition
between HRG and QGP asymptotics well in the narrow range of temperatures
$150~{}{\rm MeV}\lesssim T\lesssim 190~{}{\rm MeV}$ is a nontrivial result. In
previous effective approaches to describe the finite-temperature lattice QCD
thermodynamics results based on a spectral broadening of the HRG states, the
transition to the QGP asymptotics occurred at a much higher temperature
$250~{}{\rm MeV}\lesssim T\lesssim 400~{}{\rm MeV}$ Biro:2014sfa  or never
Turko:2011gw ; Turko:2013taa ; Turko:2014jta . In the latter case, the QGP
asymptotic behavior is mimicked by an appropriate number of unaffected low-
lying hadronic degrees of freedom.

Figure 8: The dimensionless ratio of quark number density to quark number
susceptibility $R_{12}(T)=n_{q}(T)/(\mu_{q}\chi_{q}(T))|_{\mu_{q}=0}$ as a
function of temperature for $\mu_{q}/T=0.4$ (red solid line) and
$\mu_{q}/T=0.8$ (blue dash-dotted line) compared to the lattice QCD data
PhysRevD.71.054508  $\mu_{q}/T=0.4$ (red band), $\mu_{q}/T=0.8$ (blue band).
For details, see text.

##  V Discussion and Conclusions

The main result of the present work is a unified approach to the thermodynamic
potential of hadron–quark–gluon matter at finite temperatures that is in
excellent agreement with lattice QCD thermodynamics on the temperature axis of
the QCD phase diagram. The key aspect to this approach is the quark cluster
decomposition of the thermodynamic potential within the Beth–Uhlenbeck
approach Bastian:2018wfl , which allowed us to implement the effect of Mott
dissociation to the hadron resonance gas phase of low-temperature/low-density
QCD. The MHRG model description includes, in principle, the information about
the spectral properties of all hadronic channels with their discrete and
continuous part of the spectrum, encoded in the hadronic phase shifts. Instead
of solving the equations of motion with a coupled hierarchy of Schwinger–Dyson
equations in the one, two and many-quark channels self-consistently (a
formidable task of finite-temperature quantum field theory), we applied a
schematic model for the in-medium phase shifts that was in accordance with the
Levinson theorem and sufficiently general to be applicable for all multiquark
cluster channels. This phase shift model requires only the knowledge of the
vacuum mass spectrum which can come from the particle data group tables or
from relativistic quark models and the medium dependence of the multi-quark
continuum threshold.

The latter requires the knowledge of the quark mass (i.e., the chiral
condensate) with its medium dependence as an order parameter of the chiral
symmetry breaking and restoration. Since a quark mean field model of the
(P)NJL type is not sufficient as it lacks the back-reaction from the hadron
resonance gas on the quark propagator properties, we employ here the chiral
condensate measured in continuum-extrapolated, full lattice QCD with physical
current quark masses as an input. This procedure restricts the applicability
of the present model to small chemical potentials only, where lattice QCD data
for the chiral condensate are available. In a further development of the
model, a beyond-mean-field derivation of the quark self-energy will be given.
Furthermore, at the same level of approximation, the corresponding sunset-type
diagrams for the $\Phi$ functional of the 2PI approach should be derived and
evaluated. This would allow us to calculate the generalized polarization-loop
integrals which determine the analytic properties of the multi-quark states.
These can be equivalently encoded in the corresponding medium-dependent phase
shifts of the generalized Beth–Uhlenbeck approach, as has been demonstrated in
particular examples for pions, diquarks Blaschke:2013zaa ; Blaschke:2014zsa
and nucleons Blaschke:2015sla  within the Polyakov-loop generalized NJL model.

Another important aspect of the present approach is that it leads to a
relativistic density functional theory for QCD matter in the QCD phase
diagram, with the known limits of the HRG and pQCD manifestly implemented.
Such an approach allows us to predict the existence and location of critical
endpoints in the QCD phase diagram, as has been demonstrated, e.g., in
Bastian:2018mmc , where a dependence on a free parameter could have—besides
the critical endpoint of the liquid–gas transition in the nuclear matter
phase—another endpoint for the deconfinement transition or none. This
“crossover all over” case of the QCD phase diagram is impossible to address
with two-phase approaches that use a Maxwell construction for the phase
transition. Other models that are in use for the analysis of the critical
behavior of QCD (see, e.g., Plumberg:2018fxo ; Parotto:2018pwx ) impose this
by assuming a so-called “switch function” between HRG and QGP phases. They are
valuable tools but do not have predictive power.

###### Acknowledgements.

D.B. acknowledges discussions with Konstantin Maslov, Krzysztof Redlich and
Ludwik Turko about this work.

Author contributions: Conceptualization, D.B.; Data curation, K.A.D.; Formal
analysis, D.B. and K.A.D.; Funding acquisition, D.B.; Investigation, K.A.D.;
Software, K.A.D.; Supervision, D.B. and O.K.; Validation, O.K.; Writing
original draft, D.B.; Writing review and editing, D.B. and O.K. All authors
have read and agreed to the published version of the manuscript.

Funding: The research of D.B. was supported by the Russian Fund for Basic
research under grant number 18-02-40137 and from the National Research Nuclear
University (MEPhI) in the framework of the Russian Academic Excellence Project
under contract number 02.a03.21.0005.

Conflicts of interest: The authors declare no conflict of interest.

## References

  * (1) Borsányi, S.; Fodor, Z.; Hoelbling, C.; Katz, S.D.; Krieg, S.; Ratti, C.; Szabó, K.K.; Wuppertal-Budapest Collaboration. Is there still any Tc mystery in lattice QCD? Results with physical masses in the continuum limit III. _J. High Energy Phys._ 2010, _1009_ , 73. doi:10.1007/JHEP09(2010)073.
  * (2) Bazavov, A.; Bhattacharya, T.; Cheng, M.; DeTar, C.; Ding, H.T.; Gottlieb, S.; Gupta, R.; Hegde, P.; Heller, U.M.; Karsch, F.; et al. The chiral and deconfinement aspects of the QCD transition. _Phys. Rev. D_ 2012, _85_ , 054503. doi:10.1103/PhysRevD.85.054503.
  * (3) Borsanyi, S.; Fodor, Z.; Hoelbling, C.; Katz, S.D.; Krieg, S.; Szabo, K.K. Full result for the QCD equation of state with 2+1 flavors. _Phys. Lett. B_ 2014, _730_ , 99. doi:10.1016/j.physletb.2014.01.007.
  * (4) Bazavov, A.; Bhattacharya, T.; DeTar, C.; Ding, H.T.; Gottlieb, S.; Gupta, R.; Hegde, P.; Heller, U.M.; Karsch, F.; HotQCD Collaboration; et al. Equation of state in (2 + 1)-flavor QCD. _Phys. Rev. D_ 2014, _90_ , 094503. doi:10.1103/PhysRevD.90.094503.
  * (5) Roberts, C.D.; Schmidt, S.M. Dyson-Schwinger equations: Density, temperature and continuum strong QCD. _Prog. Part. Nucl. Phys._ 2000, _45_ , S1–S103. doi:10.1016/S0146-6410(00)90011-5.
  * (6) Fischer, C.S. QCD at finite temperature and chemical potential from Dyson–Schwinger equations. _Prog. Part. Nucl. Phys._ 2019, _105_ , 1–60. doi:10.1016/j.ppnp.2019.01.002.
  * (7) Tripolt, R.A.; Strodthoff, N.; von Smekal, L.; Wambach, J. Spectral Functions for the Quark-Meson Model Phase Diagram from the Functional Renormalization Group. _Phys. Rev. D_ 2014, _89_ , 034010. doi:10.1103/PhysRevD.89.034010.
  * (8) Gao, F.; Pawlowski, J.M. QCD phase structure from functional methods. _Phys. Rev. D_ 2020, _102_ , 034027. doi:10.1103/PhysRevD.10 2.034027.
  * (9) Gerber, P.; Leutwyler, H. Hadrons Below the Chiral Phase Transition. _Nucl. Phys. B_ 1989, _321_ , 387–429. doi:10.1016/0550-3213(89)90349-0.
  * (10) Alkofer, R.; Maas, A.; Mian, W.A.; Mitter, M.; Paris-Lopez, J.; Pawlowski, J.M.; Wink, N. Bound state properties from the functional renormalization group. _Phys. Rev. D_ 2019, _99_ , 054029. doi:10.1103/PhysRevD.99.054029.
  * (11) Herbst, T.K.; Mitter, M.; Pawlowski, J.M.; Schaefer, B.J.; Stiele, R. Thermodynamics of QCD at vanishing density. _Phys. Lett. B_ 2014, _731_ , 248–256. doi:10.1016/j.physletb.2014.02.045.
  * (12) Hüfner, J.; Klevansky, S.P.; Zhuang, P.; Voss, H. Thermodynamics of a quark plasma beyond the mean field: A generalized Beth-Uhlenbeck approach. _Ann. Phys_. 1994, _234_ , 225–244. doi:10.1006/aphy.1994.1080.
  * (13) Zhuang, P.; Hüfner, J.; Klevansky, S.P. Thermodynamics of a quark - meson plasma in the Nambu-Jona-Lasinio model. _Nucl. Phys. A_ 1994, _576_ , 525–552. doi:10.1016/0375-9474(94)90743-9.
  * (14) Wergieluk, A.; Blaschke, D.; Kalinovsky, Y.L.; Friesen, A. Pion dissociation and Levinson’s theorem in hot PNJL quark matter. _Phys. Part. Nucl. Lett._ 2013, _10_ , 660–668. doi:10.1134/S1547477113070169.
  * (15) Yamazaki, K.; Matsui, T. Quark-Hadron Phase Transition in the PNJL model for interacting quarks. _Nucl. Phys. A_ 2013, _913_ , 19–50. doi:10.1016/j.nuclphysa.2013.05.018.
  * (16) Blaschke, D.; Buballa, M.; Dubinin, A.; Röpke, G.; Zablocki, D. Generalized Beth–Uhlenbeck approach to mesons and diquarks in hot, dense quark matter. _Ann. Phys._ 2014, _348_ , 228. doi:10.1016/j.aop.2014.06.002.
  * (17) Yamazaki, K.; Matsui, T. Quark-hadron phase transition in a three flavor PNJL model for interacting quarks. _Nucl. Phys. A_ 2014, _922_ , 237–261. doi:10.1016/j.nuclphysa.2013.12.010.
  * (18) Torres-Rincon, J.M.; Aichelin, J. Equation of state of a quark-Polyakov loop-meson mixture in the PNJL model at finite temperature. _arXiv_ 2016, arXiv:1601.01706.
  * (19) Torres-Rincon, J.M.; Aichelin, J. Equation of state of a quark-meson mixture in the improved Polyakov–Nambu–Jona-Lasinio model at finite chemical potential. _Phys. Rev. C_ 2017, _96_ , 045205. doi:10.1103/PhysRevC.96.045205.
  * (20) Blaschke, D.; Dubinin, A.; Radzhabov, A.; Wergieluk, A. Mott dissociation of pions and kaons in hot, dense quark matter. _Phys. Rev. D_ 2017, _96_ , 094008. doi:10.1103/PhysRevD.96.094008.
  * (21) Blaschke, D.; Dubinin, A.; Turko, L. Mott-hadron resonance gas and lattice QCD thermodynamics. _arXiv_ 2016, arXiv:1611.09845.
  * (22) Blaschke, D.; Dubinin, A.; Turko, L. Generalized Beth-Uhlenbeck approach to the equation of state for quark-hadron matter. _Acta Phys. Polon. Supp._ 2017, _10_ , 473–480. doi:10.5506/APhysPolBSupp.10.473.
  * (23) Dubinin, A.; Blaschke, D.; Kalinovsky, Y.L. Pion and sigma meson dissociation in a modified NJL model at finite temperature. _Acta Phys. Polon. Supp._ 2014, _7_ , 215–223. doi:10.5506/APhysPolBSupp.7.215.
  * (24) Cheng, M.; Christ, N.H.; Datta, S.; Heide, J.V.; Jung, C.; Karsch, F.; Kaczmarek, O.; Laermann, E.; Mawhinney, R.D.; Miao, C.; et al. The QCD equation of state with almost physical quark masses. _Phys. Rev. D_ 2008, _77_ , 014511. doi:10.1103/PhysRevD.77.014511.
  * (25) Bazavov, A.; Ding, H.T.; Hegde, P.; Kaczmarek, O.; Karsch, F.; Karthik, N.; Laermann, E.; Lahiri, A.; Larsen, R.; HotQCD; et al. Chiral crossover in QCD at zero and non-zero chemical potentials. _Phys. Lett. B_ 2019, _795_ , 15–21. doi:10.1016/j.physletb.2019.05.013.
  * (26) Turko, L.; Blaschke, D.; Prorok, D.; Berdermann, J. Mott-Hagedorn Resonance Gas and Lattice QCD Results. _Acta Phys. Polon. Supp._ 2012, 5, 485. doi:10.5506/APhysPolBSupp.5.485.
  * (27) Ratti, C.; Thaler, M.A.; Weise, W. Phases of QCD: Lattice thermodynamics and a field theoretical model. _Phys. Rev. D_ 2006, 73, 014019. doi:10.1103/PhysRevD.73.014019.
  * (28) Kapusta, J.I. _Finite Temperature Field Theory_ ; Cambridge University Press: Cambridge, UK, 1989.
  * (29) Peshier, A.; Kämpfer, B.; Soff, G. The QCD equation of state and quark star properties. _arXiv_ 2003, arXiv:hep-ph/0312080.
  * (30) Shirkov, D.V. On the Fourier transformation of renormalization invariant coupling. _Theor. Math. Phys._ 2003, _136_ , 893. doi:10.1023/A:1024560420051.
  * (31) Blaschke, D.; Kaczmarek, O.; Laermann, E.; Yudichev, V. Heavy quark potential and quarkonia dissociation rates. _Eur. Phys. J. C_ 2005, _43_ , 81. doi:10.1140/epjc/s2005-02327-5.
  * (32) Bazavov, A.; Brambilla, N.; Ding, H.T.; Petreczky, P.; Schadler, H.P.; Vairo, A.; Weber, J.H. Polyakov loop in 2+1 flavor QCD from low to high temperatures. _Phys. Rev. D_ 2016, _93_ , 114502. doi:10.1103/PhysRevD.93.114502.
  * (33) Bazavov, A.; Petreczky, P.; Weber, J.H. Equation of State in 2+1 Flavor QCD at High Temperatures. _Phys. Rev. D_ 2018, _97_ , 014510. doi:10.1103/PhysRevD.97.014510.
  * (34) Allton, C.R.; Döring, M.; Ejiri, S.; Hands, S.J.; Kaczmarek, O.; Karsch, F.; Laermann, E.; Redlich, K. Thermodynamics of two flavor QCD to sixth order in quark chemical potential. _Phys. Rev. D_ 2005, _71_ , 054508. doi:10.1103/PhysRevD.71.054508.
  * (35) Biro, T.S.; Jakovac, A. QCD above $T_{c}$: Hadrons, partons, and the continuum. _Phys. Rev. D_ 2014, _90_ , 094029. doi:10.1103/PhysRevD. 90.094029.
  * (36) Turko, L.; Blaschke, D.; Prorok, D.; Berdermann, J. An effective model of QCD thermodynamics. _J. Phys. Conf. Ser._ 2013, _455_ , 012056. doi:10.1088/1742-6596/455/1/012056.
  * (37) Turko, L.; Blaschke, D.; Prorok, D.; Berdermann, J. Effective degrees of freedom in QCD thermodynamics. _EPJ Web Conf._ 2014, _71_ , 00134. doi:10.1051/epjconf/20147100134.
  * (38) Bastian, N.U.F.; Blaschke, D.; Fischer, T.; Röpke, G. Towards a Unified Quark-Hadron Matter Equation of State for Applications in Astrophysics and Heavy-Ion Collisions. _Universe_ 2018, _4_ , 67. doi:10.3390/universe4060067.
  * (39) Blaschke, D.; Dubinin, A.; Buballa, M. Polyakov-loop suppression of colored states in a quark-meson-diquark plasma. _Phys. Rev. D_ 2015, _91_ , 125040. doi:10.1103/PhysRevD.91.125040.
  * (40) Blaschke, D.; Dubinin, A.S.; Zablocki, D. NJL model approach to diquarks and baryons in quark matter. _PoS BaldinISHEPPXXII_ 2015, 083. doi:10.22323/1.225.0083.
  * (41) Bastian, N.U.F.; Blaschke, D.B. A unified quark-nuclear matter equation of state from the cluster virial expansion within the generalized Beth–Uhlenbeck approach. _Eur. Phys. J. A_ 2021, _57_ , 35. doi:10.1140/epja/s10050-021-00343-9.
  * (42) Plumberg, C.J.; Welle, T.; Kapusta, J.I. QCD matter with a crossover and a first-order phase transition. _PoS CORFU2018_ 2018, 157. doi:10.22323/1.347.0157.
  * (43) Parotto, P.; Bluhm, M.; Mroczek, D.; Nahrgang, M.; Noronha-Hostler, J.; Rajagopal, K.; Ratti, C.; Schäfer, T.; Stephanov, M. QCD equation of state matched to lattice data and exhibiting a critical point singularity. _Phys. Rev. C_ 2020, _101_ , 034901. doi:10.1103/PhysRevC.101.034901.