# Massive Gauge Theory with Quasigluon for Hot $SU(N)$ Part I: Phase
Transition and Thermodynamics

Jun Guo<EMAIL_ADDRESS>College of Physics and Communication Electronics,
Jiangxi Normal University, Nanchang 330022, China  Zhaofeng Kang
<EMAIL_ADDRESS>School of physics, Huazhong University of Science and
Technology, Wuhan 430074, China  Jiang Zhu<EMAIL_ADDRESS>School
of physics, Huazhong University of Science and Technology, Wuhan 430074, China

###### Abstract

It is challenging to build a model that can correctly and unifiedly account
for the deconfinement phase transition and thermodynamics of the hot $SU(N)$
pure Yang-Mills (PYM) system, for any $N$. In this article, we slightly
generalize the massive PYM model to the situation with a quasigluon mass
$M_{g}(T)$ varying with temperature, inspired by the quasigluon model. In such
a framework, we can acquire an effective potential for the temporal gauge
field background by perturbative calculation, rather than adding by hand. The
resulting potential works well to describe the behavior of the hot PYM system
for all $N$, via the single parameter $M_{g}(T)$. Moreover, under the
assumption of unified eigenvalue distribution, the $M_{g}(T)$ fitted by
machine learning is found to follow $N$-universality.

###### pacs:

12.60.Jv, 14.70.Pw, 95.35.+d

##  I Introduction

To build a model that can describe the deconfinement phase transition of the
$SU(N)$ PYM system at finite temperature, which is hampered by the non-
perturbative effect, one should first figure out what knowledge we have about
such a system. In the very high-temperature region, it should recover the
Stefan-Boltzmann (SB) limit, following the asymptotic freedom of the non-
Abelian gauge theory. The more crucial information comes from the lattice
simulations, which provide a reliable way to deal with strong coupling, thus
furnishing the order of deconfinement phase transition and as well data of
therm-dynamic observables, such as pressure and the latent heat $L$ 111The
latent heat is not only important to describe the first order phase transition
(FOPT) but also critical in cosmology because the Gravitational-wave produced
during the deconfinement phase transition in the early universe is directly
related to this quantity Caprini:2015zlo ; Kubo:2018vdw ; Helmboldt:2019pan ;
Halverson:2020xpg ; Huang:2020crf ; Kang:2021epo .. The deconfinement phase
transition is a crossover for $N=2$ and FOPT for $N\geq 3$. Moreover,
combining the data for pressure and the latent heat allows us to extract the
following large $N$ scaling law Kang:2021epo

$\begin{split}p_{M}=\frac{M^{2}-1}{N^{2}-1}p_{N},\ \ \ \
L_{M}=\frac{M^{2}-1}{N^{2}-1}L_{N},\end{split}$ (1)

where $N$ and $M$ represent different color number. It is challenging to build
a model with strong theoretic ground that can correctly account for all of the
above aspects of the hot PYM system, for any $N$ beyond $N=3$. However, it is
very meaningful, not only in the theoretical sense but also in the application
to the new physics domain, where an (almost) pure $SU(N)$ gauge sector
receives wide interest baryon ; twin ; string ; Kang:2019izi ; Carenza:2022pjd
. Recently, the prospects of gravitational wave signals during the
deconfinement phase transition are studied based on different models
Kubo:2018vdw ; Halverson:2020xpg ; Huang:2020crf ; Kang:2021epo ;
Morgante:2022zvc .

The popular line is following the $Z_{N}$ center symmetry and the traced
Polyakov loop (PL) as order parameter, to construct effective PL models,
usually, the polynomial models Pisarski:2001pe ; Ratti:2006gh ; Huang:2020crf
; Kang:2021epo  also see review Fukushima:2017csk . Another line is underlined
by the Haar measure, which gains great success in the $SU(3)$ case, even
incorporating dynamic quarks Mocsy:2003qw ; Fukushima:2003fw ; Roessner:2006xn
; Fukushima:2008wg . For $SU(3)$ only, both types of model can describe the
deconfinement phase transition and as well the thermodynamics, at least in the
semi-QGP region. However, when we try to extend them to general $SU(N)$ cases,
we encounter some difficulties. The Haar-type model is shown to be
inconsistent with the above large $N$ scaling law Kang:2021epo  and moreover,
it can not be handled for very large $N$. The polynomial model proposed in Ref
Huang:2020crf ; Kang:2021epo  utilizes the competition among terms with
designed powers and signs to realize the deconfinement phase transition. Since
it basically respects just a $Z_{2}$ symmetry and thus works for any $N$, even
including $N=2$. The matrix models, inspired by the property of perturbation
potential, instead Meisinger:2001cq ; Meisinger:2003id ; Dumitru:2010mj ;
Dumitru:2012fw  treat the eigenvalues of the thermal Wilson line as
fundamental variables, which may provide a feasible way to understand the
behavior of the hot PYM system for all $N$. Largely speaking, these models are
phenomenally oriented, lacking a more profound basis to derive the shape of
the thermodynamic potential.

As long as only thermodynamics is concerned, the quasi-particle model (QPM) is
even more attractive. It is a statistical model where the gluons are assumed
to develop a temperature-dependent mass, due to the non-perturbative
interaction with the thermal environment. This picture is strongly supported
by the hard-thermal-loop perturbation theory at high-temperature regions HTLpt
. It can successfully explain the thermodynamics of the hot $SU(N)$ PYM system
from $T_{c}$ to the SB limit Goloviznin:1992ws ; Peshier:1994zf ;
Peshier:1995ty ; Gorenstein:1995vm . Later, taking into account the temporal
gauge field background $A_{0}$ brings a difference Meisinger:2003id  and opens
the possibility to describe both thermodynamics and the deconfinement phase
transition at the same time. But most studies of the interplay between
quasigluon and background focus on the modification to the pressure of hot
PYM, the critical non-perturbative dynamics driving the deconfinement phase
transition, says the Haar measure term, is still added by hand and external to
the QPM picture Sasaki:2012bi ; Ruggieri:2012ny ; Islam:2021qwh . This might
be contradictory to the spirit of QPM, where most of the non perturbation
interaction has already been “absorbed” into the quasiguon mass.

Although not a following study of QPM, the massive PYM model Curci:1976bt ;
Reinosa:2014ooa ; vanEgmond:2021jyx  shares a similar philosophy with QPM, and
it also assumes that the effective gluon mass parameter simply encodes the
non-perturbative effects. Then, the effective potential for the temporal
background can be derived, instead of added by hand, at one loop level or even
beyond Reinosa:2014zta ; Reinosa:2015gxn . This approach realizes the inverted
Weiss potential, attributed to the enhanced ghost contribution, as the
mechanism for the deconfinement phase transition. Surprisingly, for $N=2$, the
resulting effective potential indeed predicts a crossover instead of FOPT.

Thus, it is tempting to marry QPM with the massive PYM. The original massive
PYM model Reinosa:2014ooa  just takes a constant quasigluon mass, and now we
generalize it to have temperature dependence, which is in line with the
framework of hard-thermal-loop perturbation theory and may serve as a quantum
field basis for the QPM. We find that the resulting one-loop effective
potential indeed can successfully describe both the deconfinement phase
transition and the thermodynamics of the hot PYM system, for any color number
$N$. Our study is helpful to understand the deconfinement phase transition in
cosmology.

The paper is organized as follows: We give a short review of the QPM in
Section II and then goes to the generalized massive PYM according to the QPM
in Section III, where we derive the effective potential at one loop and
investigate the deconfinement phase transition with the assumption of uniform
eigenvalue distribution, which reduces the potential to one-dimension. In
Section IV we study thermodynamics from the critical temperature to the SB
limit, fitting the quasigluon mass by lattice data via machine learning.
Conclusions and discussions and as well as the appendix are cast in the
remaining two sections.

##  II Quasigluon: From HTL to $T_{c}$

For thermal gauge theories at high temperature, the classical solution should
not be described by the gluonic states without mass but with mass, which stems
from the plasma effects such as the screening of electric fields and Landau
damping. The hard-thermal-loop perturbation theory (HTLpt) HTLpt , which is a
reorganization of the perturbation series and can take into account the plasma
effects consistently. It is found that at NNLO, the hot gluon plasma can be
well described by weakly coupled quasigluons down to $~{}(2-3)T_{c}$
Andersen:2009tc ; Andersen:2010wu .

Within the HTLpt, the transverse quasigluon in the QCD medium follows the
dispersion equation

$w^{2}-k^{2}-\Pi_{t}^{*}(w,k)=0,$ (2)

where $\Pi_{t}^{*}(w,k)$ is the transverse self-energy for the hot gluons,
having weak momentum dependence but strong temperature dependence. At leading
order, it is given by HTLpt

$\Pi_{t}^{*}(w,k)=\frac{N}{6}g^{2}T^{2},$ (3)

with $g$ the gauge coupling. The gluon quasiparticles mainly propagate on
shell.

For even lower temperature, the magnetic/nonperturbative effects become
important. But it is tempting to pursue the possibility that even down to
$T_{c}$, the plasma can still be described by an ideal gas of “massive”
noninteracting “gluons”, where the strong interactions between gluon and the
in-medium have been “absorbed”, at least partially, into the quasigluon mass.
Following this line, the authors of Ref. Goloviznin:1992ws ; Peshier:1995ty
explained the lattice QCD thermodynamics near $T_{c}$ via a simple
quasiparticle model (QPM) inspired by the above HTL quasiparticle. And to
naturally match the HTLpt quasiparticle at high $T$, they simply consider such
a QPM with quasigluon mass squared 222The quasigluon mass is determined by the
pole of gluon self-energy in the complex momentum plane, but the exact
location is hampered by the non-perturbative effect.

$M_{g}^{2}(T)=\frac{N}{6}G^{2}(T)T^{2},\quad
G^{2}(T)=\frac{48\pi^{2}}{11N\log\left(\frac{T}{T_{c}/\lambda}+\frac{T_{s}}{T_{c}}\right)^{2}}$
(4)

which parameterizes the deviation from HTLpt quasiparticle via the parameter
$T_{s}$ and $\lambda$. Other form of $M_{g}(T)$ is possible, for instance the
one in Ref Schneider:2001nf .

Then, such a pool of ideal quasi gluon gas, assumed to respect the Bose
distribution $f_{B}$, has pressure

$p(T)=\frac{g(T)}{6\pi^{2}}\int_{0}^{\infty}f_{B}(E_{k})\frac{k^{4}}{E_{k}}dk-B(T),$
(5)

with $E_{k}=\sqrt{k^{2}+M_{g}^{2}(T)}$. Owing to the temperature dependence of
mass, the self-consistent thermodynamic relation for ideal gas, namely the
Gibbs-Duhem relation, $\epsilon+p=sT$ with $s=\partial p/\partial T$ is
violated. Including $B(T)$ can solve this problem Gorenstein:1995vm . It is
not independent and is determined by $M_{g}(T)$, up to a bag constant.
Surprisingly, this simple QPM is capable of reproducing the quenched QCD or
$SU(3)$ PYM lattice data in the whole region above $T_{c}$ Peshier:1995ty .
Study for Other $N=4,5,6$ is presented in Ref. Castorina:2011ra ;
Castorina:2011ja . It is common that, in order to reduce the contribution of
quasi gluons near the critical temperature, a very large quasi gluon mass is
usually required.

However, the original QPM is just a statistical model and thus can not explain
the order of $SU(N)$ phase transition. The latter is supposed to be understood
in the framework of Landau phase transition: Find a proper order parameter
$\eta$ and construct a (coarse-grained) Landau free energy as a function of
the order parameter, and then one can study the order of phase transition by
surveying its ground state. In studying the deconfinement phase transition of
$SU(N)$ PYM, the Polyakov loop (PL) associated with the center symmetry
$Z_{N}$ is identified with $\eta$; it is defined as $l_{N}={\rm
tr}\hat{L}_{F}/N$, the traced thermal Wilson line in the fundamental
representation

$\hat{L}_{F}={\cal P}e^{ig\int_{0}^{\beta}A^{a}_{4}(x,t)t^{a}dt},$ (6)

with ${\cal P}$ denoting path ordering and $t^{a}$ the generators of the
fundamental representation for $SU(N)$.

##  III Quasiparticles move in the PL background

So, it is a natural idea to combine quasi-particle model with PLM, to study
the deconfinement phase transition dynamics and thermodynamics simultaneously
Meisinger:2003id ; Ruggieri:2012ny  333This idea originated from an earlier
work Meisinger:2001cq , although there the authors have not introduced
quasigluon explicitly yet. . In such models, quasigluons moving in the PL
background generate thermodynamic potential which depends on the PL in the
adjoint representation $\hat{L}_{A}$ Meisinger:2003id :

$\Omega_{\rm QG}(\hat{L}_{A},T)=2T{\rm
tr}\int\frac{d^{3}\vec{p}}{(2\pi)^{2}}\log(1-\hat{L}_{A}e^{-E_{g}/T}).$ (7)

It is a phenomenological generalization to the usual Weiss potential
Weiss:1980rj  for the fundamental gluons to quasigluons, by replacing
$|\vec{p}|$ with $E_{g}=\sqrt{M_{g}(T)^{2}+p^{2}}$. Later, we will derive a
similarity grounded on the QFT, but with a remarkable difference. The
quasigluons dominate thermodynamics in the high-temperature region, where
$M_{g}\ll T$ and $\hat{L}_{A}\rightarrow 1$, explain the blackbody behavior.
At the lower temperature, typically below $2T_{c}$, the decreasing PL combined
with the increasing $M_{g}(T)$, is capable of explaining the deviation from
the blackbody spectrum towards $T_{c}$ Meisinger:2003id ; Sasaki:2012bi ;
Ruggieri:2012ny ; Alba:2014lda ; Islam:2021qwh .

But that’s all. We can’t expect this part to give the deconfinement phase
transition at the same time, which needs additional interaction, such as the
van der Monde determinant interaction Meisinger:2003id ; Sasaki:2012bi ;
Lo:2021qkw ; Islam:2021qwh . In this article, we follow another line proposed
in Ref. Reinosa:2014ooa , which enables us to study the non-perturbative PT in
the perturbative approach; in their philosophy, non-perturbation effects are
encoded in the gluon mass, in line with the QPM picture. In the following, we
will first present an effective model, which is a slight generalization to
that in Ref. Reinosa:2014ooa . Then, we reproduce the effective potential Eq.
(7) as well as the confining potential from the model through the leading
order thermal correction.

###  III.1 Effective model for quasigluon above $T_{c}$

The model in Ref. Reinosa:2014ooa  quantizes PYM in the background field gauge
formalism, including massive fluctuations. Then, the Faddeev-Popov gauge-fixed
Lagrangian reads

$\mathcal{L}=-\frac{1}{2g^{2}}{\rm
tr}(F_{\mu\nu}F^{\mu\nu})+\bar{D}_{\mu}\bar{c}^{a}D^{\mu}c^{a}+ih^{a}\bar{D}_{\mu}\hat{A}^{\mu,a}+\frac{1}{2}M^{2}_{g}(T)\hat{A}^{a}_{\mu}\hat{A}^{a,\mu},$
(8)

where $c,\bar{c}$ and $h$ are the Ghost fields, real Nakanishi-Lautrup field,
respectively. We have split the gauge field $A_{\mu}$ as
$A_{\mu}=\bar{A}_{\mu}+\hat{A}_{\mu}$ with $\hat{A}_{\mu}$ the massive
fluctuations. The background $\bar{A}_{\mu}$ is restricted to merely have the
constant temporal component, $\bar{A}_{\mu}=\bar{A}_{0}\delta_{0\mu}$, for the
sake of preserving invariance of the PYM system, under both temporal and
spatial translations and spatial rotations at finite $T$. The covariant
derivative acting on $\phi=(c,\bar{c},h,\hat{A}_{\mu})$ is defined as

$\bar{D}_{\mu}^{ab}=\partial_{\mu}\delta^{ab}+gf^{acb}\bar{A}^{c}_{\mu},$ (9)

where the gauge field is the background field. The above Lagrangian has
implemented the Landau-DeWitt gauge $\bar{D}_{\mu}\hat{A}^{\mu}=0$. This gauge
fixed PYM, including the gluon mass term, still respects the background local
$SU(N)$ symmetry, with covariant derivative defined above and treating $\phi$
as adjoint matter fields.

In the effective model specified by Eq. (8), the gluon mass is not originally
interpreted as quasigluon mass. Instead, it is regarded as a gauge fixing
parameter, to further remove the degeneracy among the Gribov copies, whose
existence may make the Faddeev-Popov procedure in the deep infrared region
invalid Gribov:1977wm . This region is associated with the nonperturbative
dynamics of PYM. Hence, people hope that $M_{g}$ at the same time can “absorb”
strong interactions, so that some non-perturbation phenomena can be studied by
the perturbation method. Such a philosophy is consistent with the QPM, and
therefore it is tempting to simply identify $M_{g}$ as the quasigluon mass,
which is reasonable at least at zero temperature. If such a formalism is
consistent with the Hamilton approach which establishes a QFT basis for
quasiparticle Heffner:2012sx , is open.

However, to explain thermodynamics, we need a temperature-dependent quasigluon
mass, $M_{g}(T)$. This may be odd with the usual understanding of thermal mass
origin in perturbative thermal QFT: The underlying Lagrangian is the same as
that of $T=0$ and does not include the temperature-dependent quantity, and
this kind of dependence originates from thermal correction. However, it is not
strange that the Lagrangian includes a temperature dependent quantity. In
fact, the HTL resummation scheme based on quasi-particle picture is just based
on the effective Lagrangian including thermal mass, which gives rise to the
modified propagator for the calculation of thermal corrections. Since we are
extending the quasi particle picture down to near $T_{c}$, we should naturally
include the temperature-dependent quasigluon mass term.

Therefore, as a slight generalization to the model in Ref. Reinosa:2014ooa ,
the effective Lagrangian Eq. (8) is supposed to furnish a phenomenological
framework to perturbatively study deconfinement phase transition along with
full thermodynamics above $T_{c}$.

###  III.2 The Thermodynamic Potential for Quasi-Particle Model: pure gluonic
part

In this subsection, we will calculate the thermodynamic potential for the
fundamental PL in a general PYM with gauge group $SU(N)$, following the
textbook approach. That is to integrate out all fluctuations
$\hat{A}_{\mu}=A_{\mu}-\bar{A}_{\mu}$ over the temporal background
$\bar{A}_{\mu}=\bar{A}_{0}\delta_{0\mu}$, in the $3+1$ Eucleadian QFT. For a
homogeneous background, one can always make $\bar{A}_{0}$ diagonal via some
global $SU(N)$ rotation. Therefore, we can expand $\bar{A}_{0}$ in the $su(N)$
Cartan space, which is spanned by the diagonal subgroup $\\{H^{i}\\}$
$(i=1,2,...N-1)$ with $[H^{i},H^{j}]=0$, and then
$\bar{A}_{0}=\bar{A}_{0}^{i}H^{i}$ with $\bar{A}_{0}^{i}$ is the Cartan
coordinates.

Let us first deal with the pure gluonic part of Eq. (8), from which one can
get the quadratic Lagrangian of the fluctuation field $\hat{A}_{\mu}$

$\mathcal{L}^{(2)}=-\frac{1}{2}\hat{A}_{\mu}^{a}[\delta_{ab}g^{\mu\nu}\partial^{2}-f_{abc}(\partial^{\nu}\bar{A}^{\mu,c}+2g^{\mu\nu}\bar{A}^{c}_{\rho}\partial^{\rho})+f_{acd}f_{cbe}g^{\mu\nu}\bar{A}^{d}_{\rho}\bar{A}^{\rho,e}+2f_{abc}\bar{F}^{\mu\nu,c}]\hat{A}^{b}_{\nu}.$
(10)

It can be written as the following

$\mathcal{L}^{(2)}=\frac{1}{2}\hat{A}_{\mu}^{a}(D^{-1})_{ab}\hat{A}^{\mu,b},$
(11)

with the operator defined as

$(D^{-1})_{ab}=\delta_{ab}(p^{2}+M_{g}^{2})+2i\sum_{i}f_{abi}\bar{A}^{i}_{0}p_{0}-\sum_{i,j}f_{aci}f_{cbj}\bar{A}^{i}_{0}\bar{A}^{j}_{0}.$
(12)

The temporal background contributes masses to the fluctuations, and hence
obtains a thermodynamic potential from the plasma of quasigluons calculable as
usual.

Before we calculate this potential, let’s deal with the propagator,
diagonalizing the fluctuations in the color space through a unitary
transformation. Then, the diagonal propagators takes the form of

$\tilde{D}^{-1}_{aa}(p)=(p_{0}-A_{a})^{2}+|\vec{p}|^{2}+M_{g}^{2}.$ (13)

Following the standard approach of path integral, one can get the generating
function ($Z$ below should be understood as $Z^{I}$, the gluonic part
contribution, but for the sake of simplicity, we ignore the superscript, which
we believe will not cause ambiguity)

$\log Z=\frac{1}{2}\log\det\left[-\frac{\delta^{2}\mathcal{L}^{2}}{\delta
A\delta
A}\right]=\frac{1}{2}\log\det\left[D^{-1}\right]=\frac{1}{2}\log\det\left[\tilde{D}^{-1}\right].$
(14)

where we have used the property that unitary transformation does not change
determinate. Then, using the trick that $\log\det A={\rm Tr}\log A$ we get

$\log\det[\tilde{D}^{-1}]={\rm Tr}\log[\tilde{D}^{-1}].$ (15)

“Tr” is the trace over the functional propagator operator, and can be split
into two parts: a function trace over momentum space and a color space trace
denoted by “${\rm tr_{c}}$”, explicitly,

${\rm Tr}\log[\tilde{D}^{-1}]={\rm
tr_{c}}\int\frac{d^{4}p}{(2\pi)^{4}}\log[\tilde{D}^{-1}(p)].$ (16)

In order to get the finite temperature potential, one can discretize the
energy by $p_{0}\rightarrow\omega_{n}=2i\pi nT$ and transform
$\bar{A}_{0}\rightarrow-i\bar{A}_{4}$, obtaning

$\log Z=2V{\rm
tr_{c}}\int\frac{d^{3}\vec{p}}{(2\pi)^{2}}\sum_{n=-\infty}^{\infty}\log[\tilde{D}^{-1}_{aa}(\omega_{n},|\vec{p}|)],$
(17)

where $2=\frac{1}{2}\times 4$ with 4 denoting the multiplicity from the four
components of $A_{\mu}$. From Eq. (13), the structure of the propagator
$\tilde{D}^{-1}(\omega_{n},|\vec{p}|)$ takes the form of

$\tilde{D}^{-1}_{aa}(\omega_{n},|\vec{p}|)=(\omega_{n}-A_{a})^{2}+|\vec{p}|^{2}+M_{g}^{2},$
(18)

where $A_{a}$ is a linear combination of the background $\bar{A}^{i}_{4}$,
with coefficients determined by the structure constant, but we do not find a
general expression for any $N$ yet. As a matter of fact, the concrete
expression is not important in our discussion, since later we will switch to a
parameterization of the background which is independent of $A_{a}$. Anyway, in
Appendix. B, we present the details of our calculation for $SU(4)$, and the
procedure applies to other values of $N$.

The summation of the thermal excitation modes $n$ can be done explicitly using
a trick in Appendix. C. And finally, the generating function can be compactly
written as

$\log Z=4V{\rm
tr_{c}}\int\frac{d^{3}\vec{p}}{(2\pi)^{2}}\log(1-\hat{L}_{A}e^{-E_{g}/T}),$
(19)

where $\hat{L}_{A}$ is expressed in terms of background field $\bar{A}_{\mu}$,
and it is nothing but the PL in the adjoint representation. For instance, in
$SU(3)$ it is given by

$\begin{split}\hat{L}_{A}=\operatorname{diag}[&1,1,e^{i\bar{A}^{3}_{4}/T},e^{-i\bar{A}^{3}_{4}/T},e^{i(\bar{A}^{3}_{4}+\sqrt{3}\bar{A}^{8}_{4})/2T},\\\
&e^{-i(\bar{A}^{3}_{4}+\sqrt{3}\bar{A}^{8}_{4})/2T},e^{i(\bar{A}^{3}_{4}-\sqrt{3}\bar{A}^{8}_{4})/2T},e^{-i(\bar{A}^{3}_{4}-\sqrt{3}\bar{A}^{8}_{4})/2T}],\end{split}$
(20)

where we have written it in terms of the original background. It is seen that
the temporal background behaves as an imaginary chemical potential. Eq. (19)
yields the effective potential ${\cal V}_{eff}=\frac{T}{V}\log Z$ which almost
recovers the generalized Weiss potential given in Eq. (7), up to the
coefficient. But the ghost contribution, which will be included in the
following subsection, will result in a substantial deviation related to the
deconfinement potential.

Eq. (20) demonstrates the general structure of thermal Wilson line in the
adjoint representation, i.e., its elements are organized such that it can be
rewritten in terms of the eigenphases of the fundamental thermal Wilson line
Lo:2021qkw

$\hat{L}_{F}={\cal
P}e^{ig\int_{0}^{\beta}A_{4}(x,t)dt}\rightarrow\operatorname{diag}[e^{i2\pi
q_{1}},e^{i2\pi q_{2}},...,e^{i2\pi q_{N}}],$ (21)

by virtue of the parameterization of background $\bar{A}_{4}=2\pi/(g\beta){\rm
diag}(q_{1},q_{2},...,q_{N})$, with the real $q_{i}$ satisfying the constraint
$\sum_{i=1}^{N}q_{i}=0$. As a phase factor, it is sufficient to work in the
interval $0\leq q_{i}\leq 1$. And now,

$\hat{L}_{A}=\operatorname{diag}[1,1,...,1,e^{i2\pi q_{ij}},...,e^{-i2\pi
q_{ij}}],$ (22)

where the $N-1$ “1” corresponds to the Cartan part, while the $N(N-1)/2$ pairs
of $q_{ij}\equiv q_{i}-q_{j}$ with $N\geq i>j\geq 1$ corresponds to the non-
Cartan part. The above form is more convenient and will be adopted hereafter.
Then
$1-\hat{L}_{A}e^{-E_{g}/T}=\operatorname{diag}(1-e^{-E_{g}/T},...,1-e^{i2\pi
q_{ij}-E_{g}/T},...,1-e^{-i2\pi q_{ij}-E_{g}/T})$. We also define

$\displaystyle\Omega(q_{ij},M_{g})$ $\displaystyle\equiv
2T\int\frac{d^{3}\vec{p}}{(2\pi)^{2}}\log(\det(1-\hat{L}_{A}e^{-E_{g}/T}))$
(23)
$\displaystyle=2T\sum_{i,j=1}^{N}\left(1-\frac{\delta_{ij}}{N}\right)\int\frac{d^{3}\vec{p}}{(2\pi)^{3}}\log(1-e^{-\beta
E_{g}}e^{2\pi iq_{ij}}),$

where in the second line we use $q_{ij}=-q_{ji}$ and allow $i=j$, to write the
summation compactly 444Actually, if we instead adopt the ladder basis for
$su(N)$ and start from the eigenphase parameterization of $\bar{A}_{4}$
Fukushima:2017csk , the above expression can be explicitly obtained from the
covariant derivative whose background dependent term reads
$[\bar{A}_{4},\hat{A}_{\mu}]\sim(q_{i}-q_{j})\delta_{\mu 4}$.. In this
notation, the gluonic part contribution to the effective potential ${\cal
V}_{eff}$ is $2\Omega(q_{ij},M_{g})$.

In the following, we present two important expansions of this potential, the
low temperature and the high temperature expansion. Both will be used in the
later discussions.

####  III.2.1 Low temperature expansion

In the QPM, it is found that the fitted $M_{g}(T)/T$ is sufficiently large at
least around $T_{c}$, hence one has $E_{g}/T>M_{g}/T\gtrsim{\cal O}(1)$. We
will find this is also true in our model from a full numerical study, which
enables us to make a low temperature expansion for the effective potential
around $T_{c}$. This leads to an analytical expression, which is useful in the
phase transition analysis. First, we expand the logarithm in
$\Omega(q_{ij},M_{g})$, retaining $\hat{L}_{A}$,

$\Omega(\hat{L}_{A},M_{g})=-\frac{T}{\pi^{2}}\sum_{n=1}^{\infty}\frac{1}{n}{\rm
tr}(\hat{L}_{A})^{n}\int p^{2}e^{-nE_{g}/T}dp.$ (24)

Then we substitute $p=M_{g}\sinh t$ to get

$\displaystyle\Omega(\hat{L}_{A},M_{g})=$
$\displaystyle-\frac{TM_{g}^{3}}{\pi^{2}}\sum_{n=1}^{\infty}\frac{1}{n}{\rm
tr}(\hat{L}_{A})^{n}\int_{0}^{\infty}\cosh t\sinh^{2}te^{-n(M_{g}/T)\cosh
t}dt$ (25)

Now we use the following trick to rewrite the integral as

$\Omega(\hat{L}_{A},M_{g})=\frac{T^{2}M_{g}^{3}}{\pi^{2}}\sum_{n=1}^{\infty}\frac{1}{n^{2}}{\rm
tr}(\hat{L}_{A})^{n}\frac{d}{dM_{g}}\int_{0}^{\infty}\sinh^{2}te^{-n(M_{g}/T)\cosh
t}dt,$ (26)

where the integral can be done explicitly, $\sim K_{1}(x)$, and we finally
arrive

$\displaystyle\Omega(\hat{L}_{A},M_{g})=-\frac{T^{2}M_{g}^{2}}{\pi^{2}}\sum_{n=1}^{\infty}\frac{1}{n^{2}}{\rm
tr}(\hat{L}_{A})^{n}K_{2}(nM_{g}/T).$ (27)

with $K_{i}(x)$ the modified Bessel function of the second kind, of order $i$.
For $M_{g}/T$ moderately larger than 1, the leading order is a good
approximation.

####  III.2.2 High temperature expansion

At the high temperature limit, where $1\gg M_{g}/T$, one can find a simple
analytic expression of $\Omega(q_{ij},M_{g})$; see also Ref. Meisinger:2001fi
for the complete high temperature expansions beyond the leading term. To that
end, we again expand the logarithmic function:

$\Omega(q_{ij},M_{g})=-\frac{T^{4}}{\pi^{2}}\sum_{i,j=1}^{N}\left(1-\frac{\delta_{ij}}{N}\right)\int_{0}^{\infty}dxx^{2}\sum_{n=1}^{\infty}\frac{1}{n}e^{2\pi
niq_{ij}}e^{-\sqrt{n^{2}(x^{2}+\beta^{2}M_{g}^{2})}}.$ (28)

Expand this expression according to $\beta M_{g}$, and we get

$\begin{split}\Omega(q_{ij},M_{g})&=-\frac{T^{4}}{\pi^{2}}\sum_{i,j=1}^{N}\left(1-\frac{\delta_{ij}}{N}\right)\int_{0}^{\infty}dxx^{2}\sum_{n=1}^{\infty}\frac{1}{n}e^{2\pi
niq_{ij}}\left[e^{-nx}-\beta^{2}M_{g}^{2}\frac{ne^{-nx}}{2x}\right]+\mathcal{O}(\beta^{2}M_{g}^{2}).\end{split}$
(29)

The summation over $n$ in the first term is straightforward, while in the
second term, with one more $``n"$ factor from the Taylor expansion, can be
done as the following,

$\begin{split}\sum_{n=1}^{\infty}\frac{1}{n}ne^{-nx}=-\left(\sum_{n=1}^{\infty}\frac{1}{n}e^{-nx}\right)^{\prime}=-\frac{d}{dx}\log\left(1-e^{-x}\right).\end{split}$
(30)

Similar operations can be generalized to higher order of $n$, leading to
higher derivative to $\log\left(1-e^{-x}\right)$. Now, the first two terms are
summed to

$\displaystyle\Omega(q_{ij},M_{g})$
$\displaystyle=-\frac{T^{4}}{\pi^{2}}\sum_{i,j=1}^{N}\left(1-\frac{\delta_{ij}}{N}\right)\int_{0}^{\infty}dxx^{2}[\log\left(1-e^{-x+2\pi
iq_{ij}}\right)$ (31)
$\displaystyle+\frac{\beta^{2}M_{g}^{2}}{2x}\frac{d}{dx}\log\left(1-e^{-x+2\pi
iq_{ij}}\right)].$

The result of these integration are just two Polylogarithm function

$\begin{split}\Omega(q_{ij},M_{g})=-\frac{T^{4}}{\pi^{2}}\sum_{i,j=1}^{N}\left(1-\frac{\delta_{ij}}{N}\right)\left[2Li_{4}(e^{2\pi
iq_{ij}})-\frac{\beta^{2}M_{g}^{2}}{2}Li_{2}(e^{2\pi
iq_{ij}})\right]+\mathcal{O}(\beta^{2}M_{g}^{2}).\end{split}$ (32)

Since $q_{ij}=q_{i}-q_{j}=-q_{ji}$ we can rewrite this expression into an
analytic form by Jonquière’s inversion formula

$\begin{split}Li_{n}(e^{2\pi ix})+(-1)^{n}Li_{n}(e^{-2\pi ix})=-\frac{(2\pi
i)^{n}}{n!}B_{n}(x),\end{split}$ (33)

where $B_{n}$ is the Bernoulli polynomials. In our case, we can find that

$\begin{split}&Li_{2}(e^{2\pi iq_{ij}})+Li_{2}(e^{2\pi
iq_{ji}})=2\pi^{2}B_{2}(q_{ij}),\\\ &Li_{4}(e^{2\pi iq_{ij}})+Li_{4}(e^{2\pi
iq_{ji}})=-\frac{2\pi^{4}}{3}B_{4}(q_{ij}).\end{split}$ (34)

Substitute this formula into our expression to get

$\begin{split}\Omega(q_{ij},M_{g})=\frac{1}{\pi^{2}}\sum_{i\geq
j=1}^{N}\left(1-\frac{N-1}{N}\delta_{ij}\right)\left[\frac{4\pi^{4}}{3\beta^{4}}B_{4}(q_{ij})+\frac{M_{g}^{2}\pi^{2}}{\beta^{2}}B_{2}(q_{ij})\right]+\mathcal{O}(\beta^{2}M_{g}^{2}),\end{split}$
(35)

We can write it into a more simple form

$\begin{split}\Omega(q_{ij},M_{g})=\frac{1}{\pi^{2}}\sum_{i,j=1}^{N}\left(1-\frac{\delta_{ij}}{N}\right)\left[\frac{2\pi^{4}}{3\beta^{4}}B_{4}(|q_{ij}|)+\frac{M_{g}^{2}\pi^{2}}{2\beta^{2}}B_{2}(|q_{ij}|)\right]+\mathcal{O}(\beta^{2}M_{g}^{2}).\end{split}$
(36)

In the massless limit only $B_{4}$ is present.

###  III.3 The Thermodynamic Potential: gauge-fixed part & Phase transition

In this subsection, we include the contribution from the gauge-fixed part to
the thermodynamic potential, and then study how first order deconfinement
phase transition occurs due to the ghost contribution Reinosa:2014ooa .

####  III.3.1 Infrared ghost domination

The ghost contribution is similar to that from the gluons because it is also
in the adjoint representation. However, there are two key differences, which
enable the contribution of the thermodynamic potential from the ghost fields
to successfully trigger the deconfinement phase transition. First, the ghost
fields belong to Grassmann fields, and thus there is a minus sign relative to
the gluon contribution. Second, the ghosts are still massless since the
Lattice data does not show that the correlators of ghost develop a massive
pole. Moreover, we have to take into account the contribution from the gauge-
fixing term. To deal with this term, we should do the quadratic partition
between $\hat{A}^{a}$ and field $h^{a}$, to get two quadratic terms with the
mixing term eliminated; the details can be found in Appendix A or in the
textbook Reinosa:2020mnx . The final result of the total effective potential
is given by

$\begin{split}{\cal
V}_{eff}=\frac{3}{2}\Omega(M_{g})-\frac{1}{2}\Omega(M_{g}=0),\end{split}$ (37)

The above result is in the Landau-DeWitt gauge, and as usual, the effective
potential is gauge-dependent.

Without a quasi-gluon mass, the ghost contribution cancels the nonphysical
gluonic contribution and then the potential fails to admit a phase transition.
On the contrary, the presence of $M_{g}$ makes the enhanced ghost contribution
(relative to gluon contribution) dominate the potential at low temperature,
realizing the inverted Weiss potential as the confining mechanism.

Because we are dealing with the phase transition where $T\rightarrow T_{c}$
and $M_{g}/T\gg 1$, it is fair for us to expand the first term by low
temperature expansion Eq.(27). The second term is just the zero mass limit of
Eq.(36) and thus high temperature expansion applies. Combining the above
information, we can get the following analytic form of this effective
potential

$\begin{split}\mathcal{V}_{eff}&=3T\int\frac{d^{3}\vec{p}}{(2\pi)^{3}}{\rm
tr}\log(1-\hat{L}_{A}e^{-\frac{E_{g}}{T}})-T\int\frac{d^{3}\vec{p}}{(2\pi)^{3}}{\rm
tr}\log(1-\hat{L}_{A}e^{-\frac{|\vec{p}|}{T}})\\\
&\simeq-\frac{3T^{4}}{2\pi^{2}}(\frac{M_{g}}{T})^{2}K_{2}(M_{g}/T){\rm
tr}(\hat{L}_{A})-\frac{1}{2\pi^{2}}\sum_{i,j=1}^{N}(1-\frac{\delta_{ij}}{N})\Big{[}\frac{2\pi^{4}}{3\beta^{4}}B_{4}(q_{ij})\Big{]}.\end{split}$
(38)

The first term can be translated to a function of PL, by using the identity
${\rm tr}\hat{L}_{A}={\rm tr}\hat{L}_{F}{\rm tr}\hat{L}_{F}^{\dagger}-1$.
Nevertheless, the ghost term contains $N-1$ independent variables $q_{i}$,
rather than merely the trace part of $\hat{L}_{F}$. Hence, usually, one has to
deal with a multi-dimensional field space, case by case.

A way to reduce the potential to the one-dimension problem is assuming the
uniform eigenvalue distribution, i.e., $q_{ij}=\frac{i-j}{N}s$. It is
automatically true for $N=2,3,4$ with the number of independent eigenphases
less than 4, but it is merely an ansatz for the even higher $N$. Such an
ansatz has been adopted in Ref. Dumitru:2012fw , and is shown to work well.
The ansatz is based on the observation that the confining vacuum, which is
center symmetric, is characterized by the uniform eigenvalue distribution;
dynamically, the distribution is a result of the eigenvalue repulse from the
confining potential which involves the difference between eigenvalues $q_{i}$
Ivanov:2004gq ; Pisarski:2006hz ; Dumitru:2012fw . Furthermore, it is
conjectured that the transition from the deconfining vacuum to the confining
vacuum takes the shortest path, a straight line connecting the origin and the
confining vacuum Dumitru:2012fw . 555This is not the very precise statement.
At high $T$, the origin namely $A_{0}=0$ is the deconfining vacuum. As $T$
decreases to near $T_{c}$, it moves away from the origin to the configuration
still characterized by the uniform eigenvalue distribution. In the effective
matrix model, this is shown to be a good approximation even for the relatively
small $N$ Dumitru:2012fw . Then, we get the analytic potential for any color
number

$\begin{split}\mathcal{V}_{eff}&\simeq-\frac{3T^{4}N^{2}}{2\pi^{2}}(\frac{M_{g}}{T})^{2}K_{2}(M_{g}/T)l_{N}^{2}-\frac{2\pi^{2}T^{4}}{3}\sum_{i=1}^{N}(N-i)B_{4}(\frac{i}{N}s)+f(N,T),\end{split}$
(39)

where $f(N,T)$ is a function does not depend on the order parameter, only
relevant to thermadynamics. We can easily carry out this summation and find
that

$\begin{split}\mathcal{V}_{eff}&\simeq-\frac{N^{2}}{2}\Bigg{(}\frac{3T^{4}}{\pi^{2}}(\frac{M_{g}}{T})^{2}K_{2}(M_{g}/T)l_{N}(s)^{2}\\\
&\ \ \ \ \ \ \ \
+\frac{\pi^{2}T^{4}}{45}\Big{[}(-1+s)^{2}(-1-2s+2s^{2})-\frac{5(-1+s)^{2}s^{2}}{N^{2}}+\frac{s^{3}(-4+3s)}{N^{4}}\Big{]}\Bigg{)}+f(N,T),\end{split}$
(40)

which is somehow a hybrid of the PL model and matrix model.

Usually, $l_{N}={\rm tr}\hat{L}_{F}/N$ as a function of $s$ is complicated. By
definition, we can find that

$l_{N}(s)=\left\\{\begin{aligned}
\frac{1}{N}\bigg{(}1+2\sum_{i=1}^{\frac{N-1}{2}}\cos(2\pi\frac{i}{N}s)\bigg{)}&,&N\
{\rm is\ odd},\\\
\frac{2}{N}\bigg{(}\sum_{i=1}^{\frac{N}{2}}\cos(2\pi\frac{2i-1}{2N}s)\bigg{)}&,&N{\rm\
is\ even}.\end{aligned}\right.$ (41)

The summation can be implemented by writing $\cos(nx)={\rm Re}\exp{(inx)}$,
translating it to the geometric series, and then we obtain the simple
expression

$l_{N}(s)=\frac{1}{N}\frac{\sin(\pi s)}{\sin(\pi s/N)},$ (42)

which holds both for odd and even $N$. In particular, when the color number
approaching infinity, PL takes the limit ${\sin(\pi s)}/({\pi s})$.

####  III.3.2 Deconfinment phase transition

Now we arrive the effective potential which can be used to study deconfinement
phase transition for any color number $N$. For a sufficiently large $N$, one
can simply use the following rescaled potential

$\begin{split}\frac{\mathcal{V}_{N}(s,T)}{N^{2}/2}&\simeq-\frac{3T^{4}}{\pi^{2}}\left(\frac{M_{g}}{T}\right)^{2}K_{2}(M_{g}/T)\bigg{[}\frac{\sin(\pi
s)}{\pi
s}\bigg{]}^{2}+\frac{\pi^{2}T^{4}}{45}(s-1)^{2}(1+2s-2s^{2}),\end{split}$ (43)

which is $N$-independent. To find the vacuum position of this potential, one
should calculate the derivative of this potential with respect to $s$, to
solve the following tadpole equation

$\begin{split}&\frac{6}{\pi^{2}}\left(\frac{M_{g}}{T}\right)^{2}K_{2}(M_{g}/T)\frac{\pi\sin(\pi
s)\left(\cos(\pi s)-\frac{1}{N}\sin(\pi s)\cot\left(\frac{\pi
s}{N}\right)\right)}{N^{2}\sin^{2}(\pi s/N)}\\\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad+\frac{\pi^{2}}{45}s(s-1)\Big{[}2(4s-5)-\frac{10}{N^{2}}(2s-1)+\frac{12}{N^{4}}s\Big{]}=0.\end{split}$
(44)

It has two obvious solutions: 1) $s=0$, the deconfined vacuum position at high
temperature; 2) $s=1$, the confining vacuum position, which is consistent with
the $Z_{N}$ symmetry argument: The confining vacuum should preserve the
$Z_{N}$ symmetry and then the PL value must be $l_{N}(s=1)=0$. In our model,
this is trivially satisfied for the potential from gluons, which contributes
the quadratic term $l_{N}^{2}$. However, the inverted Weiss term is not a
polynomial of $l_{N}$, and therefore $s=1$ (namely $l_{N}=0$) being its
extremum is nontrivial. It is attributed to the eigenvalue repulse of the
potential. Eq. (44) may also admit solutions for $s\neq 1$, the candidates for
the deconfined vacuum at the lower temperature.

Figure 1: Left panel: The potential behavior for $SU(2)$ theory around the
critical temperature. Right panel: The potential shape for $SU(N\geq 3)$ at
the critical temperature.

For a given $N$, the shape of the potential Eq.(40) is solely determined by
the single dimensionless parameter $M_{g}(T)/T$. Then, we plot the shape
function ${\cal U}_{N}\equiv{\cal V}_{N}(s,T)/T^{4}$ at different values of
$M_{g}(T)/T$, to search the vacuum structure; the plots are displayed in Fig.
1. For the $N=2$ case, there is only one minimum at $s\neq 1$, and eventually,
only the minimum at $s=1$ survives as the increasing $M_{g}(T)/T$; see the
left panel of Fig. 1. So, in this case, the transition from the deconfined
phase to the confining phase is crossover. On the contrary, for the case
$N\geq 3$, ${\cal U}_{N}$ has two minimums when $M_{g}(T)/T$ approaches 2.7,
with one located at $s=1$ and the other one at $s\neq 1$. It implies that the
deconfinement phase transition is first order. Furthermore, via the degeneracy
condition, we can determine the critical temperature $T_{c}$ in the unit of
$M_{g}(T_{c})$; see the right panel of Fig. 1.

The phase transition behavior predicted by the model is consistent with the
results of lattice simulation. Hence, it is of importance to understand what
causes the qualitative difference between the shape functions for $N=2$ and
$N\geq 3$. To that end, we investigate the shape function near $s=1$, and
hence it is convenient to set $t=1-s$; then, we expand it around $t=0$,keeping
the irrelevant terms up to ${\cal O}(t^{5})$,

$\begin{split}{\cal
U}_{N}=&t^{2}\left(\frac{\pi^{2}{c_{g}}\csc^{2}\frac{\pi}{N}}{N^{2}}+\frac{6}{N^{4}}-\frac{5}{N^{2}}-1\right)+t^{3}\left(\frac{2\pi^{3}c_{g}\cot\frac{\pi}{N}\csc^{2}\frac{\pi}{N}}{N^{3}}-\frac{8}{N^{4}}+\frac{10}{N^{2}}-2\right)\\\
&+t^{4}\left[c_{g}\frac{\pi^{4}\csc^{2}\frac{\pi}{N}}{3N^{2}}\left(\frac{3}{N^{2}}+\frac{9\cot^{2}\frac{\pi}{N}}{N^{2}}-1\right)+\frac{3}{N^{4}}-\frac{5}{N^{2}}+2\right],\end{split}$
(45)

where $c_{g}\equiv(M_{g}/T)^{2}K_{2}(M_{g}/T)$. The special property of $N=2$
case is that the cubic term vanishes and therefore there is no barrier. For
$N>2$, the cubic term is present and moreover carries a positive coefficient
(attributed to the gluon potential), and as a consequence ${\cal U}_{N}$ is
able to give rise to the first order deconfinement phase transition.

We end up this section with a comment on the Haar-type model Mocsy:2003qw ;
Fukushima:2003fw ; Roessner:2006xn ; Fukushima:2008wg , whose potential is
supposed to resemble

$V=-a(T)l^{2}/2+b(T)\log H_{N}(L).$ (46)

The first term of this potential comes from the kinetic term of $SU(N)$ theory
which also exists in our model. The second term, characterzied by $H_{N}(L)$,
actually is known as the Vandermonde determinant interaction of the $SU(N)$
theory; it appears mathematically to define an integration over a continuous
group, which requires an invariant Haar group measure. Such an interaction is
non perturbative, and is argued to be consistent with the picture of ghost
dominance. Explicitly, the integrand functional of our effective potential
$-\Omega(M_{g}=0)$ in the infrared regime ($E,p\rightarrow 0$) resembles the
Vandermonde determinant interaction. In other words, roughly speaking it is a
part of the ghost contribution.

##  IV Thermodynamics

Although the model can surprisingly describe the order of deconfinement phase
transition for any $N$, it is still important to check if it is able to
correctly account for the thermodynamics above $T_{c}$, in particular in the
semi-QGP region around $1.4T_{c}$ where the nonperturbative effect is
significant. We have to rely on the temperature varying $M_{g}(T)$, with
$M_{g}(T_{c})$ fixed (traded with the critical temperature $T_{c}$), to do
this job.

We should start from fitting latent heat. Because the quasi gluon mass is
temperature dependent, the latent heat is sensitive to $dM_{g}/dT$ at $T_{c}$.
Actually, the latent heat data can fix the value of this derivative at
$T_{c}$, which is crucial to fit $M_{g}(T)$ via thermodynamics.

###  IV.1 Latent heat and determination of $dM_{g}(T)/dT$ at $T_{c}$

From the thermodynamics, it is known that the latent heat $L_{N}$ released
during the first order phase transition is the energy density difference
between two vacua, $L_{N}=\varepsilon_{d}-\varepsilon_{c}$, with subscripts
$d$ and $c$ denoting for the deconfined and confining vacuum, respectively.
Then, using the second law of thermodynamics, one can find that

$\begin{split}L_{N}=T_{c}\frac{\partial\Delta P}{\partial T}-\Delta
P=-T_{c}\frac{\partial\Delta\mathcal{V}(T)}{\partial
T}|_{T=T_{c}}+\Delta\mathcal{V}_{N}(T_{c}),\end{split}$ (47)

where $\Delta\mathcal{V}_{N}=\mathcal{V}_{d}-\mathcal{V}_{c}$ is the potential
energy difference, vanishing at $T_{c}$. As a consequence, the latent heat is
determined by the entropy part, i.e.,
$L_{N}=-T_{c}\frac{\partial\Delta\mathcal{V}(T)}{\partial T}|_{T=T_{c}}$.

Note that so far we can not guarantee the confining vacuum at $s=1$ is indeed
the absolute minimum below $T_{c}$, but it is simply an inference of the
requirement that latent heat should be positive: It means
$\partial\mathcal{V}_{d}/\partial T_{c}<\partial\mathcal{V}_{c}/\partial
T_{c}$, and moreover at $T_{c}$ two vacua is degenerate,
$\mathcal{V}_{d}=\mathcal{V}_{c}$, so, below $T_{c}$ one indeed has
$\mathcal{V}_{d}<\mathcal{V}_{c}$.

We are now in the position to calculate the latent heat in our model. The
straightforward calculation of the temperature derivative of the effective
potential gives

$\begin{split}\frac{\partial\mathcal{V}(s,T)}{\partial
T}&=-\frac{3M_{g}^{2}N^{2}}{2\pi^{2}}l_{N}^{2}\Bigg{[}4TK_{2}(M_{g}/T)+K_{1}(M_{g}/T)\bigg{(}M_{g}-T\frac{dM_{g}}{dT}\bigg{)}\Bigg{]}\\\
&-N^{2}\frac{2\pi^{2}}{45}\Bigg{[}(-1+s)^{2}(-1-2s+2s^{2})-\frac{5(-1+s)^{2}s^{2}}{N^{2}}+\frac{s^{3}(-4+3s)}{N^{4}}\Bigg{]}T^{3},\end{split}$
(48)

In our model, the confining vacuum is always located at $s=1$ or $l_{N}=0$,
and thus the above derivative vanishes in this vacuum. Then, the latent heat
is determined by the contribution from the deconfined vacuum,

$\begin{split}\frac{L_{N}}{(N^{2}-1)T_{c}^{4}}&=\frac{N^{2}}{N^{2}-1}\Bigg{(}\frac{3l_{N,d}^{2}}{2\pi^{2}}\bigg{(}\frac{M_{g}}{T_{c}}\bigg{)}^{2}\Bigg{[}4K_{2}(M_{g}/T_{c})+K_{1}(M_{g}/T_{c})\bigg{(}\frac{M_{g}}{T_{c}}-\frac{dM_{g}}{dT_{c}}\bigg{)}\Bigg{]}\\\
&+\frac{2\pi^{2}}{45}\Bigg{[}(-1+s_{d})^{2}(-1-2s_{d}+2s_{d}^{2})-\frac{5(-1+s_{d})^{2}s_{d}^{2}}{N^{2}}+\frac{s_{d}^{3}(-4+3s_{d})+1}{N^{4}}\Bigg{]}\Bigg{)},\end{split}$
(49)

where $s_{d}$ (or $l_{N,d}$) is the value of $s$ (or $l_{N}$) in the
deconfined vacuum, numerically calculated by virtue of the tadpole Eq. (44),
shown in Table. 1.

On the other hand, for $N=3,...8$, the current lattice data gives the
following behavior of latent heat Datta:2010sq

$\frac{L_{N}}{(N^{2}-1)T_{c}^{4}}\simeq 0.388-\frac{1.61}{N^{2}},$ (50)

We require the calculated latent heat Eq. (49) to fit it. For the given $N$,
Eq. (49) just contains a single parameter, $dM_{g}(T)/dT$ at $T_{c}$, and
therefore its value can be uniquely fixed. We show the results in Table. 1.
The resulting values typically are around $-10$ for all $N$, indicating a
sharp increasing of quasi-gluon mass as the temperature drops down to $T_{c}$
from above. This is a well-understood behavior since it can be regarded as a
sign of the “strongest” nonperturbative effect near $T_{c}$.

Color number |  $N=3$ |  $N=4$ |  $N=5$ |  $N=6$ |  $N=7$ |  $N=8$ |  $N\rightarrow\infty$
---|---|---|---|---|---|---|---
$s_{d}$ |  0.5605 |  0.5186 |  0.5073 |  0.5033 |  0.5016 |  0.5009 |  0.5004
$l_{d}$ |  0.5910 |  0.6300 |  0.6380 |  0.6398 |  0.6400 |  0.6396 |  0.6367
$M_{g}(T_{c})/T_{c}$ |  2.7499 |  2.7203 |  2.7126 |  2.7099 |  2.7088 |  2.7083 |  2.7077
$dM_{g}(T_{c})/dT_{c}$ |  -5.7727 |  -7.9951 |  -9.2891 |  -10.0965 |  -10.6261 |  -10.9954 |  -12.3376
$L_{N}/(N^{2}-1)T_{c}^{4}$ |  0.2091 |  0.2874 |  0.3236 |  0.3433 |  0.3551 |  0.3628 |  0.3880
Table 1: Effective mass, the position of deconfined vacuum and, latent heat
around the critical temperature for different color number $N$.

###  IV.2 Fit $M_{g}$ with the thermal quantity using machine learning

According to the original idea of QPM, the proper temperature dependence
beyond $T_{c}$ of quasi-gluon mass is supposed to successfully explain the
thermodynamics of the hot PYM system up to the high $T$ region. Here, the main
thermodynamic observables of interest are the pressure $p$, the energy density
$\epsilon$, and the entropy density $s$. Actually, they are not independent
quantities. In particular, if one has $p$, then $\epsilon$ and $s$ can be
calculated by the second law of thermodynamics

$\begin{split}\epsilon=T\frac{dp}{dT}-p,\ \ \ \
s=\frac{\epsilon+p}{T}.\end{split}$ (51)

The one loop calculation leads to $p=-{\cal V}_{eff}$, given in Eq. (40).
Currently, their lattice data is available only for $N=3,4,6$ Datta:2010sq .
However, as stated in the introduction, the lattice data demonstrates $N$
scaling property, which means $P_{M}=\frac{M^{2}-1}{N^{2}-1}P_{N}$ and the
latent Eq. (50), and thus we also have “data” for other $N$ values by simple
extrapolation, for instance to $N=5,8$ used later.

In the QPM, it is known that the SB limit can be trivially recovered. The most
challenging range is the so-called semi-QGP region $T\in(T_{c},3T_{c})$, where
the deviation to the blackbody behavior becomes more and more remarkable as
$T$ approaching $T_{c}$. In the previous discussion, we have used effective
potential Eq. (38), which is based on the high and low temperature expansion,
to analyze the phase transition at $T_{c}$. Nevertheless, we do not have such
a simple analytic expression to analyze thermadynamics. It is well expected
that the low temperature expansion just holds very near $T_{c}$ and soon
becomes not reliable in the higher temperature region. Hence, we should use
its complete expression:

$\mathcal{V}_{eff}\begin{aligned}
&=\frac{T^{4}}{2\pi^{2}}\int_{0}^{\infty}dxx^{2}\Bigg{(}(N-1)\log[1-e^{-\hat{E}(x,M_{g},T)}]\\\
&+\sum_{i=1}^{N}(N-i)\log[1+e^{-2\hat{E}(x,M_{g},T)}-2e^{-\hat{E}(x,M_{g},T)}\cos(2\pi\frac{is}{N})]\Bigg{)}+...,\end{aligned}$
(52)

where $\hat{E}(x,M_{g},T)=\sqrt{x^{2}+(M_{g}/T)^{2}}$ and dots denote for the
remaining term that does not need summation.

Then, we try to obtain the interpolation function of the fitted effective
gluon mass $M_{g}(T)$ for $N=3,4,5,6,8$, through the method of machine
learning. Physical Information Neural Network Raissi:2019jcp  provides us with
a flexible and accurate method for the fitting task. It treats functions of
any complexity under fitting as a neural network, and the training goal is
making the neural network satisfying the required partial differential
relationships (such as partial differential equations and boundary conditions)
and the given data points values. In our work, we use two separate deep neural
networks $M_{g}(T)$ and $s(T)$ for the fitting task, our training goal is
making $M_{g}(T)$ and $s(T)$ to satisfy:

  * •

the extreme condition for the deconfined vacuum at $T=T_{c}$;

  * •

the degeneracy between two the deconfined vacuum and the confining vacuum;

  * •

mass parameter relationship in Table. 1 and

  * •

the lattice data for thermodynamics.

We implement the task by using TensorFlow2.0 10.5555/3026877.3026899 , both
$M_{g}(T)$ and $s(T)$ contain the same hidden layers size. For the complexity
of our problem, a two-step training is used: We pretrain $M_{g}(T)$ and $s(T)$
to fit the lattice data first, and then fine adjust $M_{g}(T)$ and $s(T)$ to
satisfy all the fitting requirements above. For more details, please check the
code in Github 666https://github.com/JGuoHep/QuasiParticle. The fitted
$M_{g}(T)$ is shown in the first panel of Fig. 2, and the perfect fitting of
pressure above $T_{c}$ is displayed in other panels of Fig. 2.

The $M_{g}(T)$ is supposed to depend on $N$: Although $M_{g}(T_{c})/T_{c}$ is
almost universal determined by the condition of degeneracy, $dM_{g}(T)/dT$
takes different values at $T_{c}$ for different color number for the sake of
correct latent heat, see Table. 1. However, it is found that the fitted
$M_{g}(T)$ are almost the same, which leads us to conjecture that this is an
universal behavior for all $N$ 777For large $N$ this is trivial, because the
$N$ dependence of the observables in our model is scaled out, well consistent
with the lattice data. But it is not trivial that it is true also for $N=3,4$.
By contrast, in the polynomial model Kang:2021epo , the fitting parameters in
the small $N$ cases are very different than those in the large $N$ cases. We
guess it is attributed to the exponential dependence of the fitting parameter
$M_{g}(T)$.. By the way, one can check the invalidation of low temperature
expansion in the region $T\gtrsim 1.4T_{c}$: The ratio $M_{g}/T$ drops to
$\approx 1.5$ as $T$ increases to $1.4T_{c}$, and then from Eq. (27) one can
see that the next leading order is only suppressed by a factor $K_{2}(2\times
1.5)/K_{2}(1.5)\sim 0.1$.

Figure 2: The first panel: The fitted $M_{g}/T$ as a function of temperature for $N=3,4,5,6,8$. The second to the sixth panels: Fitting $p/T^{4}$ in our model for various color number. $\alpha$ | $\beta$ | $\gamma$ | RMSD
---|---|---|---
0.029534 | 1.130884 | 1.541299 | 0.015707
0 | 1.186505 | 1.570699 | 0.016224
void | $\lambda$ | $T_{s}$ | RMSD
void | 10.843298 | -8.336149 | 0.081482
Table 2: $\alpha$, $\beta$ and $\gamma$ are fitting parameters in the
qausigluon mass ansatz Eq.(53), while $\lambda$ and $T_{s}$ are fitting
parameters in the conventional ansatz Eq.(4).

Actually, the $N$ universal behavior of quasigluon mass is encoded in the
quasigluon mass in the HTLpt; see the formula Eq. (4) where $N$ cancels. At
this point, our model is consistent with the HTLpt effective mass. So, it is
anticipated that the interpolation function can be fitted by the $M_{g}(T)$
with the function given in Eq.(4), with two parameters $\lambda$ and $T_{s}$.
We also try another function with three parameters

$M_{g}(T)=\alpha T+\beta T/\log(\gamma T/T_{c}).$ (53)

which is recently adopted in Ref. Islam:2021qwh . Note that unlike the
conventional $M_{g}(T)$ ansatz, which simply goes to the HTLpt quasigluon mass
in the high $T$ region, Eq. (53) does not. The fitted parameters for both
functions of $M_{g}(T)$ are shown in Table. 2. The latter has better quality,
which can be seen from the comparison in two panels of Fig. 3. This may raise
the issue of well consistence between our model with the HTLpt in the higher
$T$ region, and we will come back to this point in the Section of conclusion
and discussion. Besides, for the function Eq. (53), from Table. 2 one can see
that the values of the $\alpha$ parameter are far smaller than the other two
parameters, which means that it is almost irrelevant to fitting. So, we tried
the fitting with the vanishing $\alpha$, to find that it works equally well.

Figure 3: The fitting for the quasigluon mass ansatz Eq.(53) (left panel) and
Eq. (4) (right panel), where the fitting paramers are given in Table. 2.

##  V Conclusion and discussion

The HTL resummation in the quasi-particle picture reveals that QGP is a pool
of weakly interacting quasigluons for $T\gtrsim 2T_{c}$. Such a picture is
further used in the QPM to describe QCD thermodynamics down to $T_{c}$ and
works fairly well. The crucial idea is that the quasigluon mass could “absorb”
strong interaction and merely leaves weak interactions on quasigluons. In this
work we attempt to embed this idea to the massive PYM Reinosa:2014ooa ,
introducing a temperature-dependent quasigluon mass in the effective $SU(N)$
PYM Lagaragian Eq. (8). Via the standard perturbative calculation, we obtain
an effective model that can successfully explain the critical behavior for any
$N$, not also the first order deconfinement phase transition for $N>2$ but
also the crossover for $N=2$. Moreover, the lattice data of thermodynamics can
be fitted via the single parameter $M_{g}(T)$, which is found to demonstrate
the $N$-universal behavior, based on the available case $N=3,4,6$. This is
supported by the HTLpt quasigluon mass, but now is extended to the semi-QGP
region, and might convey some secrets of the non-perturbative effects. We look
forward to the future lattice data for other $N$, in particular, $N=5,8$ whose
“lattice data” is obtained by extrapolation via the $N$-scaling law, to test
the universal quasigluon mass conjecture.

Fitting $M_{g}(T)$ via a function that well matches with the HTLpt quasigluon
mass does not have a very good quality, and it may be improved by considering
the dressing propagator of the gluons Braun:2007bx . Then, the modified model
contains more parameter and have the potential to deal with more detailed
problems.

We are capable of conducting a unified analysis of all $N$, depending on the
assumption of uniform eigenvalue distribution of the temporal background,
which reduces the effective potential to the one-dimensional case. But it is
based on the eigenvalue repulsion and a more solid argument may be necessary.

Acknowledgements

This work is supported in part by the National Key Research and Development
Program of China Grant No. 2020YFC2201504 and in part by the National Science
Foundation of China (11775086).

Note added

Right before the submission of this work, the work by Fu-Peng Li etc.
Li:2022ozl  appeared on arxiv. They also utilize the machine learning to
reconstruct QCD equation of state in the QPM picture, which may have partial
overlap with our work.

##  Appendix A Derivation of the generating function in the Landau-DeWitt
gauge

The complete Faddeev-Popov Lagrangian in the Landau-DeWitt gauge reads

$\mathcal{L}=-\frac{1}{2g^{2}}{\rm
tr}(F_{\mu\nu}F^{\mu\nu})+\bar{D}_{\mu}\bar{c}^{a}D^{\mu}c^{a}+ih^{a}\bar{D}_{\mu}\hat{A}^{\mu,a}+\frac{1}{2}M_{g}(T)A^{a}_{tr,\mu}A^{a,\mu}_{tr},$
(54)

We are considering the constant background
$\bar{A}^{a}_{\mu}=\bar{A}^{a}_{0}\delta_{\mu 0}$ and keep only the quadratic
terms. Then the action can be split into two parts $S_{A,h}$ and $S_{c}$,
where

$\begin{split}&S_{A,h}=\int
d^{4}x\left[\frac{1}{2}\hat{A}_{\mu}^{a}(D^{-1})_{ab}\hat{A}^{\mu,b}+ih^{a}\bar{D}_{\mu}\hat{A}^{\mu,a}\right]\\\
&S_{c}=\int d^{4}x[\bar{D}_{\mu}\bar{c}^{a}\bar{D}^{\mu}c^{a}],\end{split}$
(55)

where $D^{-1}$ is given in Eq. (12).

Now we come to deal with the first part of the action. We can rewrite this
action in its color diagonalization basis $\tilde{A}^{a}_{\mu}$ and
$\tilde{h}^{a}$, to get

$\begin{split}S_{A,h}&=\int
d^{4}x\left[\frac{1}{2}\tilde{A}_{\mu}^{a}(\tilde{D}^{-1})_{a}\tilde{A}^{\mu,a}+i\tilde{h}^{a}\tilde{D}_{\mu}^{a}\tilde{A}^{\mu,a}\right]\\\
&=\int\frac{d^{4}p}{(2\pi)^{4}}\left[\frac{1}{2}\tilde{A}_{\mu}^{a}(\tilde{D}^{-1})_{a}\tilde{A}^{\mu,a}+i\tilde{h}^{a}\tilde{D}_{\mu}^{a}\tilde{A}^{\mu,a}\right],\end{split}$
(56)

where $(\tilde{D}^{-1})_{a}$ and $\tilde{D}_{\mu}^{a}$ in the momentum space
are respectively given by

$\begin{split}&\tilde{D}^{-1}_{a}(p)\equiv
D(M_{g})_{a}^{-1}=(p_{0}-A_{a})^{2}+|\vec{p}|^{2}+M_{g}^{2},\\\
&\tilde{D}_{\mu}^{a}\tilde{D}^{\mu,a}\equiv
D(0)_{a}^{-1}=(p_{0}-A_{a})^{2}+|\vec{p}|^{2}.\end{split}$ (57)

To integrate this action through path integration, we must do the quadratic
partition between $\tilde{A}^{a}$ and $\tilde{h}^{a}$. After a tedious
quadratic partition, the action takes the form of

$\displaystyle
S_{A,h}=\int\frac{d^{4}p}{(2\pi)^{4}}\frac{1}{2}\left[(\tilde{D}^{-1})_{a}(\tilde{A}_{\mu}^{a}+i\frac{\tilde{D}^{a}_{\mu}}{\sqrt{2\tilde{D}^{-1}}_{a}}\tilde{h}^{a})(\tilde{A}^{\mu,a}+i\frac{\tilde{D}^{\mu,a}}{\sqrt{2\tilde{D}^{-1}}_{a}}\tilde{h}^{a})+\frac{\tilde{D}_{\mu}^{a}\tilde{D}^{\mu,a}}{\tilde{D}^{-1}_{a}}\tilde{h}^{a}\tilde{h}^{a}\right].$
(58)

Now define a new field
$\mathcal{A}^{a}_{\mu}=\tilde{A}_{\mu}^{a}+i\frac{\tilde{D}^{a}_{\mu}}{\sqrt{2\tilde{D}^{-1}}_{a}}\tilde{h}^{a}$,
and one can rewrite the original mixed action as

$S_{A,h}=\int\frac{d^{4}p}{(2\pi)^{4}}\frac{1}{2}\left[\mathcal{A}^{a}_{\mu}D(M_{g})^{-1}_{a}\mathcal{A}^{\mu,a}+\frac{D(0)^{-1}_{a}}{D(M_{g})^{-1}_{a}}\tilde{h}^{a}\tilde{h}^{a}\right].$
(59)

The redefined Nakanishi-Lautrup field $\tilde{h}^{a}$ now gains a mass, and
its propagator is a combination of the masssive and massless propagators,
which is a result of the Landau-Dewitt gauge. Then the 1-loop effective action
is given by

$\begin{split}\log{Z}_{A,h}&=\frac{1}{2}\log\det\left[\frac{\delta^{2}S_{h,A}}{\delta\phi_{i}^{a}\phi_{j}^{a}}\right]\\\
&=2\log\det(D(M_{g})^{-1})-\frac{1}{2}\log\det(D(M_{g})^{-1})+\frac{1}{2}\log\det(D(0)^{-1})\\\
&=\frac{3}{2}\log\det(D(M_{g})^{-1})+\frac{1}{2}\log\det(D(0)^{-1}).\end{split}$
(60)

where $\phi^{a}_{i}$ represent $\\{\mathcal{A}^{a}_{\mu},\tilde{h}^{a}\\}$.
Note that there is a overall factor 4 for the $\mathcal{A}^{a}_{\mu}$
contribution, denoting for four massive modes. But the Nakanishi-Lautrup field
cancels one massive mode and effectively just leaves one massless mode.

The massless ghost contribution, taking into account its statistics, is simply
given by

$\log{Z}_{c}=-\log\det(D(0)^{-1}).$ (61)

Its contribution is halved due to the massless mode of the Nakanishi-Lautrup
field. Finally, the total effective action is

$\begin{split}\log{Z}=\frac{3}{2}\log\det(D(M_{g})^{-1})-\frac{1}{2}\log\det(D(0)^{-1}).\end{split}$
(62)

##  Appendix B Calculating the pure gluonic generating function: the $SU(4)$
sample

In this appendix we present the details of calculating the pure glunoic part,
i.e., the first term of the second line of Eq. (60), specified to $SU(4)$. Its
Cartan generators are

$T^{3}=\frac{1}{2}\left(\begin{matrix}1&0&0&0\\\ 0&-1&0&0\\\ 0&0&0&0\\\
0&0&0&0\end{matrix}\right),\ \
T^{8}=\frac{1}{2\sqrt{3}}\left(\begin{matrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&-2&0\\\
0&0&0&0\end{matrix}\right),\ \
T^{15}=\frac{1}{2\sqrt{6}}\left(\begin{matrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
0&0&0&-3\end{matrix}\right).$ (63)

Now the propagators take the form of (the quasigluon mass can be trivially
included)

$(D^{-1})_{ab}=\delta_{ab}p^{2}+2i\sum_{i=3,8,15}f_{abi}\bar{A}^{i}_{0}p_{0}-\sum_{i,j=3,8,15}f_{aci}f_{cbj}\bar{A}^{i}_{0}\bar{A}^{j}_{0}.$
(64)

After a careful calculation, one can get all the non-zero propagators

$\begin{split}&(D^{-1})_{1,1}=p^{2}+(\bar{A}^{3}_{0})^{2}\ \
(D^{-1})_{1,2}=2ip_{0}\bar{A}^{3}_{0}\\\
&(D^{-1})_{2,2}=p^{2}+(\bar{A}^{3}_{0})^{2}\ \
(D^{-1})_{1,2}=-2ip_{0}\bar{A}^{3}_{0}\\\ &(D^{-1})_{3,3}=p^{2}\\\
&(D^{-1})_{4,4}=p^{2}+\frac{1}{4}(\bar{A}^{3}_{0})^{2}+\frac{4}{3}(\bar{A}^{8}_{0})^{2}+\frac{\sqrt{3}}{2}\bar{A}^{3}_{0}\bar{A}^{8}_{0}\
\ (D^{-1})_{4,5}=ip_{0}\bar{A}^{3}_{0}+\sqrt{3}ip_{0}\bar{A}^{8}_{0}\\\
&(D^{-1})_{5,5}=p^{2}+\frac{1}{4}(\bar{A}^{3}_{0})^{2}+\frac{4}{3}(\bar{A}^{8}_{0})^{2}+\frac{\sqrt{3}}{2}\bar{A}^{3}_{0}\bar{A}^{8}_{0}\
\ (D^{-1})_{5,4}=-ip_{0}\bar{A}^{3}_{0}-\sqrt{3}ip_{0}\bar{A}^{8}_{0}\\\
&(D^{-1})_{6,6}=p^{2}+\frac{1}{4}(\bar{A}^{3}_{0})^{2}+\frac{4}{3}(\bar{A}^{8}_{0})^{2}-\frac{\sqrt{3}}{2}\bar{A}^{3}_{0}\bar{A}^{8}_{0}\
\ (D^{-1})_{6,7}=-ip_{0}\bar{A}^{3}_{0}+\sqrt{3}ip_{0}\bar{A}^{8}_{0}\\\
&(D^{-1})_{7,7}=p^{2}+\frac{1}{4}(\bar{A}^{3}_{0})^{2}+\frac{4}{3}(\bar{A}^{8}_{0})^{2}-\frac{\sqrt{3}}{2}\bar{A}^{3}_{0}\bar{A}^{8}_{0}\
\ (D^{-1})_{7,6}=ip_{0}\bar{A}^{3}_{0}-\sqrt{3}ip_{0}\bar{A}^{8}_{0}\\\
&(D^{-1})_{8,8}=p^{2}\\\
&(D^{-1})_{9,9}=p^{2}+\frac{1}{4}(\bar{A}^{3}_{0})^{2}+\frac{1}{12}(\bar{A}^{8}_{0})^{2}+\frac{1}{12}(\bar{A}^{15}_{0})^{2}+\frac{1}{2\sqrt{3}}\bar{A}^{3}_{0}\bar{A}^{8}_{0}+\frac{1}{2\sqrt{3}}\bar{A}^{3}_{0}\bar{A}^{15}_{0}+\frac{1}{6}\bar{A}^{8}_{0}\bar{A}^{15}_{0}\\\
&(D^{-1})_{9,10}=ip_{0}\bar{A}_{0}^{3}+\frac{i}{\sqrt{3}}p_{0}+\bar{A}_{0}^{8}+\frac{i}{\sqrt{3}}p_{0}\bar{A}_{0}^{15}\\\
&(D^{-1})_{10,10}=p^{2}+\frac{1}{4}(\bar{A}^{3}_{0})^{2}+\frac{1}{12}(\bar{A}^{8}_{0})^{2}+\frac{1}{12}(\bar{A}^{15}_{0})^{2}+\frac{1}{2\sqrt{3}}\bar{A}^{3}_{0}\bar{A}^{8}_{0}+\frac{1}{2\sqrt{3}}\bar{A}^{3}_{0}\bar{A}^{15}_{0}+\frac{1}{6}\bar{A}^{8}_{0}\bar{A}^{15}_{0}\\\
&(D^{-1})_{10,9}=-ip_{0}\bar{A}_{0}^{3}-\frac{i}{\sqrt{3}}p_{0}+\bar{A}_{0}^{8}-\frac{i}{\sqrt{3}}p_{0}\bar{A}_{0}^{15}\\\
&(D^{-1})_{11,11}=p^{2}+\frac{1}{4}(\bar{A}^{3}_{0})^{2}+\frac{1}{12}(\bar{A}^{8}_{0})^{2}+\frac{1}{12}(\bar{A}^{15}_{0})^{2}-\frac{1}{2\sqrt{3}}\bar{A}^{3}_{0}\bar{A}^{8}_{0}-\frac{1}{2\sqrt{3}}\bar{A}^{3}_{0}\bar{A}^{15}_{0}+\frac{1}{6}\bar{A}^{8}_{0}\bar{A}^{15}_{0}\\\
&(D^{-1})_{11,12}=-ip_{0}\bar{A}_{0}^{3}+\frac{i}{\sqrt{3}}p_{0}+\bar{A}_{0}^{8}+\frac{i}{\sqrt{3}}p_{0}\bar{A}_{0}^{15}\\\
&(D^{-1})_{12,12}=p^{2}+\frac{1}{4}(\bar{A}^{3}_{0})^{2}+\frac{1}{12}(\bar{A}^{8}_{0})^{2}+\frac{1}{12}(\bar{A}^{15}_{0})^{2}-\frac{1}{2\sqrt{3}}\bar{A}^{3}_{0}\bar{A}^{8}_{0}-\frac{1}{2\sqrt{3}}\bar{A}^{3}_{0}\bar{A}^{15}_{0}+\frac{1}{6}\bar{A}^{8}_{0}\bar{A}^{15}_{0}\\\
&(D^{-1})_{12,11}=ip_{0}\bar{A}_{0}^{3}-\frac{i}{\sqrt{3}}p_{0}+\bar{A}_{0}^{8}-\frac{i}{\sqrt{3}}p_{0}\bar{A}_{0}^{15}\\\
&(D^{-1})_{13,13}=p^{2}+\frac{1}{12}(\bar{A}^{8}_{0})^{2}+\frac{1}{12}(\bar{A}^{15}_{0})^{2}-\frac{1}{6}\bar{A}^{8}_{0}\bar{A}^{15}_{0}\
\
(D^{-1})_{13,14}=-ip_{0}\bar{A}^{8}_{0}+\frac{i}{\sqrt{3}}p_{0}\bar{A}^{15}_{0}\\\
&(D^{-1})_{14,14}=p^{2}+\frac{1}{12}(\bar{A}^{8}_{0})^{2}+\frac{1}{12}(\bar{A}^{15}_{0})^{2}-\frac{1}{6}\bar{A}^{8}_{0}\bar{A}^{15}_{0}\
\
(D^{-1})_{14,13}=ip_{0}\bar{A}^{8}_{0}-\frac{i}{\sqrt{3}}p_{0}\bar{A}^{15}_{0}\\\
&(D^{-1})_{15,15}=p^{2}\end{split}$ (65)

It is observed that the $15\times 15$ propagator matrix in the color space
$D^{-1}_{a,b}$ is a block diagonal matrix, consisting of three diagonal
elements $p^{2}$ corresponding to the Cartan part and six $2\times 2$
submatrices corresponding to the non-Cartan parts. Concretely, these six
matrices are

$M_{i}=\begin{bmatrix}p_{0}^{2}+\widetilde{A}_{i}^{2}-|\vec{p}|^{2}&ip_{0}\widetilde{A}_{i}\\\
-ip_{0}\widetilde{A}_{i}&p_{0}^{2}+\widetilde{A}_{i}^{2}-|\vec{p}|^{2}\end{bmatrix},$
(66)

where $\widetilde{A}_{i}$ is a combination of $A_{0}$. From this expression we
can see that the eigenvalues of $M_{i}$ must be
$(p_{0}+A_{i})^{2}-|\vec{p}|^{2}$ and $(p_{0}+A_{i})^{2}+|\vec{p}|^{2}$. After
a unitary diagonalization we can get $\widetilde{D}^{-1}\equiv
U^{\dagger}D^{-1}U$ as

$\begin{split}&(\tilde{D}^{-1})_{1,1}=(p_{0}-\bar{A}_{0}^{3})^{2}-|\vec{p}|^{2},\
\ (\tilde{D}^{-1})_{2,2}=(p_{0}+\bar{A}_{0}^{3})^{2}-|\vec{p}|^{2}\\\
&(\tilde{D}^{-1})_{3,3}=(\tilde{D}^{-1})_{8,8}=(\tilde{D}^{-1})_{15,15}=p^{2}\\\
&(\tilde{D}^{-1})_{4,4}=[p_{0}-\frac{1}{2}(\bar{A}_{0}^{3}+\sqrt{3}\bar{A}_{0}^{8})]^{2}-|\vec{p}|^{2},\
\
(\tilde{D}^{-1})_{5,5}=[p_{0}+\frac{1}{2}(\bar{A}_{0}^{3}+\sqrt{3}\bar{A}_{0}^{8})]^{2}-|\vec{p}|^{2}\\\
&(\tilde{D}^{-1})_{6,6}=[p_{0}+\frac{1}{2}(\bar{A}_{0}^{3}-\sqrt{3}\bar{A}_{0}^{8})]^{2}-|\vec{p}|^{2},\
\
(\tilde{D}^{-1})_{7,7}=[p_{0}-\frac{1}{2}(\bar{A}_{0}^{3}-\sqrt{3}\bar{A}_{0}^{8})]^{2}-|\vec{p}|^{2}\\\
&(\tilde{D}^{-1})_{9,9}=[p_{0}-\frac{1}{2}(\bar{A}_{0}^{3}+\frac{1}{\sqrt{3}}\bar{A}_{0}^{8}+\frac{1}{\sqrt{3}}\bar{A}_{0}^{15})]^{2}-|\vec{p}|^{2}\\\
&(\tilde{D}^{-1})_{10,10}=[p_{0}+\frac{1}{2}(\bar{A}_{0}^{3}+\frac{1}{\sqrt{3}}\bar{A}_{0}^{8}+\frac{1}{\sqrt{3}}\bar{A}_{0}^{15})]^{2}-|\vec{p}|^{2}\\\
&(\tilde{D}^{-1})_{11,11}=[p_{0}+\frac{1}{2}(\bar{A}_{0}^{3}-\frac{1}{\sqrt{3}}\bar{A}_{0}^{8}-\frac{1}{\sqrt{3}}\bar{A}_{0}^{15})]^{2}-|\vec{p}|^{2}\\\
&(\tilde{D}^{-1})_{12,12}=[p_{0}-\frac{1}{2}(\bar{A}_{0}^{3}-\frac{1}{\sqrt{3}}\bar{A}_{0}^{8}-\frac{1}{\sqrt{3}}\bar{A}_{0}^{15})]^{2}-|\vec{p}|^{2}\\\
&(\tilde{D}^{-1})_{13,13}=(p_{0}-\frac{1}{2}\frac{\bar{A}_{0}^{8}-\bar{A}_{0}^{15}}{\sqrt{3}})^{2}-|\vec{p}|^{2},\
\
(\tilde{D}^{-1})_{14,14}=(p_{0}+\frac{1}{2}\frac{\bar{A}_{0}^{8}-\bar{A}_{0}^{15}}{\sqrt{3}})^{2}-|\vec{p}|^{2}.\end{split}$
(67)

This leads to the quadratic Lagrangian written as

$\mathcal{L}=-\frac{1}{2}\widetilde{A}_{\mu}^{a}(\widetilde{D}^{-1})_{a}\widetilde{A}^{\mu,a},$
(68)

with $\widetilde{D}^{-1}={\rm
diag}((p_{0}+A_{1})^{2}-|\vec{p}|^{2},(p_{0}+A_{2})^{2}-|\vec{p}|^{2},...,(p_{0}+A_{15})^{2}-|\vec{p}|^{2})={\rm
diag}\\{(p_{0}+A_{a})^{2}-|\vec{p}|^{2}\\}$ where $A_{a}$ is zero or opposite
numbers appearing in pairs. This structure is insured by the structural
constant $f_{abc}$. One can check this structure for other $SU(N)$ theory. For
example in $SU(3)$, the diagonal propagator is the same as the first eight
propagators of $SU(4)$.

##  Appendix C Summation over the thermal modes

In this appendix we explicitly implement the summation over the thermal modes
present in Eq. (17), rewritten as

$\log Z=2V{\rm tr_{c}}\int\frac{d^{3}\vec{p}}{(2\pi)^{2}}\nu(E_{g}),$ (69)

where we have introduced the function

$\begin{split}\nu(E_{g})\equiv\sum_{n=-\infty}^{\infty}\log[\tilde{D}^{-1}_{aa}]=\sum_{n=-\infty}^{\infty}\log[(\omega_{n}-A_{a})^{2}+E_{g}^{2}],\end{split}$
(70)

with $E_{g}^{2}=|\vec{p}|^{2}+M_{g}^{2}$. To pull out the object to be summed
from the logarithm, we differentiate $\nu(E_{g})$ with respect to $E_{g}$,

$\begin{split}\frac{\partial\nu(E_{g})}{\partial
E_{g}}&=\sum_{n=-\infty}^{\infty}\frac{2E_{g}}{(\omega_{n}-A_{a})^{2}+E_{g}^{2}}=\frac{1}{\pi
T}\sum_{n=-\infty}^{\infty}\frac{E_{g}/2\pi T}{(n-A_{a}/2\pi
T)^{2}+(E_{g}/2\pi T)^{2}}.\end{split}$ (71)

Such a series can be summed explicitly, to get

$\begin{split}\frac{\partial\nu(E_{g})}{\partial
E_{g}}=\frac{1}{T}\frac{\sinh^{2}\left(E_{g}/2\pi
T\right)}{\sin^{2}(A_{a}/2\pi T)+\sinh^{2}(E_{g}/2\pi
T)}\left(1+\frac{e^{-E_{g}/T}}{1-e^{-E_{g}/T}}\right).\end{split}$ (72)

Then, integrating both sides over $E_{g}$, we have

$\nu(E_{g})=\frac{E_{g}}{T}+2\log\left[\sqrt{\left(1-e^{-i\frac{A_{a}}{T}}e^{\frac{-E_{g}}{T}}\right)\left(1-e^{i\frac{A_{a}}{T}}e^{\frac{-E_{g}}{T}}\right)}\right]+(E_{g}\
\rm independent~{}terms).$ (73)

Using the identity $\log M=\log\det M$ and the fact that $\nu(E_{g})$ is a
diagonal matrix in the color space thus a simple trace operation, we obtain

$\log
Z=2V\int\frac{d^{3}\vec{p}}{(2\pi)^{2}}2\log\prod_{a}\left[\sqrt{\left(1-e^{-i\frac{A_{a}}{T}}e^{\frac{-E_{g}}{T}}\right)\left(1-e^{i\frac{A_{a}}{T}}e^{\frac{-E_{g}}{T}}\right)}\right],$
(74)

where we have ignored the infinite vacuum energy and $E_{g}$ independent
terms. One should notice that each $A_{a}$ is paired with another
$A_{b}=-A_{a}$. Eventually, the generating function can be written as a more
compacted form:

$\log Z=4V{\rm
tr_{c}}\int\frac{d^{3}\vec{p}}{(2\pi)^{2}}\log\left(1-\hat{L}_{A}e^{-E_{g}/T}\right),$
(75)

with $\hat{L}_{A}={\rm
diag}(\exp(-iA_{1}/T),\exp(-iA_{2}/T),...,\exp(-iA_{N^{2}-1}/T))$ where,
again, $\pm A_{a}$ pairly appear.

## References

  * (1) Z. Kang, J. Zhu and S. Matsuzaki, JHEP 09, 060 (2021).
  * (2) H. An, S. L. Chen, R. N. Mohapatra and Y. Zhang, JHEP 1003, 124 (2010).
  * (3) Z. Chacko, H. S. Goh and R. Harnik, Phys. Rev. Lett. 96, 231802 (2006).
  * (4) R. Blumenhagen, M. Cvetic, P. Langacker and G. Shiu, Ann. Rev. Nucl. Part. Sci. 55, 71 (2005).
  * (5) Z. Kang, Phys. Lett. B 801, 135149 (2020).
  * (6) P. Carenza, R. Pasechnik, G. Salinas and Z. W. Wang, [arXiv:2207.13716 [hep-ph]].
  * (7) J. Halverson, C. Long, A. Maiti, B. Nelson and G. Salinas, JHEP 05, 154 (2021).
  * (8) J. Kubo and M. Yamada, JHEP 10, 003 (2018).
  * (9) W. C. Huang, M. Reichert, F. Sannino and Z. W. Wang, Phys. Rev. D 104, no.3, 035005 (2021).
  * (10) E. Morgante, N. Ramberg and P. Schwaller, [arXiv:2210.11821 [hep-ph]].
  * (11) C. Ratti, M. A. Thaler and W. Weise, [arXiv:nucl-th/0604025 [nucl-th]].
  * (12) R. D. Pisarski, Nucl. Phys. A 702, 151-158 (2002).
  * (13) K. Fukushima and V. Skokov, Prog. Part. Nucl. Phys. 96, 154-199 (2017).
  * (14) S. Roessner, C. Ratti and W. Weise, Phys. Rev. D 75, 034007 (2007).
  * (15) K. Fukushima, Phys. Lett. B 591, 277-284 (2004).
  * (16) K. Fukushima, Phys. Rev. D 77, 114028 (2008) [erratum: Phys. Rev. D 78, 039902 (2008)].
  * (17) A. Mocsy, F. Sannino and K. Tuominen, Phys. Rev. Lett. 92, 182302 (2004).
  * (18) P. N. Meisinger, T. R. Miller and M. C. Ogilvie, Phys. Rev. D 65, 034009 (2002).
  * (19) P. N. Meisinger, M. C. Ogilvie and T. R. Miller, Phys. Lett. B 585, 149-154 (2004).
  * (20) A. Dumitru, Y. Guo, Y. Hidaka, C. P. K. Altes and R. D. Pisarski, Phys. Rev. D 83, 034022 (2011).
  * (21) A. Dumitru, Y. Guo, Y. Hidaka, C. P. K. Altes and R. D. Pisarski, Phys. Rev. D 86, 105017 (2012).
  * (22) J.O. Andersen, E. Braaten, and M. Strickland, Phys. Rev. Lett. 83, 2139 (1999); Phys. Rev. D 62, 045004 (2000); J.O. Andersen, E. Braaten, E. Petitgirard, and M. Strickland, Phys. Rev. D 66, 085016 (2002).
  * (23) V. Goloviznin and H. Satz, Z. Phys. C 57, 671-676 (1993).
  * (24) A. Peshier, B. Kampfer, O. P. Pavlenko and G. Soff, Phys. Lett. B 337, 235-239 (1994).
  * (25) A. Peshier, B. Kampfer, O. P. Pavlenko and G. Soff, Phys. Rev. D 54, 2399-2402 (1996).
  * (26) M. I. Gorenstein and S. N. Yang, Phys. Rev. D 52, 5206-5212 (1995).
  * (27) A. J. Helmboldt, J. Kubo and S. van der Woude, Phys. Rev. D 100, no.5, 055025 (2019).
  * (28) C. Caprini, M. Hindmarsh, S. Huber, T. Konstandin, J. Kozaczuk, G. Nardini, J. M. No, A. Petiteau, P. Schwaller and G. Servant, et al. JCAP 04, 001 (2016).
  * (29) C. Sasaki and K. Redlich, Phys. Rev. D 86, 014007 (2012).
  * (30) M. Ruggieri, P. Alba, P. Castorina, S. Plumari, C. Ratti and V. Greco, Phys. Rev. D 86, 054007 (2012).
  * (31) C. A. Islam, M. G. Mustafa, R. Ray and P. Singha, [arXiv:2109.13321 [hep-ph]].
  * (32) G. Curci and R. Ferrari, Nuovo Cim. A 32, 151-168 (1976).
  * (33) U. Reinosa, J. Serreau, M. Tissier and N. Wschebor, Phys. Lett. B 742, 61-68 (2015).
  * (34) D. M. van Egmond, U. Reinosa, J. Serreau and M. Tissier, SciPost Phys. 12, no.3, 087 (2022).
  * (35) U. Reinosa, J. Serreau, M. Tissier and N. Wschebor, Phys. Rev. D 91, 045035 (2015).
  * (36) U. Reinosa, J. Serreau, M. Tissier and N. Wschebor, Phys. Rev. D 93, no.10, 105002 (2016).
  * (37) J. O. Andersen, M. Strickland and N. Su, Phys. Rev. Lett. 104, 122003 (2010).
  * (38) J. O. Andersen, L. E. Leganger, M. Strickland and N. Su, Phys. Lett. B 696, 468-472 (2011).
  * (39) R. A. Schneider and W. Weise, Phys. Rev. C 64, 055201 (2001).
  * (40) P. Castorina, D. E. Miller and H. Satz, Eur. Phys. J. C 71, 1673 (2011).
  * (41) P. Castorina, V. Greco, D. Jaccarino and D. Zappala, Eur. Phys. J. C 71, 1826 (2011).
  * (42) N. Weiss, Phys. Rev. D 24, 475 (1981).
  * (43) P. Alba, W. Alberico, M. Bluhm, V. Greco, C. Ratti and M. Ruggieri, Nucl. Phys. A 934, 41-51 (2014).
  * (44) P. M. Lo, K. Redlich and C. Sasaki, Phys. Rev. D 103, no.7, 074026 (2021).
  * (45) P. N. Meisinger and M. C. Ogilvie, Phys. Rev. D 65, 056013 (2002).
  * (46) V. N. Gribov, Nucl. Phys. B 139, 1 (1978) doi:10.1016/0550-3213(78)90175-X
  * (47) J. Heffner, H. Reinhardt and D. R. Campagnari, Phys. Rev. D 85, 125029 (2012).
  * (48) U. Reinosa, [arXiv:2009.04933 [hep-th]].
  * (49) Y. B. Ivanov, V. V. Skokov and V. D. Toneev, Phys. Rev. D 71, 014005 (2005).
  * (50) R. D. Pisarski, Phys. Rev. D 74, 121703 (2006).
  * (51) S. Datta and S. Gupta, Phys. Rev. D 82, 114505 (2010).
  * (52) M. Raissi and P. Perdikaris and G.E. Karniadakis, Journal of Computational Physics 378, 686-707 (2019).
  * (53) M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean et al., _Tensorflow: A system for large-scale machine learning_ , in _Proceedings of the 12th USENIX Conference on Operating Systems Design and Implementation_ , OSDI’16, (USA), p. 265-283, USENIX Association, 2016.
  * (54) J. Braun, H. Gies and J. M. Pawlowski, Phys. Lett. B 684, 262-267 (2010).
  * (55) F. P. Li, H. L. Lü, L. G. Pang and G. Y. Qin, [arXiv:2211.07994 [hep-ph]].